Belle II Software  release-05-01-25
GeneratedVertexDisplacerModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2019 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Sascha Dreyer, Savino Longo *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <generators/modules/GeneratedVertexDisplacerModule.h>
12 #include <framework/core/ModuleParam.templateDetails.h>
13 
14 
15 
16 using namespace Belle2;
17 using namespace std;
18 
19 //-----------------------------------------------------------------
20 // Register the Module
21 //-----------------------------------------------------------------
22 REG_MODULE(GeneratedVertexDisplacer)
23 
24 //-----------------------------------------------------------------
25 // Implementation
26 //-----------------------------------------------------------------
27 
29 {
30  // Set module properties
31  setDescription(
32  R"DOC(""Takes a list of PDG values and lifetime paramters to displaces the vertex of MCParticles with matching PDG value corresponding to the given lifetime parameter. Can be used betwenerator and the detector simulation.)DOC");
33 
34  // Parameter definitions
35  addParam("lifetimeOption", m_lifetimeOption,
36  "Set the lifetime option, either 'fixed', 'flat' or 'exponential'. Option is set globally.", string("fixed"));
37  addParam("lifetime", m_lifetime,
38  "Set the numerical value of the lifetime c*tau in units of [cm]. Example: (lifetime_1, lifetime_2, ... ,lifetime_N). Lifetime_X is set for PDG_X.",
39  std::vector<float> {0});
40  addParam("pdgVal", m_pdgVals,
41  "PDG values of MCParticles that should be displaced. Subsequent daughters of these will be displaced accordingly. Example: (PDG_1, PDG_2, ... , PDG_N). PDG_X corresponds to lifetime_X.",
42  std::vector<int> { -999});
43  addParam("maxDecayTime", m_maxDecayTime,
44  "Set the maximal allowed decay time in c*tau [cm] for the option 'flat'. The 'exponential' option is not bound. Default 300cm.",
45  static_cast<float>(300));
46  addParam("ctau", m_ctau,
47  "Input unit option. True (default): module expects lifetime as ctau [cm]. False: lifetime as tau [ns].",
48  true);
49 
50 }
51 
52 
54 {
55 
56  B2DEBUG(0, "Initialize GeneratedVertexDisplacerModule");
57  m_mcparticles.isRequired(m_particleList);
58 
59  if (m_pdgVals.size() != m_lifetime.size()) {
60  B2FATAL("List of PDG values and lifetime parameters have different sizes. Specify a lifetime for each PDG value.");
61  }
62 
63  if ((m_lifetimeOption.compare("fixed") != 0) && (m_lifetimeOption.compare("flat") != 0)
64  && (m_lifetimeOption.compare("exponential") != 0)) {
65  B2FATAL("Lifetime option must be 0 (fixed), 1 (flat) or 2 (exponential).");
66  }
67 }
68 
69 
71 {
72 
73  if (m_mcparticles.getEntries() < 1) {
74  B2WARNING("MC particle container empty. Calling next event.");
75  return;
76  }
77 
78  for (int mcp_index = 0; mcp_index < m_mcparticles.getEntries(); mcp_index++) {
79 
80  MCParticle& mcp = *m_mcparticles[mcp_index];
81 
83 
84  int mcp_pdg = mcp.getPDG();
85 
86  for (unsigned int param_index = 0; param_index < m_pdgVals.size(); param_index++) {
87  if (m_pdgVals.at(param_index) == mcp_pdg) {
88  B2DEBUG(0, "Displacing Vertex of particle with pdg_id: " << mcp_pdg << " with lifetime: " << m_lifetime.at(
89  param_index) << " and option: " << m_lifetimeOption);
90  displace(mcp, m_lifetime.at(param_index));
91  B2DEBUG(0, "new vertex (x,y,z) at:" << "(" << mcp.getDecayVertex().X() << "," << mcp.getDecayVertex().Y() << "," <<
92  mcp.getDecayVertex().Z() << ")" << " with decaylength=" << std::sqrt(std::pow(mcp.getDecayVertex().X(),
93  2) + std::pow(mcp.getDecayVertex().Y(), 2) + std::pow(mcp.getDecayVertex().Z(), 2))) ;
94  }
95  }
96  }
97 }
98 
99 
101 {
102  TLorentzVector* displacementVector = new TLorentzVector();
103  getDisplacement(particle, lifetime, *displacementVector);
104 
105  TVector3 newDecayVertex = particle.getProductionVertex();
106  newDecayVertex.SetX(newDecayVertex.X() + displacementVector->X());
107  newDecayVertex.SetY(newDecayVertex.Y() + displacementVector->Y());
108  newDecayVertex.SetZ(newDecayVertex.Z() + displacementVector->Z());
109 
110  particle.setDecayVertex(newDecayVertex);
111  particle.setDecayTime(particle.getProductionTime() + displacementVector->T());
112  particle.setValidVertex(true);
113 
114  // displace direct daughters
115  if (particle.getNDaughters()) {
116  displaceDaughter(*displacementVector, particle.getDaughters());
117  }
118 
119  delete displacementVector;
120 }
121 
122 
123 void GeneratedVertexDisplacerModule::displaceDaughter(TLorentzVector& motherDisplacementVector, std::vector<MCParticle*> daughters)
124 {
125  for (unsigned int daughter_index = 0; daughter_index < daughters.size(); daughter_index++) {
126 
127  MCParticle* daughter_mcp = daughters.at(daughter_index);
128  int daughter_mcpArrayIndex = daughter_mcp->getArrayIndex();
129 
130  // getDaughters returns a copied list, need to change parameters of the original particle in the MCParticle StoreArray.
131  MCParticle& mcp = *m_mcparticles[daughter_mcpArrayIndex];
132  mcp.setProductionVertex(mcp.getProductionVertex() + motherDisplacementVector.Vect());
133  mcp.setProductionTime(mcp.getProductionTime() + motherDisplacementVector.T());
134  mcp.setValidVertex(true);
135 
136  // Displace subsequent daughters
137  if (daughter_mcp->getNDaughters()) {
138  displaceDaughter(motherDisplacementVector, daughter_mcp->getDaughters());
139  }
140  }
141 }
142 
143 
145  const MCParticle& particle, float lifetime, TLorentzVector& displacement)
146 {
147  TLorentzVector fourVector_mcp = particle.get4Vector();
148  float decayLength_mcp = 0;
149  if (!m_ctau) lifetime *= Const::speedOfLight;
150 
151  if (m_lifetimeOption.compare("fixed") == 0) decayLength_mcp = lifetime;
152  else if (m_lifetimeOption.compare("flat") == 0) decayLength_mcp = gRandom->Uniform(0, m_maxDecayTime);
153  else {
154  decayLength_mcp = -1 * std::log(gRandom->Uniform(0, 1.)) * lifetime * fourVector_mcp.Gamma() * fourVector_mcp.Beta();
155 
156  if (!particle.getMass()) {
157  B2WARNING("Displacing a particle with zero mass. Forcing Gamma=1 for the decay time.");
158  decayLength_mcp = -1 * std::log(gRandom->Uniform(0, 1.)) * lifetime;
159  }
160  }
161 
162  float pMag = fourVector_mcp.P();
163  displacement.SetX(decayLength_mcp * fourVector_mcp.X() / pMag);
164  displacement.SetY(decayLength_mcp * fourVector_mcp.Y() / pMag);
165  displacement.SetZ(decayLength_mcp * fourVector_mcp.Z() / pMag);
166  displacement.SetT(decayLength_mcp / (Const::speedOfLight * fourVector_mcp.Beta()));
167 }
168 
169 
171 {
172 
173 }
174 
175 
Belle2::GeneratedVertexDisplacerModule::displaceDaughter
void displaceDaughter(TLorentzVector &motherDisplacementVector, std::vector< MCParticle * > daughters)
Helper function to loop over subsequent daughters and displaces their vertices corresponding to their...
Definition: GeneratedVertexDisplacerModule.cc:123
Belle2::GeneratedVertexDisplacerModule
"Takes a list of PDG values and lifetime paramters to displaces the vertex of MCParticles with matchi...
Definition: GeneratedVertexDisplacerModule.h:46
Belle2::GeneratedVertexDisplacerModule::displace
void displace(MCParticle &particle, float lifetime)
Helper function to displace the mother particles (corresponding to given pdf values).
Definition: GeneratedVertexDisplacerModule.cc:100
Belle2::GeneratedVertexDisplacerModule::event
virtual void event() override
Method is called for each event.
Definition: GeneratedVertexDisplacerModule.cc:70
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::MCParticle::getArrayIndex
int getArrayIndex() const
Get 0-based index of the particle in the corresponding MCParticle list.
Definition: MCParticle.h:255
Belle2::MCParticle::getDecayVertex
TVector3 getDecayVertex() const
Return decay vertex.
Definition: MCParticle.h:230
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2::MCParticle::hasStatus
bool hasStatus(unsigned short int bitmask) const
Return if specific status bit is set.
Definition: MCParticle.h:140
Belle2::Const::speedOfLight
static const double speedOfLight
[cm/ns]
Definition: Const.h:568
Belle2::MCParticle::setProductionTime
void setProductionTime(float time)
Set production time.
Definition: MCParticle.h:392
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::MCParticle::c_Initial
@ c_Initial
bit 5: Particle is initial such as e+ or e- and not going to Geant4
Definition: MCParticle.h:68
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::MCParticle::setProductionVertex
void setProductionVertex(const TVector3 &vertex)
Set production vertex position.
Definition: MCParticle.h:404
Belle2::GeneratedVertexDisplacerModule::getDisplacement
void getDisplacement(const MCParticle &particle, float lifetime, TLorentzVector &displacement)
Helper function to calculate the numerical value of the vertex displacement (x,y,z,...
Definition: GeneratedVertexDisplacerModule.cc:144
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
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:43
Belle2::MCParticle::getNDaughters
int getNDaughters() const
Return number of daughter MCParticles.
Definition: MCParticle.cc:66
Belle2::MCParticle::setValidVertex
void setValidVertex(bool valid)
Set indication wether vertex and time information is valid or just default.
Definition: MCParticle.h:386
Belle2::GeneratedVertexDisplacerModule::terminate
virtual void terminate() override
Terminates the module.
Definition: GeneratedVertexDisplacerModule.cc:170
Belle2::MCParticle::c_PrimaryParticle
@ c_PrimaryParticle
bit 0: Particle is primary particle.
Definition: MCParticle.h:58
Belle2::GeneratedVertexDisplacerModule::initialize
virtual void initialize() override
Register input and output data, initialises the module.
Definition: GeneratedVertexDisplacerModule.cc:53
Belle2::MCParticle::getProductionTime
float getProductionTime() const
Return production time in ns.
Definition: MCParticle.h:170