LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
larg4::MCTruthEventActionService Class Reference

#include "MCTruthEventAction_service.h"

Inheritance diagram for larg4::MCTruthEventActionService:
artg4tk::PrimaryGeneratorActionBase artg4tk::ActionBase

Public Member Functions

 MCTruthEventActionService (fhicl::ParameterSet const &)
 
 ~MCTruthEventActionService ()
 
void setInputCollections (std::vector< art::Handle< std::vector< simb::MCTruth >>> const &mclists)
 
std::string const & myName () const
 
virtual void initialize ()
 

Private Member Functions

void generatePrimaries (G4Event *anEvent) override
 

Private Attributes

std::vector< art::Handle< std::vector< simb::MCTruth > > > const * fMCLists
 MCTruthCollection input lists. More...
 
std::map< G4int, G4int > fUnknownPDG
 map of unknown PDG codes to instances More...
 
std::map< G4int, G4int > fNon1StatusPDG
 PDG codes skipped because not status 1. More...
 
std::map< G4int, G4int > fProcessedPDG
 PDG codes processed. More...
 

Static Private Attributes

static G4ParticleTable * fParticleTable = nullptr
 Geant4's table of particle definitions. More...
 

Detailed Description

Definition at line 40 of file MCTruthEventAction_service.h.

Constructor & Destructor Documentation

larg4::MCTruthEventActionService::MCTruthEventActionService ( fhicl::ParameterSet const &  p)

Definition at line 42 of file MCTruthEventAction.cc.

43  : PrimaryGeneratorActionBase(p.get<string>("name", "MCTruthEventActionService"))
44 {}
larg4::MCTruthEventActionService::~MCTruthEventActionService ( )

Definition at line 46 of file MCTruthEventAction.cc.

References fNon1StatusPDG, fProcessedPDG, fUnknownPDG, and n.

47 {
48  // Print out a list of "unknown" PDG codes we saw in the input.
49  if (!fUnknownPDG.empty()) {
50  std::ostringstream badtxt;
51  for (auto const [pdg, n] : fUnknownPDG) {
52  badtxt << "\n Unknown PDG code = " << pdg << ", seen " << n << " times.";
53  if (genie_specific(pdg)) { badtxt << " (GENIE specific)"; }
54  }
55 
56  mf::LogWarning("MCTruthEventAction") << "The following unknown PDG codes were present in the "
57  << "simb::MCTruth input.\n"
58  << "They were not processed by Geant4." << badtxt.str();
59  }
60 
61  // Print out a list of PDG codes we saw in the input that weren't status=1
62  if (!fNon1StatusPDG.empty()) {
63  std::ostringstream non1txt;
64  for (auto const [pdg, n] : fNon1StatusPDG) {
65  non1txt << "\n PDG code = " << pdg << ", seen " << n << " times.";
66  if (genie_specific(pdg)) { non1txt << " (GENIE specific)"; }
67  }
68 
69  mf::LogDebug("MCTruthEventAction") << "The following PDG codes were present in the "
70  << "simb::MCTruth input with Status != 1.\n"
71  << "They were not processed by Geant4." << non1txt.str();
72  }
73 
74  // Print out a list of PDG codes that were processed
75  if (!fProcessedPDG.empty()) {
76  std::ostringstream goodtxt;
77  for (auto const [pdg, n] : fProcessedPDG) {
78  goodtxt << "\n PDG code = " << pdg << ", seen " << n << " times.";
79  }
80  mf::LogDebug("MCTruthEventAction") << "The following PDG codes were present in the "
81  << "simb::MCTruth input with Status = 1 "
82  << "and were processed by Geant4" << goodtxt.str();
83  }
84 }
std::map< G4int, G4int > fProcessedPDG
PDG codes processed.
std::map< G4int, G4int > fNon1StatusPDG
PDG codes skipped because not status 1.
std::map< G4int, G4int > fUnknownPDG
map of unknown PDG codes to instances
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Char_t n[5]

