File indexing completed on 2024-05-10 02:21:24
0001 #include "SimG4Core/CheckSecondary/interface/CheckSecondary.h"
0002 #include "SimG4Core/CheckSecondary/interface/TreatSecondary.h"
0003 #include "SimG4Core/Physics/interface/G4ProcessTypeEnumerator.h"
0004
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/Framework/interface/EventSetup.h"
0007
0008 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0009 #include "SimG4Core/Notification/interface/BeginOfTrack.h"
0010 #include "SimG4Core/Notification/interface/EndOfEvent.h"
0011
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013
0014 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0015 #include <CLHEP/Units/SystemOfUnits.h>
0016 #include "G4HCofThisEvent.hh"
0017 #include "G4Step.hh"
0018 #include "G4Track.hh"
0019 #include "G4VProcess.hh"
0020
0021 #include <cmath>
0022 #include <iomanip>
0023 #include <iostream>
0024
0025 CheckSecondary::CheckSecondary(const edm::ParameterSet &p)
0026 : treatSecondary(nullptr),
0027 typeEnumerator(nullptr),
0028 nsec(nullptr),
0029 procids(nullptr),
0030 px(nullptr),
0031 py(nullptr),
0032 pz(nullptr),
0033 mass(nullptr),
0034 deltae(nullptr),
0035 procs(nullptr),
0036 file(nullptr),
0037 tree(nullptr) {
0038 edm::ParameterSet m_p = p.getParameter<edm::ParameterSet>("CheckSecondary");
0039 std::string saveFile = m_p.getUntrackedParameter<std::string>("SaveInFile", "None");
0040 treatSecondary = new TreatSecondary(m_p);
0041 typeEnumerator = new G4ProcessTypeEnumerator();
0042
0043 nsec = new std::vector<int>();
0044 px = new std::vector<double>();
0045 py = new std::vector<double>();
0046 pz = new std::vector<double>();
0047 mass = new std::vector<double>();
0048 deltae = new std::vector<double>();
0049 procids = new std::vector<int>();
0050 procs = new std::vector<std::string>();
0051
0052 if (saveFile != "None") {
0053 saveToTree = true;
0054 tree = bookTree(saveFile);
0055 edm::LogInfo("CheckSecondary") << "Instantiate CheckSecondary with first"
0056 << " hadronic interaction information"
0057 << " to be saved in file " << saveFile;
0058 } else {
0059 saveToTree = false;
0060 edm::LogInfo("CheckSecondary") << "Instantiate CheckSecondary with first"
0061 << " hadronic interaction information"
0062 << " not saved";
0063 }
0064 }
0065
0066 CheckSecondary::~CheckSecondary() {
0067 if (saveToTree)
0068 endTree();
0069 if (nsec)
0070 delete nsec;
0071 if (px)
0072 delete px;
0073 if (py)
0074 delete py;
0075 if (pz)
0076 delete pz;
0077 if (mass)
0078 delete mass;
0079 if (deltae)
0080 delete deltae;
0081 if (procs)
0082 delete procs;
0083 if (procids)
0084 delete procids;
0085 if (typeEnumerator)
0086 delete typeEnumerator;
0087 if (treatSecondary)
0088 delete treatSecondary;
0089 }
0090
0091 TTree *CheckSecondary::bookTree(std::string fileName) {
0092 file = new TFile(fileName.c_str(), "RECREATE");
0093 file->cd();
0094
0095 TTree *t1 = new TTree("T1", "Secondary Particle Information");
0096 t1->Branch("SecondaryPx", "std::vector<double>", &px);
0097 t1->Branch("SecondaryPy", "std::vector<double>", &py);
0098 t1->Branch("SecondaryPz", "std::vector<double>", &pz);
0099 t1->Branch("SecondaryMass", "std::vector<double>", &mass);
0100 t1->Branch("NumberSecondaries", "std::vector<int>", &nsec);
0101 t1->Branch("DeltaEinInteract", "std::vector<double>", &deltae);
0102 t1->Branch("ProcessID", "std::vector<int>", &procids);
0103 t1->Branch("ProcessNames", "std::vector<std::string>", &procs);
0104 return t1;
0105 }
0106
0107 void CheckSecondary::endTree() {
0108 edm::LogInfo("CheckSecondary") << "Save the Secondary Tree " << tree->GetName() << " (" << tree << ") in file "
0109 << file->GetName() << " (" << file << ")";
0110 file->cd();
0111 tree->Write();
0112 file->Close();
0113 delete file;
0114 }
0115
0116 void CheckSecondary::update(const BeginOfEvent *evt) {
0117 int iev = (*evt)()->GetEventID();
0118 LogDebug("CheckSecondary") << "CheckSecondary::=====> Begin of event = " << iev;
0119
0120 (*nsec).clear();
0121 (*procs).clear();
0122 (*procids).clear();
0123 (*deltae).clear();
0124 (*px).clear();
0125 (*py).clear();
0126 (*pz).clear();
0127 (*mass).clear();
0128 }
0129
0130 void CheckSecondary::update(const BeginOfTrack *trk) {
0131 const G4Track *thTk = (*trk)();
0132 treatSecondary->initTrack(thTk);
0133 if (thTk->GetParentID() <= 0)
0134 storeIt = true;
0135 else
0136 storeIt = false;
0137 nHad = 0;
0138 LogDebug("CheckSecondary") << "CheckSecondary:: Track " << thTk->GetTrackID() << " Parent " << thTk->GetParentID()
0139 << " Flag " << storeIt;
0140 }
0141
0142 void CheckSecondary::update(const G4Step *aStep) {
0143 std::string name;
0144 int procID;
0145 bool hadrInt;
0146 double deltaE;
0147 std::vector<int> charge;
0148 std::vector<math::XYZTLorentzVector> tracks = treatSecondary->tracks(aStep, name, procID, hadrInt, deltaE, charge);
0149 if (storeIt && hadrInt) {
0150 double pInit = (aStep->GetPreStepPoint()->GetMomentum()).mag();
0151 double pEnd = (aStep->GetPostStepPoint()->GetMomentum()).mag();
0152 nHad++;
0153 int sec = (int)(tracks.size());
0154 LogDebug("CheckSecondary") << "CheckSecondary:: Hadronic Interaction " << nHad << " of type " << name << " ID "
0155 << procID << " with " << sec << " secondaries "
0156 << " and Momentum (MeV/c) at start " << pInit << " and at end " << pEnd;
0157 (*nsec).push_back(sec);
0158 (*procs).push_back(name);
0159 (*procids).push_back(procID);
0160 (*deltae).push_back(deltaE);
0161 if (nHad == 1) {
0162 for (int i = 0; i < sec; i++) {
0163 double m = tracks[i].M();
0164 if (charge[i] < 0)
0165 m = -m;
0166 (*px).push_back(tracks[i].Px());
0167 (*py).push_back(tracks[i].Py());
0168 (*pz).push_back(tracks[i].Pz());
0169 (*mass).push_back(m);
0170 }
0171 }
0172 }
0173 }
0174
0175 void CheckSecondary::update(const EndOfEvent *evt) {
0176 LogDebug("CheckSecondary") << "CheckSecondary::EndofEvent =====> Event " << (*evt)()->GetEventID() << " with "
0177 << (*nsec).size() << " hadronic collisions with"
0178 << " secondaries produced in each step";
0179 for (unsigned int i = 0; i < (*nsec).size(); i++)
0180 LogDebug("CheckSecondary") << " " << (*nsec)[i] << " from " << (*procs)[i] << " ID " << (*procids)[i] << " ("
0181 << typeEnumerator->processG4Name((*procids)[i]) << ") deltaE = " << (*deltae)[i]
0182 << " MeV";
0183 LogDebug("CheckSecondary") << "And " << (*mass).size() << " secondaries "
0184 << "produced in the first interactions";
0185 for (unsigned int i = 0; i < (*mass).size(); i++)
0186 LogDebug("CheckSecondary") << "Secondary " << i << " (" << (*px)[i] << ", " << (*py)[i] << ", " << (*pz)[i] << ", "
0187 << (*mass)[i] << ")";
0188
0189 if (saveToTree)
0190 tree->Fill();
0191 }
