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