Belle II Software development
MCParticleGenerator Class Reference

This class creates Geant4 primary particles from MCParticle list The generator assumes that each particle in the list has exactly one mother. More...

#include <MCParticleGenerator.h>

Inheritance diagram for MCParticleGenerator:

Public Member Functions

 MCParticleGenerator (const std::string &mcCollectionName, MCParticleGraph &mcParticleGraph)
 The constructor of the MCParticleGenerator class.
 
virtual ~MCParticleGenerator ()
 The destructor of the MCParticleGenerator class.
 
virtual void GeneratePrimaryVertex (G4Event *event)
 Create G4 primary particles from MCParticle list.
 

Protected Member Functions

void addParticle (const MCParticle &mcParticle, G4Event *event, G4PrimaryParticle *lastG4Mother, int motherIndex, bool useTime)
 Takes a MCParticle and creates a primary particle for Geant4.
 
G4PrimaryVertex * determineVertex (const MCParticleGraph::GraphParticle &p, double &productionTimeShift)
 Create a simulation vertex for the given particle.
 

Protected Attributes

G4VSolid * m_topSolid {nullptr}
 Pointer to the top level solid to check if particles are inside the simulation volume.
 
std::string m_mcCollectionName
 Name of the MCParticle collection.
 
MCParticleGraphm_mcParticleGraph
 Reference to the MCParticle graph.
 

Detailed Description

This class creates Geant4 primary particles from MCParticle list The generator assumes that each particle in the list has exactly one mother.

Definition at line 33 of file MCParticleGenerator.h.

Constructor & Destructor Documentation

◆ MCParticleGenerator()

MCParticleGenerator ( const std::string &  mcCollectionName,
MCParticleGraph mcParticleGraph 
)

The constructor of the MCParticleGenerator class.

Parameters
mcCollectionNameMCParticle collection from which MCParticles are read.
mcParticleGraphReference to the MCParticle graph that will be filled

Definition at line 32 of file MCParticleGenerator.cc.

32 :
33 G4VPrimaryGenerator(), m_mcCollectionName(mcCollectionName), m_mcParticleGraph(mcParticleGraph)
34{
35 G4VPhysicalVolume* topVolume = geometry::GeometryManager::getInstance().getTopVolume();
36 if (topVolume) {
37 m_topSolid = topVolume->GetLogicalVolume()->GetSolid();
38 }
39}
std::string m_mcCollectionName
Name of the MCParticle collection.
MCParticleGraph & m_mcParticleGraph
Reference to the MCParticle graph.
G4VSolid * m_topSolid
Pointer to the top level solid to check if particles are inside the simulation volume.
G4VPhysicalVolume * getTopVolume()
Return a pointer to the top volume.
static GeometryManager & getInstance()
Return a reference to the instance.

◆ ~MCParticleGenerator()

~MCParticleGenerator ( )
virtual

The destructor of the MCParticleGenerator class.

Definition at line 42 of file MCParticleGenerator.cc.

43{
44
45}

Member Function Documentation

◆ addParticle()

void addParticle ( const MCParticle mcParticle,
G4Event *  event,
G4PrimaryParticle *  lastG4Mother,
int  motherIndex,
bool  useTime 
)
protected

Takes a MCParticle and creates a primary particle for Geant4.

The daughters of the specified MCParticle are added recursively.

Parameters
mcParticlea reference to MCParticle
eventa pointer to Geant4 event
lastG4Mothera pointer to G4 mother particle
motherIndexmother index from MCParticle graph
useTimeuse MCParticle decay time (as given by the generator)

Definition at line 112 of file MCParticleGenerator.cc.

