9#include <simulation/kernel/MCParticleGenerator.h>
11#include <mdst/dataobjects/MCParticle.h>
12#include <framework/gearbox/Unit.h>
13#include <framework/logging/Logger.h>
14#include <framework/datastore/StoreArray.h>
15#include <simulation/kernel/UserInfo.h>
16#include <geometry/GeometryManager.h>
18#include <G4ParticleTable.hh>
19#include <G4VUserPrimaryParticleInformation.hh>
20#include <G4ParticleDefinition.hh>
21#include <G4ParticleTypes.hh>
22#include <G4VPhysicalVolume.hh>
23#include <G4LogicalVolume.hh>
25#include <CLHEP/Vector/LorentzVector.h>
26#include <G4PrimaryParticle.hh>
31using namespace Belle2::Simulation;
39 m_topSolid = topVolume->GetLogicalVolume()->GetSolid();
67 for (
int iPart = 0; iPart < nPart; iPart++) {
69 if (currParticle->
getMother() != NULL)
continue;
80 ROOT::Math::XYZVector mcProdVtx = p.getProductionVertex();
81 G4ThreeVector pos(mcProdVtx.X() /
Unit::mm * CLHEP::mm,
82 mcProdVtx.Y() /
Unit::mm * CLHEP::mm,
83 mcProdVtx.Z() /
Unit::mm * CLHEP::mm);
84 productionTimeShift = 0;
89 ROOT::Math::PxPyPzEVector mcPartMom4 = p.get4Vector();
90 CLHEP::HepLorentzVector mom(mcPartMom4.X() /
Unit::MeV * CLHEP::MeV,
94 const auto dir = mom.vect().unit();
95 double distance =
m_topSolid->DistanceToIn(pos, dir);
97 if (distance == kInfinity) {
98 B2DEBUG(10,
"Particle starts outside simulation volume and isn't intersecting, skipping");
102 const double speed = (p.getMass() == 0 ? mom.beta() : 1) * CLHEP::c_light;
103 const double flightLength = speed * p.getLifetime();
104 if (distance > flightLength) {
105 B2DEBUG(10,
"Particle starts outside simulation volume and doesn't reach it before decaying, skipping");
108 pos += dir * distance;
109 productionTimeShift = (distance / speed);
111 return new G4PrimaryVertex(pos.x(), pos.y(), pos.z(), p.getProductionTime() + productionTimeShift);
116 G4PrimaryParticle* lastG4Mother,
120 G4PrimaryParticle* g4Mother = lastG4Mother;
131 G4ParticleDefinition* pdef = NULL;
133 pdef = G4OpticalPhoton::OpticalPhotonDefinition();
135 pdef = G4ParticleTable::GetParticleTable()->FindParticle(mcParticle.
getPDG());
136 if (pdef == NULL && neutral_llp) {
137 pdef = G4ParticleTable::GetParticleTable()->FindParticle(
"LongLivedNeutralParticle");
144 B2WARNING(
"PDG code " << mcParticle.
getPDG() <<
" unknown to Geant4 - particle skipped.");
167 G4PrimaryParticle* newPart =
nullptr;
169 G4PrimaryVertex* vertex =
nullptr;
172 double productionTimeShift{0};
176 if (addToG4 and lastG4Mother ==
nullptr) {
179 if (not vertex) addToG4 =
false;
184 const ROOT::Math::PxPyPzEVector mcPartMom4 = graphParticle.
get4Vector();
185 newPart =
new G4PrimaryParticle(pdef,
189 mcPartMom4.E() /
Unit::MeV * CLHEP::MeV);
192 ROOT::Math::XYZVector polarization = graphParticle.
getDecayVertex();
193 newPart->SetPolarization(polarization.X(), polarization.Y(), polarization.Z());
195 newPart->SetUserInformation(
new ParticleInfo(graphParticle));
202 double propertime = (graphParticle.
getLifetime() - productionTimeShift) / mcPartMom4.Gamma();
203 newPart->SetProperTime(propertime);
207 lastG4Mother->SetDaughter(newPart);
213 vertex->SetPrimary(newPart);
214 event->AddPrimaryVertex(vertex);
223 B2DEBUG(100,
"The particle " << mcParticle.
getIndex() <<
" (PDG " << mcParticle.
getPDG() <<
") was not added to Geant4");
229 addParticle(*daughter, event, g4Mother, currMotherIndex, useTime);
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.
Class to build, validate and sort a particle decay chain.
A Class to store the Monte Carlo particle information.
float getEnergy() const
Return particle energy in GeV.
void setDecayTime(float time)
Set decay time.
void setMass(float mass)
Set particle mass.
void setDecayVertex(const ROOT::Math::XYZVector &vertex)
Set decay vertex.
int getIndex() const
Get 1-based index of the particle in the corresponding MCParticle list.
ROOT::Math::XYZVector getDecayVertex() const
Return decay vertex.
float getLifetime() const
Return the lifetime in ns.
void setEnergy(float energy)
Set energy.
std::vector< Belle2::MCParticle * > getDaughters() const
Get vector of all daughter particles, empty vector if none.
float getMass() const
Return the particle mass in GeV.
ROOT::Math::XYZVector getProductionVertex() const
Return production vertex position.
float getDecayTime() const
Return the decay time in ns.
void setValidVertex(bool valid)
Set indication whether vertex and time information is valid or just default.
void setProductionVertex(const ROOT::Math::XYZVector &vertex)
Set production vertex position.
bool hasValidVertex() const
Indication whether vertex and time information is useful or just default.
unsigned int getStatus(unsigned short int bitmask=USHRT_MAX) const
Return status code of particle.
int getLastDaughter() const
Get 1-based index of last daughter, 0 if no daughters.
float getCharge() const
Return the particle charge defined in TDatabasePDG.
ROOT::Math::PxPyPzEVector get4Vector() const
Return 4Vector of particle.
void setPDG(int pdg)
Set PDG code of the particle.
int getPDG() const
Return PDG code of particle.
float getProductionTime() const
Return production time in ns.
ROOT::Math::XYZVector getMomentum() const
Return momentum.
void setMomentum(const ROOT::Math::XYZVector &momentum)
Set particle momentum.
void setStatus(unsigned short int status)
Set Status code for the particle.
int getFirstDaughter() const
Get 1-based index of first daughter, 0 if no daughters.
void setProductionTime(float time)
Set production time.
std::string m_mcCollectionName
Name of the MCParticle collection.
virtual ~MCParticleGenerator()
The destructor of the MCParticleGenerator class.
MCParticleGraph & m_mcParticleGraph
Reference to the MCParticle graph.
void addParticle(const MCParticle &mcParticle, G4Event *event, G4PrimaryParticle *lastG4Mother, int motherIndex, bool useTime)
Takes a MCParticle and creates a primary particle for Geant4.
virtual void GeneratePrimaryVertex(G4Event *event)
Create G4 primary particles from MCParticle list.
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.
MCParticleGenerator(const std::string &mcCollectionName, MCParticleGraph &mcParticleGraph)
The constructor of the MCParticleGenerator class.
Accessor to arrays stored in the data store.
int getEntries() const
Get the number of objects in the array.
static const double mm
[millimeters]
static const double eV
[electronvolt]
static const double MeV
[megaelectronvolt]
G4VPhysicalVolume * getTopVolume()
Return a pointer to the top volume.
static GeometryManager & getInstance()
Return a reference to the instance.
bool isInitial() const
Check if particle is an initial particle such as ISR.
MCParticle * getMother() const
Returns a pointer to the mother particle.
bool isVirtual() const
Check if particle is virtual.
Abstract base class for different kinds of events.