11 #include <generators/modules/GeneratedVertexDisplacerModule.h>
12 #include <framework/core/ModuleParam.templateDetails.h>
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");
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].",
56 B2DEBUG(0,
"Initialize GeneratedVertexDisplacerModule");
57 m_mcparticles.isRequired(m_particleList);
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.");
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).");
73 if (m_mcparticles.getEntries() < 1) {
74 B2WARNING(
"MC particle container empty. Calling next event.");
78 for (
int mcp_index = 0; mcp_index < m_mcparticles.getEntries(); mcp_index++) {
84 int mcp_pdg = mcp.
getPDG();
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));
102 TLorentzVector* displacementVector =
new TLorentzVector();
103 getDisplacement(particle, lifetime, *displacementVector);
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());
110 particle.setDecayVertex(newDecayVertex);
111 particle.setDecayTime(particle.getProductionTime() + displacementVector->T());
112 particle.setValidVertex(
true);
115 if (particle.getNDaughters()) {
116 displaceDaughter(*displacementVector, particle.getDaughters());
119 delete displacementVector;
125 for (
unsigned int daughter_index = 0; daughter_index < daughters.size(); daughter_index++) {
127 MCParticle* daughter_mcp = daughters.at(daughter_index);
131 MCParticle& mcp = *m_mcparticles[daughter_mcpArrayIndex];
138 displaceDaughter(motherDisplacementVector, daughter_mcp->
getDaughters());
145 const MCParticle& particle,
float lifetime, TLorentzVector& displacement)
147 TLorentzVector fourVector_mcp = particle.get4Vector();
148 float decayLength_mcp = 0;
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);
154 decayLength_mcp = -1 * std::log(gRandom->Uniform(0, 1.)) * lifetime * fourVector_mcp.Gamma() * fourVector_mcp.Beta();
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;
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);