117{
118 G4PrimaryParticle* g4Mother = lastG4Mother;
119
120 //Check if the particle should be added to Geant4
121 //Only add particle if its PDG code is known to Geant4 and it is not flagged as virtual or initial.
122 bool addToG4 = (!mcParticle.isVirtual()) && (!mcParticle.isInitial());
123
124 //MS: distinguish optical photon from the rest of particles
125 bool opticalPhoton = mcParticle.getPDG() == 0 && mcParticle.getEnergy() < 10.0 * Unit::eV;
126 //SD: add neutral long-lived particle
127 bool neutral_llp = mcParticle.getCharge() == 0.0 && mcParticle.getLifetime() > 0.0;
128
129 G4ParticleDefinition* pdef = NULL;
130 if (opticalPhoton) {
131 pdef = G4OpticalPhoton::OpticalPhotonDefinition();
132 } else {
133 pdef = G4ParticleTable::GetParticleTable()->FindParticle(mcParticle.getPDG());
134 if (pdef == NULL && neutral_llp) {
135 pdef = G4ParticleTable::GetParticleTable()->FindParticle("LongLivedNeutralParticle");
136 }
137 }
138
139 if (pdef == NULL) {
140 // RG, DK when the skipped particle has a long lived and it is not flagged as virtual
141 if (addToG4 && mcParticle.getLifetime() != 0.0)
142 B2WARNING("PDG code " << mcParticle.getPDG() << " unknown to Geant4 - particle skipped.");
143 addToG4 = false;
144 }
145
146 //Add particle to MCParticle graph
148 graphParticle.setPDG(mcParticle.getPDG());
149 graphParticle.setStatus(mcParticle.getStatus());
150 graphParticle.setMass(mcParticle.getMass());
151 graphParticle.setEnergy(mcParticle.getEnergy());
152 graphParticle.setValidVertex(mcParticle.hasValidVertex());
153 graphParticle.setProductionTime(mcParticle.getProductionTime());
154 graphParticle.setDecayTime(mcParticle.getDecayTime());
155 graphParticle.setProductionVertex(mcParticle.getProductionVertex());
156 graphParticle.setMomentum(mcParticle.getMomentum());
157 graphParticle.setDecayVertex(mcParticle.getDecayVertex());
158 graphParticle.setFirstDaughter(mcParticle.getFirstDaughter());
159 graphParticle.setLastDaughter(mcParticle.getLastDaughter());
160 //Add decay to MCParticle graph
161 if (motherIndex > 0) graphParticle.comesFrom(m_mcParticleGraph[motherIndex - 1]);
162
163 //Create a new Geant4 Primary particle and
164 //store the link to the GraphMCParticle object as user info.
165 G4PrimaryParticle* newPart = nullptr;
166 // And maybe we need to generate a new primary Geant4 vertex for this particle.
167 G4PrimaryVertex* vertex = nullptr;
168 // And in the rare case that the particle started outside the simulation volume
169 // we need an adjusted production time
170 double productionTimeShift{0};
171
172 // Ok, we want to add the particle and it doesn't have a mother so lets see
173 // if we can find a proper vertex
174 if (addToG4 and lastG4Mother == nullptr) {
175 vertex = determineVertex(graphParticle, productionTimeShift);
176 // we don't have a mother and no vertex so skip this particle
177 if (not vertex) addToG4 = false;
178 }
179 // So now we have a particle to add and either a mother or a vertex.
180 // Let's create the particle in Geant4
181 if (addToG4) {
182 const ROOT::Math::PxPyPzEVector mcPartMom4 = graphParticle.get4Vector();
183 newPart = new G4PrimaryParticle(pdef,
184 mcPartMom4.X() / Unit::MeV * CLHEP::MeV,
185 mcPartMom4.Y() / Unit::MeV * CLHEP::MeV,
186 mcPartMom4.Z() / Unit::MeV * CLHEP::MeV,
187 mcPartMom4.E() / Unit::MeV * CLHEP::MeV);
188 newPart->SetMass(graphParticle.getMass() / Unit::MeV * CLHEP::MeV);
189 if (opticalPhoton) {
190 ROOT::Math::XYZVector polarization = graphParticle.getDecayVertex(); // temporary stored here
191 newPart->SetPolarization(polarization.X(), polarization.Y(), polarization.Z());
192 }
193 newPart->SetUserInformation(new ParticleInfo(graphParticle));
194
195 //Set propagation time only if useTime is true, the MCparticle has a valid vertex and has children.
196 useTime &= graphParticle.hasValidVertex() && mcParticle.getFirstDaughter() > 0;
197 if (useTime) {
198 //ProperTime is in particle eigentime, so convert lab lifetime to eigentime
199 //correct for any possible shift between simulation vertex and generation vertex
200 double propertime = (graphParticle.getLifetime() - productionTimeShift) / mcPartMom4.Gamma();
201 newPart->SetProperTime(propertime);
202 }
203
204 if (lastG4Mother) {
205 lastG4Mother->SetDaughter(newPart);
206 // This particle will be produced at the decay vertex of its mother so we **want** Geant4 to
207 // overwrite production time/vertex, so as convention having NaN here means we'll
208 // take over the values from Geant4 instead of the ones from the generator
209 graphParticle.setProductionTime(std::numeric_limits<float>::quiet_NaN());
210 } else {
211 vertex->SetPrimary(newPart);
212 event->AddPrimaryVertex(vertex);
213 }
214
215 g4Mother = newPart;
216
217 //Do not store the generator info in the MCParticles block unless Geant4 creates a track in the detector.
218 graphParticle.setIgnore();
219
220 } else {
221 B2DEBUG(100, "The particle " << mcParticle.getIndex() << " (PDG " << mcParticle.getPDG() << ") was not added to Geant4");
222 }
223
224 //Add all children
225 int currMotherIndex = m_mcParticleGraph.size();
226 for (MCParticle* daughter : mcParticle.getDaughters()) {
227 addParticle(*daughter, event, g4Mother, currMotherIndex, useTime);
228 }
229}
Class to represent Particle data in graph.
void comesFrom(GraphParticle &mother)
Tells the graph that this particle is a decay product of mother.
void setFirstDaughter(int daughter)
Set the 1-based index of the first daughter, 0 means no daughters.
void setLastDaughter(int daughter)
Set the 1-based index of the last daughter, 0 means no daughters.
void setIgnore(bool ignore=true)
Set or remove the ignore flag.
size_t size() const
Return the number of particles in the graph.
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:32
float getEnergy() const
Return particle energy in GeV.
Definition: MCParticle.h:147
void setDecayTime(float time)
Set decay time.
Definition: MCParticle.h:390
void setMass(float mass)
Set particle mass.
Definition: MCParticle.h:366
void setDecayVertex(const ROOT::Math::XYZVector &vertex)
Set decay vertex.
Definition: MCParticle.h:447
int getIndex() const
Get 1-based index of the particle in the corresponding MCParticle list.
Definition: MCParticle.h:230
ROOT::Math::XYZVector getDecayVertex() const
Return decay vertex.
Definition: MCParticle.h:219
float getLifetime() const
Return the lifetime in ns.
Definition: MCParticle.h:177
void setEnergy(float energy)
Set energy.
Definition: MCParticle.h:372
std::vector< Belle2::MCParticle * > getDaughters() const
Get vector of all daughter particles, empty vector if none.
Definition: MCParticle.cc:52
float getMass() const
Return the particle mass in GeV.
Definition: MCParticle.h:135
ROOT::Math::XYZVector getProductionVertex() const
Return production vertex position.
Definition: MCParticle.h:189
float getDecayTime() const
Return the decay time in ns.
Definition: MCParticle.h:168
void setValidVertex(bool valid)
Set indication wether vertex and time information is valid or just default.
Definition: MCParticle.h:378
void setProductionVertex(const ROOT::Math::XYZVector &vertex)
Set production vertex position.
Definition: MCParticle.h:396
bool hasValidVertex() const
Indication whether vertex and time information is useful or just default.
Definition: MCParticle.h:153
unsigned int getStatus(unsigned short int bitmask=USHRT_MAX) const
Return status code of particle.
Definition: MCParticle.h:122
int getLastDaughter() const
Get 1-based index of last daughter, 0 if no daughters.
Definition: MCParticle.h:259
float getCharge() const
Return the particle charge defined in TDatabasePDG.
Definition: MCParticle.cc:36
ROOT::Math::PxPyPzEVector get4Vector() const
Return 4Vector of particle.
Definition: MCParticle.h:207
void setPDG(int pdg)
Set PDG code of the particle.
Definition: MCParticle.h:335
int getPDG() const
Return PDG code of particle.
Definition: MCParticle.h:112
float getProductionTime() const
Return production time in ns.
Definition: MCParticle.h:159
ROOT::Math::XYZVector getMomentum() const
Return momentum.
Definition: MCParticle.h:198
void setMomentum(const ROOT::Math::XYZVector &momentum)
Set particle momentum.
Definition: MCParticle.h:417
void setStatus(unsigned short int status)
Set Status code for the particle.
Definition: MCParticle.h:346
int getFirstDaughter() const
Get 1-based index of first daughter, 0 if no daughters.
Definition: MCParticle.h:251
void setProductionTime(float time)
Set production time.
Definition: MCParticle.h:384
void addParticle(const MCParticle &mcParticle, G4Event *event, G4PrimaryParticle *lastG4Mother, int motherIndex, bool useTime)
Takes a MCParticle and creates a primary particle for Geant4.
G4PrimaryVertex * determineVertex(const MCParticleGraph::GraphParticle &p, double &productionTimeShift)
Create a simulation vertex for the given particle.
UserInfo class which is used to attach additional information to Geant4 particles and tracks.
Definition: UserInfo.h:36
static const double eV
[electronvolt]
Definition: Unit.h:112
static const double MeV
[megaelectronvolt]
Definition: Unit.h:114
bool isInitial() const
Check if particle is an initial particle such as ISR.
Definition: MCParticle.h:590
GraphParticle & addParticle()
Add new particle to the graph.
bool isVirtual() const
Check if particle is virtual.
Definition: MCParticle.h:575

