Belle II Software development
PrintMCParticlesModule.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <analysis/modules/PrintMCParticles/PrintMCParticlesModule.h>
10#include <analysis/dataobjects/StringWrapper.h>
11
12#include <mdst/dataobjects/MCParticle.h>
13
14#include <framework/logging/LogConnectionConsole.h>
15
16#include <boost/format.hpp>
17#include <boost/algorithm/string.hpp>
18
19#include <TDatabasePDG.h>
20#include <Math/Vector3D.h>
21
22using namespace Belle2;
23
24//-----------------------------------------------------------------
25// Register the Module
26//-----------------------------------------------------------------
27REG_MODULE(PrintMCParticles);
28
29//-----------------------------------------------------------------
30// Implementation
31//-----------------------------------------------------------------
32
33namespace {
35 std::string statusToString(const MCParticle& mc)
36 {
37 static std::map<int, std::string> names {
38 {MCParticle::c_PrimaryParticle, "PrimaryParticle"},
39 {MCParticle::c_StableInGenerator, "StableInGenerator"},
40 {MCParticle::c_LeftDetector, "LeftDetector"},
41 {MCParticle::c_StoppedInDetector, "StoppedInDetector"},
42 {MCParticle::c_IsVirtual, "IsVirtual"},
43 {MCParticle::c_Initial, "Initial"},
44 {MCParticle::c_IsISRPhoton, "IsISRPhoton"},
45 {MCParticle::c_IsFSRPhoton, "IsFSRPhoton"},
46 {MCParticle::c_IsPHOTOSPhoton, "IsPHOTOSPhoton"},
47 };
48 std::vector<std::string> set;
49 for (auto&& [status, name] : names) {
50 if (mc.hasStatus(status)) set.emplace_back(name);
51 }
52 return boost::join(set, ", ");
53 }
54
56 std::string g4ProcessToName(int process)
57 {
58 static std::map<int, std::string> translation{
59 {1, "CoulombScattering"},
60 {2, "Ionisation"},
61 {3, "Bremsstrahlung"},
62 {4, "PairProdByCharged"},
63 {5, "Annihilation"},
64 {6, "AnnihilationToMuMu"},
65 {7, "AnnihilationToHadrons"},
66 {8, "NuclearStopping"},
67 {9, "ElectronGeneralProcess"},
68 {10, "MultipleScattering"},
69 {11, "Rayleigh"},
70 {12, "PhotoElectricEffect"},
71 {13, "ComptonScattering"},
72 {14, "GammaConversion"},
73 {15, "GammaConversionToMuMu"},
74 {16, "GammaGeneralProcess"},
75 {17, "PositronGeneralProcess"},
76 {18, "AnnihilationToTauTau"},
77 {21, "Cerenkov"},
78 {22, "Scintillation"},
79 {23, "SynchrotronRadiation"},
80 {24, "TransitionRadiation"},
81 {25, "SurfaceReflection"},
82 {40, "DarkBremsstrahlung"},
83 {49, "MuonPairProdByCharged"},
84 {111, "HadronElastic"},
85 {116, "NeutronGeneral"},
86 {121, "HadronInelastic"},
87 {131, "Capture"},
88 {132, "MuAtomicCapture"},
89 {141, "Fission"},
90 {151, "HadronAtRest"},
91 {152, "LeptonAtRest"},
92 {161, "ChargeExchange"},
93 // {210, "RadioactiveDecay"},
94 {201, "Decay"},
95 {202, "Decay_WithSpin"},
96 {203, "Decay_PionMakeSpin"},
97 {210, "Decay_Radioactive"},
98 {211, "Decay_Unknown"},
99 {221, "Decay_MuAtom "},
100 {231, "Decay_External"},
101 {310, "EMDissociation"},
102 };
103 if (auto it = translation.find(process); it != translation.end()) {
104 return it->second;
105 }
106 return "[Unknown process]";
107 }
108}
109
111{
112 //Set module properties
113 setDescription("Print an MCParticle List");
115
116 //Parameter definition
117 addParam("storeName", m_particleList, "Name of the StoreArray to print", m_particleList);
118 addParam("onlyPrimaries", m_onlyPrimaries, "Show only primary particles", true);
119 addParam("maxLevel", m_maxLevel, "Show only up to specified depth level, -1 means no limit", -1);
120 addParam("showVertices", m_showVertices, "Show also the particle production vertices and times", false);
121 addParam("showMomenta", m_showMomenta, "Show also the particle momenta", false);
122 addParam("showProperties", m_showProperties, "Show the basic particle properties", false);
123 addParam("showStatus", m_showStatus, "Show extendend status information of the particle", false);
124 addParam("suppressPrint", m_suppressPrint, "Suppress print the information", false);
125}
126
128{
129 m_mcparticles.isRequired(m_particleList);
130 m_stringWrapper.registerInDataStore("MCDecayString");
131}
132
134{
135 // clear any previous outputs
136 m_output.str(std::string());
137
138 if (not m_particleList.empty()) {
139 m_output << "Content from MCParticle list '" << m_mcparticles.getName() << "'" << std::endl;
140 } else {
141 m_output << "Content of MCParticle list" << std::endl;
142 }
143
144 //Loop over the top level particles (no mother particle exists)
145 std::vector<MCParticle*> first_gen;
146 for (MCParticle& mc : m_mcparticles) {
147 if (mc.getMother() != nullptr) continue;
148 first_gen.emplace_back(&mc);
149 }
150 filterPrimaryOnly(first_gen);
151 printTree(first_gen);
152
153 if (!m_suppressPrint)
154 B2INFO(m_output.str());
155
156 if (not m_stringWrapper.isValid())
157 m_stringWrapper.create();
158
159 m_stringWrapper->setString(std::string(m_output.str()));
160}
161
162void PrintMCParticlesModule::filterPrimaryOnly(std::vector<MCParticle*>& particles) const
163{
164 if (!m_onlyPrimaries or particles.empty()) return;
165 particles.erase(std::remove_if(particles.begin(), particles.end(), [](MCParticle * mc) {
166 return not mc->hasStatus(MCParticle::c_PrimaryParticle);
167 }), particles.end());
168}
169
170void PrintMCParticlesModule::printTree(const std::vector<MCParticle*>& particles, int level, const std::string& indent)
171{
172 //If we show extra content make the particle name and pdg code bold if supported
173 //And if we also show secondaries make those red to distinguish
174 const bool useColor = LogConnectionConsole::terminalSupportsColors(STDOUT_FILENO);
175 const bool anyExtraInfo = m_showProperties or m_showMomenta or m_showVertices or m_showStatus;
176 std::string colorStart[] = {"", ""};
177 std::string colorEnd = "";
178 if (useColor) {
179 // make it bold if we have extra info, otherwise no need
180 colorStart[0] = anyExtraInfo ? "\x1b[31;1m" : "\x1b[31m";
181 colorStart[1] = anyExtraInfo ? "\x1b[1m" : "";
182 colorEnd = "\x1b[m";
183 }
184
185 TDatabasePDG* pdb = TDatabasePDG::Instance();
186
187 // ok, go over all particles and print them
188 const auto num_particles = particles.size();
189 size_t particle_index{0};
190 for (const auto* mc : particles) {
191 // are we the last in the list? if so we need a different box char
192 const bool last = ++particle_index == num_particles;
193 m_output << indent << (last ? "╰" : "├");
194 m_output << (mc->hasStatus(MCParticle::c_PrimaryParticle) ? "── " : "╶╶ ");
195 // now print the name and pdg code. Optionally with color
196 TParticlePDG* pdef = pdb->GetParticle(mc->getPDG());
197 m_output << colorStart[mc->hasStatus(MCParticle::c_PrimaryParticle)] << (pdef ? pdef->GetTitle() : "[UNKNOWN]");
198 m_output << " (" << mc->getPDG() << ")" << colorEnd;
199
200 // Particle itself is done but we might need to show more things like properties and daughters so ...
201 auto daughters = mc->getDaughters();
202 filterPrimaryOnly(daughters);
203 const bool showDaughters = (not daughters.empty()) and (m_maxLevel <= 0 or level < m_maxLevel);
204
205 // new indent level for any information related to this particle.
206 const std::string newIndent = indent + (last ? " " : "│") + " ";
207 // and the indent level for any properties associated to this particle
208 const std::string propIndent = "\n" + newIndent + (showDaughters ? "│ " : "");
209 // limited by level restriction? add an indicator that something is skipped
210 if ((not showDaughters) and (not daughters.empty())) {
211 m_output << " → …";
212 }
213 // Now show whatever the user wants
214 if (m_showProperties) {
215 m_output << propIndent;
216 m_output << boost::format("mass=%.3g energy=%.3g charge=%d") % mc->getMass() % mc->getEnergy() % mc->getCharge();
217 if (not mc->hasStatus(MCParticle::c_LeftDetector)) m_output << boost::format(" lifetime=%.3g") % mc->getLifetime();
218 }
219 if (m_showMomenta) {
220 const ROOT::Math::XYZVector& p = mc->getMomentum();
221 m_output << propIndent;
222 m_output << boost::format("p=(%.3g, %.3g, %.3g) |p|=%.3g") % p.X() % p.Y() % p.Z() % p.R();
223 }
224 if (m_showVertices) {
225 const ROOT::Math::XYZVector& v = mc->getVertex();
226 m_output << propIndent;
227 m_output << boost::format("production vertex=(%.3g, %.3g, %.3g), time=%.3g") % v.X() % v.Y() % v.Z() % mc->getProductionTime();
228 }
229 if (m_showStatus) {
230 m_output << propIndent;
231 m_output << "status flags=" << statusToString(*mc);
232 if (not mc->hasStatus(MCParticle::c_PrimaryParticle)) {
233 m_output << propIndent;
234 m_output << "creation process=" << g4ProcessToName(mc->getSecondaryPhysicsProcess());
235 }
236 m_output << propIndent;
237 m_output << "list index=" << mc->getIndex();
238 }
239 // and if we have shown any extra info and have daughters leave a blank line
240 if (showDaughters and anyExtraInfo) m_output << propIndent;
241 // now we're done with the particle ...
242 m_output << std::endl;
243 // but maybe we need to show the daughters
244 if (showDaughters) {
245 printTree(daughters, level + 1, newIndent);
246 }
247 // and if we show any extra info also leave blank line to sibling
248 if (anyExtraInfo and not last) m_output << newIndent << std::endl;
249 }
250}
static bool terminalSupportsColors(int fileDescriptor)
Returns true if the given file descriptor is a tty and supports colors.
A Class to store the Monte Carlo particle information.
Definition MCParticle.h:32
@ c_IsFSRPhoton
bit 7: Particle is from final state radiation
Definition MCParticle.h:61
@ c_Initial
bit 5: Particle is initial such as e+ or e- and not going to Geant4
Definition MCParticle.h:57
@ c_IsPHOTOSPhoton
bit 8: Particle is an radiative photon from PHOTOS
Definition MCParticle.h:63
@ c_PrimaryParticle
bit 0: Particle is primary particle.
Definition MCParticle.h:47
@ c_LeftDetector
bit 2: Particle left the detector (the simulation volume).
Definition MCParticle.h:51
@ c_IsVirtual
bit 4: Particle is virtual and not going to Geant4.
Definition MCParticle.h:55
@ c_StableInGenerator
bit 1: Particle is stable, i.e., not decaying in the generator.
Definition MCParticle.h:49
@ c_StoppedInDetector
bit 3: Particle was stopped in the detector (the simulation volume).
Definition MCParticle.h:53
@ c_IsISRPhoton
bit 6: Particle is from initial state radiation
Definition MCParticle.h:59
void setDescription(const std::string &description)
Sets the description of the module.
Definition Module.cc:214
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
Definition Module.cc:208
Module()
Constructor.
Definition Module.cc:30
@ c_ParallelProcessingCertified
This module can be run in parallel processing mode safely (All I/O must be done through the data stor...
Definition Module.h:80
@ c_TerminateInAllProcesses
When using parallel processing, call this module's terminate() function in all processes().
Definition Module.h:83
std::stringstream m_output
Buffer to keep all the output while building the tree.
virtual void initialize() override
init.
bool m_showProperties
Show remaining properties.
virtual void event() override
Method is called for each event.
StoreArray< MCParticle > m_mcparticles
store array for the MCParticles
bool m_showStatus
Show extended status information.
int m_maxLevel
Show only up to specified depth level.
bool m_showVertices
Show particle production vertices.
void filterPrimaryOnly(std::vector< MCParticle * > &particles) const
if m_onlyPrimary is set remove all non primary particles from the given vector inplace
StoreObjPtr< StringWrapper > m_stringWrapper
string wrapper to store the MCDecayString
void printTree(const std::vector< MCParticle * > &particles, int level=1, const std::string &indent="")
Loops recursively over the MCParticle list and prints information about each particle.
bool m_showMomenta
Show particle momenta.
std::string m_particleList
The name of the MCParticle collection.
bool m_suppressPrint
Suppress print the information.
bool m_onlyPrimaries
Print only primary particles.
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition Module.h:559
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition Module.h:649
Abstract base class for different kinds of events.