Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-26 05:07:16 

0001 // -*- C++ -*-
0002 //
0003 // Package:     ShowerLibraryProducer
0004 // Class  :     CastorShowerLibraryMaker
0005 //
0006 // Implementation:
0007 //     <Notes on implementation>
0008 //
0009 // Original Author: P. Katsas
0010 //         Created: 02/2007
0011 //
0012 // Adapted by W. Carvalho , 02/2009
0013 //
0014 //////////////////////////////////////////////////////////////
0015 
0016 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0017 #include "DataFormats/Math/interface/Point3D.h"
0018 
0019 #include "FWCore/Utilities/interface/EDMException.h"
0020 #include "FWCore/Utilities/interface/Exception.h"
0021 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0023 #include "FWCore/ServiceRegistry/interface/Service.h"
0024 
0025 // Classes for shower library Root file
0026 #include "SimDataFormats/CaloHit/interface/CastorShowerLibraryInfo.h"
0027 #include "SimDataFormats/CaloHit/interface/CastorShowerEvent.h"
0028 
0029 #include "SimG4Core/Notification/interface/SimG4Exception.h"
0030 #include "SimG4Core/Notification/interface/BeginOfJob.h"
0031 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0032 #include "SimG4Core/Notification/interface/EndOfRun.h"
0033 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0034 #include "SimG4Core/Notification/interface/EndOfEvent.h"
0035 #include "SimG4Core/Notification/interface/Observer.h"
0036 #include "SimG4Core/Watcher/interface/SimWatcher.h"
0037 
0038 #include "SimG4CMS/Calo/interface/CaloG4Hit.h"
0039 #include "SimG4CMS/Calo/interface/CaloG4HitCollection.h"
0040 
0041 #include "SimG4CMS/Forward/interface/CastorNumberingScheme.h"
0042 
0043 #include "G4RunManager.hh"
0044 #include "G4SDManager.hh"
0045 #include "G4Step.hh"
0046 #include "G4Track.hh"
0047 #include "G4Event.hh"
0048 #include "G4PrimaryVertex.hh"
0049 #include "G4VProcess.hh"
0050 #include "G4HCofThisEvent.hh"
0051 #include "G4UserEventAction.hh"
0052 
0053 #include <CLHEP/Random/Randomize.h>
0054 #include <CLHEP/Units/SystemOfUnits.h>
0055 #include <CLHEP/Units/PhysicalConstants.h>
0056 
0057 #include "TROOT.h"
0058 #include "TFile.h"
0059 #include "TH1.h"
0060 #include "TH2.h"
0061 #include "TProfile.h"
0062 #include "TNtuple.h"
0063 #include "TRandom.h"
0064 #include "TLorentzVector.h"
0065 #include "TUnixSystem.h"
0066 #include "TSystem.h"
0067 #include "TMath.h"
0068 #include "TF1.h"
0069 
0070 #include <cassert>
0071 #include <cmath>
0072 #include <cstdlib>
0073 #include <iomanip>
0074 #include <iostream>
0075 #include <map>
0076 #include <set>
0077 #include <sstream>
0078 #include <string>
0079 #include <vector>
0080 
0081 typedef std::vector<std::vector<std::vector<std::vector<CastorShowerEvent> > > > SLBin3D;  // bin in energy, eta and phi
0082 
0083 class CastorShowerLibraryMaker : public SimWatcher,
0084                                  public Observer<const BeginOfJob*>,
0085                                  public Observer<const BeginOfRun*>,
0086                                  public Observer<const EndOfRun*>,
0087                                  public Observer<const BeginOfEvent*>,
0088                                  public Observer<const EndOfEvent*>,
0089                                  public Observer<const G4Step*> {
0090 public:
0091   CastorShowerLibraryMaker(const edm::ParameterSet& p);
0092   ~CastorShowerLibraryMaker() override;
0093 
0094 private:
0095   typedef int ebin;
0096   typedef int etabin;
0097   typedef int phibin;
0098   // private structures
0099   struct ShowerLib {
0100     CastorShowerLibraryInfo SLInfo;  // the info
0101     SLBin3D SLCollection;            // the showers
0102     std::vector<double> SLEnergyBins;
0103     std::vector<double> SLEtaBins;
0104     std::vector<double> SLPhiBins;
0105     unsigned int nEvtPerBinE;
0106     unsigned int nEvtPerBinEta;
0107     unsigned int nEvtPerBinPhi;
0108     std::vector<int> nEvtInBinE;
0109     std::vector<std::vector<int> > nEvtInBinEta;
0110     std::vector<std::vector<std::vector<int> > > nEvtInBinPhi;
0111   };
0112 
0113   // observer classes
0114   void update(const BeginOfJob* run) override;
0115   void update(const BeginOfRun* run) override;
0116   void update(const EndOfRun* run) override;
0117   void update(const BeginOfEvent* evt) override;
0118   void update(const EndOfEvent* evt) override;
0119   void update(const G4Step* step) override;
0120 
0121 private:
0122   void Finish();
0123 
0124   // Job general parameters
0125   int verbosity;
0126   std::string eventNtFileName;
0127 
0128   unsigned int NPGParticle;                  // number of particles requested to Particle Gun
0129   std::vector<int> PGParticleIDs;            //p. gun particle IDs
0130   bool DoHadSL;                              // true if hadronic SL should be produced
0131   bool DoEmSL;                               // true if electromag. SL should be produced
0132   bool InsideCastor;                         // true if particle step inside CASTOR
0133   bool DeActivatePhysicsProcess;             //cfg parameter: True if phys. proc. should be off from IP to Castor
0134   std::vector<G4PrimaryParticle*> thePrims;  // list of primaries for this event
0135 
0136   // Pointers for user defined class objects to be stored to Root file
0137   CastorShowerLibraryInfo* emInfo;
0138   CastorShowerLibraryInfo* hadInfo;
0139   CastorShowerEvent* emShower;
0140   CastorShowerEvent* hadShower;
0141   ShowerLib emSLHolder;
0142   ShowerLib hadSLHolder;
0143   ShowerLib* SLShowerptr;                          // pointer to the current shower collection (above)
0144   std::map<int, std::set<int> > MapOfSecondaries;  // map to hold all secondaries ID keyed by
0145                                                    // the PDG code of the primary
0146 
0147   std::map<int, G4ThreeVector> PrimaryMomentum;
0148   std::map<int, G4ThreeVector> PrimaryPosition;
0149   double MaxEta;  // limits the eta region, the lower limit is given by the SL bins
0150   double MaxPhi;  // limits the phi region, the lower limit is given by the SL bins
0151                   // private methods
0152   int FindEnergyBin(double e);
0153   int FindEtaBin(double eta);
0154   int FindPhiBin(double phi);
0155   bool SLacceptEvent(int, int, int);
0156   bool IsSLReady();
0157   void GetKinematics(G4PrimaryParticle*, double& px, double& py, double& pz, double& pInit, double& eta, double& phi);
0158   void GetKinematics(int, double& px, double& py, double& pz, double& pInit, double& eta, double& phi);
0159 
0160   std::vector<G4PrimaryParticle*> GetPrimary(const G4Event*);
0161   bool FillShowerEvent(CaloG4HitCollection*, CastorShowerEvent*, int);
0162   void InitSLHolder(ShowerLib&);
0163 
0164   void printSLstatus(int, int, int);
0165   int& SLnEvtInBinE(int ebin);
0166   int& SLnEvtInBinEta(int ebin, int etabin);
0167   int& SLnEvtInBinPhi(int ebin, int etabin, int phibin);
0168   bool SLisEBinFilled(int ebin);
0169   bool SLisEtaBinFilled(int ebin, int etabin);
0170   bool SLisPhiBinFilled(int ebin, int etabin, int phibin);
0171   void KillSecondaries(const G4Step* step);
0172   void GetMissingEnergy(CaloG4HitCollection*, double&, double&);
0173 
0174   // Root pointers
0175   TFile* theFile;
0176   TTree* theTree;
0177 
0178   int eventIndex;
0179   int stepIndex;  // ignore, please
0180 };
0181 
0182 CastorShowerLibraryMaker::CastorShowerLibraryMaker(const edm::ParameterSet& p)
0183     : NPGParticle(0),
0184       DoHadSL(false),
0185       DoEmSL(false),
0186       DeActivatePhysicsProcess(false),
0187       emShower(nullptr),
0188       hadShower(nullptr) {
0189   MapOfSecondaries.clear();
0190   hadInfo = nullptr;
0191   emInfo = nullptr;
0192   edm::ParameterSet p_SLM = p.getParameter<edm::ParameterSet>("CastorShowerLibraryMaker");
0193   verbosity = p_SLM.getParameter<int>("Verbosity");
0194   eventNtFileName = p_SLM.getParameter<std::string>("EventNtupleFileName");
0195   hadSLHolder.nEvtPerBinPhi = p_SLM.getParameter<int>("nhadEvents");
0196   emSLHolder.nEvtPerBinPhi = p_SLM.getParameter<int>("nemEvents");
0197   hadSLHolder.SLEnergyBins = p_SLM.getParameter<std::vector<double> >("SLhadEnergyBins");
0198   hadSLHolder.SLEtaBins = p_SLM.getParameter<std::vector<double> >("SLhadEtaBins");
0199   hadSLHolder.SLPhiBins = p_SLM.getParameter<std::vector<double> >("SLhadPhiBins");
0200   emSLHolder.SLEnergyBins = p_SLM.getParameter<std::vector<double> >("SLemEnergyBins");
0201   emSLHolder.SLEtaBins = p_SLM.getParameter<std::vector<double> >("SLemEtaBins");
0202   emSLHolder.SLPhiBins = p_SLM.getParameter<std::vector<double> >("SLemPhiBins");
0203   PGParticleIDs = p_SLM.getParameter<std::vector<int> >("PartID");
0204   DeActivatePhysicsProcess = p_SLM.getParameter<bool>("DeActivatePhysicsProcess");
0205   MaxPhi = p_SLM.getParameter<double>("SLMaxPhi");
0206   MaxEta = p_SLM.getParameter<double>("SLMaxEta");
0207   //
0208   NPGParticle = PGParticleIDs.size();
0209   for (unsigned int i = 0; i < PGParticleIDs.size(); i++) {
0210     switch (std::abs(PGParticleIDs.at(i))) {
0211       case 11:
0212       case 22:
0213         DoEmSL = true;
0214         break;
0215       default:
0216         DoHadSL = true;
0217     }
0218   }
0219   hadSLHolder.nEvtPerBinEta = (hadSLHolder.nEvtPerBinPhi) * (hadSLHolder.SLPhiBins.size());
0220   hadSLHolder.nEvtPerBinE = (hadSLHolder.nEvtPerBinEta) * (hadSLHolder.SLEtaBins.size());
0221   emSLHolder.nEvtPerBinEta = (emSLHolder.nEvtPerBinPhi) * (emSLHolder.SLPhiBins.size());
0222   emSLHolder.nEvtPerBinE = (emSLHolder.nEvtPerBinEta) * (emSLHolder.SLEtaBins.size());
0223 
0224   edm::LogVerbatim("HcalSim") << "============================================================================";
0225   edm::LogVerbatim("HcalSim") << "CastorShowerLibraryMaker:: Initialized as observer";
0226   edm::LogVerbatim("HcalSim") << " Event Ntuple will be created";
0227   edm::LogVerbatim("HcalSim") << " Event Ntuple file: " << eventNtFileName;
0228   edm::LogVerbatim("HcalSim") << " Number of Hadronic events in E   bins: " << hadSLHolder.nEvtPerBinE;
0229   edm::LogVerbatim("HcalSim") << " Number of Hadronic events in Eta bins: " << hadSLHolder.nEvtPerBinEta;
0230   edm::LogVerbatim("HcalSim") << " Number of Hadronic events in Phi bins: " << hadSLHolder.nEvtPerBinPhi;
0231   edm::LogVerbatim("HcalSim") << " Number of Electromag. events in E   bins: " << emSLHolder.nEvtPerBinE;
0232   edm::LogVerbatim("HcalSim") << " Number of Electromag. events in Eta bins: " << emSLHolder.nEvtPerBinEta;
0233   edm::LogVerbatim("HcalSim") << " Number of Electromag. events in Phi bins: " << emSLHolder.nEvtPerBinPhi;
0234   edm::LogVerbatim("HcalSim") << "============================================================================\n";
0235 
0236   // Initializing the SL collections
0237   InitSLHolder(hadSLHolder);
0238   InitSLHolder(emSLHolder);
0239 }
0240 void CastorShowerLibraryMaker::InitSLHolder(ShowerLib& showerholder) {
0241   int nBinsE, nBinsEta, nBinsPhi, nEvtPerBinPhi;
0242   nBinsE = showerholder.SLEnergyBins.size();
0243   nBinsEta = showerholder.SLEtaBins.size();
0244   nBinsPhi = showerholder.SLPhiBins.size();
0245   nEvtPerBinPhi = showerholder.nEvtPerBinPhi;
0246   //
0247   // Info
0248   //
0249   showerholder.SLInfo.Energy.setNEvts(nEvtPerBinPhi * nBinsPhi * nBinsEta * nBinsE);
0250   showerholder.SLInfo.Energy.setNEvtPerBin(nEvtPerBinPhi * nBinsPhi * nBinsEta);
0251   showerholder.SLInfo.Energy.setNBins(nBinsE);
0252   showerholder.SLInfo.Energy.setBin(showerholder.SLEnergyBins);
0253   //
0254   showerholder.SLInfo.Eta.setNEvts(nEvtPerBinPhi * nBinsPhi * nBinsEta);
0255   showerholder.SLInfo.Eta.setNEvtPerBin(nEvtPerBinPhi * nBinsPhi);
0256   showerholder.SLInfo.Eta.setNBins(nBinsEta);
0257   showerholder.SLInfo.Eta.setBin(showerholder.SLEtaBins);
0258   //
0259   showerholder.SLInfo.Phi.setNEvts(nEvtPerBinPhi * nBinsPhi);
0260   showerholder.SLInfo.Phi.setNEvtPerBin(nEvtPerBinPhi);
0261   showerholder.SLInfo.Phi.setNBins(nBinsPhi);
0262   showerholder.SLInfo.Phi.setBin(showerholder.SLPhiBins);
0263   //
0264   // Shower
0265   showerholder.SLCollection.assign(nBinsE, std::vector<std::vector<std::vector<CastorShowerEvent> > >());
0266   showerholder.nEvtInBinE.assign(nBinsE, 0);
0267   showerholder.nEvtInBinEta.assign(nBinsE, std::vector<int>(0));
0268   showerholder.nEvtInBinPhi.assign(nBinsE, std::vector<std::vector<int> >());
0269   for (int i = 0; i < nBinsE; i++) {
0270     showerholder.SLCollection.at(i).assign(nBinsEta, std::vector<std::vector<CastorShowerEvent> >());
0271     showerholder.nEvtInBinEta.at(i).assign(nBinsEta, 0);
0272     showerholder.nEvtInBinPhi.at(i).assign(nBinsEta, std::vector<int>(0));
0273     for (int j = 0; j < nBinsEta; j++) {
0274       showerholder.SLCollection.at(i).at(j).assign(nBinsPhi, std::vector<CastorShowerEvent>());
0275       showerholder.nEvtInBinPhi.at(i).at(j).assign(nBinsPhi, 0);
0276       for (int k = 0; k < nBinsPhi; k++)
0277         showerholder.SLCollection.at(i).at(j).at(k).assign(nEvtPerBinPhi, CastorShowerEvent());
0278     }
0279   }
0280 }
0281 
0282 //===============================================================================================
0283 
0284 CastorShowerLibraryMaker::~CastorShowerLibraryMaker() {
0285   Finish();
0286 
0287   edm::LogVerbatim("HcalSim") << "CastorShowerLibraryMaker: End of process";
0288 }
0289 
0290 //=================================================================== per EVENT
0291 void CastorShowerLibraryMaker::update(const BeginOfJob* job) {
0292   edm::LogVerbatim("HcalSim") << " CastorShowerLibraryMaker::Starting new job ";
0293 }
0294 
0295 //==================================================================== per RUN
0296 void CastorShowerLibraryMaker::update(const BeginOfRun* run) {
0297   edm::LogVerbatim("HcalSim") << "\nCastorShowerLibraryMaker: Starting Run";
0298 
0299   edm::LogVerbatim("HcalSim") << "CastorShowerLibraryMaker: output event root file created";
0300 
0301   TString eventfilename = eventNtFileName;
0302   theFile = new TFile(eventfilename, "RECREATE");
0303   theTree = new TTree("CastorCherenkovPhotons", "Cherenkov Photons");
0304 
0305   Int_t split = 1;
0306   Int_t bsize = 64000;
0307   emInfo = new CastorShowerLibraryInfo();
0308   emShower = new CastorShowerEvent();
0309   hadInfo = new CastorShowerLibraryInfo();
0310   hadShower = new CastorShowerEvent();
0311   // Create Branchs
0312   theTree->Branch("emShowerLibInfo.", "CastorShowerLibraryInfo", &emInfo, bsize, split);
0313   theTree->Branch("emParticles.", "CastorShowerEvent", &emShower, bsize, split);
0314   theTree->Branch("hadShowerLibInfo.", "CastorShowerLibraryInfo", &hadInfo, bsize, split);
0315   theTree->Branch("hadParticles.", "CastorShowerEvent", &hadShower, bsize, split);
0316 
0317   // set the Info for electromagnetic shower
0318   // set the energy bins info
0319   emInfo->Energy.setNEvts(emSLHolder.nEvtPerBinE * emSLHolder.SLEnergyBins.size());
0320   emInfo->Energy.setNBins(emSLHolder.SLEnergyBins.size());
0321   emInfo->Energy.setNEvtPerBin(emSLHolder.nEvtPerBinE);
0322   emInfo->Energy.setBin(emSLHolder.SLEnergyBins);
0323   // set the eta bins info
0324   emInfo->Eta.setNEvts(emSLHolder.nEvtPerBinEta * emSLHolder.SLEtaBins.size());
0325   emInfo->Eta.setNBins(emSLHolder.SLEtaBins.size());
0326   emInfo->Eta.setNEvtPerBin(emSLHolder.nEvtPerBinEta);
0327   emInfo->Eta.setBin(emSLHolder.SLEtaBins);
0328   // set the eta bins info
0329   emInfo->Phi.setNEvts(emSLHolder.nEvtPerBinPhi * emSLHolder.SLPhiBins.size());
0330   emInfo->Phi.setNBins(emSLHolder.SLPhiBins.size());
0331   emInfo->Phi.setNEvtPerBin(emSLHolder.nEvtPerBinPhi);
0332   emInfo->Phi.setBin(emSLHolder.SLPhiBins);
0333   // The same for the hadronic shower
0334   // set the energy bins info
0335   hadInfo->Energy.setNEvts(hadSLHolder.nEvtPerBinE * hadSLHolder.SLEnergyBins.size());
0336   hadInfo->Energy.setNBins(hadSLHolder.SLEnergyBins.size());
0337   hadInfo->Energy.setNEvtPerBin(hadSLHolder.nEvtPerBinE);
0338   hadInfo->Energy.setBin(hadSLHolder.SLEnergyBins);
0339   // set the eta bins info
0340   hadInfo->Eta.setNEvts(hadSLHolder.nEvtPerBinEta * hadSLHolder.SLEtaBins.size());
0341   hadInfo->Eta.setNBins(hadSLHolder.SLEtaBins.size());
0342   hadInfo->Eta.setNEvtPerBin(hadSLHolder.nEvtPerBinEta);
0343   hadInfo->Eta.setBin(hadSLHolder.SLEtaBins);
0344   // set the eta bins info
0345   hadInfo->Phi.setNEvts(hadSLHolder.nEvtPerBinPhi * hadSLHolder.SLPhiBins.size());
0346   hadInfo->Phi.setNBins(hadSLHolder.SLPhiBins.size());
0347   hadInfo->Phi.setNEvtPerBin(hadSLHolder.nEvtPerBinPhi);
0348   hadInfo->Phi.setBin(hadSLHolder.SLPhiBins);
0349   // int flag = theTree->GetBranch("CastorShowerLibInfo")->Fill();
0350   // Loop on all leaves of this branch to fill Basket buffer.
0351   // The function returns the number of bytes committed to the memory basket.
0352   // If a write error occurs, the number of bytes returned is -1.
0353   // If no data are written, because e.g. the branch is disabled,
0354   // the number of bytes returned is 0.
0355   // if(flag==-1) {
0356   //    edm::LogVerbatim("CastorAnalyzer") << " WARNING: Error writing to Branch \"CastorShowerLibInfo\" \n" ;
0357   // } else
0358   // if(flag==0) {
0359   //    edm::LogVerbatim("CastorAnalyzer") << " WARNING: No data written to Branch \"CastorShowerLibInfo\" \n" ;
0360   // }
0361 
0362   // Initialize "accounting" variables
0363 
0364   eventIndex = 0;
0365 }
0366 
0367 //=================================================================== per EVENT
0368 void CastorShowerLibraryMaker::update(const BeginOfEvent* evt) {
0369   eventIndex++;
0370   stepIndex = 0;
0371   InsideCastor = false;
0372   PrimaryMomentum.clear();
0373   PrimaryPosition.clear();
0374   int NAccepted = 0;
0375   // reset the pointers to the shower objects
0376   SLShowerptr = nullptr;
0377   MapOfSecondaries.clear();
0378   thePrims.clear();
0379   G4EventManager* e_mgr = G4EventManager::GetEventManager();
0380   if (IsSLReady()) {
0381     printSLstatus(-1, -1, -1);
0382     update((EndOfRun*)nullptr);
0383     return;
0384   }
0385 
0386   thePrims = GetPrimary((*evt)());
0387   for (unsigned int i = 0; i < thePrims.size(); i++) {
0388     G4PrimaryParticle* thePrim = thePrims.at(i);
0389     int particleType = thePrim->GetPDGcode();
0390 
0391     std::string SLType("");
0392     if (particleType == 11) {
0393       SLShowerptr = &emSLHolder;
0394       SLType = "Electromagnetic";
0395     } else {
0396       SLShowerptr = &hadSLHolder;
0397       SLType = "Hadronic";
0398     }
0399     double px = 0., py = 0., pz = 0., pInit = 0., eta = 0., phi = 0.;
0400     GetKinematics(thePrim, px, py, pz, pInit, eta, phi);
0401     int ebin = FindEnergyBin(pInit);
0402     int etabin = FindEtaBin(eta);
0403     int phibin = FindPhiBin(phi);
0404     if (verbosity)
0405       edm::LogVerbatim("HcalSim") << "\n========================================================================"
0406                                   << "\nBeginOfEvent: E   : " << pInit << "\t  bin : " << ebin
0407                                   << "\n              Eta : " << eta << "\t  bin : " << etabin
0408                                   << "\n              Phi : " << phi << "\t  bin : " << phibin
0409                                   << "\n========================================================================";
0410 
0411     if (ebin < 0 || etabin < 0 || phibin < 0)
0412       continue;
0413     bool accept = false;
0414     if (!(SLacceptEvent(ebin, etabin, phibin))) {
0415       /*
0416 // To increase the chance of a particle arriving at CASTOR inside a not full bin,
0417 // check if there is available phase space in the neighboring bins
0418         unsigned int ebin_min = std::max(0,ebin-3);
0419         unsigned int eta_bin_min = std::max(0,etabin-2);
0420         unsigned int eta_bin_max = std::min(etabin,etabin+2);
0421         unsigned int phi_bin_min = std::max(0,phibin-2);
0422         unsigned int phi_bin_max = std::min(phibin,phibin+2);
0423         for(unsigned int i_ebin=ebin_min;i_ebin<=(unsigned int)ebin;i_ebin++) {
0424            for (unsigned int i_etabin=eta_bin_min;i_etabin<=eta_bin_max;i_etabin++) {
0425                for (unsigned int i_phibin=phi_bin_min;i_phibin<=phi_bin_max;i_phibin++) {
0426                    if (SLacceptEvent((int)i_ebin,(int)i_etabin,(int)i_phibin)) {accept=true;break;}
0427                }
0428                if (accept) break;
0429            }
0430            if (accept) break;
0431         }
0432 */
0433       if (!accept)
0434         edm::LogVerbatim("CastorShowerLibraryMaker")
0435             << "Event not accepted for ebin=" << ebin << ",etabin=" << etabin << ",phibin=" << phibin;
0436     } else {
0437       accept = true;
0438     }
0439     if (accept)
0440       NAccepted++;
0441   }
0442 
0443   if (NAccepted == 0) {
0444     const_cast<G4Event*>((*evt)())->SetEventAborted();
0445     const_cast<G4Event*>((*evt)())->KeepTheEvent((G4bool) false);
0446     e_mgr->AbortCurrentEvent();
0447   }
0448   SLShowerptr = nullptr;
0449   //
0450   edm::LogVerbatim("HcalSim") << "CastorShowerLibraryMaker: Processing Event Number: " << eventIndex;
0451 }
0452 
0453 //=================================================================== per STEP
0454 void CastorShowerLibraryMaker::update(const G4Step* aStep) {
0455   static thread_local int CurrentPrimary = 0;
0456   G4Track* trk = aStep->GetTrack();
0457   int pvec_size;
0458   if (trk->GetCurrentStepNumber() == 1) {
0459     if (trk->GetParentID() == 0) {
0460       CurrentPrimary = (int)trk->GetDynamicParticle()->GetPDGcode();
0461       if (CurrentPrimary == 0)
0462         SimG4Exception("CastorShowerLibraryMaker::update(G4Step) -> Primary particle undefined");
0463       InsideCastor = false;
0464       // Deactivate the physics process
0465       if (DeActivatePhysicsProcess) {
0466         G4ProcessManager* p_mgr = trk->GetDefinition()->GetProcessManager();
0467         G4ProcessVector* pvec = p_mgr->GetProcessList();
0468         pvec_size = pvec->size();
0469         for (int i = 0; i < pvec_size; i++) {
0470           G4VProcess* proc = (*pvec)(i);
0471           if (proc->GetProcessName() != "Transportation" && proc->GetProcessName() != "Decay") {
0472             edm::LogVerbatim("HcalSim") << "DeActivating process: " << proc->GetProcessName();
0473             p_mgr->SetProcessActivation(proc, false);
0474           }
0475         }
0476       }
0477       // move track to z of CASTOR
0478       G4ThreeVector pos;
0479       pos.setZ(-14390);
0480       double t = std::abs((pos.z() - trk->GetPosition().z())) / trk->GetVelocity();
0481       double r = (pos.z() - trk->GetPosition().z()) / trk->GetMomentum().cosTheta();
0482       pos.setX(r * sin(trk->GetMomentum().theta()) * cos(trk->GetMomentum().phi()) + trk->GetPosition().x());
0483       pos.setY(r * sin(trk->GetMomentum().theta()) * sin(trk->GetMomentum().phi()) + trk->GetPosition().y());
0484       trk->SetPosition(pos);
0485       trk->SetGlobalTime(trk->GetGlobalTime() + t);
0486       trk->AddTrackLength(r);
0487     } else if (!InsideCastor) {
0488       edm::LogVerbatim("HcalSim") << "CastorShowerLibraryMaker::update(G4Step) -> Killing spurious track";
0489       trk->SetTrackStatus(fKillTrackAndSecondaries);
0490       return;
0491     }
0492     MapOfSecondaries[CurrentPrimary].insert((int)trk->GetTrackID());
0493   }
0494   // Checks if primary already inside CASTOR
0495   std::string CurVolume = trk->GetVolume()->GetName();
0496   if (!InsideCastor && (
0497                            //CurVolume=="C3EF"||CurVolume=="C4EF"||CurVolume=="CAEL"||
0498                            //CurVolume=="CAHL"||CurVolume=="C3HF"||CurVolume=="C4HF")) {
0499                            //CurVolume=="CastorB"||
0500                            CurVolume == "CAST")) {
0501     //CurVolume=="CAIR")) {
0502     InsideCastor = true;
0503     // Activate the physics process
0504     if (trk->GetParentID() == 0 && DeActivatePhysicsProcess) {
0505       G4ProcessManager* p_mgr = trk->GetDefinition()->GetProcessManager();
0506       G4ProcessVector* pvec = p_mgr->GetProcessList();
0507       pvec_size = pvec->size();
0508       for (int i = 0; i < pvec_size; i++) {
0509         G4VProcess* proc = (*pvec)(i);
0510         if (proc->GetProcessName() != "Transportation" && proc->GetProcessName() != "Decay") {
0511           edm::LogVerbatim("HcalSim") << "  Activating process: " << proc->GetProcessName();
0512           p_mgr->SetProcessActivation(proc, true);
0513         }
0514       }
0515     }
0516     //PrimaryMomentum[CurrentPrimary]=aStep->GetPreStepPoint()->GetMomentum();
0517     // check fiducial eta and phi
0518     if (trk->GetMomentum().phi() > MaxPhi || trk->GetMomentum().eta() > MaxEta) {
0519       trk->SetTrackStatus(fKillTrackAndSecondaries);
0520       InsideCastor = false;
0521       return;
0522     }
0523     PrimaryMomentum[CurrentPrimary] = trk->GetMomentum();
0524     PrimaryPosition[CurrentPrimary] = trk->GetPosition();
0525     KillSecondaries(aStep);
0526     return;
0527   }
0528   // Kill the secondaries if they have been produced before entering castor
0529   if (CurrentPrimary != 0 && trk->GetParentID() == 0 && !InsideCastor) {
0530     KillSecondaries(aStep);
0531     if (verbosity) {
0532       double pre_phi = aStep->GetPreStepPoint()->GetMomentum().phi();
0533       double cur_phi = trk->GetMomentum().phi();
0534       if (pre_phi != cur_phi) {
0535         edm::LogVerbatim("HcalSim") << "Primary track phi :  " << pre_phi << " changed in current step: " << cur_phi
0536                                     << " by processes:";
0537         const G4VProcess* proc = aStep->GetPreStepPoint()->GetProcessDefinedStep();
0538         edm::LogVerbatim("HcalSim") << "           " << proc->GetProcessName() << "  In volume " << CurVolume;
0539       }
0540     }
0541   }
0542 
0543   //==============================================
0544   /*
0545 */
0546   /*
0547   if(aStep->IsFirstStepInVolume()) { 
0548     edm::LogVerbatim("CastorShowerLibraryMaker") << "CastorShowerLibraryMaker::update(const G4Step * aStep):"
0549                                              << "\n IsFirstStepInVolume , " 
0550                                              << "time = " << aStep->GetTrack()->GetGlobalTime() ; 
0551   }
0552   stepIndex++;
0553 */
0554 }
0555 
0556 //================= End of EVENT ===============
0557 void CastorShowerLibraryMaker::update(const EndOfEvent* evt) {
0558   // check if the job is done!
0559   if ((*evt)()->IsAborted()) {
0560     edm::LogVerbatim("HcalSim") << "\n========================================================================"
0561                                 << "\nEndOfEvent: EVENT ABORTED"
0562                                 << "\n========================================================================";
0563     return;
0564   }
0565   //DynamicRangeFlatRandomEGunProducer* pgun = edm::DynamicRangeFlatRandomEGunKernel::get_instance();
0566   //edm::LogVerbatim("HcalSim") << pgun->EGunMaxE();
0567   /*
0568   edm::LogVerbatim("HcalSim") << "Minimum energy in Particle Gun : " << pgun->EGunMinE() << "\nMaximum energy in Particle Gun : " << pgun->EGunMaxE();
0569 */
0570   if (verbosity)
0571     edm::LogVerbatim("HcalSim") << "CastorShowerLibraryMaker: End of Event: " << eventIndex;
0572   // Get the pointer to the primary particle
0573   if (thePrims.empty()) {
0574     edm::LogVerbatim("CastorShowerLibraryMaker") << "No valid primary particle found. Skipping event";
0575     return;
0576   }
0577   // access to the G4 hit collections
0578   G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
0579   int CAFIid = G4SDManager::GetSDMpointer()->GetCollectionID("CastorFI");
0580   CaloG4HitCollection* theCAFI = (CaloG4HitCollection*)allHC->GetHC(CAFIid);
0581   if (verbosity)
0582     edm::LogVerbatim("CastorShowerLibraryMaker") << " update(*evt) --> accessed all HC ";
0583   edm::LogVerbatim("CastorShowerLibraryMaker") << "Found " << theCAFI->entries() << " hits in G4HitCollection";
0584   if (theCAFI->entries() == 0) {
0585     edm::LogVerbatim("CastorShowerLibraryMaker") << "\n Empty G4HitCollection";
0586     return;
0587   }
0588 
0589   // Loop over primaries
0590   int NEvtAccepted = 0;
0591   int NHitInEvent = 0;
0592   for (unsigned int i = 0; i < thePrims.size(); i++) {
0593     G4PrimaryParticle* thePrim = thePrims.at(i);
0594     if (!thePrim) {
0595       edm::LogVerbatim("CastorShowerLibraryMaker") << "nullptr Pointer to the primary";
0596       continue;
0597     }
0598     // Check primary particle type
0599     int particleType = thePrim->GetPDGcode();
0600 
0601     // set the pointer to the shower collection
0602     std::string SLType("");
0603     if (particleType == 11) {
0604       SLShowerptr = &emSLHolder;
0605       SLType = "Electromagnetic";
0606     } else {
0607       SLShowerptr = &hadSLHolder;
0608       SLType = "Hadronic";
0609     }
0610     edm::LogVerbatim("CastorShowerLibraryMaker") << "\n Primary (thePrim) trackID is " << thePrim->GetTrackID() << "\n";
0611 
0612     // Obtain primary particle's initial momentum (pInit)
0613     double px = 0., py = 0., pz = 0., pInit = 0., eta = 0., phi = 0.;
0614     GetKinematics(particleType, px, py, pz, pInit, eta, phi);
0615     // Check if current event falls into any bin
0616     // first: energy
0617     if (pInit == 0) {
0618       edm::LogVerbatim("CastorShowerLibraryMaker") << "Primary did not hit CASTOR";
0619       continue;
0620     }
0621     int ebin = FindEnergyBin(pInit);
0622     int etabin = FindEtaBin(eta);
0623     int phibin = FindPhiBin(phi);
0624     edm::LogVerbatim("HcalSim") << SLType;
0625     printSLstatus(ebin, etabin, phibin);
0626     if (!SLacceptEvent(ebin, etabin, phibin)) {
0627       edm::LogVerbatim("CastorShowerLibraryMaker")
0628           << "Event not accepted for ebin=" << ebin << ",etabin=" << etabin << ",phibin=" << phibin << "(" << pInit
0629           << "," << eta << "," << phi << ")";
0630       continue;
0631     }
0632     //
0633     // event passed. Fill the vector accordingly
0634     //
0635     // Look for the Hit Collection
0636     edm::LogVerbatim("CastorShowerLibraryMaker")
0637         << "\n CastorShowerLibraryMaker::update(EndOfEvent * evt) - event #" << (*evt)()->GetEventID();
0638 
0639     /*
0640      edm::LogVerbatim("HcalSim") << "Number of collections : " << allHC->GetNumberOfCollections();
0641      for(int ii = 0;ii<allHC->GetNumberOfCollections();ii++) 
0642         edm::LogVerbatim("HcalSim") << "Name of collection " << ii << " : " << allHC->GetHC(ii)->GetName();
0643 */
0644 
0645     CastorShowerEvent* shower = nullptr;
0646     int cur_evt_idx = SLShowerptr->nEvtInBinPhi.at(ebin).at(etabin).at(phibin);
0647     shower = &(SLShowerptr->SLCollection.at(ebin).at(etabin).at(phibin).at(cur_evt_idx));
0648 
0649     // Get Hit information
0650     if (FillShowerEvent(theCAFI, shower, particleType)) {
0651       //  Primary particle information
0652       /*
0653         edm::LogVerbatim("CastorShowerLibraryMaker") << "New SL event: Primary = " << particleType << "; Energy = " << pInit << "; Eta = " << eta << "; Phi = " << phi << "; Nhits = " << shower->getNhit();
0654 */
0655       shower->setPrimE(pInit);
0656       shower->setPrimEta(eta);
0657       shower->setPrimPhi(phi);
0658       shower->setPrimX(PrimaryPosition[particleType].x());
0659       shower->setPrimY(PrimaryPosition[particleType].y());
0660       shower->setPrimZ(PrimaryPosition[particleType].z());
0661       SLnEvtInBinE(ebin)++;
0662       SLnEvtInBinEta(ebin, etabin)++;
0663       SLnEvtInBinPhi(ebin, etabin, phibin)++;
0664       NHitInEvent += shower->getNhit();
0665       NEvtAccepted++;
0666     } else {
0667       shower->Clear();
0668     }
0669   }
0670   // Check for unassociated energy
0671   int thecafi_entries = theCAFI->entries();
0672   if (NEvtAccepted == int(thePrims.size()) && thecafi_entries != NHitInEvent) {
0673     edm::LogWarning("HcalSim") << "WARNING: Inconsistent Number of Hits -> Hits in collection: " << theCAFI->entries()
0674                                << "   Hits in the showers: " << NHitInEvent;
0675     double miss_energy = 0;
0676     double tot_energy = 0;
0677     GetMissingEnergy(theCAFI, miss_energy, tot_energy);
0678     if (miss_energy > 0) {
0679       edm::LogVerbatim("HcalSim") << "Total missing energy: " << miss_energy
0680                                   << " for an incident energy: " << tot_energy;
0681     }
0682   }
0683 }
0684 
0685 //========================= End of RUN ======================
0686 void CastorShowerLibraryMaker::update(const EndOfRun* run) {
0687   // Fill the tree with the collected objects
0688   if (!IsSLReady())
0689     SimG4Exception("\n\nShower Library  NOT READY.\n\n");
0690 
0691   unsigned int ibine, ibineta, ibinphi, ievt;  // indexes for em shower
0692   unsigned int jbine, jbineta, jbinphi, jevt;  // indexes for had shower
0693 
0694   ibine = ibineta = ibinphi = ievt = jbine = jbineta = jbinphi = jevt = 0;
0695 
0696   int nEvtInTree = 0;
0697   int nEMevt = emSLHolder.nEvtPerBinE * emSLHolder.SLEnergyBins.size();
0698   int nHadevt = hadSLHolder.nEvtPerBinE * hadSLHolder.SLEnergyBins.size();
0699   int maxEvtInTree = std::max(nEMevt, nHadevt);
0700 
0701   emInfo = &emSLHolder.SLInfo;
0702   hadInfo = &hadSLHolder.SLInfo;
0703 
0704   while (nEvtInTree < maxEvtInTree) {
0705     if (emShower)
0706       emShower->Clear();
0707     if (hadShower)
0708       hadShower->Clear();
0709     while (ibine < emSLHolder.SLEnergyBins.size() && nEMevt > 0) {
0710       emShower = &(emSLHolder.SLCollection.at(ibine).at(ibineta).at(ibinphi).at(ievt));
0711       ievt++;
0712       if (ievt == emSLHolder.nEvtPerBinPhi) {
0713         ievt = 0;
0714         ibinphi++;
0715       }
0716       if (ibinphi == emSLHolder.SLPhiBins.size()) {
0717         ibinphi = 0;
0718         ibineta++;
0719       }
0720       if (ibineta == emSLHolder.SLEtaBins.size()) {
0721         ibineta = 0;
0722         ibine++;
0723       }
0724       break;
0725     }
0726     while (jbine < hadSLHolder.SLEnergyBins.size() && nHadevt > 0) {
0727       hadShower = &(hadSLHolder.SLCollection.at(jbine).at(jbineta).at(jbinphi).at(jevt));
0728       jevt++;
0729       if (jevt == hadSLHolder.nEvtPerBinPhi) {
0730         jevt = 0;
0731         jbinphi++;
0732       }
0733       if (jbinphi == hadSLHolder.SLPhiBins.size()) {
0734         jbinphi = 0;
0735         jbineta++;
0736       }
0737       if (jbineta == hadSLHolder.SLEtaBins.size()) {
0738         jbineta = 0;
0739         jbine++;
0740       }
0741       break;
0742     }
0743     theTree->Fill();
0744     nEvtInTree++;
0745     if (nEvtInTree == 1) {
0746       theTree->SetBranchStatus("emShowerLibInfo.", false);
0747       theTree->SetBranchStatus("hadShowerLibInfo.", false);
0748     }
0749   }
0750   // check if run is nullptr and exit
0751   if (run == nullptr)
0752     throw SimG4Exception("\n\nNumber of needed trigger events reached in CastorShowerLibraryMaker\n\n");
0753 }
0754 
0755 //============================================================
0756 void CastorShowerLibraryMaker::Finish() {
0757   // if (doNTcastorevent) {
0758 
0759   theFile->cd();
0760   theTree->Write("", TObject::kOverwrite);
0761   edm::LogVerbatim("HcalSim") << "CastorShowerLibraryMaker: Ntuple event written";
0762   theFile->Close();
0763   edm::LogVerbatim("HcalSim") << "CastorShowerLibraryMaker: Event file closed";
0764 
0765   // Delete pointers to objects, now that TTree has been written and TFile closed
0766   //  delete      info;
0767   //  delete  emShower;
0768   //  delete hadShower;
0769   // }
0770 }
0771 int CastorShowerLibraryMaker::FindEnergyBin(double energy) {
0772   //
0773   // returns the integer index of the energy bin, taken from SLenergies vector
0774   // returns -1 if ouside valid range
0775   //
0776   if (!SLShowerptr) {
0777     edm::LogVerbatim("CastorShowerLibraryMaker") << "\n\nFindEnergyBin can be called only after BeginOfEvent\n\n";
0778     throw SimG4Exception("\n\nnullptr Pointer to the shower library.\n\n");
0779   }
0780   const std::vector<double>& SLenergies = SLShowerptr->SLEnergyBins;
0781   if (energy >= SLenergies.back())
0782     return SLenergies.size() - 1;
0783 
0784   unsigned int i = 0;
0785   for (; i < SLenergies.size() - 1; i++)
0786     if (energy >= SLenergies.at(i) && energy < SLenergies.at(i + 1))
0787       return (int)i;
0788 
0789   // now i points to the last but 1 bin
0790   if (energy >= SLenergies.at(i))
0791     return (int)i;
0792   // energy outside bin range
0793   return -1;
0794 }
0795 int CastorShowerLibraryMaker::FindEtaBin(double eta) {
0796   //
0797   // returns the integer index of the eta bin, taken from SLetas vector
0798   // returns -1 if ouside valid range
0799   //
0800   if (!SLShowerptr) {
0801     edm::LogVerbatim("CastorShowerLibraryMaker") << "\n\nFindEtaBin can be called only after BeginOfEvent\n\n";
0802     throw SimG4Exception("\n\nnullptr Pointer to the shower library.\n\n");
0803   }
0804   const std::vector<double>& SLetas = SLShowerptr->SLEtaBins;
0805   if (eta >= SLetas.back())
0806     return SLetas.size() - 1;
0807   unsigned int i = 0;
0808   for (; i < SLetas.size() - 1; i++)
0809     if (eta >= SLetas.at(i) && eta < SLetas.at(i + 1))
0810       return (int)i;
0811   // now i points to the last but 1 bin
0812   if (eta >= SLetas.at(i))
0813     return (int)i;
0814   // eta outside bin range
0815   return -1;
0816 }
0817 int CastorShowerLibraryMaker::FindPhiBin(double phi) {
0818   //
0819   // returns the integer index of the phi bin, taken from SLphis vector
0820   // returns -1 if ouside valid range
0821   //
0822   // needs protection in case phi is outside range -pi,pi
0823   //
0824   if (!SLShowerptr) {
0825     edm::LogVerbatim("CastorShowerLibraryMaker") << "\n\nFindPhiBin can be called only after BeginOfEvent\n\n";
0826     throw SimG4Exception("\n\nnullptr Pointer to the shower library.\n\n");
0827   }
0828   const std::vector<double>& SLphis = SLShowerptr->SLPhiBins;
0829   if (phi >= SLphis.back())
0830     return SLphis.size() - 1;
0831   unsigned int i = 0;
0832   for (; i < SLphis.size() - 1; i++)
0833     if (phi >= SLphis.at(i) && phi < SLphis.at(i + 1))
0834       return (int)i;
0835   // now i points to the last but 1 bin
0836   if (phi >= SLphis.at(i))
0837     return (int)i;
0838   // phi outside bin range
0839   return -1;
0840 }
0841 bool CastorShowerLibraryMaker::IsSLReady() {
0842   // at this point, the pointer to the shower library should be nullptr
0843   if (SLShowerptr) {
0844     edm::LogVerbatim("CastorShowerLibraryMaker") << "\n\nIsSLReady must be called when a new event starts.\n\n";
0845     throw SimG4Exception("\n\nNOT nullptr Pointer to the shower library.\n\n");
0846   }
0847   // it is enough to check if all the energy bin is filled
0848   if (DoEmSL) {
0849     SLShowerptr = &emSLHolder;
0850     for (unsigned int i = 0; i < SLShowerptr->SLEnergyBins.size(); i++) {
0851       if (!SLisEBinFilled(i)) {
0852         SLShowerptr = nullptr;
0853         return false;
0854       }
0855     }
0856   }
0857   if (DoHadSL) {
0858     SLShowerptr = &hadSLHolder;
0859     for (unsigned int i = 0; i < SLShowerptr->SLEnergyBins.size(); i++) {
0860       if (!SLisEBinFilled(i)) {
0861         SLShowerptr = nullptr;
0862         return false;
0863       }
0864     }
0865   }
0866   SLShowerptr = nullptr;
0867   return true;
0868 }
0869 void CastorShowerLibraryMaker::GetKinematics(
0870     int thePrim, double& px, double& py, double& pz, double& pInit, double& eta, double& phi) {
0871   if (thePrim == 0)
0872     return;
0873   if (PrimaryMomentum.find(thePrim) == PrimaryMomentum.end())
0874     return;
0875   px = PrimaryMomentum[thePrim].x() / GeV;
0876   py = PrimaryMomentum[thePrim].y() / GeV;
0877   pz = PrimaryMomentum[thePrim].z() / GeV;
0878   pInit = PrimaryMomentum[thePrim].mag() / GeV;
0879   if (pInit == 0)
0880     return;
0881   double costheta = pz / pInit;
0882   double theta = acos(std::min(std::max(costheta, double(-1.)), double(1.)));
0883   eta = -log(tan(theta / 2.0));
0884   phi = (px == 0 && py == 0) ? 0 : atan2(py, px);  // the recommended way of calculating phi
0885   phi = PrimaryMomentum[thePrim].phi();
0886 }
0887 void CastorShowerLibraryMaker::GetKinematics(
0888     G4PrimaryParticle* thePrim, double& px, double& py, double& pz, double& pInit, double& eta, double& phi) {
0889   px = py = pz = phi = eta = 0.0;
0890   if (thePrim == nullptr)
0891     return;
0892   px = thePrim->GetMomentum().x() / GeV;
0893   py = thePrim->GetMomentum().y() / GeV;
0894   pz = thePrim->GetMomentum().z() / GeV;
0895   pInit = thePrim->GetMomentum().mag() / GeV;
0896   //pInit = sqrt(pow(px,2.)+pow(py,2.)+pow(pz,2.));
0897   if (pInit == 0)
0898     return;
0899   double costheta = pz / pInit;
0900   double theta = acos(std::min(std::max(costheta, double(-1.)), double(1.)));
0901   eta = -log(tan(theta / 2.0));
0902   phi = (px == 0 && py == 0) ? 0 : atan2(py, px);  // the recommended way of calculating phi
0903   phi = thePrim->GetMomentum().phi();
0904   //if (px!=0) phi=atan(py/px);
0905 }
0906 std::vector<G4PrimaryParticle*> CastorShowerLibraryMaker::GetPrimary(const G4Event* evt) {
0907   // Find Primary info:
0908   int trackID = 0;
0909   std::vector<G4PrimaryParticle*> thePrims;
0910   G4PrimaryParticle* thePrim = nullptr;
0911   G4int nvertex = evt->GetNumberOfPrimaryVertex();
0912   edm::LogVerbatim("CastorShowerLibraryMaker") << "Event has " << nvertex << " vertex";
0913   if (nvertex != 1) {
0914     edm::LogVerbatim("CastorShowerLibraryMaker") << "CastorShowerLibraryMaker::GetPrimary ERROR: no vertex";
0915     return thePrims;
0916   }
0917 
0918   for (int i = 0; i < nvertex; i++) {
0919     G4PrimaryVertex* avertex = evt->GetPrimaryVertex(i);
0920     if (avertex == nullptr) {
0921       edm::LogVerbatim("CastorShowerLibraryMaker")
0922           << "CastorShowerLibraryMaker::GetPrimary ERROR: pointer to vertex = 0";
0923       continue;
0924     }
0925     unsigned int npart = avertex->GetNumberOfParticle();
0926     if (npart != NPGParticle)
0927       continue;
0928     for (unsigned int j = 0; j < npart; j++) {
0929       unsigned int k = 0;
0930       //int test_pID = 0;
0931       trackID = j;
0932       thePrim = avertex->GetPrimary(trackID);
0933       while (k < NPGParticle && PGParticleIDs.at(k++) != thePrim->GetPDGcode()) {
0934         ;
0935       };
0936       if (k > NPGParticle)
0937         continue;  // ID not found in the requested particles
0938       thePrims.push_back(thePrim);
0939     }
0940   }
0941   return thePrims;
0942 }
0943 void CastorShowerLibraryMaker::printSLstatus(int ebin, int etabin, int phibin) {
0944   if (!SLShowerptr) {
0945     edm::LogVerbatim("CastorShowerLibraryInfo") << "nullptr shower pointer. Printing both";
0946     edm::LogVerbatim("HcalSim") << "Electromagnetic";
0947     SLShowerptr = &emSLHolder;
0948     this->printSLstatus(ebin, etabin, phibin);
0949     edm::LogVerbatim("HcalSim") << "Hadronic";
0950     SLShowerptr = &hadSLHolder;
0951     this->printSLstatus(ebin, etabin, phibin);
0952     SLShowerptr = nullptr;
0953     return;
0954   }
0955   int nBinsE = SLShowerptr->SLEnergyBins.size();
0956   int nBinsEta = SLShowerptr->SLEtaBins.size();
0957   int nBinsPhi = SLShowerptr->SLPhiBins.size();
0958   std::vector<double> SLenergies = SLShowerptr->SLEnergyBins;
0959   std::ostringstream st1;
0960   for (int n = 0; n < 11 + (nBinsEta * nBinsPhi); n++)
0961     st1 << "=";
0962   edm::LogVerbatim("HcalSim") << st1.str();
0963   for (int i = 0; i < nBinsE; i++) {
0964     std::ostringstream st1;
0965     st1 << "E bin " << std::setw(6) << SLenergies.at(i) << " : ";
0966     for (int j = 0; j < nBinsEta; j++) {
0967       for (int k = 0; k < nBinsPhi; k++) {
0968         (SLisPhiBinFilled(i, j, k)) ? st1 << "1" : st1 << "-";
0969       }
0970       if (j < nBinsEta - 1)
0971         st1 << "|";
0972     }
0973     st1 << " (" << SLnEvtInBinE(i) << " events)";
0974     edm::LogVerbatim("HcalSim") << st1.str();
0975     if (ebin != i)
0976       continue;
0977     std::ostringstream st2;
0978     st2 << "               ";
0979     for (int j = 0; j < nBinsEta; j++) {
0980       for (int k = 0; k < nBinsPhi; k++) {
0981         (ebin == i && etabin == j && phibin == k) ? st2 << "^" : st2 << " ";
0982       }
0983       if (j < nBinsEta - 1)
0984         st2 << " ";
0985     }
0986     edm::LogVerbatim("HcalSim") << st2.str();
0987   }
0988   std::ostringstream st2;
0989   for (int n = 0; n < 11 + (nBinsEta * nBinsPhi); n++)
0990     st2 << "=";
0991   edm::LogVerbatim("HcalSim") << st2.str();
0992 }
0993 bool CastorShowerLibraryMaker::SLacceptEvent(int ebin, int etabin, int phibin) {
0994   if (SLShowerptr == nullptr) {
0995     edm::LogVerbatim("CastorShowerLibraryMaker::SLacceptEvent:") << "Error. nullptr pointer to CastorShowerEvent";
0996     return false;
0997   }
0998   if (ebin < 0 || ebin >= int(SLShowerptr->SLEnergyBins.size()))
0999     return false;
1000   if (SLisEBinFilled(ebin))
1001     return false;
1002 
1003   if (etabin < 0 || etabin >= int(SLShowerptr->SLEtaBins.size()))
1004     return false;
1005   if (SLisEtaBinFilled(ebin, etabin))
1006     return false;
1007 
1008   if (phibin < 0 || phibin >= int(SLShowerptr->SLPhiBins.size()))
1009     return false;
1010   if (SLisPhiBinFilled(ebin, etabin, phibin))
1011     return false;
1012   return true;
1013 }
1014 bool CastorShowerLibraryMaker::FillShowerEvent(CaloG4HitCollection* theCAFI, CastorShowerEvent* shower, int ipart) {
1015   unsigned int volumeID = 0;
1016   double en_in_fi = 0.;
1017   //double totalEnergy = 0;
1018 
1019   int nentries = theCAFI->entries();
1020 
1021   // Compute Total Energy in CastorFI volume
1022   /*
1023      for(int ihit = 0; ihit < nentries; ihit++) {
1024        CaloG4Hit* aHit = (*theCAFI)[ihit];
1025        totalEnergy += aHit->getEnergyDeposit();
1026      }
1027 */
1028   if (!shower) {
1029     edm::LogVerbatim("CastorShowerLibraryMaker") << "Error. nullptr pointer to CastorShowerEvent";
1030     return false;
1031   }
1032 
1033   CastorNumberingScheme* theCastorNumScheme = new CastorNumberingScheme();
1034   // Hit position
1035   math::XYZPoint entry;
1036   math::XYZPoint position;
1037   int nHits;
1038   nHits = 0;
1039   for (int ihit = 0; ihit < nentries; ihit++) {
1040     CaloG4Hit* aHit = (*theCAFI)[ihit];
1041     int hit_particleID = aHit->getTrackID();
1042     if (MapOfSecondaries[ipart].find(hit_particleID) == MapOfSecondaries[ipart].end()) {
1043       if (verbosity)
1044         edm::LogVerbatim("CastorShowerLibraryMaker") << "Skipping hit from trackID " << hit_particleID;
1045       continue;
1046     }
1047     volumeID = aHit->getUnitID();
1048     double hitEnergy = aHit->getEnergyDeposit();
1049     en_in_fi += aHit->getEnergyDeposit();
1050     float time = aHit->getTimeSlice();
1051     int zside, sector, zmodule;
1052     theCastorNumScheme->unpackIndex(volumeID, zside, sector, zmodule);
1053     entry = aHit->getEntry();
1054     position = aHit->getPosition();
1055     if (verbosity)
1056       edm::LogVerbatim("CastorShowerLibraryMaker") << "\n side , sector , module = " << zside << " , " << sector
1057                                                    << " , " << zmodule << "\n nphotons = " << hitEnergy;
1058 
1059     if (verbosity)
1060       edm::LogVerbatim("CastorShowerLibraryMaker")
1061           << "\n packIndex = " << theCastorNumScheme->packIndex(zside, sector, zmodule);
1062 
1063     if (verbosity && time > 100.) {
1064       edm::LogVerbatim("CastorShowerLibraryMaker")
1065           << "\n nentries = " << nentries << "\n     time[" << ihit << "] = " << time << "\n  trackID[" << ihit
1066           << "] = " << aHit->getTrackID() << "\n volumeID[" << ihit << "] = " << volumeID << "\n nphotons[" << ihit
1067           << "] = " << hitEnergy << "\n side, sector, module  = " << zside << ", " << sector << ", " << zmodule
1068           << "\n packIndex " << theCastorNumScheme->packIndex(zside, sector, zmodule) << "\n X,Y,Z = " << entry.x()
1069           << "," << entry.y() << "," << entry.z();
1070     }
1071     if (verbosity)
1072       edm::LogVerbatim("CastorShowerLibraryMaker") << "\n    Incident Energy = " << aHit->getIncidentEnergy() << " \n";
1073 
1074     //  CaloG4Hit information
1075     shower->setDetID(volumeID);
1076     shower->setHitPosition(position);
1077     shower->setNphotons(hitEnergy);
1078     shower->setTime(time);
1079     nHits++;
1080   }
1081   // Write number of hits to CastorShowerEvent instance
1082   if (nHits == 0) {
1083     edm::LogVerbatim("CastorShowerLibraryMaker") << "No hits found for this track (trackID=" << ipart << ").";
1084     if (theCastorNumScheme)
1085       delete theCastorNumScheme;
1086     return false;
1087   }
1088   shower->setNhit(nHits);
1089 
1090   // update the event counters
1091   if (theCastorNumScheme)
1092     delete theCastorNumScheme;
1093   return true;
1094 }
1095 int& CastorShowerLibraryMaker::SLnEvtInBinE(int ebin) {
1096   if (!SLShowerptr) {
1097     edm::LogVerbatim("CastorShowerLibraryMaker") << "\n\nSLnEvtInBinE can be called only after BeginOfEvent\n\n";
1098     throw SimG4Exception("\n\nnullptr Pointer to the shower library.");
1099   }
1100   return SLShowerptr->nEvtInBinE.at(ebin);
1101 }
1102 
1103 int& CastorShowerLibraryMaker::SLnEvtInBinEta(int ebin, int etabin) {
1104   if (!SLShowerptr) {
1105     edm::LogVerbatim("CastorShowerLibraryMaker") << "\n\nSLnEvtInBinEta can be called only after BeginOfEvent\n\n";
1106     throw SimG4Exception("\n\nnullptr Pointer to the shower library.");
1107   }
1108   return SLShowerptr->nEvtInBinEta.at(ebin).at(etabin);
1109 }
1110 
1111 int& CastorShowerLibraryMaker::SLnEvtInBinPhi(int ebin, int etabin, int phibin) {
1112   if (!SLShowerptr) {
1113     edm::LogVerbatim("CastorShowerLibraryMaker") << "\n\nSLnEvtInBinPhi can be called only after BeginOfEvent\n\n";
1114     throw SimG4Exception("\n\nnullptr Pointer to the shower library.");
1115   }
1116   return SLShowerptr->nEvtInBinPhi.at(ebin).at(etabin).at(phibin);
1117 }
1118 bool CastorShowerLibraryMaker::SLisEBinFilled(int ebin) {
1119   if (!SLShowerptr) {
1120     edm::LogVerbatim("CastorShowerLibraryMaker") << "\n\nSLisEBinFilled can be called only after BeginOfEvent\n\n";
1121     throw SimG4Exception("\n\nnullptr Pointer to the shower library.");
1122   }
1123   if (SLShowerptr->nEvtInBinE.at(ebin) < (int)SLShowerptr->nEvtPerBinE)
1124     return false;
1125   return true;
1126 }
1127 bool CastorShowerLibraryMaker::SLisEtaBinFilled(int ebin, int etabin) {
1128   if (!SLShowerptr) {
1129     edm::LogVerbatim("CastorShowerLibraryMaker") << "\n\nSLisEtaBinFilled can be called only after BeginOfEvent\n\n";
1130     throw SimG4Exception("\n\nnullptr Pointer to the shower library.");
1131   }
1132   if (SLShowerptr->nEvtInBinEta.at(ebin).at(etabin) < (int)SLShowerptr->nEvtPerBinEta)
1133     return false;
1134   return true;
1135 }
1136 bool CastorShowerLibraryMaker::SLisPhiBinFilled(int ebin, int etabin, int phibin) {
1137   if (!SLShowerptr) {
1138     edm::LogVerbatim("CastorShowerLibraryMaker") << "\n\nSLisPhiBinFilled can be called only after BeginOfEvent\n\n";
1139     throw SimG4Exception("\n\nnullptr Pointer to the shower library.");
1140   }
1141   if (SLShowerptr->nEvtInBinPhi.at(ebin).at(etabin).at(phibin) < (int)SLShowerptr->nEvtPerBinPhi)
1142     return false;
1143   return true;
1144 }
1145 void CastorShowerLibraryMaker::KillSecondaries(const G4Step* aStep) {
1146   const G4TrackVector* p_sec = aStep->GetSecondary();
1147   for (int i = 0; i < int(p_sec->size()); i++) {
1148     edm::LogVerbatim("HcalSim") << "Killing track ID " << p_sec->at(i)->GetTrackID()
1149                                 << " and its secondaries Produced by Process "
1150                                 << p_sec->at(i)->GetCreatorProcess()->GetProcessName() << " in the volume "
1151                                 << aStep->GetTrack()->GetVolume()->GetName();
1152     p_sec->at(i)->SetTrackStatus(fKillTrackAndSecondaries);
1153   }
1154 }
1155 
1156 void CastorShowerLibraryMaker::GetMissingEnergy(CaloG4HitCollection* theCAFI, double& miss_energy, double& tot_energy) {
1157   // Get the total deposited energy and the one from hit not associated to a primary
1158   miss_energy = 0;
1159   tot_energy = 0;
1160   int nhits = theCAFI->entries();
1161   for (int ihit = 0; ihit < nhits; ihit++) {
1162     int id = (*theCAFI)[ihit]->getTrackID();
1163     tot_energy += (*theCAFI)[ihit]->getEnergyDeposit();
1164     int hit_prim = 0;
1165     for (unsigned int i = 0; i < thePrims.size(); i++) {
1166       int particleType = thePrims.at(i)->GetPDGcode();
1167       if (MapOfSecondaries[particleType].find(id) != MapOfSecondaries[particleType].end())
1168         hit_prim = particleType;
1169     }
1170     if (hit_prim == 0) {
1171       edm::LogVerbatim("HcalSim") << "Track ID " << id << " produced a hit which is not associated with a primary.";
1172       miss_energy += (*theCAFI)[ihit]->getEnergyDeposit();
1173     }
1174   }
1175 }
1176 
1177 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
1178 #include "FWCore/PluginManager/interface/ModuleDef.h"
1179 
1180 DEFINE_SIMWATCHER(CastorShowerLibraryMaker);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!