MpGeometry.cc

Go to the documentation of this file.

#include "MuonProphet.h"
 
#include "GaudiKernel/Point3DTypes.h"
#include "GaudiKernel/Vector3DTypes.h"

StatusCode MuonProphet::geometryCalculationInit()
{
  /*
  double width_near = 10.0*CLHEP::meter ; // width of near pool                          
  double width_far  = 16.0*CLHEP::meter ; // width of far pool                           
  double length     = 16.0*CLHEP::meter ; // length of near and far pool                 
  double depth      = 10.0*CLHEP::meter ; // depth of near and far pool 
  double PoolDeadThickness = 84*CLHEP::mm;

  // OWS box is constructed.
  if(m_site == Site::kDayaBay || m_site == Site::kLingAo) {
    
    m_waterPool = new RectangularBox(length/2, -length/2,
                                     width_near/2, -width_near/2,
                                     depth/2, -depth/2+PoolDeadThickness);
  } else if ( m_site == Site::kFar ) {
    m_waterPool = new RectangularBox(length/2, -length/2,
                                     width_far/2, -width_far/2,
                                     depth/2, -depth/2+PoolDeadThickness);
  }
  
  info()<<"Water pool geometry "<<*m_waterPool<<endreq;
  info()<<"TrkLengthInWaterThres= "<<m_trkLengthInWaterThres/CLHEP::cm<<"cm"<<endreq;


  double height_RPC = 0.50*CLHEP::meter ; // RPCs are nominally 50cm above  surface of shield volume (includes 20cm of air above water)
  double thick_RPC  = (0.056+0.030)*CLHEP::meter ; // 1.4cm/layer*4layers + 1cm spacing*3spaces between layers
  double air_thick  =  0.2*CLHEP::meter ; // 20cm of air above water
  double extend_RPC = 1*CLHEP::meter;
  if(m_site == Site::kDayaBay || m_site == Site::kLingAo) {

    m_rpcPlane = new RpcPlane(length/2+extend_RPC, -length/2-extend_RPC,
                              width_near/2+extend_RPC, -width_near/2-extend_RPC,
                              depth/2+air_thick+height_RPC+thick_RPC/2);
  } else if ( m_site == Site::kFar ) {

    m_rpcPlane = new RpcPlane(length/2+extend_RPC, -length/2-extend_RPC,
                              width_far/2+extend_RPC, -width_far/2-extend_RPC,
                              depth/2+air_thick+height_RPC+thick_RPC/2);
  }
  info()<<"RPC geometry "<<*m_rpcPlane<<endreq;
  */

  // 1. DetectorDataSvc
  IDataProviderSvc* detSvc = 0;
  StatusCode sc = service("DetectorDataSvc",detSvc,true);
  if (sc.isFailure()) return sc;
 
  // 2. Top Detector Element
  debug()<<"DetectorElement input: "<<m_topDetectorElementName<<endreq;
  SmartDataPtr<IDetectorElement> topDE(detSvc,m_topDetectorElementName);
  if (!topDE) return StatusCode::FAILURE;
  debug()<<"DetectorElement output: "<<topDE->name()<<endreq;
 
  // 3. GeometryInfo
  m_topGeoInfo = topDE->geometry();
  if(!m_topGeoInfo) return StatusCode::FAILURE;
  debug()<<"Top Geometry "<<m_topGeoInfo->lvolumeName()<<endreq;

  // 4. get top LVolume
  const ILVolume* cILV = m_topGeoInfo->lvolume();
  m_topLVolume = const_cast<ILVolume*>(cILV);
  if(!m_topLVolume) return StatusCode::FAILURE;
  info()<<"Got top logic volume: "<<m_topLVolume->name()<<endreq;

  // 5. get the pointer to all charactoristic materials
  m_mixGas = 0;    
  m_owsWater = 0;  
  m_iwsWater = 0;  
  m_mineralOil=0;  
  m_rock = 0;      

  DataObject* pObj;

  sc = detSvc->retrieveObject("/dd/Materials/MixGas", pObj);
  if(sc.isFailure()) return sc;
  m_mixGas = dynamic_cast<Material *> (pObj);
  //SmartDataPtr<Material> m_mixGas(detSvc,"/dd/Materials/MixGas");
  if(!m_mixGas) return StatusCode::FAILURE;
  debug()<<"Got charactoristic material for RPC: "<<m_mixGas->name()<<endreq;



  sc = detSvc->retrieveObject("/dd/Materials/OwsWater", pObj);
  if(sc.isFailure()) return sc;
  m_owsWater = dynamic_cast<Material *> (pObj);
  //SmartDataPtr<Material> m_mixGas(detSvc,"/dd/Materials/OwsWater");
  if(!m_owsWater) return StatusCode::FAILURE;
  debug()<<"Got charactoristic material for OWS: "<<m_owsWater->name()<<endreq;

  sc = detSvc->retrieveObject("/dd/Materials/IwsWater", pObj);
  if(sc.isFailure()) return sc;
  m_iwsWater = dynamic_cast<Material *> (pObj);
  //SmartDataPtr<Material> m_mixGas(detSvc,"/dd/Materials/IwsWater");
  if(!m_iwsWater) return StatusCode::FAILURE;
  debug()<<"Got charactoristic material for IWS: "<<m_iwsWater->name()<<endreq;

  // Rock
  sc = detSvc->retrieveObject("/dd/Materials/Rock", pObj);
  if(sc.isFailure()) return sc;
  m_rock = dynamic_cast<Material *> (pObj);
  //SmartDataPtr<Material> m_rock(detSvc,"/dd/Materials/Rock");
  if(!m_rock) return StatusCode::FAILURE;
  debug()<<"Got charactoristic material for rock: "<<m_rock->name()<<endreq;

  // MineralOil
  sc = detSvc->retrieveObject("/dd/Materials/MineralOil", pObj);
  if(sc.isFailure()) return sc;
  m_mineralOil = dynamic_cast<Material *> (pObj);
  //SmartDataPtr<Material> m_rock(detSvc,"/dd/Materials/MineralOil");
  if(!m_mineralOil) return StatusCode::FAILURE;
  debug()<<"Got charactoristic material for AD: "<<m_mineralOil->name()<<endreq;

  //  debug()<<*m_mixGas<<endreq;
  //  debug()<<int(m_mixGas)<<" "<<int(m_owsWater)<<" "<<int(m_rock)<<endreq;

  vector<string> ADNames;
  
  if(m_site == Site::kDayaBay) 
    {
      ADNames.push_back("/dd/Structure/AD/db-oil1");
      ADNames.push_back("/dd/Structure/AD/db-oil2");
    } 
  if(m_site == Site::kLingAo)
    {
      ADNames.push_back("/dd/Structure/AD/la-oil1");
      ADNames.push_back("/dd/Structure/AD/la-oil2");
    }
  if(m_site == Site::kFar)  
    {
      ADNames.push_back("/dd/Structure/AD/far-oil1");
      ADNames.push_back("/dd/Structure/AD/far-oil2");
      ADNames.push_back("/dd/Structure/AD/far-oil3");
      ADNames.push_back("/dd/Structure/AD/far-oil4");
    }

  IGeometryInfo* pGI;
  for(unsigned int idx = 0; idx < ADNames.size(); ++idx)
    {
      // Get Detector Element
      SmartDataPtr<IDetectorElement> pDE(detSvc,ADNames[idx]);
      if (!pDE) return StatusCode::FAILURE;
      debug()<<"Got DetectorElement: "<<pDE->name()<<endreq;
 
      // Then geometryInfo
      pGI = pDE->geometry();
      if(!pGI) return StatusCode::FAILURE;
      debug()<<"Got GeometryInfo: "<<pGI->lvolumeName()<<endreq;
      
      m_ADs.push_back(pGI);
    }

  return StatusCode::SUCCESS;

}


