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