◆ determineVertex()

G4PrimaryVertex * determineVertex ( const MCParticleGraph::GraphParticle p,
double &  productionTimeShift 
)
protected

Create a simulation vertex for the given particle.

This checks if the particle starts inside the simulation volume. If it's inside trivially just create a vertex at the given position and time. If not check if the particle intersects with the simulation volume.

If it does intersect and will survive long enough to get there we create a vertex at the simulation volume boundary and also set productionTimeShift to the flight time to get to the simulation volume.

Parameters
pParticle to create the vertex for
productionTimeShiftreturn the flight time of the particle before reaching the simulation volume
Returns
a new G4PrimaryVertex or nullptr if the particle didn't intersect with the simulation volume in time.

Definition at line 74 of file MCParticleGenerator.cc.

75{
76 // We want to determine the simulation vertex for the given particle.
77 // So lets see where it is
78 ROOT::Math::XYZVector mcProdVtx = p.getProductionVertex();
79 G4ThreeVector pos(mcProdVtx.X() / Unit::mm * CLHEP::mm,
80 mcProdVtx.Y() / Unit::mm * CLHEP::mm,
81 mcProdVtx.Z() / Unit::mm * CLHEP::mm);
82 productionTimeShift = 0;
83 // and make sure its inside the simulation volume
84 if (m_topSolid && m_topSolid->Inside(pos) == kOutside) {
85 // Ok, we're outside the simulation volume. Check if we can actually
86 // get to the simulation volume in time, no fields outside the volume
87 ROOT::Math::PxPyPzEVector mcPartMom4 = p.get4Vector();
88 CLHEP::HepLorentzVector mom(mcPartMom4.X() / Unit::MeV * CLHEP::MeV,
89 mcPartMom4.Y() / Unit::MeV * CLHEP::MeV,
90 mcPartMom4.Z() / Unit::MeV * CLHEP::MeV,
91 mcPartMom4.E() / Unit::MeV * CLHEP::MeV);
92 const auto dir = mom.vect().unit();
93 double distance = m_topSolid->DistanceToIn(pos, dir);
94
95 if (distance == kInfinity) {
96 B2DEBUG(10, "Particle starts outside simulation volume and isn't intersecting, skipping");
97 return nullptr;
98 }
99
100 const double speed = (p.getMass() == 0 ? mom.beta() : 1) * CLHEP::c_light;
101 const double flightLength = speed * p.getLifetime();
102 if (distance > flightLength) {
103 B2DEBUG(10, "Particle starts outside simulation volume and doesn't reach it before decaying, skipping");
104 return nullptr;
105 }
106 pos += dir * distance;
107 productionTimeShift = (distance / speed);
108 }
109 return new G4PrimaryVertex(pos.x(), pos.y(), pos.z(), p.getProductionTime() + productionTimeShift);
110}
static const double mm
[millimeters]
Definition: Unit.h:70