// calculate all needed geometry parameters
StatusCode MuonProphet::geometryCalculation(HepMC::GenParticle * pMuon)
{
  /*
  HepGeom::Point3D<double>  aPoint(pMuon->production_vertex()->position().x(),
                                   pMuon->production_vertex()->position().y(),
                                   pMuon->production_vertex()->position().z());
  HepGeom::Point3D<double> aSlope(pMuon->momentum().px(),
                                  pMuon->momentum().py(),
                                  pMuon->momentum().pz());
 
  aSlope = aSlope.unit();

  m_crossWater = m_waterPool->intersect(aPoint,aSlope,m_crossWaterA,m_crossWaterB);
  
  if(m_crossWater) {
    debug() << "Intersect with water OWS at " << m_crossWaterA <<" and "<< m_crossWaterB <<endreq;
    
    m_trkLengthInWater = m_crossWaterA.distance(m_crossWaterB);

    debug() << "Track length in water is "<<m_trkLengthInWater/CLHEP::cm<<"cm"<<endreq;
  }

  
  m_crossRpc = m_rpcPlane->intersect(aPoint,aSlope,m_rpcIntersect);
  
  if(m_crossRpc) {
    debug()<<"Intersect with RPC at "<<m_rpcIntersect <<endreq;
  }
  */

  //
  // muon track is expressed as x(t)=m_point+m_unit*t

  m_point.SetXYZ(pMuon->production_vertex()->position().x(),
                 pMuon->production_vertex()->position().y(),
                 pMuon->production_vertex()->position().z());
  m_vec.SetXYZ(pMuon->momentum().px(),
               pMuon->momentum().py(),
               pMuon->momentum().pz());
  m_unit = m_vec.Unit();

  //Gaudi::XYZPoint point(0,0,0);
  debug()<<"Point  "<<m_point<<endreq;
  string path = m_topGeoInfo->belongsToPath(m_point,-1);
  debug()<<"Belong to path "<<path<<endreq;
  
  //Gaudi::XYZPoint point(0,0,0);
  //Gaudi::XYZVector vec(1,0,0);

  unsigned int N = m_topLVolume->intersectLine(m_point, 
                                               m_unit, 
                                               m_intersections, 
                                               0,100*CLHEP::meter,
                                               -1);
  
  debug()<<"Found "<<N<<" intersections"<<endreq;
  unsigned int idx;
  for (idx=0; idx<N; ++idx) {
    debug()<<"intersection: "<<m_intersections[idx].first<<" "      
          <<m_intersections[idx].second->name()<<"  "
          <<m_intersections[idx].second->density()<<endreq;
    //<<" "<<int(m_intersections[idx].second)<<endreq;
  }

  // Stop muon. Stopping power = 2 [MeVcm2/g]
  // Begin with an approximation
  double ke = pMuon->momentum().e() - pMuon->momentum().m();
  double sp = 2 * CLHEP::MeV * CLHEP::cm2 / CLHEP::g;
  double trkLength = ke/sp;
  //debug()<< CLHEP::MeV << "  " << CLHEP::cm2 << "  " << CLHEP::g <<endreq;
  debug()<<"ke= "<<ke<<"  sp="<<sp<<"  trkLength="<<trkLength<<endreq;
  
  for (idx=0; idx<N; ++idx) {
    double Lpred = m_intersections[idx].first.second - m_intersections[idx].first.first;
    debug()<<"Length in this volume "<< Lpred*m_intersections[idx].second->density()<<endreq;
    if( Lpred*m_intersections[idx].second->density() < trkLength )  {
      trkLength -= Lpred*m_intersections[idx].second->density();
    } else { // end of track
      m_intersections[idx].first.second = 
        m_intersections[idx].first.first + 
        trkLength / m_intersections[idx].second->density();
      break;
    }
  }
  
  if( idx<N-1 )  {  // erase [ first, last ) 
    m_intersections.erase( m_intersections.begin()+idx+1, m_intersections.end() );
  }
  
  N = m_intersections.size();
  debug()<<"Left "<<N<<" intersections"<<endreq;
  for (idx=0; idx<N; ++idx) {
    debug()<<"intersection: "<<m_intersections[idx].first<<" "
          <<m_intersections[idx].second->name()<<"  "
          <<m_intersections[idx].second->density()<<endreq;
    //<<" "<<int(m_intersections[idx].second)<<endreq;
  }



  for (idx=0; idx<N; ++idx) {
    if(m_intersections[idx].second == m_mixGas) {
      m_crossRpc = true;
      break;
    }
  }
  
  bool first = true;
  m_tickWaterMin = 0;
  m_tickWaterMax = 0;

  for (idx=0; idx<N; ++idx) {
    if(m_intersections[idx].second == m_owsWater ||
       m_intersections[idx].second == m_iwsWater ) {
      // first time find OwsWater
      if(first) {
        m_tickWaterMin = m_intersections[idx].first.first;
        m_tickWaterMax = m_intersections[idx].first.second;
        first = false;
      } else {
        if(m_intersections[idx].first.first < m_tickWaterMin ) {
          m_tickWaterMin = m_intersections[idx].first.first;
        }
        if(m_intersections[idx].first.second > m_tickWaterMax ) {
          m_tickWaterMax = m_intersections[idx].first.second;
        }
      }
    }
  }

   
  m_trkLengthInWater = m_tickWaterMax - m_tickWaterMin;
  debug()<<"Track length in water "<< m_trkLengthInWater/CLHEP::m <<"m"<<endreq;

  if( m_trkLengthInWater>0 ) m_crossWater = true;

  m_crucialSegments.clear();
  
  for (idx=0; idx<N; ++idx) {
    if(m_intersections[idx].second == m_rock ||
       m_intersections[idx].second == m_owsWater ||
       m_intersections[idx].second == m_iwsWater ||
       m_intersections[idx].second == m_mineralOil ) 
      {
        if(m_crucialSegments.empty())    
          {
            ILVolume::Intersection * seg = new ILVolume::Intersection(m_intersections[idx]);
            m_crucialSegments.push_back(*seg);
          } else {   
            // check whether last one has the same material as current volume
            // iwsWater and owsWater are considered as the same
            if((m_crucialSegments.back().second == m_intersections[idx].second) ||
               (m_intersections[idx].second == m_iwsWater && m_crucialSegments.back().second == m_owsWater) ||
               (m_intersections[idx].second == m_owsWater && m_crucialSegments.back().second == m_iwsWater) )
              { // If yes, combine them. Assume LS and oil are same
                m_crucialSegments.back().first.second = m_intersections[idx].first.second;
              } else {
                ILVolume::Intersection * seg = new ILVolume::Intersection(m_intersections[idx]);
                m_crucialSegments.push_back(*seg);
              }
          }
      }
  }

  unsigned int NJiaodian = m_crucialSegments.size();
  for (idx=0; idx<NJiaodian; ++idx) {
    debug()<<"intersection: "<<m_crucialSegments[idx].first<<" "
           <<m_crucialSegments[idx].second->name()<<endreq;
    //           <<" "<<int(m_crucialSegments[idx].second)<<endreq;
  }

  return StatusCode::SUCCESS;
  
}