File indexing completed on 2023-03-17 11:24:40
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
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
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;
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
0101 struct ShowerLib {
0102 CastorShowerLibraryInfo SLInfo;
0103 SLBin3D SLCollection;
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
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
0127 int verbosity;
0128 std::string eventNtFileName;
0129
0130 unsigned int NPGParticle;
0131 std::vector<int> PGParticleIDs;
0132 bool DoHadSL;
0133 bool DoEmSL;
0134 bool InsideCastor;
0135 bool DeActivatePhysicsProcess;
0136 std::vector<G4PrimaryParticle*> thePrims;
0137
0138
0139 CastorShowerLibraryInfo* emInfo;
0140 CastorShowerLibraryInfo* hadInfo;
0141 CastorShowerEvent* emShower;
0142 CastorShowerEvent* hadShower;
0143 ShowerLib emSLHolder;
0144 ShowerLib hadSLHolder;
0145 ShowerLib* SLShowerptr;
0146 std::map<int, std::set<int> > MapOfSecondaries;
0147
0148
0149 std::map<int, G4ThreeVector> PrimaryMomentum;
0150 std::map<int, G4ThreeVector> PrimaryPosition;
0151 double MaxEta;
0152 double MaxPhi;
0153
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
0177 TFile* theFile;
0178 TTree* theTree;
0179
0180 int eventIndex;
0181 int stepIndex;
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
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
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
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
0293 void CastorShowerLibraryMaker::update(const BeginOfJob* job) {
0294 edm::LogVerbatim("HcalSim") << " CastorShowerLibraryMaker::Starting new job ";
0295 }
0296
0297
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
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
0320
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
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
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
0336
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
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
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
0352
0353
0354
0355
0356
0357
0358
0359
0360
0361
0362
0363
0364
0365
0366 eventIndex = 0;
0367 }
0368
0369
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
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
0419
0420
0421
0422
0423
0424
0425
0426
0427
0428
0429
0430
0431
0432
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
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
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
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
0497 std::string CurVolume = trk->GetVolume()->GetName();
0498 if (!InsideCastor && (
0499
0500
0501
0502 CurVolume == "CAST")) {
0503
0504 InsideCastor = true;
0505
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
0519
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
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
0550
0551
0552
0553
0554
0555
0556 }
0557
0558
0559 void CastorShowerLibraryMaker::update(const EndOfEvent* evt) {
0560
0561 if ((*evt)()->IsAborted()) {
0562 edm::LogVerbatim("HcalSim") << "\n========================================================================"
0563 << "\nEndOfEvent: EVENT ABORTED"
0564 << "\n========================================================================";
0565 return;
0566 }
0567
0568
0569
0570
0571
0572 if (verbosity)
0573 edm::LogVerbatim("HcalSim") << "CastorShowerLibraryMaker: End of Event: " << eventIndex;
0574
0575 if (thePrims.empty()) {
0576 edm::LogVerbatim("CastorShowerLibraryMaker") << "No valid primary particle found. Skipping event";
0577 return;
0578 }
0579
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
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
0601 int particleType = thePrim->GetPDGcode();
0602
0603
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
0615 double px = 0., py = 0., pz = 0., pInit = 0., eta = 0., phi = 0.;
0616 GetKinematics(particleType, px, py, pz, pInit, eta, phi);
0617
0618
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
0636
0637
0638 edm::LogVerbatim("CastorShowerLibraryMaker")
0639 << "\n CastorShowerLibraryMaker::update(EndOfEvent * evt) - event #" << (*evt)()->GetEventID();
0640
0641
0642
0643
0644
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
0652 if (FillShowerEvent(theCAFI, shower, particleType)) {
0653
0654
0655
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
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
0688
0689
0690
0691
0692
0693
0694
0695
0696
0697
0698
0699
0700
0701
0702
0703
0704 }
0705
0706
0707 void CastorShowerLibraryMaker::update(const EndOfRun* run) {
0708
0709 if (!IsSLReady())
0710 SimG4Exception("\n\nShower Library NOT READY.\n\n");
0711
0712 unsigned int ibine, ibineta, ibinphi, ievt;
0713 unsigned int jbine, jbineta, jbinphi, jevt;
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
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
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
0787
0788
0789
0790
0791 }
0792 int CastorShowerLibraryMaker::FindEnergyBin(double energy) {
0793
0794
0795
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
0811 if (energy >= SLenergies.at(i))
0812 return (int)i;
0813
0814 return -1;
0815 }
0816 int CastorShowerLibraryMaker::FindEtaBin(double eta) {
0817
0818
0819
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
0833 if (eta >= SLetas.at(i))
0834 return (int)i;
0835
0836 return -1;
0837 }
0838 int CastorShowerLibraryMaker::FindPhiBin(double phi) {
0839
0840
0841
0842
0843
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
0857 if (phi >= SLphis.at(i))
0858 return (int)i;
0859
0860 return -1;
0861 }
0862 bool CastorShowerLibraryMaker::IsSLReady() {
0863
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
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);
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
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);
0924 phi = thePrim->GetMomentum().phi();
0925
0926 }
0927 std::vector<G4PrimaryParticle*> CastorShowerLibraryMaker::GetPrimary(const G4Event* evt) {
0928
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
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;
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
1039
1040 int nentries = theCAFI->entries();
1041
1042
1043
1044
1045
1046
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
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
1096 shower->setDetID(volumeID);
1097 shower->setHitPosition(position);
1098 shower->setNphotons(hitEnergy);
1099 shower->setTime(time);
1100 nHits++;
1101 }
1102
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
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
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);