CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/SimG4CMS/Forward/src/CastorShowerLibrary.cc

Go to the documentation of this file.
00001 
00002 // File: CastorShowerLibrary.cc
00003 // Description: Shower library for CASTOR calorimeter
00004 //              Adapted from HFShowerLibrary
00005 //
00006 //  Wagner Carvalho
00007 //
00009 
00010 #include "SimG4CMS/Forward/interface/CastorShowerLibrary.h"
00011 #include "FWCore/Utilities/interface/Exception.h"
00012 
00013 #include "G4VPhysicalVolume.hh"
00014 #include "G4Step.hh"
00015 #include "G4Track.hh"
00016 #include "G4ParticleTable.hh"
00017 #include "Randomize.hh"
00018 #include "CLHEP/Units/GlobalSystemOfUnits.h"
00019 
00020 //#define DebugLog
00021 
00022 CastorShowerLibrary::CastorShowerLibrary(std::string & name, edm::ParameterSet const & p) 
00023                                           : hf(0), evtInfo(0), emBranch(0), hadBranch(0),
00024                                             nMomBin(0), totEvents(0), evtPerBin(0),
00025                                             nBinsE(0),nBinsEta(0),nBinsPhi(0),
00026                                             nEvtPerBinE(0),nEvtPerBinEta(0),nEvtPerBinPhi(0),
00027                                             SLenergies(),SLetas(),SLphis() {
00028   
00029   initFile(p);
00030   
00031 }
00032 
00033 //=============================================================================================
00034 
00035 
00036 CastorShowerLibrary::~CastorShowerLibrary() {
00037   if (hf)     hf->Close();
00038 }
00039 
00040 
00041 //=============================================================================================
00042 
00043 void CastorShowerLibrary::initFile(edm::ParameterSet const & p) {
00045 //
00046 //  Init TFile and associated TBranch's of CASTOR Root file 
00047 //  holding library events 
00048 //
00050 
00051   //
00052   //  Read PSet for Castor shower library
00053   //
00054 
00055   edm::ParameterSet m_CS   = p.getParameter<edm::ParameterSet>("CastorShowerLibrary");
00056   edm::FileInPath fp       = m_CS.getParameter<edm::FileInPath>("FileName");
00057   std::string pTreeName    = fp.fullPath();
00058   std::string branchEvInfo = m_CS.getUntrackedParameter<std::string>("BranchEvt");
00059   std::string branchEM     = m_CS.getUntrackedParameter<std::string>("BranchEM");
00060   std::string branchHAD    = m_CS.getUntrackedParameter<std::string>("BranchHAD");
00061   verbose                  = m_CS.getUntrackedParameter<bool>("Verbosity",false);
00062 
00063   // Open TFile 
00064   if (pTreeName.find(".") == 0) pTreeName.erase(0,2);
00065   const char* nTree = pTreeName.c_str();
00066   hf                = TFile::Open(nTree);
00067 
00068   // Check that TFile has been successfully opened
00069   if (!hf->IsOpen()) { 
00070      edm::LogError("CastorShower") << "CastorShowerLibrary: opening " << nTree << " failed";
00071      throw cms::Exception("Unknown", "CastorShowerLibrary") 
00072                                    << "Opening of " << pTreeName << " fails\n";
00073   } else {
00074      edm::LogInfo("CastorShower")  << "CastorShowerLibrary: opening " << nTree << " successfully"; 
00075   }
00076 
00077   // Check for the TBranch holding EventInfo in "Events" TTree
00078   TTree* event = (TTree *) hf ->Get("CastorCherenkovPhotons");
00079   if (event) {
00080      evtInfo = (TBranchObject *) event->GetBranch(branchEvInfo.c_str());
00081      if (evtInfo) {
00082         loadEventInfo(evtInfo);
00083      } else {
00084         edm::LogError("CastorShower") << "CastorShowerLibrary: " << branchEvInfo.c_str()
00085                                       << " Branch does not exit in Event";
00086         throw cms::Exception("Unknown", "CastorShowerLibrary") << "Event information absent\n";
00087      }
00088   } else {
00089      edm::LogError("CastorShower") << "CastorShowerLibrary: Events Tree does not exist";
00090      throw cms::Exception("Unknown", "CastorShowerLibrary") << "Events tree absent\n";
00091   }
00092   
00093   // Get EM and HAD Branchs
00094   emBranch         = (TBranchObject *) event->GetBranch(branchEM.c_str());
00095   if (verbose) emBranch->Print();
00096   hadBranch        = (TBranchObject *) event->GetBranch(branchHAD.c_str());
00097   if (verbose) hadBranch->Print();
00098   edm::LogInfo("CastorShower") << "CastorShowerLibrary: Branch " << branchEM 
00099                                << " has " << emBranch->GetEntries() 
00100                                << " entries and Branch " << branchHAD 
00101                                << " has " << hadBranch->GetEntries() 
00102                                << " entries";
00103 
00104 }
00105 
00106 //=============================================================================================
00107 
00108 void CastorShowerLibrary::loadEventInfo(TBranchObject* branch) {
00110 //
00111 //  Get EventInfo from the "TBranch* branch" of Root file 
00112 //  holding library events 
00113 //
00114 //  Based on HFShowerLibrary::loadEventInfo
00115 //
00117 
00118   eventInfo = new CastorShowerLibraryInfo();
00119   branch->SetAddress(&eventInfo);
00120   branch->GetEntry(0);
00121   // Initialize shower library general parameters
00122 
00123   totEvents   = eventInfo->Energy.getNEvts();
00124 //  nMomBin     = eventInfo->Energy.getNBins();
00125 //  evtPerBin   = eventInfo->Energy.getNEvtPerBin();
00126 //  pmom        = eventInfo->Energy.getBin();
00127   nBinsE        = eventInfo->Energy.getNBins();
00128   nEvtPerBinE   = eventInfo->Energy.getNEvtPerBin();
00129   SLenergies    = eventInfo->Energy.getBin();
00130   nBinsEta      = eventInfo->Eta.getNBins();
00131   nEvtPerBinEta = eventInfo->Eta.getNEvtPerBin();
00132   SLetas        = eventInfo->Eta.getBin();
00133   nBinsPhi      = eventInfo->Phi.getNBins();
00134   nEvtPerBinPhi = eventInfo->Phi.getNEvtPerBin();
00135   SLphis        = eventInfo->Phi.getBin();
00136   
00137   // Convert from GeV to MeV
00138   for (unsigned int i=0; i<SLenergies.size(); i++) SLenergies[i] *= GeV;
00139   
00140   edm::LogInfo("CastorShower") << " CastorShowerLibrary::loadEventInfo : " 
00141                                << "\n \n Total number of events     :  " << totEvents 
00142                                <<    "\n   Number of bins  (E)       :  " << nBinsE
00143                                <<    "\n   Number of events/bin (E)  :  " << nEvtPerBinE
00144                                <<    "\n   Number of bins  (Eta)       :  " << nBinsEta
00145                                <<    "\n   Number of events/bin (Eta)  :  " << nEvtPerBinEta 
00146                                <<    "\n   Number of bins  (Phi)       :  " << nBinsPhi
00147                                <<    "\n   Number of events/bin (Phi)  :  " << nEvtPerBinPhi << "\n";
00148   
00149   for (unsigned int i=0; i<nBinsE; i++)
00150      edm::LogInfo("CastorShower") << "CastorShowerLibrary: SLenergies[" << i << "] = "
00151                                   << SLenergies[i]/GeV << " GeV";
00152   for (unsigned int i=0; i<nBinsEta; i++)
00153      edm::LogInfo("CastorShower") << "CastorShowerLibrary: SLetas[" << i << "] = "
00154                                   << SLetas[i];
00155   for (unsigned int i=0; i<nBinsPhi; i++)
00156      edm::LogInfo("CastorShower") << "CastorShowerLibrary: SLphis[" << i << "] = "
00157                                   << SLphis[i] << " rad";
00158 }
00159 
00160 //=============================================================================================
00161 
00162 void CastorShowerLibrary::initParticleTable(G4ParticleTable * theParticleTable) {
00164 //
00165 //  Set particle codes according to PDG encoding 
00166 //
00167 //  Based on HFShowerLibrary::initRun
00168 //
00170 
00171   G4String particleName;
00172 
00173   edm::LogInfo("CastorShower") << "CastorShowerLibrary::initParticleTable"
00174                                << " ***  Accessing PDGEncoding  ***" ;
00175 
00176   emPDG       = theParticleTable->FindParticle(particleName="e-")->GetPDGEncoding();
00177   epPDG       = theParticleTable->FindParticle(particleName="e+")->GetPDGEncoding();
00178   gammaPDG    = theParticleTable->FindParticle(particleName="gamma")->GetPDGEncoding();
00179   pi0PDG      = theParticleTable->FindParticle(particleName="pi0")->GetPDGEncoding();
00180   etaPDG      = theParticleTable->FindParticle(particleName="eta")->GetPDGEncoding();
00181   nuePDG      = theParticleTable->FindParticle(particleName="nu_e")->GetPDGEncoding();
00182   numuPDG     = theParticleTable->FindParticle(particleName="nu_mu")->GetPDGEncoding();
00183   nutauPDG    = theParticleTable->FindParticle(particleName="nu_tau")->GetPDGEncoding();
00184   anuePDG     = theParticleTable->FindParticle(particleName="anti_nu_e")->GetPDGEncoding();
00185   anumuPDG    = theParticleTable->FindParticle(particleName="anti_nu_mu")->GetPDGEncoding();
00186   anutauPDG   = theParticleTable->FindParticle(particleName="anti_nu_tau")->GetPDGEncoding();
00187   geantinoPDG = theParticleTable->FindParticle(particleName="geantino")->GetPDGEncoding();
00188   mumPDG      = theParticleTable->FindParticle(particleName="mu-")->GetPDGEncoding();
00189   mupPDG      = theParticleTable->FindParticle(particleName="mu+")->GetPDGEncoding();
00190   
00191 //#ifdef DebugLog
00192   LogDebug("CastorShower") << "CastorShowerLibrary: Particle codes for e- = " << emPDG
00193                            << ", e+ = " << epPDG << ", gamma = " << gammaPDG 
00194                            << ", pi0 = " << pi0PDG << ", eta = " << etaPDG
00195                            << ", geantino = " << geantinoPDG << "\n        nu_e = "
00196                            << nuePDG << ", nu_mu = " << numuPDG << ", nu_tau = "
00197                            << nutauPDG << ", anti_nu_e = " << anuePDG
00198                            << ", anti_nu_mu = " << anumuPDG << ", anti_nu_tau = "
00199                            << anutauPDG << ", mu- = " << mumPDG << ", mu+ = "
00200                            << mupPDG;
00201 //#endif
00202 
00203     edm::LogInfo("CastorShower") << "  *****   Successfully called:  " 
00204                                  << "CastorShowerLibrary::initParticleTable()   ***** " ;
00205 
00206 }
00207 
00208 
00209 //=============================================================================================
00210 
00211 CastorShowerEvent CastorShowerLibrary::getShowerHits(G4Step * aStep, bool & ok) {
00212 
00213   G4StepPoint * preStepPoint  = aStep->GetPreStepPoint(); 
00214   // G4StepPoint * postStepPoint = aStep->GetPostStepPoint(); 
00215   G4Track *     track         = aStep->GetTrack();
00216   // Get Z-direction 
00217   const G4DynamicParticle *aParticle = track->GetDynamicParticle();
00218   G4ThreeVector               momDir = aParticle->GetMomentumDirection();
00219   //  double mom = aParticle->GetTotalMomentum();
00220 
00221   G4ThreeVector hitPoint = preStepPoint->GetPosition();   
00222   G4String      partType = track->GetDefinition()->GetParticleName();
00223   int           parCode  = track->GetDefinition()->GetPDGEncoding();
00224 
00225   CastorShowerEvent hit;
00226   hit.Clear();
00227   
00228   ok = false;
00229   if (parCode == pi0PDG   || parCode == etaPDG    || parCode == nuePDG  ||
00230       parCode == numuPDG  || parCode == nutauPDG  || parCode == anuePDG ||
00231       parCode == anumuPDG || parCode == anutauPDG || parCode == geantinoPDG) 
00232     return hit;
00233   ok = true;
00234 
00235   double pin    = preStepPoint->GetTotalEnergy();
00236   double etain  = momDir.getEta();
00237   double phiin  = momDir.getPhi();
00238   
00239   double zint = hitPoint.z();
00240   double R=sqrt(hitPoint.x()*hitPoint.x() + hitPoint.y()*hitPoint.y());
00241   double theta = atan2(R,std::abs(zint));
00242   phiin = atan2(hitPoint.y(),hitPoint.x());
00243   etain = -1*(std::log(std::tan((M_PI-theta)/2)));
00244 
00245   // Replace "interpolation/extrapolation" by new method "select" that just randomly 
00246   // selects a record from the appropriate energy bin and fills its content to  
00247   // "showerEvent" instance of class CastorShowerEvent
00248   
00249   if (parCode == emPDG || parCode == epPDG || parCode == gammaPDG ) {
00250     select(0, pin, etain, phiin);
00251     // if (pin<SLenergies[nBinsE-1]) {
00252     //   interpolate(0, pin);
00253     // } else {
00254     //   extrapolate(0, pin);
00255     // }
00256   } else {
00257     select(1, pin, etain, phiin);
00258     // if (pin<SLenergies[nBinsE-1]) {
00259     //   interpolate(1, pin);
00260     // } else {
00261     //   extrapolate(1, pin);
00262     // }
00263   }
00264     
00265   hit = (*showerEvent);
00266   return hit;
00267 
00268 }
00269 
00270 //=============================================================================================
00271 
00272 void CastorShowerLibrary::getRecord(int type, int record) {
00274 //
00275 //  Retrieve event # "record" from the library and stores it   
00276 //  into  vector<CastorShowerHit> showerHit;
00277 //
00278 //  Based on HFShowerLibrary::getRecord
00279 //
00280 //  Modified 02/02/09 by W. Carvalho
00281 //
00283 
00284 #ifdef DebugLog
00285   LogDebug("CastorShower") << "CastorShowerLibrary::getRecord: ";
00286 #endif  
00287   int nrc  = record;
00288   showerEvent = new CastorShowerEvent();
00289   if (type > 0) {
00290     hadBranch->SetAddress(&showerEvent);
00291     hadBranch->GetEntry(nrc);
00292   } else {
00293     emBranch->SetAddress(&showerEvent);
00294     emBranch->GetEntry(nrc);
00295   }
00296 
00297 #ifdef DebugLog
00298   int nHit = showerEvent->getNhit();
00299   LogDebug("CastorShower") << "CastorShowerLibrary::getRecord: Record " << record
00300                            << " of type " << type << " with " << nHit 
00301                            << " CastorShowerHits";
00302 
00303 #endif
00304 
00305 }
00306 
00307 //=======================================================================================
00308 
00309 void CastorShowerLibrary::select(int type, double pin, double etain, double phiin) {
00311 //
00312 //  Selects an event from the library based on 
00313 //
00314 //    type:   0 --> em
00315 //           >0 --> had
00316 //    pin  :  momentum
00317 //    etain:  eta (if not given, disregard the eta binning
00318 //    phiin:  phi (if not given, disregard the phi binning
00319 //
00320 //  Created 30/01/09 by W. Carvalho
00321 //
00323 
00324   int irec;                         // to hold record number
00325   double r = G4UniformRand();       // to randomly select within an energy bin (r=[0,1])
00326 
00327   // Randomly select a record from the right energy bin in the library, 
00328   // based on track momentum (pin) 
00329    
00330   //   pin < SLenergies[MIN]
00331 /*
00332   if (pin<SLenergies[0]) {
00333      edm::LogWarning("CastorShower") << "CastorShowerLibrary: Warning, pin = " << pin 
00334                                      << " less than minimum SLenergies " << SLenergies[0] << " in library." 
00335                                      << " For the moment, selecting hit from the first bin" ;
00336      irec = int(nEvtPerBinE*r) + 1;
00337   //   pin > SLenergies[MAX]
00338   } else if (pin>SLenergies[nBinsE]) {
00339 
00340 // This part needs rethinking because the last element of SLenergies is no longer the upper limit of the bins
00341     edm::LogWarning("CastorShower") << "CastorShowerLibrary: Warning, pin = " << pin 
00342                                     << " greater than maximum SLenergies " << SLenergies[nBinsE] << " in library." 
00343                                     << " For the moment, selecting hit from the last bin";
00344     irec = (nBinsE-1)*nEvtPerBinE + int(evtPerBin*r) + 1;
00345   //   SLenergies[MIN] < pin < SLenergies[MAX]
00346   } else {
00347      for (unsigned int j=0; j<nBinsE; j++) {
00348         if (pin >= SLenergies[j] && pin < SLenergies[j+1]) {
00349            irec = j*evtPerBin + int(evtPerBin*r) + 1 ;
00350            if (irec<0) {
00351               edm::LogWarning("CastorShower") << "CastorShowerLibrary:: Illegal irec = "
00352                                               << irec << " now set to 0";
00353               irec = 0;
00354            } else if (irec > totEvents) {
00355               edm::LogWarning("CastorShower") << "CastorShowerLibrary:: Illegal irec = "
00356                                               << irec << " now set to "<< totEvents;
00357               irec = totEvents;
00358            }
00359         }
00360      }
00361   }
00362 */
00363   int ienergy = FindEnergyBin(pin);
00364   int ieta    = FindEtaBin(etain);
00365 #ifdef DebugLog
00366   if (verbose) edm::LogInfo("CastorShower") << " ienergy = " << ienergy ;
00367   if (verbose) edm::LogInfo("CastorShower") << " ieta = " << ieta;
00368 #endif
00369 
00370   int iphi;
00371   double phiMin = 0. ;
00372   double phiMax = M_PI/4 ;     // 45 * (pi/180)  rad 
00373   if(phiin < phiMin) phiin = phiin + M_PI ;
00374   if(phiin >= phiMin && phiin < phiMax) {
00375      iphi = FindPhiBin(phiin) ;
00376   } else {
00377      double remainder = fmod(phiin , M_PI/4) ;
00378      phiin = phiMin + remainder ;
00379      iphi = FindPhiBin(phiin) ;
00380   }
00381 #ifdef DebugLog
00382   if (verbose) edm::LogInfo("CastorShower") << " iphi = " << iphi;
00383 #endif
00384   if (ienergy==-1) ienergy=0;    // if pin < first bin, choose an event in the first one
00385   if (ieta==-1)    ieta=0; // idem for eta
00386   if (iphi!=-1) irec = int(nEvtPerBinE*ienergy+nEvtPerBinEta*ieta+nEvtPerBinPhi*(iphi+r));
00387   else irec = int(nEvtPerBinE*(ienergy+r));
00388 
00389 #ifdef DebugLog
00390   edm::LogInfo("CastorShower") << "CastorShowerLibrary:: Select record " << irec 
00391                                << " of type " << type ; 
00392 #endif
00393 
00394   //  Retrieve record number "irec" from the library
00395   getRecord (type, irec);
00396   
00397 }
00398 int CastorShowerLibrary::FindEnergyBin(double energy) {
00399   //
00400   // returns the integer index of the energy bin, taken from SLenergies vector
00401   // returns -1 if ouside valid range
00402   //
00403   if (energy >= SLenergies.back()) return SLenergies.size()-1;
00404 
00405   unsigned int i = 0;
00406   for(;i<SLenergies.size()-1;i++)
00407     if (energy >= SLenergies.at(i) && energy < SLenergies.at(i+1)) return (int)i;
00408 
00409   // now i points to the last but 1 bin
00410   if (energy>=SLenergies.at(i)) return (int)i;
00411   // energy outside bin range
00412   return -1;
00413 }
00414 int CastorShowerLibrary::FindEtaBin(double eta) {
00415   //
00416   // returns the integer index of the eta bin, taken from SLetas vector
00417   // returns -1 if ouside valid range
00418   //
00419   if (eta>=SLetas.back()) return SLetas.size()-1;
00420   unsigned int i = 0;
00421   for(;i<SLetas.size()-1;i++)
00422      if (eta >= SLetas.at(i) && eta < SLetas.at(i+1)) return (int)i;
00423   // now i points to the last but 1 bin
00424   if (eta>=SLetas.at(i)) return (int)i;
00425   // eta outside bin range
00426   return -1;
00427 }
00428 int CastorShowerLibrary::FindPhiBin(double phi) {
00429   //
00430   // returns the integer index of the phi bin, taken from SLphis vector
00431   // returns -1 if ouside valid range
00432   //
00433   // needs protection in case phi is outside range -pi,pi
00434   //
00435   if (phi>=SLphis.back()) return SLphis.size()-1;
00436   unsigned int i = 0;
00437   for(;i<SLphis.size()-1;i++)
00438      if (phi >= SLphis.at(i) && phi < SLphis.at(i+1)) return (int)i;
00439   // now i points to the last but 1 bin
00440   if (phi>=SLphis.at(i)) return (int)i;
00441   // phi outside bin range
00442   return -1;
00443 }