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

In This Package:

MuonProphet.cc
Go to the documentation of this file.
00001 
00016 #include "MuonProphet.h"
00017 // 
00018 MuonProphet::MuonProphet(const std::string& type,
00019                          const std::string& name,
00020                          const IInterface* parent)
00021   : GaudiTool(type,name,parent),
00022     m_muon(0),
00023     m_topGeoInfo(0),
00024     m_topLVolume(0)
00025 {
00026   declareInterface<IHepMCEventMutator>(this);
00027   //
00028   declareProperty("Active",            m_active=true, "Use false to turn off this module");
00029   //
00030   declareProperty("GenTools",          m_genToolNames,"Tools to generate HepMC::GenEvents");
00031   declareProperty("GenYields",         m_genYields,   "Yield for each generator");
00032   declareProperty("GenYieldMeasuredAt",m_genYieldMeasuredAt,"The energy at which the yield is measured");
00033   declareProperty("GenLifetimes",      m_genLifetimes,"Lifetime of each generator");
00034   //
00035   declareProperty("ActNeutron",        m_actNeutron=true,"Turn on or off neutron production");
00036   declareProperty("NeutronYields",     m_neutronYield=-1,"Neutron yield. See MuonProphet.h for description");
00037   //
00038   declareProperty("Site",              m_siteName="DayaBay", "Site");
00039 
00041   declareProperty("TrkLengthInWaterThres", m_trkLengthInWaterThres, 
00042                   "If a track length in water is longer than this, it has high trigger rate.");
00043 
00044   declareProperty("WaterPoolTriggerEff",   m_waterPoolTriggerEff=-1,
00045                   "Trigger efficiency in high trigger rate region, not in use!");
00046 
00048   declareProperty("TopDetectorElementName",
00049                   m_topDetectorElementName="/dd/Structure/DayaBay",
00050                   "Name of the DetectorElement that holds all others");
00051 }
00052 
00053 
00054 MuonProphet::~MuonProphet()
00055 {}
00056 
00057   
00058 StatusCode MuonProphet::initialize()
00059 {
00060   info() << "Initializing MuonProphet" << endreq;
00061   
00062   size_t s=m_genToolNames.size();
00063 
00064   if(s != m_genYields.size() ||
00065      s != m_genLifetimes.size() ||
00066      s != m_genYieldMeasuredAt.size() )  {
00067     error()<<"Imcomplete parameter set"<<endreq;
00068     return StatusCode::FAILURE;
00069   }
00070 
00071   // print all user configuration
00072   info()<<"MuonProphet is active: "<<m_active<<endreq;
00073   m_site=Site::FromString(m_siteName.c_str());
00074   info()<<"MuonProphet is set to site: "<<m_siteName<<" "<<m_site<<endreq;
00075   info()<<"Name  ||  Yield  ||  Yield measured at  ||  Lifetime"<<endreq;
00076   for (size_t idx=0; idx<s; ++idx) {
00077     info() << m_genToolNames[idx]<< " "
00078            << m_genYields[idx]/(CLHEP::cm2/CLHEP::gram) <<"cm2/g "
00079            << m_genYieldMeasuredAt[idx]/CLHEP::GeV <<"GeV "
00080            << m_genLifetimes[idx]/CLHEP::second<<"s" <<endreq;
00081   }
00082 
00083   // get all spallation background gentool
00084   string gtn;
00085   for (size_t idx=0; idx<s; ++idx) {
00086     gtn = m_genToolNames[idx];
00087     try {
00088       m_genTools.push_back(tool<IHepMCEventMutator>(gtn));
00089     }
00090     catch(const GaudiException& exg) {
00091       fatal() << "Failed to get generator: \"" << gtn << "\"" << endreq;
00092       return StatusCode::FAILURE;
00093     }
00094     info () << "Added spallation background gentool " << gtn << endreq;
00095   }
00096 
00099   if(geometryCalculationInit().isFailure()) return StatusCode::FAILURE;
00101   // Spallation neutron related distribution initialization
00102   if(spallationNeutronsInit().isFailure()) return StatusCode::FAILURE;
00103 
00104 
00105   return StatusCode::SUCCESS;
00106 }
00107 
00108 
00109 StatusCode MuonProphet::finalize()
00110 {
00111   return StatusCode::SUCCESS;
00112 }
00113   
00114 // HepMCEventMutator interface
00115 StatusCode MuonProphet::mutate(HepMC::GenEvent& event)
00116 {
00117   // Not active. Don't bother to run it.
00118   if(!m_active) return StatusCode::SUCCESS;
00119   debug()<<"Processing a GenEvent"<<endreq;
00120 
00121   printout(event);
00122 
00123   m_crossRpc = false;
00124   m_crossWater = false;
00125 
00126   // Predict the fate of muon and change it.
00127   MpMuonFate fate = predictMuonFate(event);
00128   debug()<<fate<<endreq;
00129 
00130   if (fate.getCode() == MpMuonFate::kNeedSim) {
00131     // do nothing, return
00132     // to do:
00133     // But if geant4 can't produce these background, I still need to do it here.
00134     // For example, Li9 and He8.
00135     // Then it will be consistent within whole simulation scheme.
00136     // Right now no kNeedSim decision has been made.
00137   } else if (fate.getCode() == MpMuonFate::kUnknown) {
00138     // do nothing, return. Let geant4 to handle it.
00139   } else { 
00143     int ntrk=0;
00144     size_t s=m_genTools.size();
00145 
00146     for (size_t idx=0; idx<s; ++idx) {
00147       
00148       ntrk=spallationProducts(event, idx);
00149       
00150       if(ntrk<0) {
00151         error()<<"failed to run one spallation background production"<<endreq;
00152         return StatusCode::FAILURE;
00153       }
00154     }
00155     
00156     // spallation neutron
00157     // Since it will be triggered, won't be bothered with gamma, 
00158     // because it can only cause prompt background and it is buried in muon signal.
00159     if( m_actNeutron ) {
00160       ntrk=spallationNeutrons(event);    
00161       if(ntrk<0) {
00162         error()<<"failed to run spallation neutron background production"<<endreq;
00163         return StatusCode::FAILURE;
00164       }
00165     }
00166   }
00167 
00170   HepMC::GenParticle* particle = new HepMC::GenParticle(HepMC::FourVector(0.01*CLHEP::eV,
00171                                                                           0,
00172                                                                           0,
00173                                                                           0.01*CLHEP::eV),
00174                                                         22,
00175                                                         1/*=status*/);
00177   (*(event.vertices_begin()))->add_particle_out(particle);
00178   
00180   printout(event);
00181 
00182   return StatusCode::SUCCESS;
00183 }
00184   
00185 
00186 // Predict the fate of a muon, alter the status of the muon track
00187 // This will affect all the following simulation behaviour, FAST or FULL?
00188 MpMuonFate MuonProphet::predictMuonFate(HepMC::GenEvent& event)
00189 {
00190   
00191   HepMC::GenEvent::particle_iterator p;
00192   HepMC::GenEvent::particle_iterator particleEnd=event.particles_end();
00193   
00194   int nMuon=0;
00195   
00196   m_muon=0;
00199   for ( p = event.particles_begin(); p!=particleEnd; ++p) {
00200 
00201     // muon pdg id is 13 or -13
00202     if( (*p)->pdg_id() == 13 || (*p)->pdg_id() == -13 ) {
00203       // Get the Muon pointer
00204       m_muon=*p;
00205       // In case of two or more than two muon tracks, don't work.
00206       ++nMuon;      
00207       if(nMuon>=2) {
00208         return MpMuonFate::kUnknown;
00209       }
00210     }
00211   }
00212 
00213   if(!m_muon) {
00214     info()<<"No muon track found"<<endreq;
00215     return MpMuonFate::kUnknown;
00216   }
00217 
00218   // calculate all needed geometry parameters
00219   geometryCalculation(m_muon);
00220 
00221   // get the trigger status of RPC and waterpool.  
00222   //   Note: No MpTriggerStat::kNeedSim will be issued.
00223   m_rpcTriggerStatus = triggeredByRpc(m_muon);
00224   //   Note: No MpTriggerStat::kNeedSim will be issued.
00225   m_poolTriggerStatus = triggeredByWaterPool(m_muon);
00226 
00227   debug()<<"RPC trigger status:  "<<m_rpcTriggerStatus<<endreq;
00228   debug()<<"Pool trigger status: "<<m_poolTriggerStatus<<endreq;
00229 
00230   
00231   MpMuonFate fate(m_rpcTriggerStatus.getCode(), m_poolTriggerStatus.getCode());
00232   m_muon->set_status(fate.getCode());
00233 
00234   return fate;
00235 
00236   return MpMuonFate::kUnknown;
00237 }
00238 
00239 
00240 StatusCode MuonProphet::printout(HepMC::GenEvent& event)
00241 {
00242   // every vertex
00243   HepMC::GenEvent::vertex_iterator vtx, done = event.vertices_end();
00244   for (vtx = event.vertices_begin(); vtx != done; ++vtx) {
00245     debug()<<*(*vtx)<<endreq;
00246     // every particle
00247     HepMC::GenVertex::particles_in_const_iterator p;
00248     HepMC::GenVertex::particles_in_const_iterator inEnd=(*vtx)->particles_in_const_end();
00249     for ( p = (*vtx)->particles_in_const_begin(); p!=inEnd; ++p) {
00250       debug()<<*(*p)<<endreq;
00251     }
00252     HepMC::GenVertex::particles_out_const_iterator outEnd=(*vtx)->particles_out_const_end();
00253     for ( p = (*vtx)->particles_out_const_begin(); p!=outEnd; ++p) {
00254       debug()<<*(*p)<<endreq;
00255     }
00256   }
00257   return StatusCode::SUCCESS;
00258 }
00259 
| 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