Forward/src/CastorShowerLibrary.cc

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: CastorShowerLibrary.cc
0003 // Description: Shower library for CASTOR calorimeter
0004 //              Adapted from HFShowerLibrary
0005 //
0006 //  Wagner Carvalho
0007 //
0008 ///////////////////////////////////////////////////////////////////////////////
0009
0010 #include "SimG4CMS/Forward/interface/CastorShowerLibrary.h"
0011 #include "SimG4Core/Notification/interface/G4TrackToParticleID.h"
0012 #include "FWCore/Utilities/interface/Exception.h"
0013
0014 #include "G4VPhysicalVolume.hh"
0015 #include "G4Step.hh"
0016 #include "G4Track.hh"
0017 #include "G4ParticleTable.hh"
0018 #include "Randomize.hh"
0019 #include <CLHEP/Units/PhysicalConstants.h>
0020 #include <CLHEP/Units/SystemOfUnits.h>
0021
0022 //ROOT
0023 #include "TROOT.h"
0024 #include "TFile.h"
0025 #include "TTree.h"
0026 #include "TBranch.h"
0027
0028 //#define EDM_ML_DEBUG
0029
0030 CastorShowerLibrary::CastorShowerLibrary(const std::string& name, edm::ParameterSet const& p)
0031     : hf(nullptr),
0032       emBranch(nullptr),
0033       hadBranch(nullptr),
0034       verbose(false),
0035       nMomBin(0),
0036       totEvents(0),
0037       evtPerBin(0),
0038       nBinsE(0),
0039       nBinsEta(0),
0040       nBinsPhi(0),
0041       nEvtPerBinE(0),
0042       nEvtPerBinEta(0),
0043       nEvtPerBinPhi(0),
0044       SLenergies(),
0045       SLetas(),
0046       SLphis() {
0047   initFile(p);
0048 }
0049
0050 //=============================================================================================
0051
0052 CastorShowerLibrary::~CastorShowerLibrary() {
0053   if (hf)
0054     hf->Close();
0055 }
0056
0057 //=============================================================================================
0058
0059 void CastorShowerLibrary::initFile(edm::ParameterSet const& p) {
0060   //////////////////////////////////////////////////////////
0061   //
0062   //  Init TFile and associated TBranch's of CASTOR Root file
0063   //  holding library events
0064   //
0065   //////////////////////////////////////////////////////////
0066
0067   //
0068   //  Read PSet for Castor shower library
0069   //
0070
0071   edm::ParameterSet m_CS = p.getParameter<edm::ParameterSet>("CastorShowerLibrary");
0072   edm::FileInPath fp = m_CS.getParameter<edm::FileInPath>("FileName");
0073   std::string pTreeName = fp.fullPath();
0074   std::string branchEvInfo = m_CS.getUntrackedParameter<std::string>("BranchEvt");
0075   std::string branchEM = m_CS.getUntrackedParameter<std::string>("BranchEM");
0076   std::string branchHAD = m_CS.getUntrackedParameter<std::string>("BranchHAD");
0077   verbose = m_CS.getUntrackedParameter<bool>("Verbosity", false);
0078
0079   // Open TFile
0080   if (pTreeName.find('.') == 0)
0081     pTreeName.erase(0, 2);
0082   const char* nTree = pTreeName.c_str();
0083   hf = TFile::Open(nTree);
0084
0085   // Check that TFile has been successfully opened
0086   if (!hf->IsOpen()) {
0087     edm::LogError("CastorShower") << "CastorShowerLibrary: opening " << nTree << " failed";
0088     throw cms::Exception("Unknown", "CastorShowerLibrary") << "Opening of " << pTreeName << " fails\n";
0089   } else {
0090     edm::LogVerbatim("CastorShower") << "CastorShowerLibrary: opening " << nTree << " successfully";
0091   }
0092
0093   // Check for the TBranch holding EventVerbatim in "Events" TTree
0094   TTree* event = hf->Get<TTree>("CastorCherenkovPhotons");
0095   if (event) {
0096     auto evtInfo = event->GetBranch(branchEvInfo.c_str());
0097     if (evtInfo) {
0098       loadEventInfo(evtInfo);
0099     } else {
0100       edm::LogError("CastorShower") << "CastorShowerLibrary: " << branchEvInfo.c_str()
0101                                     << " Branch does not exit in Event";
0102       throw cms::Exception("Unknown", "CastorShowerLibrary") << "Event information absent\n";
0103     }
0104   } else {
0105     edm::LogError("CastorShower") << "CastorShowerLibrary: Events Tree does not exist";
0106     throw cms::Exception("Unknown", "CastorShowerLibrary") << "Events tree absent\n";
0107   }
0108
0109   // Get EM and HAD Branchs
0110   emBranch = event->GetBranch(branchEM.c_str());
0111   if (verbose)
0112     emBranch->Print();
0113   hadBranch = event->GetBranch(branchHAD.c_str());
0114   if (verbose)
0115     hadBranch->Print();
0116   edm::LogVerbatim("CastorShower") << "CastorShowerLibrary: Branch " << branchEM << " has " << emBranch->GetEntries()
0117                                    << " entries and Branch " << branchHAD << " has " << hadBranch->GetEntries()
0118                                    << " entries";
0119 }
0120
0121 //=============================================================================================
0122
0123 void CastorShowerLibrary::loadEventInfo(TBranch* branch) {
0124   //////////////////////////////////////////////////////////
0125   //
0126   //  Get EventInfo from the "TBranch* branch" of Root file
0127   //  holding library events
0128   //
0129   //  Based on HFShowerLibrary::loadEventInfo
0130   //
0131   //////////////////////////////////////////////////////////
0132
0133   CastorShowerLibraryInfo tempInfo;
0134   auto* eventInfo = &tempInfo;
0135   branch->SetAddress(&eventInfo);
0136   branch->GetEntry(0);
0137   // Initialize shower library general parameters
0138
0139   totEvents = eventInfo->Energy.getNEvts();
0140   nBinsE = eventInfo->Energy.getNBins();
0141   nEvtPerBinE = eventInfo->Energy.getNEvtPerBin();
0142   SLenergies = eventInfo->Energy.getBin();
0143   nBinsEta = eventInfo->Eta.getNBins();
0144   nEvtPerBinEta = eventInfo->Eta.getNEvtPerBin();
0145   SLetas = eventInfo->Eta.getBin();
0146   nBinsPhi = eventInfo->Phi.getNBins();
0147   nEvtPerBinPhi = eventInfo->Phi.getNEvtPerBin();
0148   SLphis = eventInfo->Phi.getBin();
0149
0150   // Convert from GeV to MeV
0151   for (unsigned int i = 0; i < SLenergies.size(); i++) {
0152     SLenergies[i] *= CLHEP::GeV;
0153   }
0154
0155   edm::LogVerbatim("CastorShower") << " CastorShowerLibrary::loadEventInfo : "
0156                                    << "\n \n Total number of events     :  " << totEvents
0157                                    << "\n   Number of bins  (E)       :  " << nBinsE
0158                                    << "\n   Number of events/bin (E)  :  " << nEvtPerBinE
0159                                    << "\n   Number of bins  (Eta)       :  " << nBinsEta
0160                                    << "\n   Number of events/bin (Eta)  :  " << nEvtPerBinEta
0161                                    << "\n   Number of bins  (Phi)       :  " << nBinsPhi
0162                                    << "\n   Number of events/bin (Phi)  :  " << nEvtPerBinPhi << "\n";
0163
0164   std::stringstream ss1;
0165   ss1 << "CastorShowerLibrary: energies in GeV:\n";
0166   for (unsigned int i = 0; i < nBinsE; ++i) {
0167     if (i > 0 && i / 10 * 10 == i) {
0168       ss1 << "\n";
0169     }
0170     ss1 << " " << SLenergies[i] / CLHEP::GeV;
0171   }
0172   ss1 << "\nCastorShowerLibrary: etas:\n";
0173   for (unsigned int i = 0; i < nBinsEta; ++i) {
0174     if (i > 0 && i / 10 * 10 == i) {
0175       ss1 << "\n";
0176     }
0177     ss1 << " " << SLetas[i];
0178   }
0179   ss1 << "\nCastorShowerLibrary: phis:\n";
0180   for (unsigned int i = 0; i < nBinsPhi; ++i) {
0181     if (i > 0 && i / 10 * 10 == i) {
0182       ss1 << "\n";
0183     }
0184     ss1 << " " << SLphis[i];
0185   }
0186   edm::LogVerbatim("CastorShower") << ss1.str();
0187 }
0188
0189 //=============================================================================================
0190
0191 CastorShowerEvent CastorShowerLibrary::getShowerHits(const G4Step* aStep, bool& ok) {
0192   G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0193   G4Track* track = aStep->GetTrack();
0194
0195   const G4ThreeVector& hitPoint = preStepPoint->GetPosition();
0196
0197   CastorShowerEvent hit;
0198   hit.Clear();
0199
0200   ok = false;
0201   bool isEM = G4TrackToParticleID::isGammaElectronPositron(track);
0202   if (!isEM && !G4TrackToParticleID::isStableHadronIon(track)) {
0203     return hit;
0204   }
0205   ok = true;
0206
0207   double pin = preStepPoint->GetTotalEnergy();
0208   double zint = hitPoint.z();
0209   double R = sqrt(hitPoint.x() * hitPoint.x() + hitPoint.y() * hitPoint.y());
0210   double theta = atan2(R, std::abs(zint));
0211   double phiin = atan2(hitPoint.y(), hitPoint.x());
0212   double etain = -std::log(std::tan((CLHEP::pi - theta) * 0.5));
0213
0214   // Replace "interpolation/extrapolation" by new method "select" that just randomly
0215   // selects a record from the appropriate energy bin and fills its content to
0216   // "showerEvent" instance of class CastorShowerEvent
0217
0218   if (isEM) {
0219     hit = select(0, pin, etain, phiin);
0220   } else {
0221     hit = select(1, pin, etain, phiin);
0222   }
0223
0224   return hit;
0225 }
0226
0227 //=============================================================================================
0228
0229 CastorShowerEvent CastorShowerLibrary::getRecord(int type, int record) {
0230   //////////////////////////////////////////////////////////////
0231   //
0232   //  Retrieve event # "record" from the library and stores it
0233   //  into  vector<CastorShowerHit> showerHit;
0234   //
0235   //  Based on HFShowerLibrary::getRecord
0236   //
0237   //  Modified 02/02/09 by W. Carvalho
0238   //
0239   //////////////////////////////////////////////////////////////
0240
0241 #ifdef EDM_ML_DEBUG
0242   LogDebug("CastorShower") << "CastorShowerLibrary::getRecord: ";
0243 #endif
0244   int nrc = record;
0245   CastorShowerEvent retValue;
0246   CastorShowerEvent* showerEvent = &retValue;
0247   if (type > 0) {
0248     hadBranch->SetAddress(&showerEvent);
0249     hadBranch->GetEntry(nrc);
0250   } else {
0251     emBranch->SetAddress(&showerEvent);
0252     emBranch->GetEntry(nrc);
0253   }
0254
0255 #ifdef EDM_ML_DEBUG
0256   int nHit = showerEvent->getNhit();
0257   LogDebug("CastorShower") << "CastorShowerLibrary::getRecord: Record " << record << " of type " << type << " with "
0258                            << nHit << " CastorShowerHits";
0259
0260 #endif
0261   return retValue;
0262 }
0263
0264 //=======================================================================================
0265
0266 CastorShowerEvent CastorShowerLibrary::select(int type, double pin, double etain, double phiin) {
0267   ////////////////////////////////////////////////////////
0268   //
0269   //  Selects an event from the library based on
0270   //
0271   //    type:   0 --> em
0272   //           >0 --> had
0273   //    pin  :  momentum
0274   //    etain:  eta (if not given, disregard the eta binning
0275   //    phiin:  phi (if not given, disregard the phi binning
0276   //
0277   //  Created 30/01/09 by W. Carvalho
0278   //
0279   ////////////////////////////////////////////////////////
0280
0281   int irec;                    // to hold record number
0282   double r = G4UniformRand();  // to randomly select within an energy bin (r=[0,1])
0283
0284   // Randomly select a record from the right energy bin in the library,
0285   // based on track momentum (pin)
0286
0287   int ienergy = FindEnergyBin(pin);
0288   int ieta = FindEtaBin(etain);
0289 #ifdef EDM_ML_DEBUG
0290   if (verbose)
0291     edm::LogVerbatim("CastorShower") << " ienergy = " << ienergy;
0292   if (verbose)
0293     edm::LogVerbatim("CastorShower") << " ieta = " << ieta;
0294 #endif
0295
0296   int iphi;
0297   const double phiMin = 0.;
0298   const double phiMax = M_PI / 4.;  // 45 * (pi/180)  rad
0299   if (phiin < phiMin)
0300     phiin = phiin + M_PI;
0301   if (phiin >= phiMin && phiin < phiMax) {
0302     iphi = FindPhiBin(phiin);
0303   } else {
0304     double remainder = fmod(phiin, phiMax);
0305     phiin = phiMin + remainder;
0306     iphi = FindPhiBin(phiin);
0307   }
0308 #ifdef EDM_ML_DEBUG
0309   if (verbose)
0310     edm::LogVerbatim("CastorShower") << " iphi = " << iphi;
0311 #endif
0312   if (ienergy == -1)
0313     ienergy = 0;  // if pin < first bin, choose an event in the first one
0314   if (ieta == -1)
0315     ieta = 0;  // idem for eta
0316   if (iphi != -1)
0317     irec = int(nEvtPerBinE * ienergy + nEvtPerBinEta * ieta + nEvtPerBinPhi * (iphi + r));
0318   else
0319     irec = int(nEvtPerBinE * (ienergy + r));
0320
0321 #ifdef EDM_ML_DEBUG
0322   edm::LogVerbatim("CastorShower") << "CastorShowerLibrary:: Select record " << irec << " of type " << type;
0323 #endif
0324
0325   //  Retrieve record number "irec" from the library
0326   return getRecord(type, irec);
0327 }
0328
0329 //=======================================================================================
0330
0331 int CastorShowerLibrary::FindEnergyBin(double energy) {
0332   //
0333   // returns the integer index of the energy bin, taken from SLenergies vector
0334   // returns -1 if ouside valid range
0335   //
0336   if (energy >= SLenergies.back())
0337     return SLenergies.size() - 1;
0338
0339   unsigned int i = 0;
0340   for (; i < SLenergies.size() - 1; i++)
0341     if (energy >= SLenergies.at(i) && energy < SLenergies.at(i + 1))
0342       return (int)i;
0343
0344   // now i points to the last but 1 bin
0345   if (energy >= SLenergies.at(i))
0346     return (int)i;
0347   // energy outside bin range
0348   return -1;
0349 }
0350
0351 //=======================================================================================
0352
0353 int CastorShowerLibrary::FindEtaBin(double eta) {
0354   //
0355   // returns the integer index of the eta bin, taken from SLetas vector
0356   // returns -1 if ouside valid range
0357   //
0358   if (eta >= SLetas.back())
0359     return SLetas.size() - 1;
0360   unsigned int i = 0;
0361   for (; i < SLetas.size() - 1; i++)
0362     if (eta >= SLetas.at(i) && eta < SLetas.at(i + 1))
0363       return (int)i;
0364   // now i points to the last but 1 bin
0365   if (eta >= SLetas.at(i))
0366     return (int)i;
0367   // eta outside bin range
0368   return -1;
0369 }
0370
0371 //=======================================================================================
0372
0373 int CastorShowerLibrary::FindPhiBin(double phi) {
0374   //
0375   // returns the integer index of the phi bin, taken from SLphis vector
0376   // returns -1 if ouside valid range
0377   //
0378   // needs protection in case phi is outside range -pi,pi
0379   //
0380   if (phi >= SLphis.back())
0381     return SLphis.size() - 1;
0382   unsigned int i = 0;
0383   for (; i < SLphis.size() - 1; i++)
0384     if (phi >= SLphis.at(i) && phi < SLphis.at(i + 1))
0385       return (int)i;
0386   // now i points to the last but 1 bin
0387   if (phi >= SLphis.at(i))
0388     return (int)i;
0389   // phi outside bin range
0390   return -1;
0391 }
0392
0393 //=======================================================================================