Member Function Documentation

void larg4::MCTruthEventActionService::generatePrimaries ( G4Event *  anEvent)
overrideprivatevirtual

Reimplemented from artg4tk::PrimaryGeneratorActionBase.

Definition at line 88 of file MCTruthEventAction.cc.

References fMCLists, fNon1StatusPDG, fParticleTable, fProcessedPDG, fUnknownPDG, art::Ptr< T >::get(), simb::MCTruth::GetParticle(), MF_LOG_DEBUG, MF_LOG_INFO, simb::MCParticle::Momentum(), simb::MCTruth::NParticles(), simb::MCParticle::PdgCode(), simb::MCParticle::Polarization(), g4b::PrimaryParticleInformation::SetMCTruth(), simb::MCParticle::StatusCode(), simb::MCParticle::T(), simb::MCParticle::TrackId(), simb::MCParticle::Vx(), simb::MCParticle::Vy(), simb::MCParticle::Vz(), x, y, Z, and z.

89 {
90  // For each MCTruth (probably only one, but you never know):
91  // index keeps track of which MCTruth object you are using
92  auto const& mclistHandles = *fMCLists;
93 
94  size_t mclSize = mclistHandles.size(); // -- should match the number of generators
95  mf::LogDebug("generatePrimaries") << "MCTruth Handles Size: " << mclSize;
96 
97  // -- Loop over MCTruth Handle List
98  size_t index = 0;
99  std::map<CLHEP::HepLorentzVector, G4PrimaryVertex*> vertexMap;
100  for (size_t mcl = 0; mcl < mclSize; ++mcl) {
101  mf::LogDebug("generatePrimaries")
102  << "MCTruth Handle Number: " << (mcl + 1) << " of " << mclSize;
103  art::Handle<std::vector<simb::MCTruth>> mclistHandle = mclistHandles[mcl];
104  // Loop over all MCTruth handle entries for a given generator,
105  // usually only one, but you never know
106  for (size_t i = 0; i < mclistHandle->size(); ++i) {
107  art::Ptr<simb::MCTruth> mclist(mclistHandle, i);
108 
109  mf::LogDebug("generatePrimaries")
110  << "art::Ptr number " << (i + 1) << " of " << mclistHandle->size() << ", Ptr: " << mclist;
111  int nPart = mclist->NParticles();
112  MF_LOG_INFO("generatePrimaries") << "Generating " << nPart << " particles";
113 
114  // -- Loop over all particles in MCTruth Object
115  for (int m = 0; m != nPart; ++m) {
116  simb::MCParticle const& particle = mclist->GetParticle(m);
117  G4int const pdgCode = particle.PdgCode();
118 
119  if (particle.StatusCode() != 1) {
120  MF_LOG_DEBUG("generatePrimaries")
121  << "Status code != 1, skipping particle number with MCTruth index = "
122  << "(" << mcl << "," << i << ")"
123  << " and particle index = " << m;
124  fNon1StatusPDG[pdgCode] += 1;
125  continue;
126  }
127 
128  if (((m + 1) % nPart) < 2) // -- only first and last will satisfy this
129  {
130  mf::LogDebug("generatePrimaries") << "Particle Number: " << (m + 1) << " of " << nPart;
131  mf::LogDebug("generatePrimaries") << "TrackID: " << particle.TrackId();
132  mf::LogDebug("generatePrimaries") << "index: " << index;
133  }
134 
135  G4double x = particle.Vx() * CLHEP::cm;
136  G4double y = particle.Vy() * CLHEP::cm;
137  G4double z = particle.Vz() * CLHEP::cm;
138  G4double t = particle.T() * CLHEP::ns;
139 
140  // Create a CLHEP four-vector from the particle's vertex.
141  CLHEP::HepLorentzVector fourpos(x, y, z, t);
142 
143  // Is this vertex already in our map?
144  G4PrimaryVertex* vertex = nullptr;
145  auto result = vertexMap.find(fourpos);
146  if (result == vertexMap.end()) {
147  // No, it's not, so create a new vertex and add it to the map.
148  vertex = new G4PrimaryVertex(x, y, z, t);
149  vertexMap[fourpos] = vertex;
150 
151  // Add the vertex to the G4Event.
152  anEvent->AddPrimaryVertex(vertex);
153  }
154  else {
155  // Yes, it is, so use the existing vertex.
156  vertex = result->second;
157  }
158 
159  // Get additional particle information.
160  TLorentzVector const& momentum = particle.Momentum(); // (px,py,pz,E)
161  TVector3 const& polarization = particle.Polarization();
162 
163  // Get the particle table if necessary. (Note: we're doing
164  // this "late" because I'm not sure at what point the G4
165  // particle table is initialized in the loading process.)
166  if (fParticleTable == nullptr) { fParticleTable = G4ParticleTable::GetParticleTable(); }
167 
168  // Get Geant4's definition of the particle.
169  G4ParticleDefinition* particleDefinition = nullptr;
170 
171  if (pdgCode == 0) { particleDefinition = fParticleTable->FindParticle("opticalphoton"); }
172  else {
173  particleDefinition = fParticleTable->FindParticle(pdgCode);
174  }
175 
176  if (pdgCode >= 2'000'000'000) { // If the particle is generator-specific
177  mf::LogDebug("ConvertPrimaryToGeant4")
178  << ": %%% Will skip particle with generator-specific PDG code = " << pdgCode;
179  fUnknownPDG[pdgCode] += 1;
180  continue;
181  }
182  if (pdgCode > 1'000'000'000) { // If the particle is a nucleus
183  mf::LogDebug("ConvertPrimaryToGeant4")
184  << ": %%% Nuclear PDG code = " << pdgCode << " (x,y,z,t)=(" << x << "," << y << "," << z
185  << "," << t << ")"
186  << " P=" << momentum.P() << ", E=" << momentum.E();
187  // If the particle table doesn't have a definition yet, ask
188  // the ion table for one. This will create a new ion
189  // definition as needed.
190  if (!particleDefinition) {
191  int const Z = (pdgCode % 10'000'000) / 10'000; // atomic number
192  int const A = (pdgCode % 10'000) / 10; // mass number
193  particleDefinition = fParticleTable->GetIonTable()->GetIon(Z, A, 0.);
194  }
195  }
196 
197  // What if the PDG code is unknown? This has been a known
198  // issue with GENIE.
199  if (particleDefinition == nullptr) {
200  mf::LogDebug("ConvertPrimaryToGeant4") << ": %%% Code not found = " << pdgCode;
201  fUnknownPDG[pdgCode] += 1;
202  continue;
203  }
204 
205  fProcessedPDG[pdgCode] += 1;
206 
207  // Create a Geant4 particle to add to the vertex.
208  auto* g4particle = new G4PrimaryParticle(particleDefinition,
209  momentum.Px() * CLHEP::GeV,
210  momentum.Py() * CLHEP::GeV,
211  momentum.Pz() * CLHEP::GeV);
212 
213  // Add more particle information the Geant4 particle.
214  g4particle->SetCharge(particleDefinition->GetPDGCharge());
215  g4particle->SetPolarization(polarization.x(), polarization.y(), polarization.z());
216 
217  // Add the particle to the vertex.
218  vertex->SetPrimary(g4particle);
219 
220  // Create a PrimaryParticleInformation object, and save the
221  // MCTruth pointer in it. This will allow the
222  // ParticleActionList class to access MCTruth information
223  // during Geant4's tracking.
224 
225  auto* primaryParticleInfo = new g4b::PrimaryParticleInformation;
226  primaryParticleInfo->SetMCTruth(mclist.get(), index, m);
227 
228  // Save the PrimaryParticleInformation in the
229  // G4PrimaryParticle for access during tracking.
230  g4particle->SetUserInformation(primaryParticleInfo);
231  } // -- for each particle in MCTruth
232  ++index;
233  } // -- for each MCTruth entry
234  } // -- for each MCTruth handle
235 } // -- generatePrimaries()
Float_t x
Definition: compare.C:6
const TVector3 & Polarization() const
Definition: MCParticle.h:215
int PdgCode() const
Definition: MCParticle.h:213
static G4ParticleTable * fParticleTable
Geant4&#39;s table of particle definitions.
Float_t y
Definition: compare.C:6
Double_t z
Definition: plot.C:276
int StatusCode() const
Definition: MCParticle.h:212
int TrackId() const
Definition: MCParticle.h:211
Float_t Z
Definition: plot.C:37
std::vector< art::Handle< std::vector< simb::MCTruth > > > const * fMCLists
MCTruthCollection input lists.
double T(const int i=0) const
Definition: MCParticle.h:225
std::map< G4int, G4int > fProcessedPDG
PDG codes processed.
#define MF_LOG_INFO(category)
std::map< G4int, G4int > fNon1StatusPDG
PDG codes skipped because not status 1.
double Vx(const int i=0) const
Definition: MCParticle.h:222
std::map< G4int, G4int > fUnknownPDG
map of unknown PDG codes to instances
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:221
#define MF_LOG_DEBUG(id)
double Vz(const int i=0) const
Definition: MCParticle.h:224
void SetMCTruth(const simb::MCTruth *m, size_t idx=0, GeneratedParticleIndex_t indexInTruth=simb::NoGeneratedParticleIndex)
double Vy(const int i=0) const
Definition: MCParticle.h:223
vertex reconstruction
virtual void artg4tk::ActionBase::initialize ( )
inlinevirtualinherited

Reimplemented in artg4tk::myParticleGunActionService, and artg4tk::HepevtInputActionService.

Definition at line 36 of file ActionBase.hh.

37  {}
std::string const& artg4tk::ActionBase::myName ( ) const
inlineinherited

Definition at line 25 of file ActionBase.hh.

References artg4tk::ActionBase::myName_.

26  {
27  return myName_;
28  }
std::string myName_
Definition: ActionBase.hh:41
void larg4::MCTruthEventActionService::setInputCollections ( std::vector< art::Handle< std::vector< simb::MCTruth >>> const &  mclists)
inline

Definition at line 45 of file MCTruthEventAction_service.h.

46  {
47  fMCLists = &mclists;
48  }
std::vector< art::Handle< std::vector< simb::MCTruth > > > const * fMCLists
MCTruthCollection input lists.

Member Data Documentation

std::vector<art::Handle<std::vector<simb::MCTruth> > > const* larg4::MCTruthEventActionService::fMCLists
private

MCTruthCollection input lists.

Definition at line 58 of file MCTruthEventAction_service.h.

Referenced by generatePrimaries().

std::map<G4int, G4int> larg4::MCTruthEventActionService::fNon1StatusPDG
private

PDG codes skipped because not status 1.

Definition at line 60 of file MCTruthEventAction_service.h.

Referenced by generatePrimaries(), and ~MCTruthEventActionService().

G4ParticleTable * larg4::MCTruthEventActionService::fParticleTable = nullptr
staticprivate

Geant4's table of particle definitions.

Definition at line 56 of file MCTruthEventAction_service.h.

Referenced by generatePrimaries().

std::map<G4int, G4int> larg4::MCTruthEventActionService::fProcessedPDG
private

PDG codes processed.

Definition at line 61 of file MCTruthEventAction_service.h.

Referenced by generatePrimaries(), and ~MCTruthEventActionService().

std::map<G4int, G4int> larg4::MCTruthEventActionService::fUnknownPDG
private

map of unknown PDG codes to instances

Definition at line 59 of file MCTruthEventAction_service.h.

Referenced by generatePrimaries(), and ~MCTruthEventActionService().


The documentation for this class was generated from the following files: