Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:29:35

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/GlobalSystemOfUnits.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 }