Belle II Software  release-05-01-25
MCParticleGenerator.cc
1 /*****************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010-2011 Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Andreas Moll, Martin Ritter, Doris Kim, Romulus Godang, *
7  * Marko Staric *
8  * *
9  * This software is provided "as is" without any warranty. *
10  * *
11  * 7/28/13 Modified (addToG4) to (addToG4 && mcParticle.getLifetime() != 0.0)*
12  * to deal with the resonance particles that enter the detector in GEANT4 *
13  * by R. Godang and Doris Kim *
14  *****************************************************************************/
15 
16 #include <simulation/kernel/MCParticleGenerator.h>
17 
18 #include <mdst/dataobjects/MCParticle.h>
19 #include <framework/gearbox/Unit.h>
20 #include <framework/logging/Logger.h>
21 #include <framework/datastore/StoreArray.h>
22 #include <simulation/kernel/UserInfo.h>
23 #include <geometry/GeometryManager.h>
24 
25 #include <G4ParticleTable.hh>
26 #include <G4VUserPrimaryParticleInformation.hh>
27 #include <G4ParticleDefinition.hh>
28 #include <G4ParticleTypes.hh>
29 #include <G4VPhysicalVolume.hh>
30 #include <G4LogicalVolume.hh>
31 #include <G4VSolid.hh>
32 #include <CLHEP/Vector/LorentzVector.h>
33 
34 #include <TLorentzVector.h>
35 
36 using namespace std;
37 using namespace Belle2;
38 using namespace Belle2::Simulation;
39 
40 
41 MCParticleGenerator::MCParticleGenerator(const string& mcCollectionName, MCParticleGraph& mcParticleGraph) :
42  G4VPrimaryGenerator(), m_mcCollectionName(mcCollectionName), m_mcParticleGraph(mcParticleGraph)
43 {
44  G4VPhysicalVolume* topVolume = geometry::GeometryManager::getInstance().getTopVolume();
45  if (topVolume) {
46  m_topSolid = topVolume->GetLogicalVolume()->GetSolid();
47  }
48 }
49 
50 
52 {
53 
54 }
55 
56 
58 {
59  //Get the MCParticle collection from the DataStore
61 
62  //Prepare the MCParticle graph
64 
65  //Loop over the primary particles. The MCParticle collection has to be
66  //sorted breadth first: primary particles come first and then the daughters
67  // ... at least that would be the case if there is only one generator, but if
68  // more than one generator is used that is not necessarily true. So ignore
69  // that last statement.
70  //
71  // Let's add all top level particles, i.e. without mother, and addParticle
72  // will be called recursively for all daughters
73  int nPart = mcParticles.getEntries();
74  for (int iPart = 0; iPart < nPart; iPart++) {
75  MCParticle* currParticle = mcParticles[iPart];
76  if (currParticle->getMother() != NULL) continue;
77 
78  //Add primary particle (+ daughters) and the vertex
79  addParticle(*currParticle, event, NULL, 0, true);
80  }
81 }
82 
83 G4PrimaryVertex* MCParticleGenerator::determineVertex(const MCParticleGraph::GraphParticle& p, double& productionTimeShift)
84 {
85  // We want to determine the simulation vertex for the given particle.
86  // So lets see where it is
87  TVector3 mcProdVtx = p.getProductionVertex();
88  G4ThreeVector pos(mcProdVtx.X() / Unit::mm * CLHEP::mm,
89  mcProdVtx.Y() / Unit::mm * CLHEP::mm,
90  mcProdVtx.Z() / Unit::mm * CLHEP::mm);
91  productionTimeShift = 0;
92  // and make sure its inside the simulation volume
93  if (m_topSolid && m_topSolid->Inside(pos) == kOutside) {
94  // Ok, we're outside the simulation volume. Check if we can actually
95  // get to the simulation volume in time, no fields outside the volume
96  TLorentzVector mcPartMom4 = p.get4Vector();
97  CLHEP::HepLorentzVector mom(mcPartMom4.X() / Unit::MeV * CLHEP::MeV,
98  mcPartMom4.Y() / Unit::MeV * CLHEP::MeV,
99  mcPartMom4.Z() / Unit::MeV * CLHEP::MeV,
100  mcPartMom4.E() / Unit::MeV * CLHEP::MeV);
101  const auto dir = mom.vect().unit();
102  double distance = m_topSolid->DistanceToIn(pos, dir);
103 
104  if (distance == kInfinity) {
105  B2DEBUG(10, "Particle starts outside simulation volume and isn't intersecting, skipping");
106  return nullptr;
107  }
108 
109  const double speed = (p.getMass() == 0 ? mom.beta() : 1) * CLHEP::c_light;
110  const double flightLength = speed * p.getLifetime();
111  if (distance > flightLength) {
112  B2DEBUG(10, "Particle starts outside simulation volume and doesn't reach it before decaying, skipping");
113  return nullptr;
114  }
115  pos += dir * distance;
116  productionTimeShift = (distance / speed);
117  }
118  return new G4PrimaryVertex(pos.x(), pos.y(), pos.z(), p.getProductionTime() + productionTimeShift);
119 }
120 
122  G4Event* event,
123  G4PrimaryParticle* lastG4Mother,
124  int motherIndex,
125  bool useTime)
126 {
127  G4PrimaryParticle* g4Mother = lastG4Mother;
128 
129  //Check if the particle should be added to Geant4
130  //Only add particle if its PDG code is known to Geant4 and it is not flagged as virtual or initial.
131  bool addToG4 = (!mcParticle.isVirtual()) && (!mcParticle.isInitial());
132 
133  //MS: distinguish optical photon from the rest of particles
134  bool opticalPhoton = mcParticle.getPDG() == 0 && mcParticle.getEnergy() < 10.0 * Unit::eV;
135  //SD: add neutral long-lived particle
136  bool neutral_llp = mcParticle.getCharge() == 0.0 && mcParticle.getLifetime() > 0.0;
137 
138  G4ParticleDefinition* pdef = NULL;
139  if (opticalPhoton) {
140  pdef = G4OpticalPhoton::OpticalPhotonDefinition();
141  } else {
142  pdef = G4ParticleTable::GetParticleTable()->FindParticle(mcParticle.getPDG());
143  if (pdef == NULL && neutral_llp) {
144  pdef = G4ParticleTable::GetParticleTable()->FindParticle("LongLivedNeutralParticle");
145  }
146  }
147 
148  if (pdef == NULL) {
149  // RG, DK when the skipped particle has a long lived and it is not flagged as virtual
150  if (addToG4 && mcParticle.getLifetime() != 0.0)
151  B2WARNING("PDG code " << mcParticle.getPDG() << " unknown to Geant4 - particle skipped.");
152  addToG4 = false;
153  }
154 
155  //Add particle to MCParticle graph
157  graphParticle.setPDG(mcParticle.getPDG());
158  graphParticle.setStatus(mcParticle.getStatus());
159  graphParticle.setMass(mcParticle.getMass());
160  graphParticle.setEnergy(mcParticle.getEnergy());
161  graphParticle.setValidVertex(mcParticle.hasValidVertex());
162  graphParticle.setProductionTime(mcParticle.getProductionTime());
163  graphParticle.setDecayTime(mcParticle.getDecayTime());
164  graphParticle.setProductionVertex(mcParticle.getProductionVertex());
165  graphParticle.setMomentum(mcParticle.getMomentum());
166  graphParticle.setDecayVertex(mcParticle.getDecayVertex());
167  graphParticle.setFirstDaughter(mcParticle.getFirstDaughter());
168  graphParticle.setLastDaughter(mcParticle.getLastDaughter());
169  //Add decay to MCParticle graph
170  if (motherIndex > 0) graphParticle.comesFrom(m_mcParticleGraph[motherIndex - 1]);
171 
172  //Create a new Geant4 Primary particle and
173  //store the link to the GraphMCParticle object as user info.
174  G4PrimaryParticle* newPart = nullptr;
175  // And maybe we need to generate a new primary Geant4 vertex for this particle.
176  G4PrimaryVertex* vertex = nullptr;
177  // And in the rare case that the particle started outside the simulation volume
178  // we need an adjusted production time
179  double productionTimeShift{0};
180 
181  // Ok, we want to add the particle and it doesn't have a mother so lets see
182  // if we can find a proper vertex
183  if (addToG4 and lastG4Mother == nullptr) {
184  vertex = determineVertex(graphParticle, productionTimeShift);
185  // we don't have a mother and no vertex so skip this particle
186  if (not vertex) addToG4 = false;
187  }
188  // So now we have a particle to add and either a mother or a vertex.
189  // Let's create the particle in Geant4
190  if (addToG4) {
191  const TLorentzVector mcPartMom4 = graphParticle.get4Vector();
192  newPart = new G4PrimaryParticle(pdef,
193  mcPartMom4.X() / Unit::MeV * CLHEP::MeV,
194  mcPartMom4.Y() / Unit::MeV * CLHEP::MeV,
195  mcPartMom4.Z() / Unit::MeV * CLHEP::MeV,
196  mcPartMom4.E() / Unit::MeV * CLHEP::MeV);
197  newPart->SetMass(graphParticle.getMass() / Unit::MeV * CLHEP::MeV);
198  if (opticalPhoton) {
199  TVector3 polarization = graphParticle.getDecayVertex(); // temporary stored here
200  newPart->SetPolarization(polarization.X(), polarization.Y(), polarization.Z());
201  }
202  newPart->SetUserInformation(new ParticleInfo(graphParticle));
203 
204  //Set propagation time only if useTime is true, the MCparticle has a valid vertex and has children.
205  useTime &= graphParticle.hasValidVertex() && mcParticle.getFirstDaughter() > 0;
206  if (useTime) {
207  //ProperTime is in particle eigentime, so convert lab lifetime to eigentime
208  //correct for any possible shift between simulation vertex and generation vertex
209  double propertime = (graphParticle.getLifetime() - productionTimeShift) / mcPartMom4.Gamma();
210  newPart->SetProperTime(propertime);
211  }
212 
213  if (lastG4Mother) {
214  lastG4Mother->SetDaughter(newPart);
215  // This particle will be produced at the decay vertex of its mother so we **want** Geant4 to
216  // overwrite production time/vertex, so as convention having NaN here means we'll
217  // take over the values from Geant4 instead of the ones from the generator
218  graphParticle.setProductionTime(std::numeric_limits<float>::quiet_NaN());
219  } else {
220  vertex->SetPrimary(newPart);
221  event->AddPrimaryVertex(vertex);
222  }
223 
224  g4Mother = newPart;
225 
226  //Do not store the generator info in the MCParticles block unless Geant4 creates a track in the detector.
227  graphParticle.setIgnore();
228 
229  } else {
230  B2DEBUG(100, "The particle " << mcParticle.getIndex() << " (PDG " << mcParticle.getPDG() << ") was not added to Geant4");
231  }
232 
233  //Add all children
234  int currMotherIndex = m_mcParticleGraph.size();
235  for (MCParticle* daughter : mcParticle.getDaughters()) {
236  addParticle(*daughter, event, g4Mother, currMotherIndex, useTime);
237  }
238 }
239 
240 
241 //===================================================================
242 // Protected methods
243 //===================================================================
Belle2::MCParticle::getFirstDaughter
int getFirstDaughter() const
Get 1-based index of first daughter, 0 if no daughters.
Definition: MCParticle.h:262
Belle2::MCParticle::getIndex
int getIndex() const
Get 1-based index of the particle in the corresponding MCParticle list.
Definition: MCParticle.h:241
Belle2::MCParticle::getEnergy
float getEnergy() const
Return particle energy in GeV.
Definition: MCParticle.h:158
Belle2::MCParticle::getCharge
float getCharge() const
Return the particle charge defined in TDatabasePDG.
Definition: MCParticle.cc:35
Belle2::MCParticleGraph::size
size_t size() const
Return the number of particles in the graph.
Definition: MCParticleGraph.h:253
Belle2::MCParticle::getDecayTime
float getDecayTime() const
Return the decay time in ns.
Definition: MCParticle.h:179
Belle2::Simulation::MCParticleGenerator::~MCParticleGenerator
virtual ~MCParticleGenerator()
The destructor of the MCParticleGenerator class.
Definition: MCParticleGenerator.cc:51
Belle2::MCParticleGraph::GraphParticle::setIgnore
void setIgnore(bool ignore=true)
Set or remove the ignore flag.
Definition: MCParticleGraph.h:143
Belle2::MCParticleGraph
Class to build, validate and sort a particle decay chain.
Definition: MCParticleGraph.h:48
Belle2::Simulation::MCParticleGenerator::m_mcParticleGraph
MCParticleGraph & m_mcParticleGraph
Reference to the MCParticle graph.
Definition: MCParticleGenerator.h:62
Belle2::Simulation::MCParticleGenerator::m_mcCollectionName
std::string m_mcCollectionName
Name of the MCParticle collection.
Definition: MCParticleGenerator.h:61
Belle2::Simulation::MCParticleGenerator::m_topSolid
G4VSolid * m_topSolid
Pointer to the top level solid to check if particles are inside the simulation volume.
Definition: MCParticleGenerator.h:60
Belle2::Simulation::MCParticleGenerator::determineVertex
G4PrimaryVertex * determineVertex(const MCParticleGraph::GraphParticle &p, double &productionTimeShift)
Create a simulation vertex for the given particle.
Definition: MCParticleGenerator.cc:83
Belle2::Simulation::UserInfo
UserInfo class which is used to attach additional information to Geant4 particles and tracks.
Definition: UserInfo.h:46
Belle2::geometry::GeometryManager::getTopVolume
G4VPhysicalVolume * getTopVolume()
Return a pointer to the top volume.
Definition: GeometryManager.h:59
Belle2::MCParticle::setDecayVertex
void setDecayVertex(const TVector3 &vertex)
Set decay vertex.
Definition: MCParticle.h:452
Belle2::MCParticle::setStatus
void setStatus(unsigned short int status)
Set Status code for the particle.
Definition: MCParticle.h:354
Belle2::geometry::GeometryManager::getInstance
static GeometryManager & getInstance()
Return a reference to the instance.
Definition: GeometryManager.cc:98
Belle2::Simulation::MCParticleGenerator::GeneratePrimaryVertex
virtual void GeneratePrimaryVertex(G4Event *event)
Create G4 primary particles from MCParticle list.
Definition: MCParticleGenerator.cc:57
Belle2::MCParticle::isInitial
bool isInitial() const
Check if particle is an initial particle such as ISR.
Definition: MCParticle.h:583
Belle2::MCParticle::getDecayVertex
TVector3 getDecayVertex() const
Return decay vertex.
Definition: MCParticle.h:230
Belle2::MCParticle::hasValidVertex
bool hasValidVertex() const
Indication whether vertex and time information is useful or just default.
Definition: MCParticle.h:164
Belle2::Unit::MeV
static const double MeV
[megaelectronvolt]
Definition: Unit.h:124
Belle2::Unit::eV
static const double eV
[electronvolt]
Definition: Unit.h:122
Belle2::MCParticle::getStatus
unsigned int getStatus(unsigned short int bitmask=USHRT_MAX) const
Return status code of particle.
Definition: MCParticle.h:133
Belle2::MCParticle::setProductionTime
void setProductionTime(float time)
Set production time.
Definition: MCParticle.h:392
Belle2::MCParticleGraph::GraphParticle::setLastDaughter
void setLastDaughter(int daughter)
Set the 1-based index of the last daughter, 0 means no daughters.
Definition: MCParticleGraph.h:133
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::MCParticle::setPDG
void setPDG(int pdg)
Set PDG code of the particle.
Definition: MCParticle.h:343
Belle2::MCParticle::getProductionVertex
TVector3 getProductionVertex() const
Return production vertex position.
Definition: MCParticle.h:200
Belle2::MCParticle::getPDG
int getPDG() const
Return PDG code of particle.
Definition: MCParticle.h:123
Belle2::MCParticleGraph::addParticle
GraphParticle & addParticle()
Add new particle to the graph.
Definition: MCParticleGraph.h:305
Belle2::MCParticle::setProductionVertex
void setProductionVertex(const TVector3 &vertex)
Set production vertex position.
Definition: MCParticle.h:404
Belle2::MCParticle::getMass
float getMass() const
Return the particle mass in GeV.
Definition: MCParticle.h:146
Belle2::MCParticle::getMother
MCParticle * getMother() const
Returns a pointer to the mother particle.
Definition: MCParticle.h:593
Belle2::MCParticle::getLastDaughter
int getLastDaughter() const
Get 1-based index of last daughter, 0 if no daughters.
Definition: MCParticle.h:270
Belle2::MCParticle::isVirtual
bool isVirtual() const
Check if particle is virtual.
Definition: MCParticle.h:568
Belle2::MCParticle::getDaughters
std::vector< Belle2::MCParticle * > getDaughters() const
Get vector of all daughter particles, empty vector if none.
Definition: MCParticle.cc:51
Belle2::MCParticle::setEnergy
void setEnergy(float energy)
Set energy.
Definition: MCParticle.h:380
Belle2::Unit::mm
static const double mm
[millimeters]
Definition: Unit.h:80
Belle2::MCParticle::get4Vector
TLorentzVector get4Vector() const
Return 4Vector of particle.
Definition: MCParticle.h:218
Belle2::MCParticle::getLifetime
float getLifetime() const
Return the lifetime in ns.
Definition: MCParticle.h:188
Belle2::MCParticle::getMomentum
TVector3 getMomentum() const
Return momentum.
Definition: MCParticle.h:209
Belle2::MCParticle::setMass
void setMass(float mass)
Set particle mass.
Definition: MCParticle.h:374
Belle2::MCParticleGraph::GraphParticle::comesFrom
void comesFrom(GraphParticle &mother)
Tells the graph that this particle is a decay product of mother.
Definition: MCParticleGraph.h:116
Belle2::MCParticle
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:43
Belle2::StoreArray< MCParticle >
Belle2::MCParticle::setDecayTime
void setDecayTime(float time)
Set decay time.
Definition: MCParticle.h:398
Belle2::MCParticleGraph::clear
void clear()
Reset particles and decay information to make the class reusable.
Definition: MCParticleGraph.h:298
Belle2::MCParticle::setMomentum
void setMomentum(const TVector3 &momentum)
Set particle momentum.
Definition: MCParticle.h:425
Belle2::MCParticle::setValidVertex
void setValidVertex(bool valid)
Set indication wether vertex and time information is valid or just default.
Definition: MCParticle.h:386
Belle2::StoreArray::getEntries
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:226
Belle2::Simulation::MCParticleGenerator::addParticle
void addParticle(const MCParticle &mcParticle, G4Event *event, G4PrimaryParticle *lastG4Mother, int motherIndex, bool useTime)
Takes a MCParticle and creates a primary particle for Geant4.
Definition: MCParticleGenerator.cc:121
Belle2::MCParticleGraph::GraphParticle
Class to represent Particle data in graph.
Definition: MCParticleGraph.h:86
Belle2::MCParticle::getProductionTime
float getProductionTime() const
Return production time in ns.
Definition: MCParticle.h:170
Belle2::MCParticleGraph::GraphParticle::setFirstDaughter
void setFirstDaughter(int daughter)
Set the 1-based index of the first daughter, 0 means no daughters.
Definition: MCParticleGraph.h:125