Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:24:40 

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