/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

MpGeometry.cc
Go to the documentation of this file.
00001 
00012 #include "MuonProphet.h"
00013  
00014 #include "GaudiKernel/Point3DTypes.h"
00015 #include "GaudiKernel/Vector3DTypes.h"
00016 
00017 StatusCode MuonProphet::geometryCalculationInit()
00018 {
00019   /*
00021   double width_near = 10.0*CLHEP::meter ; // width of near pool                          
00022   double width_far  = 16.0*CLHEP::meter ; // width of far pool                           
00023   double length     = 16.0*CLHEP::meter ; // length of near and far pool                 
00024   double depth      = 10.0*CLHEP::meter ; // depth of near and far pool 
00025   double PoolDeadThickness = 84*CLHEP::mm;
00026 
00027   // OWS box is constructed.
00028   if(m_site == Site::kDayaBay || m_site == Site::kLingAo) {
00029     
00030     m_waterPool = new RectangularBox(length/2, -length/2,
00031                                      width_near/2, -width_near/2,
00032                                      depth/2, -depth/2+PoolDeadThickness);
00033   } else if ( m_site == Site::kFar ) {
00034     m_waterPool = new RectangularBox(length/2, -length/2,
00035                                      width_far/2, -width_far/2,
00036                                      depth/2, -depth/2+PoolDeadThickness);
00037   }
00038   
00039   info()<<"Water pool geometry "<<*m_waterPool<<endreq;
00040   info()<<"TrkLengthInWaterThres= "<<m_trkLengthInWaterThres/CLHEP::cm<<"cm"<<endreq;
00041 
00042 
00044   double height_RPC = 0.50*CLHEP::meter ; // RPCs are nominally 50cm above  surface of shield volume (includes 20cm of air above water)
00045   double thick_RPC  = (0.056+0.030)*CLHEP::meter ; // 1.4cm/layer*4layers + 1cm spacing*3spaces between layers
00046   double air_thick  =  0.2*CLHEP::meter ; // 20cm of air above water
00047   double extend_RPC = 1*CLHEP::meter;
00048   if(m_site == Site::kDayaBay || m_site == Site::kLingAo) {
00049 
00050     m_rpcPlane = new RpcPlane(length/2+extend_RPC, -length/2-extend_RPC,
00051                               width_near/2+extend_RPC, -width_near/2-extend_RPC,
00052                               depth/2+air_thick+height_RPC+thick_RPC/2);
00053   } else if ( m_site == Site::kFar ) {
00054 
00055     m_rpcPlane = new RpcPlane(length/2+extend_RPC, -length/2-extend_RPC,
00056                               width_far/2+extend_RPC, -width_far/2-extend_RPC,
00057                               depth/2+air_thick+height_RPC+thick_RPC/2);
00058   }
00059   info()<<"RPC geometry "<<*m_rpcPlane<<endreq;
00060   */
00061 
00063   // 1. DetectorDataSvc
00064   IDataProviderSvc* detSvc = 0;
00065   StatusCode sc = service("DetectorDataSvc",detSvc,true);
00066   if (sc.isFailure()) return sc;
00067  
00068   // 2. Top Detector Element
00069   debug()<<"DetectorElement input: "<<m_topDetectorElementName<<endreq;
00070   SmartDataPtr<IDetectorElement> topDE(detSvc,m_topDetectorElementName);
00071   if (!topDE) return StatusCode::FAILURE;
00072   debug()<<"DetectorElement output: "<<topDE->name()<<endreq;
00073  
00074   // 3. GeometryInfo
00075   m_topGeoInfo = topDE->geometry();
00076   if(!m_topGeoInfo) return StatusCode::FAILURE;
00077   debug()<<"Top Geometry "<<m_topGeoInfo->lvolumeName()<<endreq;
00078 
00079   // 4. get top LVolume
00080   const ILVolume* cILV = m_topGeoInfo->lvolume();
00081   m_topLVolume = const_cast<ILVolume*>(cILV);
00082   if(!m_topLVolume) return StatusCode::FAILURE;
00083   info()<<"Got top logic volume: "<<m_topLVolume->name()<<endreq;
00084 
00085   // 5. get the pointer to all charactoristic materials
00086   m_mixGas = 0;    
00087   m_owsWater = 0;  
00088   m_iwsWater = 0;  
00089   m_mineralOil=0;  
00090   m_rock = 0;      
00091 
00092   DataObject* pObj;
00093 
00097   sc = detSvc->retrieveObject("/dd/Materials/MixGas", pObj);
00098   if(sc.isFailure()) return sc;
00099   m_mixGas = dynamic_cast<Material *> (pObj);
00100   //SmartDataPtr<Material> m_mixGas(detSvc,"/dd/Materials/MixGas");
00101   if(!m_mixGas) return StatusCode::FAILURE;
00102   debug()<<"Got charactoristic material for RPC: "<<m_mixGas->name()<<endreq;
00103 
00104 
00105 
00107   sc = detSvc->retrieveObject("/dd/Materials/OwsWater", pObj);
00108   if(sc.isFailure()) return sc;
00109   m_owsWater = dynamic_cast<Material *> (pObj);
00110   //SmartDataPtr<Material> m_mixGas(detSvc,"/dd/Materials/OwsWater");
00111   if(!m_owsWater) return StatusCode::FAILURE;
00112   debug()<<"Got charactoristic material for OWS: "<<m_owsWater->name()<<endreq;
00113 
00115   sc = detSvc->retrieveObject("/dd/Materials/IwsWater", pObj);
00116   if(sc.isFailure()) return sc;
00117   m_iwsWater = dynamic_cast<Material *> (pObj);
00118   //SmartDataPtr<Material> m_mixGas(detSvc,"/dd/Materials/IwsWater");
00119   if(!m_iwsWater) return StatusCode::FAILURE;
00120   debug()<<"Got charactoristic material for IWS: "<<m_iwsWater->name()<<endreq;
00121 
00122   // Rock
00123   sc = detSvc->retrieveObject("/dd/Materials/Rock", pObj);
00124   if(sc.isFailure()) return sc;
00125   m_rock = dynamic_cast<Material *> (pObj);
00126   //SmartDataPtr<Material> m_rock(detSvc,"/dd/Materials/Rock");
00127   if(!m_rock) return StatusCode::FAILURE;
00128   debug()<<"Got charactoristic material for rock: "<<m_rock->name()<<endreq;
00129 
00130   // MineralOil
00131   sc = detSvc->retrieveObject("/dd/Materials/MineralOil", pObj);
00132   if(sc.isFailure()) return sc;
00133   m_mineralOil = dynamic_cast<Material *> (pObj);
00134   //SmartDataPtr<Material> m_rock(detSvc,"/dd/Materials/MineralOil");
00135   if(!m_mineralOil) return StatusCode::FAILURE;
00136   debug()<<"Got charactoristic material for AD: "<<m_mineralOil->name()<<endreq;
00137 
00138   //  debug()<<*m_mixGas<<endreq;
00139   //  debug()<<int(m_mixGas)<<" "<<int(m_owsWater)<<" "<<int(m_rock)<<endreq;
00140 
00142   vector<string> ADNames;
00143   
00144   if(m_site == Site::kDayaBay) 
00145     {
00146       ADNames.push_back("/dd/Structure/AD/db-oil1");
00147       ADNames.push_back("/dd/Structure/AD/db-oil2");
00148     } 
00149   if(m_site == Site::kLingAo)
00150     {
00151       ADNames.push_back("/dd/Structure/AD/la-oil1");
00152       ADNames.push_back("/dd/Structure/AD/la-oil2");
00153     }
00154   if(m_site == Site::kFar)  
00155     {
00156       ADNames.push_back("/dd/Structure/AD/far-oil1");
00157       ADNames.push_back("/dd/Structure/AD/far-oil2");
00158       ADNames.push_back("/dd/Structure/AD/far-oil3");
00159       ADNames.push_back("/dd/Structure/AD/far-oil4");
00160     }
00161 
00162   IGeometryInfo* pGI;
00163   for(unsigned int idx = 0; idx < ADNames.size(); ++idx)
00164     {
00165       // Get Detector Element
00166       SmartDataPtr<IDetectorElement> pDE(detSvc,ADNames[idx]);
00167       if (!pDE) return StatusCode::FAILURE;
00168       debug()<<"Got DetectorElement: "<<pDE->name()<<endreq;
00169  
00170       // Then geometryInfo
00171       pGI = pDE->geometry();
00172       if(!pGI) return StatusCode::FAILURE;
00173       debug()<<"Got GeometryInfo: "<<pGI->lvolumeName()<<endreq;
00174       
00175       m_ADs.push_back(pGI);
00176     }
00177 
00178   return StatusCode::SUCCESS;
00179 
00180 }
00181 
00182 
00183 // calculate all needed geometry parameters
00184 StatusCode MuonProphet::geometryCalculation(HepMC::GenParticle * pMuon)
00185 {
00187   /*
00188   HepGeom::Point3D<double>  aPoint(pMuon->production_vertex()->position().x(),
00189                                    pMuon->production_vertex()->position().y(),
00190                                    pMuon->production_vertex()->position().z());
00191   HepGeom::Point3D<double> aSlope(pMuon->momentum().px(),
00192                                   pMuon->momentum().py(),
00193                                   pMuon->momentum().pz());
00194  
00195   aSlope = aSlope.unit();
00196 
00197   m_crossWater = m_waterPool->intersect(aPoint,aSlope,m_crossWaterA,m_crossWaterB);
00198   
00199   if(m_crossWater) {
00200     debug() << "Intersect with water OWS at " << m_crossWaterA <<" and "<< m_crossWaterB <<endreq;
00201     
00202     m_trkLengthInWater = m_crossWaterA.distance(m_crossWaterB);
00203 
00204     debug() << "Track length in water is "<<m_trkLengthInWater/CLHEP::cm<<"cm"<<endreq;
00205   }
00206 
00207   
00208   m_crossRpc = m_rpcPlane->intersect(aPoint,aSlope,m_rpcIntersect);
00209   
00210   if(m_crossRpc) {
00211     debug()<<"Intersect with RPC at "<<m_rpcIntersect <<endreq;
00212   }
00213   */
00214 
00216   //
00217   // muon track is expressed as x(t)=m_point+m_unit*t
00218 
00219   m_point.SetXYZ(pMuon->production_vertex()->position().x(),
00220                  pMuon->production_vertex()->position().y(),
00221                  pMuon->production_vertex()->position().z());
00222   m_vec.SetXYZ(pMuon->momentum().px(),
00223                pMuon->momentum().py(),
00224                pMuon->momentum().pz());
00225   m_unit = m_vec.Unit();
00226 
00228   //Gaudi::XYZPoint point(0,0,0);
00229   debug()<<"Point  "<<m_point<<endreq;
00230   string path = m_topGeoInfo->belongsToPath(m_point,-1);
00231   debug()<<"Belong to path "<<path<<endreq;
00232   
00233   //Gaudi::XYZPoint point(0,0,0);
00234   //Gaudi::XYZVector vec(1,0,0);
00235 
00236   unsigned int N = m_topLVolume->intersectLine(m_point, 
00237                                                m_unit, 
00238                                                m_intersections, 
00239                                                0,100*CLHEP::meter,
00240                                                -1);
00241   
00242   debug()<<"Found "<<N<<" intersections"<<endreq;
00243   unsigned int idx;
00244   for (idx=0; idx<N; ++idx) {
00245     debug()<<"intersection: "<<m_intersections[idx].first<<" "      
00246           <<m_intersections[idx].second->name()<<"  "
00247           <<m_intersections[idx].second->density()<<endreq;
00248     //<<" "<<int(m_intersections[idx].second)<<endreq;
00249   }
00250 
00251   // Stop muon. Stopping power = 2 [MeVcm2/g]
00252   // Begin with an approximation
00253   double ke = pMuon->momentum().e() - pMuon->momentum().m();
00254   double sp = 2 * CLHEP::MeV * CLHEP::cm2 / CLHEP::g;
00255   double trkLength = ke/sp;
00256   //debug()<< CLHEP::MeV << "  " << CLHEP::cm2 << "  " << CLHEP::g <<endreq;
00257   debug()<<"ke= "<<ke<<"  sp="<<sp<<"  trkLength="<<trkLength<<endreq;
00258   
00259   for (idx=0; idx<N; ++idx) {
00260     double Lpred = m_intersections[idx].first.second - m_intersections[idx].first.first;
00261     debug()<<"Length in this volume "<< Lpred*m_intersections[idx].second->density()<<endreq;
00262     if( Lpred*m_intersections[idx].second->density() < trkLength )  {
00263       trkLength -= Lpred*m_intersections[idx].second->density();
00264     } else { // end of track
00265       m_intersections[idx].first.second = 
00266         m_intersections[idx].first.first + 
00267         trkLength / m_intersections[idx].second->density();
00268       break;
00269     }
00270   }
00271   
00272   if( idx<N-1 )  {  // erase [ first, last ) 
00273     m_intersections.erase( m_intersections.begin()+idx+1, m_intersections.end() );
00274   }
00275   
00276   N = m_intersections.size();
00277   debug()<<"Left "<<N<<" intersections"<<endreq;
00278   for (idx=0; idx<N; ++idx) {
00279     debug()<<"intersection: "<<m_intersections[idx].first<<" "
00280           <<m_intersections[idx].second->name()<<"  "
00281           <<m_intersections[idx].second->density()<<endreq;
00282     //<<" "<<int(m_intersections[idx].second)<<endreq;
00283   }
00284 
00285 
00286 
00288   for (idx=0; idx<N; ++idx) {
00289     if(m_intersections[idx].second == m_mixGas) {
00290       m_crossRpc = true;
00291       break;
00292     }
00293   }
00294   
00296   bool first = true;
00297   m_tickWaterMin = 0;
00298   m_tickWaterMax = 0;
00299 
00300   for (idx=0; idx<N; ++idx) {
00301     if(m_intersections[idx].second == m_owsWater ||
00302        m_intersections[idx].second == m_iwsWater ) {
00303       // first time find OwsWater
00304       if(first) {
00305         m_tickWaterMin = m_intersections[idx].first.first;
00306         m_tickWaterMax = m_intersections[idx].first.second;
00307         first = false;
00308       } else {
00309         if(m_intersections[idx].first.first < m_tickWaterMin ) {
00310           m_tickWaterMin = m_intersections[idx].first.first;
00311         }
00312         if(m_intersections[idx].first.second > m_tickWaterMax ) {
00313           m_tickWaterMax = m_intersections[idx].first.second;
00314         }
00315       }
00316     }
00317   }
00318 
00319    
00320   m_trkLengthInWater = m_tickWaterMax - m_tickWaterMin;
00321   debug()<<"Track length in water "<< m_trkLengthInWater/CLHEP::m <<"m"<<endreq;
00322 
00324   if( m_trkLengthInWater>0 ) m_crossWater = true;
00326 
00330   m_crucialSegments.clear();
00331   
00332   for (idx=0; idx<N; ++idx) {
00333     if(m_intersections[idx].second == m_rock ||
00334        m_intersections[idx].second == m_owsWater ||
00335        m_intersections[idx].second == m_iwsWater ||
00336        m_intersections[idx].second == m_mineralOil ) 
00337       {
00338         if(m_crucialSegments.empty())    
00339           {
00340             ILVolume::Intersection * seg = new ILVolume::Intersection(m_intersections[idx]);
00341             m_crucialSegments.push_back(*seg);
00342           } else {   
00343             // check whether last one has the same material as current volume
00344             // iwsWater and owsWater are considered as the same
00345             if((m_crucialSegments.back().second == m_intersections[idx].second) ||
00346                (m_intersections[idx].second == m_iwsWater && m_crucialSegments.back().second == m_owsWater) ||
00347                (m_intersections[idx].second == m_owsWater && m_crucialSegments.back().second == m_iwsWater) )
00348               { // If yes, combine them. Assume LS and oil are same
00349                 m_crucialSegments.back().first.second = m_intersections[idx].first.second;
00350               } else {
00351                 ILVolume::Intersection * seg = new ILVolume::Intersection(m_intersections[idx]);
00352                 m_crucialSegments.push_back(*seg);
00353               }
00354           }
00355       }
00356   }
00357 
00358   unsigned int NJiaodian = m_crucialSegments.size();
00359   for (idx=0; idx<NJiaodian; ++idx) {
00360     debug()<<"intersection: "<<m_crucialSegments[idx].first<<" "
00361            <<m_crucialSegments[idx].second->name()<<endreq;
00362     //           <<" "<<int(m_crucialSegments[idx].second)<<endreq;
00363   }
00364 
00365   return StatusCode::SUCCESS;
00366   
00367 }
00368 
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:21:29 for MuonProphet by doxygen 1.7.4