ConvertMCTruthToG4.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \file ConvertMCTruthToG4.cxx
3 /// \brief Convert MCTruth to G4Event; Geant4 event generator
4 ///
5 /// \version $Id: ConvertMCTruthToG4.cxx,v 1.7 2012-09-24 15:19:29 brebel Exp $
6 /// \author seligman@nevis.columbia.edu, brebel@fnal.gov
7 ////////////////////////////////////////////////////////////////////////
8 
9 #include "nutools/G4Base/ConvertMCTruthToG4.h"
10 #include "nutools/G4Base/PrimaryParticleInformation.h"
13 
14 #include "Geant4/G4Event.hh"
15 #include "Geant4/G4IonTable.hh"
16 #include "Geant4/G4PrimaryVertex.hh"
17 #include "Geant4/G4ParticleDefinition.hh"
18 #include "Geant4/G4ParticleTable.hh"
19 
20 #include <TParticle.h>
21 
22 #include <CLHEP/Vector/LorentzVector.h>
23 
24 #include <iostream>
25 #include <vector>
26 #include <map>
27 
29 
30 namespace g4b{
31 
32  // Static variables used by class.
33 
34  // Geant4's particle table.
35  G4ParticleTable* ConvertMCTruthToG4::fParticleTable = 0;
36 
37  //-----------------------------------------------------
38  // Constructor and destructor.
40  {
41  }
42 
43  //-----------------------------------------------------
45  {
46  // Print out a list of "unknown" PDG codes we saw in the input.
47  if ( ! fUnknownPDG.empty() ){
48  std::ostringstream badtxt;
50  for ( ; i != fUnknownPDG.end(); ++i ){
51  int pdg = (*i).first;
52  badtxt << "\n Unknown PDG code = " << pdg
53  << ", seen " << (*i).second << " times.";
54 
55  const int genieLo = 2000000001;
56  const int genieHi = 2000000202;
57  if ( pdg >= genieLo && pdg <= genieHi ) {
58  badtxt << " (GENIE specific)";
59  }
60  }
61 
62  mf::LogWarning("ConvertPrimaryToGeant4")
63  << "The following unknown PDG codes were present in the "
64  << "simb::MCTruth input.\n"
65  << "They were not processed by Geant4."
66  << badtxt.str();
67  }
68 
69  }
70 
71  //-----------------------------------------------------
73  {
74  fConvertList.clear();
75  }
76 
77  //-----------------------------------------------------
79  {
80  this->Append( mct.get() );
81  }
82 
83  //-----------------------------------------------------
85  {
86  fConvertList.push_back( mct );
87  }
88 
89  //-----------------------------------------------------
91  {
92  // Take the contents of MCTruth objects and use them to
93  // initialize the G4Event.
94 
95  // A G4Event organizes its particles in terms of "vertices" and
96  // "particles", like HepMC. Unfortunately, ROOT doesn't use
97  // HepMC, so the MCTruth objects aren't organized that way.
98  // For most of the work that we'll ever do, there'll be only one
99  // vertex in the event. However, just in case there are multiple
100  // vertices (e.g., overlays, double vertex studies) I want the
101  // code to function properly.
102 
103  // So create a map of particle positions and associated
104  // G4PrimaryVertex*. Note that the map must use CLHEP's
105  // LorentzVector, and not ROOT's, since ROOT does not define an
106  // operator< for its physics vectors.
107  std::map< CLHEP::HepLorentzVector, G4PrimaryVertex* > vertexMap;
109 
110  // For each MCTruth (probably only one, but you never know):
111  // index keeps track of which MCTruth object you are using
112  size_t index = 0;
113  for( auto const* mct : fConvertList){
114 
115  // For each simb::MCParticle in the MCTruth:
116  for ( int p = 0; p != mct->NParticles(); ++p ){
117  // Implementation note: the following statement copies the
118  // MCParticle from the MCTruth, instead of getting a
119  // const reference.
120  simb::MCParticle const& particle = mct->GetParticle(p);
121 
122  // status code == 1 means "track this particle." Any
123  // other status code should be ignored by the Monte Carlo.
124  if ( particle.StatusCode() != 1 ) continue;
125 
126  // Get the Particle Data Group code for the particle.
127  G4int pdgCode = particle.PdgCode();
128 
129  // Get the vertex. Note that LArSoft/ROOT uses cm, but
130  // Geant4/CLHEP uses mm.
131  G4double x = particle.Vx() * CLHEP::cm;
132  G4double y = particle.Vy() * CLHEP::cm;
133  G4double z = particle.Vz() * CLHEP::cm;
134  G4double t = particle.T() * CLHEP::ns;
135 
136  // Create a CLHEP four-vector from the particle's vertex.
137  CLHEP::HepLorentzVector fourpos(x,y,z,t);
138 
139  // Is this vertex already in our map?
140  G4PrimaryVertex* vertex = 0;
142  if ( result == vertexMap.end() ){
143  // No, it's not, so create a new vertex and add it to the
144  // map.
145  vertex = new G4PrimaryVertex(x, y, z, t);
146  vertexMap[ fourpos ] = vertex;
147 
148  // Add the vertex to the G4Event.
149  event->AddPrimaryVertex( vertex );
150  }
151  else{
152  // Yes, it is, so use the existing vertex.
153  vertex = (*result).second;
154  }
155 
156  // Get additional particle information.
157  TLorentzVector momentum = particle.Momentum(); // (px,py,pz,E)
158  TVector3 polarization = particle.Polarization();
159 
160  // Get the particle table if necessary. (Note: we're
161  // doing this "late" because I'm not sure at what point
162  // the G4 particle table is initialized in the loading process.
163  if ( fParticleTable == 0 ){
164  fParticleTable = G4ParticleTable::GetParticleTable();
165  }
166 
167  // Get Geant4's definition of the particle.
168  G4ParticleDefinition* particleDefinition;
169 
170  if(pdgCode==0){
171  particleDefinition = fParticleTable->FindParticle("opticalphoton");
172  }
173  else
174  particleDefinition = fParticleTable->FindParticle(pdgCode);
175 
176  if ( pdgCode > 1000000000) { // If the particle is a nucleus
177  MF_LOG_DEBUG("ConvertPrimaryToGeant4") << ": %%% Nuclear PDG code = " << pdgCode
178  << " (x,y,z,t)=(" << x
179  << "," << y
180  << "," << z
181  << "," << t << ")"
182  << " P=" << momentum.P()
183  << ", E=" << momentum.E();
184  // If the particle table doesn't have a definition yet, ask the ion
185  // table for one. This will create a new ion definition as needed.
186  if (!particleDefinition) {
187  int Z = (pdgCode % 10000000) / 10000; // atomic number
188  int A = (pdgCode % 10000) / 10; // mass number
189  particleDefinition = fParticleTable->GetIonTable()->GetIon(Z, A, 0.);
190  }
191  }
192 
193  // What if the PDG code is unknown? This has been a known
194  // issue with GENIE.
195  if ( particleDefinition == 0 ){
196  MF_LOG_DEBUG("ConvertPrimaryToGeant4") << ": %%% Code not found = " << pdgCode;
197  fUnknownPDG[ pdgCode ] += 1;
198  continue;
199  }
200 
201  // Create a Geant4 particle to add to the vertex.
202  G4PrimaryParticle* g4particle = new G4PrimaryParticle( particleDefinition,
203  momentum.Px() * CLHEP::GeV,
204  momentum.Py() * CLHEP::GeV,
205  momentum.Pz() * CLHEP::GeV);
206 
207  // Add more particle information the Geant4 particle.
208  G4double charge = particleDefinition->GetPDGCharge();
209  g4particle->SetCharge( charge );
210  g4particle->SetPolarization( polarization.x(),
211  polarization.y(),
212  polarization.z() );
213 
214  // Add the particle to the vertex.
215  vertex->SetPrimary( g4particle );
216 
217  // Create a PrimaryParticleInformation object, and save
218  // the MCTruth pointer in it. This will allow the
219  // ParticleActionList class to access MCTruth
220  // information during Geant4's tracking.
221  PrimaryParticleInformation* primaryParticleInfo = new PrimaryParticleInformation;
222  primaryParticleInfo->SetMCTruth( mct, index, p );
223 
224  // Save the PrimaryParticleInformation in the
225  // G4PrimaryParticle for access during tracking.
226  g4particle->SetUserInformation( primaryParticleInfo );
227 
228  MF_LOG_DEBUG("ConvertPrimaryToGeant4") << ": %%% primary PDG=" << pdgCode
229  << ", (x,y,z,t)=(" << x
230  << "," << y
231  << "," << z
232  << "," << t << ")"
233  << " P=" << momentum.P()
234  << ", E=" << momentum.E();
235 
236  } // for each particle in MCTruth
237  ++index;
238  } // for each MCTruth
239  }// GeneratePrimaries
240 }// namespace
intermediate_table::iterator iterator
const TVector3 & Polarization() const
Definition: MCParticle.h:213
int PdgCode() const
Definition: MCParticle.h:211
virtual void GeneratePrimaries(G4Event *)
static QCString result
void Reset()
Get ready to convert a new set of MCTruth objects.
static G4ParticleTable * fParticleTable
Geant4&#39;s table of particle definitions.
int StatusCode() const
Definition: MCParticle.h:210
intermediate_table::const_iterator const_iterator
Particle class.
double y
ConvertMCTruthToG4()
Standard constructor and destructor.
double z
double T(const int i=0) const
Definition: MCParticle.h:223
basic interface to Geant4 for ART-based software
p
Definition: test.py:223
double Vx(const int i=0) const
Definition: MCParticle.h:220
std::map< G4int, G4int > fUnknownPDG
map of unknown PDG codes to instances
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:219
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
#define MF_LOG_DEBUG(id)
double Vz(const int i=0) const
Definition: MCParticle.h:222
static const double cm
Definition: Units.h:76
list x
Definition: train.py:276
void SetMCTruth(const simb::MCTruth *m, size_t idx=0, GeneratedParticleIndex_t indexInTruth=simb::NoGeneratedParticleIndex)
void Append(art::Ptr< simb::MCTruth > &mct)
Event generator information.
Definition: MCTruth.h:32
T const * get() const
Definition: Ptr.h:149
std::vector< const simb::MCTruth * > fConvertList
List of MCTruth objects to convert for this spill.
static const double GeV
Definition: Units.h:29
QAsciiDict< Entry > ns
double Vy(const int i=0) const
Definition: MCParticle.h:221
Event finding and building.
vertex reconstruction