CheckSecondary/src/TreatSecondary.cc

0001 #include "SimG4Core/CheckSecondary/interface/TreatSecondary.h"
0002 #include "SimG4Core/Physics/interface/G4ProcessTypeEnumerator.h"
0003
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005
0006 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0007 #include <CLHEP/Units/SystemOfUnits.h>
0008 #include "G4HCofThisEvent.hh"
0009 #include "G4Step.hh"
0010 #include "G4Track.hh"
0011 #include "G4VProcess.hh"
0012
0013 #include <cmath>
0014 #include <iomanip>
0015 #include <iostream>
0016
0017 using CLHEP::GeV;
0018 using CLHEP::MeV;
0019
0020 TreatSecondary::TreatSecondary(const edm::ParameterSet &p) : typeEnumerator(nullptr) {
0021   verbosity = p.getUntrackedParameter<int>("Verbosity", 0);
0022   killAfter = p.getUntrackedParameter<int>("KillAfter", -1);
0023   minDeltaE = p.getUntrackedParameter<double>("MinimumDeltaE", 10.0) * MeV;
0024
0025   edm::LogInfo("CheckSecondary") << "Instantiate CheckSecondary with Flag"
0026                                  << " for Killing track after " << killAfter
0027                                  << " hadronic interactions\nDefine inelastic"
0028                                  << " if > 1 seondary or change in KE > " << minDeltaE << " MeV\n";
0029
0030   typeEnumerator = new G4ProcessTypeEnumerator();
0031 }
0032
0033 TreatSecondary::~TreatSecondary() {
0034   if (typeEnumerator)
0035     delete typeEnumerator;
0036 }
0037
0038 void TreatSecondary::initTrack(const G4Track *thTk) {
0039   step = 0;
0040   nsecL = 0;
0041   nHad = 0;
0042   eTrack = thTk->GetKineticEnergy() / MeV;
0043   LogDebug("CheckSecondary") << "TreatSecondary::initTrack:Track: " << thTk->GetTrackID()
0044                              << " Type: " << thTk->GetDefinition()->GetParticleName() << " KE "
0045                              << thTk->GetKineticEnergy() / GeV << " GeV p " << thTk->GetMomentum().mag() / GeV
0046                              << " GeV daughter of particle " << thTk->GetParentID();
0047 }
0048
0049 std::vector<math::XYZTLorentzVector> TreatSecondary::tracks(
0050     const G4Step *aStep, std::string &name, int &procid, bool &hadrInt, double &deltaE, std::vector<int> &charges) {
0051   step++;
0052   procid = -1;
0053   name = "Unknown";
0054   hadrInt = false;
0055   deltaE = 0;
0056   std::vector<math::XYZTLorentzVector> secondaries;
0057   charges.clear();
0058
0059   if (aStep != nullptr) {
0060     const G4TrackVector *tkV = aStep->GetSecondary();
0061     G4Track *thTk = aStep->GetTrack();
0062     const G4StepPoint *preStepPoint = aStep->GetPreStepPoint();
0063     const G4StepPoint *postStepPoint = aStep->GetPostStepPoint();
0064     double eTrackNew = thTk->GetKineticEnergy() / MeV;
0065     deltaE = eTrack - eTrackNew;
0066     eTrack = eTrackNew;
0067     if (tkV != nullptr && postStepPoint != nullptr) {
0068       int nsc = (*tkV).size();
0069       const G4VProcess *proc = postStepPoint->GetProcessDefinedStep();
0070       if (proc != nullptr) {
0071         G4ProcessType type = proc->GetProcessType();
0072         procid = typeEnumerator->processIdLong(proc);
0073         name = proc->GetProcessName();
0074         int sec = nsc - nsecL;
0075         LogDebug("CheckSecondary") << sec << " secondaries in step " << thTk->GetCurrentStepNumber() << " of track "
0076                                    << thTk->GetTrackID() << " from " << name << " of type " << type << " ID " << procid
0077                                    << " (" << typeEnumerator->processG4Name(procid) << ")";
0078
0079         // hadronic interaction
0080         if (procid >= 121 && procid <= 151) {
0081           LogDebug("CheckSecondary") << "Hadronic Interaction " << nHad << " of Type " << procid << " with " << sec
0082                                      << " secondaries from process " << proc->GetProcessName() << " Delta E " << deltaE
0083                                      << " Flag " << hadrInt;
0084           math::XYZTLorentzVector secondary;
0085           for (int i = nsecL; i < nsc; ++i) {
0086             G4Track *tk = (*tkV)[i];
0087             G4ThreeVector pp = tk->GetMomentum();
0088             double ee = tk->GetTotalEnergy();
0089             secondary = math::XYZTLorentzVector(pp.x(), pp.y(), pp.z(), ee);
0090             secondaries.push_back(secondary);
0091             int charge = (int)(tk->GetDefinition()->GetPDGCharge());
0092             charges.push_back(charge);
0093           }
0094           if (verbosity > 0) {
0095             for (int i = nsecL; i < nsc; i++) {
0096               G4Track *tk = (*tkV)[i];
0097               LogDebug("CheckSecondary") << "Secondary: " << sec << " ID " << tk->GetTrackID() << " Status "
0098                                          << tk->GetTrackStatus() << " Particle "
0099                                          << tk->GetDefinition()->GetParticleName() << " Position " << tk->GetPosition()
0100                                          << " KE " << tk->GetKineticEnergy() << " Time " << tk->GetGlobalTime();
0101             }
0102           }
0103         }
0104         nsecL = nsc;
0105       }
0106     }
0107     if (verbosity > 1) {
0108       LogDebug("CheckSecondary") << "Track: " << thTk->GetTrackID() << " Status " << thTk->GetTrackStatus()
0109                                  << " Particle " << thTk->GetDefinition()->GetParticleName() << " at "
0110                                  << preStepPoint->GetPosition() << " Step: " << step << " KE "
0111                                  << thTk->GetKineticEnergy() / GeV << " GeV; p " << thTk->GetMomentum().mag() / GeV
0112                                  << " GeV/c; Step Length " << aStep->GetStepLength() << " Energy Deposit "
0113                                  << aStep->GetTotalEnergyDeposit() / MeV << " MeV; Interaction " << hadrInt;
0114     }
0115   }
0116   return secondaries;
0117 }