File indexing completed on 2024-05-10 02:21:24
0001 #include "SimG4Core/CustomPhysics/interface/RHStopTracer.h"
0002
0003 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0004 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0005 #include "SimG4Core/Notification/interface/BeginOfTrack.h"
0006 #include "SimG4Core/Notification/interface/EndOfTrack.h"
0007 #include "SimG4Core/Notification/interface/TrackInformation.h"
0008 #include "SimG4Core/Application/interface/TrackingAction.h"
0009
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012
0013 #include "G4Track.hh"
0014 #include "G4Run.hh"
0015 #include "G4Event.hh"
0016 #include <CLHEP/Units/SystemOfUnits.h>
0017 #include "G4ParticleTable.hh"
0018 #include "G4ParticleDefinition.hh"
0019
0020 RHStopTracer::RHStopTracer(edm::ParameterSet const& p) {
0021 edm::ParameterSet parameters = p.getParameter<edm::ParameterSet>("RHStopTracer");
0022 mStopRegular = parameters.getUntrackedParameter<bool>("stopRegularParticles", false);
0023 mTraceEnergy = parameters.getUntrackedParameter<double>("traceEnergy", 1.e20);
0024 mTraceParticleName = parameters.getParameter<std::string>("traceParticle");
0025 minPdgId = parameters.getUntrackedParameter<int>("minPdgId", 1000000);
0026 maxPdgId = parameters.getUntrackedParameter<int>("maxPdgId", 2000000);
0027 otherPdgId = parameters.getUntrackedParameter<int>("otherPdgId", 17);
0028 produces<std::vector<std::string> >("StoppedParticlesName");
0029 produces<std::vector<float> >("StoppedParticlesX");
0030 produces<std::vector<float> >("StoppedParticlesY");
0031 produces<std::vector<float> >("StoppedParticlesZ");
0032 produces<std::vector<float> >("StoppedParticlesTime");
0033 produces<std::vector<int> >("StoppedParticlesPdgId");
0034 produces<std::vector<float> >("StoppedParticlesMass");
0035 produces<std::vector<float> >("StoppedParticlesCharge");
0036
0037 mTraceEnergy *= CLHEP::GeV;
0038 rePartName = mTraceParticleName;
0039 setMT(true);
0040
0041 edm::LogVerbatim("SimG4CoreCustomPhysics")
0042 << "RHStopTracer::RHStopTracer " << mTraceParticleName << " Eth(GeV)= " << mTraceEnergy;
0043 }
0044
0045 void RHStopTracer::update(const BeginOfRun* fRun) {
0046 LogDebug("SimG4CoreCustomPhysics") << "RHStopTracer::update-> begin of the run " << (*fRun)()->GetRunID();
0047 }
0048
0049 void RHStopTracer::update(const BeginOfEvent* fEvent) {
0050 LogDebug("SimG4CoreCustomPhysics") << "RHStopTracer::update-> begin of the event " << (*fEvent)()->GetEventID();
0051 }
0052
0053 void RHStopTracer::update(const BeginOfTrack* fTrack) {
0054 const G4Track* track = (*fTrack)();
0055 const G4ParticleDefinition* part = track->GetDefinition();
0056 const std::string& stringPartName = part->GetParticleName();
0057 bool matched = false;
0058 int pdgid = std::abs(part->GetPDGEncoding());
0059 if ((pdgid > minPdgId && pdgid < maxPdgId) || pdgid == otherPdgId)
0060 matched = std::regex_match(stringPartName, rePartName);
0061 if (matched || track->GetKineticEnergy() > mTraceEnergy) {
0062 LogDebug("SimG4CoreCustomPhysics") << "RHStopTracer::update-> new track: ID/Name/pdgId/mass/charge/Parent: "
0063 << track->GetTrackID() << '/' << part->GetParticleName() << '/'
0064 << part->GetPDGEncoding() << '/' << part->GetPDGMass() / CLHEP::GeV << " GeV/"
0065 << part->GetPDGCharge() << '/' << track->GetParentID()
0066 << " Position: " << track->GetPosition() << ' '
0067 << " R/phi: " << track->GetPosition().perp() << '/' << track->GetPosition().phi()
0068 << " 4vec " << track->GetMomentum();
0069 } else if (mStopRegular) {
0070 const_cast<G4Track*>(track)->SetTrackStatus(fStopAndKill);
0071 }
0072 }
0073
0074 void RHStopTracer::update(const EndOfTrack* fTrack) {
0075 const G4Track* track = (*fTrack)();
0076 const G4ParticleDefinition* part = track->GetDefinition();
0077 const std::string& stringPartName = part->GetParticleName();
0078 bool matched = false;
0079 int pdgid = std::abs(part->GetPDGEncoding());
0080 if ((pdgid > minPdgId && pdgid < maxPdgId) || pdgid == otherPdgId)
0081 matched = std::regex_match(stringPartName, rePartName);
0082 if (matched || track->GetKineticEnergy() > mTraceEnergy) {
0083 LogDebug("SimG4CoreCustomPhysics") << "RHStopTracer::update-> stop track: ID/Name/pdgId/mass/charge/Parent: "
0084 << track->GetTrackID() << '/' << part->GetParticleName() << '/'
0085 << part->GetPDGEncoding() << '/' << part->GetPDGMass() / CLHEP::GeV << " GeV/"
0086 << part->GetPDGCharge() << '/' << track->GetParentID()
0087 << " Position: " << track->GetPosition() << ' '
0088 << " R/phi: " << track->GetPosition().perp() << '/' << track->GetPosition().phi()
0089 << " 4vec " << track->GetMomentum();
0090 if (track->GetMomentum().mag() < 0.001) {
0091 LogDebug("SimG4CoreCustomPhysics") << "RHStopTracer:: track has stopped, so making StopPoint";
0092 mStopPoints.push_back(StopPoint(track->GetDefinition()->GetParticleName(),
0093 track->GetPosition().x(),
0094 track->GetPosition().y(),
0095 track->GetPosition().z(),
0096 track->GetGlobalTime(),
0097 track->GetDefinition()->GetPDGEncoding(),
0098 track->GetDefinition()->GetPDGMass() / CLHEP::GeV,
0099 track->GetDefinition()->GetPDGCharge()));
0100 }
0101 }
0102 }
0103
0104 void RHStopTracer::produce(edm::Event& fEvent, const edm::EventSetup&) {
0105 LogDebug("SimG4CoreCustomPhysics") << "RHStopTracer::produce->";
0106
0107 std::unique_ptr<std::vector<std::string> > names(new std::vector<std::string>);
0108 std::unique_ptr<std::vector<float> > xs(new std::vector<float>);
0109 std::unique_ptr<std::vector<float> > ys(new std::vector<float>);
0110 std::unique_ptr<std::vector<float> > zs(new std::vector<float>);
0111 std::unique_ptr<std::vector<float> > ts(new std::vector<float>);
0112 std::unique_ptr<std::vector<int> > ids(new std::vector<int>);
0113 std::unique_ptr<std::vector<float> > masses(new std::vector<float>);
0114 std::unique_ptr<std::vector<float> > charges(new std::vector<float>);
0115
0116 std::vector<StopPoint>::const_iterator stopPoint = mStopPoints.begin();
0117 for (; stopPoint != mStopPoints.end(); ++stopPoint) {
0118 names->push_back(stopPoint->name);
0119 xs->push_back(stopPoint->x);
0120 ys->push_back(stopPoint->y);
0121 zs->push_back(stopPoint->z);
0122 ts->push_back(stopPoint->t);
0123 ids->push_back(stopPoint->id);
0124 masses->push_back(stopPoint->mass);
0125 charges->push_back(stopPoint->charge);
0126 }
0127 fEvent.put(std::move(names), "StoppedParticlesName");
0128 fEvent.put(std::move(xs), "StoppedParticlesX");
0129 fEvent.put(std::move(ys), "StoppedParticlesY");
0130 fEvent.put(std::move(zs), "StoppedParticlesZ");
0131 fEvent.put(std::move(ts), "StoppedParticlesTime");
0132 fEvent.put(std::move(ids), "StoppedParticlesPdgId");
0133 fEvent.put(std::move(masses), "StoppedParticlesMass");
0134 fEvent.put(std::move(charges), "StoppedParticlesCharge");
0135 mStopPoints.clear();
0136 }
