 |
Belle II Software
release-05-01-25
|
16 #include <simulation/kernel/MCParticleGenerator.h>
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>
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>
34 #include <TLorentzVector.h>
38 using namespace Belle2::Simulation;
41 MCParticleGenerator::MCParticleGenerator(
const string& mcCollectionName,
MCParticleGraph& mcParticleGraph) :
42 G4VPrimaryGenerator(), m_mcCollectionName(mcCollectionName), m_mcParticleGraph(mcParticleGraph)
46 m_topSolid = topVolume->GetLogicalVolume()->GetSolid();
74 for (
int iPart = 0; iPart < nPart; iPart++) {
76 if (currParticle->
getMother() != NULL)
continue;
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;
96 TLorentzVector mcPartMom4 = p.get4Vector();
97 CLHEP::HepLorentzVector mom(mcPartMom4.X() /
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);
104 if (distance == kInfinity) {
105 B2DEBUG(10,
"Particle starts outside simulation volume and isn't intersecting, skipping");
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");
115 pos += dir * distance;
116 productionTimeShift = (distance / speed);
118 return new G4PrimaryVertex(pos.x(), pos.y(), pos.z(), p.getProductionTime() + productionTimeShift);
123 G4PrimaryParticle* lastG4Mother,
127 G4PrimaryParticle* g4Mother = lastG4Mother;
138 G4ParticleDefinition* pdef = NULL;
140 pdef = G4OpticalPhoton::OpticalPhotonDefinition();
142 pdef = G4ParticleTable::GetParticleTable()->FindParticle(mcParticle.
getPDG());
143 if (pdef == NULL && neutral_llp) {
144 pdef = G4ParticleTable::GetParticleTable()->FindParticle(
"LongLivedNeutralParticle");
151 B2WARNING(
"PDG code " << mcParticle.
getPDG() <<
" unknown to Geant4 - particle skipped.");
174 G4PrimaryParticle* newPart =
nullptr;
176 G4PrimaryVertex* vertex =
nullptr;
179 double productionTimeShift{0};
183 if (addToG4 and lastG4Mother ==
nullptr) {
186 if (not vertex) addToG4 =
false;
191 const TLorentzVector mcPartMom4 = graphParticle.
get4Vector();
192 newPart =
new G4PrimaryParticle(pdef,
196 mcPartMom4.E() /
Unit::MeV * CLHEP::MeV);
200 newPart->SetPolarization(polarization.X(), polarization.Y(), polarization.Z());
202 newPart->SetUserInformation(
new ParticleInfo(graphParticle));
209 double propertime = (graphParticle.
getLifetime() - productionTimeShift) / mcPartMom4.Gamma();
210 newPart->SetProperTime(propertime);
214 lastG4Mother->SetDaughter(newPart);
220 vertex->SetPrimary(newPart);
221 event->AddPrimaryVertex(vertex);
230 B2DEBUG(100,
"The particle " << mcParticle.
getIndex() <<
" (PDG " << mcParticle.
getPDG() <<
") was not added to Geant4");
236 addParticle(*daughter, event, g4Mother, currMotherIndex, useTime);
int getFirstDaughter() const
Get 1-based index of first daughter, 0 if no daughters.
int getIndex() const
Get 1-based index of the particle in the corresponding MCParticle list.
float getEnergy() const
Return particle energy in GeV.
float getCharge() const
Return the particle charge defined in TDatabasePDG.
size_t size() const
Return the number of particles in the graph.
float getDecayTime() const
Return the decay time in ns.
virtual ~MCParticleGenerator()
The destructor of the MCParticleGenerator class.
void setIgnore(bool ignore=true)
Set or remove the ignore flag.
Class to build, validate and sort a particle decay chain.
MCParticleGraph & m_mcParticleGraph
Reference to the MCParticle graph.
std::string m_mcCollectionName
Name of the MCParticle collection.
G4VSolid * m_topSolid
Pointer to the top level solid to check if particles are inside the simulation volume.
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.
G4VPhysicalVolume * getTopVolume()
Return a pointer to the top volume.
void setDecayVertex(const TVector3 &vertex)
Set decay vertex.
void setStatus(unsigned short int status)
Set Status code for the particle.
static GeometryManager & getInstance()
Return a reference to the instance.
virtual void GeneratePrimaryVertex(G4Event *event)
Create G4 primary particles from MCParticle list.
bool isInitial() const
Check if particle is an initial particle such as ISR.
TVector3 getDecayVertex() const
Return decay vertex.
bool hasValidVertex() const
Indication whether vertex and time information is useful or just default.
static const double MeV
[megaelectronvolt]
static const double eV
[electronvolt]
unsigned int getStatus(unsigned short int bitmask=USHRT_MAX) const
Return status code of particle.
void setProductionTime(float time)
Set production time.
void setLastDaughter(int daughter)
Set the 1-based index of the last daughter, 0 means no daughters.
Abstract base class for different kinds of events.
void setPDG(int pdg)
Set PDG code of the particle.
TVector3 getProductionVertex() const
Return production vertex position.
int getPDG() const
Return PDG code of particle.
GraphParticle & addParticle()
Add new particle to the graph.
void setProductionVertex(const TVector3 &vertex)
Set production vertex position.
float getMass() const
Return the particle mass in GeV.
MCParticle * getMother() const
Returns a pointer to the mother particle.
int getLastDaughter() const
Get 1-based index of last daughter, 0 if no daughters.
bool isVirtual() const
Check if particle is virtual.
std::vector< Belle2::MCParticle * > getDaughters() const
Get vector of all daughter particles, empty vector if none.
void setEnergy(float energy)
Set energy.
static const double mm
[millimeters]
TLorentzVector get4Vector() const
Return 4Vector of particle.
float getLifetime() const
Return the lifetime in ns.
TVector3 getMomentum() const
Return momentum.
void setMass(float mass)
Set particle mass.
void comesFrom(GraphParticle &mother)
Tells the graph that this particle is a decay product of mother.
A Class to store the Monte Carlo particle information.
void setDecayTime(float time)
Set decay time.
void clear()
Reset particles and decay information to make the class reusable.
void setMomentum(const TVector3 &momentum)
Set particle momentum.
void setValidVertex(bool valid)
Set indication wether vertex and time information is valid or just default.
int getEntries() const
Get the number of objects in the array.
void addParticle(const MCParticle &mcParticle, G4Event *event, G4PrimaryParticle *lastG4Mother, int motherIndex, bool useTime)
Takes a MCParticle and creates a primary particle for Geant4.
Class to represent Particle data in graph.
float getProductionTime() const
Return production time in ns.
void setFirstDaughter(int daughter)
Set the 1-based index of the first daughter, 0 means no daughters.