Back to home page

Project CMSSW displayed by LXR

 
 

    


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

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 "G4SystemOfUnits.hh"
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 
0040   edm::LogInfo("SimG4CoreCustomPhysics") << "RHStopTracer::RHStopTracer " << mTraceParticleName
0041                                          << " Eth(GeV)= " << mTraceEnergy;
0042 }
0043 
0044 RHStopTracer::~RHStopTracer() {}
0045 
0046 void RHStopTracer::update(const BeginOfRun* fRun) {
0047   LogDebug("SimG4CoreCustomPhysics") << "RHStopTracer::update-> begin of the run " << (*fRun)()->GetRunID();
0048 }
0049 
0050 void RHStopTracer::update(const BeginOfEvent* fEvent) {
0051   LogDebug("SimG4CoreCustomPhysics") << "RHStopTracer::update-> begin of the event " << (*fEvent)()->GetEventID();
0052 }
0053 
0054 void RHStopTracer::update(const BeginOfTrack* fTrack) {
0055   const G4Track* track = (*fTrack)();
0056   const G4ParticleDefinition* part = track->GetDefinition();
0057   const std::string& stringPartName = part->GetParticleName();
0058   bool matched = false;
0059   int pdgid = std::abs(part->GetPDGEncoding());
0060   if ((pdgid > minPdgId && pdgid < maxPdgId) || pdgid == otherPdgId)
0061     matched = std::regex_match(stringPartName, rePartName);
0062   if (matched || track->GetKineticEnergy() > mTraceEnergy) {
0063     LogDebug("SimG4CoreCustomPhysics") << "RHStopTracer::update-> new track: ID/Name/pdgId/mass/charge/Parent: "
0064                                        << track->GetTrackID() << '/' << part->GetParticleName() << '/'
0065                                        << part->GetPDGEncoding() << '/' << part->GetPDGMass() / GeV << " GeV/"
0066                                        << part->GetPDGCharge() << '/' << track->GetParentID()
0067                                        << " Position: " << track->GetPosition() << ' '
0068                                        << " R/phi: " << track->GetPosition().perp() << '/' << track->GetPosition().phi()
0069                                        << "   4vec " << track->GetMomentum();
0070   } else if (mStopRegular) {  // kill regular particles
0071     const_cast<G4Track*>(track)->SetTrackStatus(fStopAndKill);
0072   }
0073 }
0074 
0075 void RHStopTracer::update(const EndOfTrack* fTrack) {
0076   const G4Track* track = (*fTrack)();
0077   const G4ParticleDefinition* part = track->GetDefinition();
0078   const std::string& stringPartName = part->GetParticleName();
0079   bool matched = false;
0080   int pdgid = std::abs(part->GetPDGEncoding());
0081   if ((pdgid > minPdgId && pdgid < maxPdgId) || pdgid == otherPdgId)
0082     matched = std::regex_match(stringPartName, rePartName);
0083   if (matched || track->GetKineticEnergy() > mTraceEnergy) {
0084     LogDebug("SimG4CoreCustomPhysics") << "RHStopTracer::update-> stop track: ID/Name/pdgId/mass/charge/Parent: "
0085                                        << track->GetTrackID() << '/' << part->GetParticleName() << '/'
0086                                        << part->GetPDGEncoding() << '/' << part->GetPDGMass() / GeV << " GeV/"
0087                                        << part->GetPDGCharge() << '/' << track->GetParentID()
0088                                        << " Position: " << track->GetPosition() << ' '
0089                                        << " R/phi: " << track->GetPosition().perp() << '/' << track->GetPosition().phi()
0090                                        << "   4vec " << track->GetMomentum();
0091     if (track->GetMomentum().mag() < 0.001) {
0092       LogDebug("SimG4CoreCustomPhysics") << "RHStopTracer:: track has stopped, so making StopPoint";
0093       mStopPoints.push_back(StopPoint(track->GetDefinition()->GetParticleName(),
0094                                       track->GetPosition().x(),
0095                                       track->GetPosition().y(),
0096                                       track->GetPosition().z(),
0097                                       track->GetGlobalTime(),
0098                                       track->GetDefinition()->GetPDGEncoding(),
0099                                       track->GetDefinition()->GetPDGMass() / GeV,
0100                                       track->GetDefinition()->GetPDGCharge()));
0101     }
0102   }
0103 }
0104 
0105 void RHStopTracer::produce(edm::Event& fEvent, const edm::EventSetup&) {
0106   LogDebug("SimG4CoreCustomPhysics") << "RHStopTracer::produce->";
0107 
0108   std::unique_ptr<std::vector<std::string> > names(new std::vector<std::string>);
0109   std::unique_ptr<std::vector<float> > xs(new std::vector<float>);
0110   std::unique_ptr<std::vector<float> > ys(new std::vector<float>);
0111   std::unique_ptr<std::vector<float> > zs(new std::vector<float>);
0112   std::unique_ptr<std::vector<float> > ts(new std::vector<float>);
0113   std::unique_ptr<std::vector<int> > ids(new std::vector<int>);
0114   std::unique_ptr<std::vector<float> > masses(new std::vector<float>);
0115   std::unique_ptr<std::vector<float> > charges(new std::vector<float>);
0116 
0117   std::vector<StopPoint>::const_iterator stopPoint = mStopPoints.begin();
0118   for (; stopPoint != mStopPoints.end(); ++stopPoint) {
0119     names->push_back(stopPoint->name);
0120     xs->push_back(stopPoint->x);
0121     ys->push_back(stopPoint->y);
0122     zs->push_back(stopPoint->z);
0123     ts->push_back(stopPoint->t);
0124     ids->push_back(stopPoint->id);
0125     masses->push_back(stopPoint->mass);
0126     charges->push_back(stopPoint->charge);
0127   }
0128   fEvent.put(std::move(names), "StoppedParticlesName");
0129   fEvent.put(std::move(xs), "StoppedParticlesX");
0130   fEvent.put(std::move(ys), "StoppedParticlesY");
0131   fEvent.put(std::move(zs), "StoppedParticlesZ");
0132   fEvent.put(std::move(ts), "StoppedParticlesTime");
0133   fEvent.put(std::move(ids), "StoppedParticlesPdgId");
0134   fEvent.put(std::move(masses), "StoppedParticlesMass");
0135   fEvent.put(std::move(charges), "StoppedParticlesCharge");
0136   mStopPoints.clear();
0137 }