◆ GeneratePrimaryVertex()

void GeneratePrimaryVertex ( G4Event *  event)
virtual

Create G4 primary particles from MCParticle list.

This method is called for each event.

Parameters
eventPointer to the G4Event.

Definition at line 48 of file MCParticleGenerator.cc.

49{
50 //Get the MCParticle collection from the DataStore
52
53 //Prepare the MCParticle graph
55
56 //Loop over the primary particles. The MCParticle collection has to be
57 //sorted breadth first: primary particles come first and then the daughters
58 // ... at least that would be the case if there is only one generator, but if
59 // more than one generator is used that is not necessarily true. So ignore
60 // that last statement.
61 //
62 // Let's add all top level particles, i.e. without mother, and addParticle
63 // will be called recursively for all daughters
64 int nPart = mcParticles.getEntries();
65 for (int iPart = 0; iPart < nPart; iPart++) {
66 MCParticle* currParticle = mcParticles[iPart];
67 if (currParticle->getMother() != NULL) continue;
68
69 //Add primary particle (+ daughters) and the vertex
70 addParticle(*currParticle, event, NULL, 0, true);
71 }
72}
Accessor to arrays stored in the data store.
Definition: StoreArray.h:113
MCParticle * getMother() const
Returns a pointer to the mother particle.
Definition: MCParticle.h:600
void clear()
Reset particles and decay information to make the class reusable.

Member Data Documentation

◆ m_mcCollectionName

std::string m_mcCollectionName
protected

Name of the MCParticle collection.

Definition at line 59 of file MCParticleGenerator.h.

◆ m_mcParticleGraph

MCParticleGraph& m_mcParticleGraph
protected

Reference to the MCParticle graph.

Definition at line 60 of file MCParticleGenerator.h.

◆ m_topSolid

G4VSolid* m_topSolid {nullptr}
protected

Pointer to the top level solid to check if particles are inside the simulation volume.

Definition at line 58 of file MCParticleGenerator.h.


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