Belle II Software light-2406-ragdoll
ParticleMCDecayStringModule.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/ParticleMCDecayString/ParticleMCDecayStringModule.h>
10
11#include <mdst/dataobjects/MCParticle.h>
12
13#include <framework/logging/Logger.h>
14#include <framework/pcore/ProcHandler.h>
15#include <framework/utilities/RootFileCreationManager.h>
16
17#include <string>
18#include <vector>
19
20#include <boost/algorithm/string.hpp>
21
22using namespace Belle2;
23
24REG_MODULE(ParticleMCDecayString);
25
26ParticleMCDecayStringModule::ParticleMCDecayStringModule() : Module(), m_tree("", DataStore::c_Persistent), m_hashset("",
27 DataStore::c_Persistent), m_decayHash(0.0), m_decayHashExtended(0.0)
28{
29 setDescription("Creates the Monte Carlo decay string of a Particle and its daughters. "
30 "The MC decay string of the particle is hashed and saved as a 32bit pattern in the extra info field decayHash of the particle. "
31 "The MC decay string of the particle + its daughters is hashed as well and saved as another 32bit pattern in the extra info field decayHashExtended of the particle. "
32 "The mapping hash <-> MC decay string in saved in a TTree by this module. "
33 "The 32bit pattern must be saved as a float (because our extra info field, variable manager and ntuple output only supports float) "
34 "but they just represent 32 bits of a hash! "
35 "The MC decay string can also be stored in an analysis ROOT file using the MCDecayString NtupleTool. "
36 "Details on the MC decay string format can be found here: `MCDecayString`");
38 addParam("listName", m_listName, "Particles from these ParticleList are used as input.");
39 addParam("fileName", m_fileName, "Filename in which the hash strings are saved, if empty the strings are not saved",
40 std::string(""));
41 addParam("treeName", m_treeName, "Tree name in which the hash strings are saved", std::string("hashtable"));
42 addParam("conciseString", m_useConciseString, "If set to true, the code will use a more concise format for the string.", false);
43 addParam("identifiers", m_identifiers, "Identifiers used to identify particles in the concise format.",
44 std::string("abcdefghijklmnopqrstuvwxyz"));
45}
46
48{
49 m_pList.isRequired(m_listName);
50
51 //This might not work for non-default names of Particle array:
53
54 m_stringWrapperArray.registerInDataStore();
56
57
58 // Initializing the output root file
59 if (m_fileName != "") {
61 if (!m_file) {
62 B2WARNING("Could not create file " << m_fileName);
63 return;
64 }
65
66 m_file->cd();
67
68 // check if TTree with that name already exists
69 if (m_file->Get(m_treeName.c_str())) {
70 B2WARNING("Tree with this name already exists: " << m_fileName);
71 return;
72 }
73
75 m_tree.construct(m_treeName.c_str(), "Decay Hash Map");
76 m_tree->get().Branch("decayHash", &m_decayHash);
77 m_tree->get().Branch("decayHashExtended", &m_decayHashExtended);
78 m_tree->get().Branch("decayString", &m_decayString);
79 m_tree->get().SetBasketSize("*", 1600);
80 m_tree->get().SetCacheSize(100000);
81 }
82
83 m_hashset.registerInDataStore(m_fileName + m_treeName + "_hashset", DataStore::c_DontWriteOut);
84 m_hashset.construct();
85
86}
87
89{
90
91 for (unsigned iParticle = 0; iParticle < m_pList->getListSize(); ++iParticle) {
92 Particle* particle = m_pList->getParticle(iParticle);
93
94 const std::string decayString = getMCDecayStringFromMCParticle(particle->getRelatedTo<MCParticle>());
95 std::string decayStringExtended = getDecayString(*particle); //removed const to allow string to be modified to a different format.
96
97 if (m_useConciseString) {convertToConciseString(decayStringExtended);}
98
99 uint32_t decayHash = m_hasher(decayString);
100 uint32_t decayHashExtended = m_hasher(decayStringExtended);
101
102 uint64_t m_decayHashFull = decayHash;
103 m_decayHashFull <<= 32;
104 m_decayHashFull += decayHashExtended;
105
106 // Convert unsigned int decay hash into a float keeping the same bit pattern
107 assert(sizeof(float) == sizeof(uint32_t));
108
109 union convert {
110 uint32_t i;
111 float f;
112 };
113 convert bitconverter;
114
115 bitconverter.i = decayHash;
116 m_decayHash = bitconverter.f;
118
119 // cppcheck doesn't like this use of union and throws warnings
120 // cppcheck-suppress redundantAssignment
121 bitconverter.i = decayHashExtended;
122 m_decayHashExtended = bitconverter.f;
124
125 m_decayString = decayStringExtended;
126
127 StringWrapper* stringWrapper = m_stringWrapperArray.appendNew();
128 particle->addRelationTo(stringWrapper);
129 stringWrapper->setString(m_decayString);
130
131 auto it = m_hashset->get().find(m_decayHashFull);
132 if (it == m_hashset->get().end()) {
133 m_hashset->get().insert(m_decayHashFull);
134
135 if (m_tree.isValid()) {
136 m_tree->get().Fill();
137 }
138 }
139
140 }
141}
142
144{
146 if (m_tree.isValid()) {
147 B2INFO("Writing NTuple " << m_treeName);
148 TDirectory::TContext directoryGuard(m_file.get());
149 m_tree->write(m_file.get());
150
151 const bool writeError = m_file->TestBit(TFile::kWriteError);
152 m_file.reset();
153 if (writeError) {
154 B2FATAL("A write error occurred while saving '" << m_fileName << "', please check if enough disk space is available.");
155 }
156 }
157 }
158}
159
160
162{
163 const MCParticle* mcPMother = mcP->getMother();
164 if (mcPMother == nullptr) {
165 return mcP;
166 } else {
167 return getInitialParticle(mcPMother);
168 }
169}
170
172bool isFSP(int pdg)
173{
174 switch (abs(pdg)) {
175 case 211: //pi^+
176 case 321: //K^+
177 case 11: //e
178 case 12: //nu_e
179 case 13: //mu
180 case 14: //nu_mu
181 case 16: //nu_tau
182 case 22: //gamma
183 case 310: //K_S
184 case 130: //K_L
185 case 2112: //n
186 case 2212: //p
187 return true;
188 default:
189 return false;
190 }
191}
192
194{
195
196 std::string output;
197 output += getDecayStringFromParticle(&p) + " | ";
198 output += getMCDecayStringFromParticle(&p);
199 return output;
200
201}
202
204{
205
206 std::string output = " ";
207
208 output += std::to_string(p->getPDGCode());
209
210 if (not isFSP(p->getPDGCode())) {
211 output += " (-->";
212 for (auto daughter : p->getDaughters()) {
213 output += getDecayStringFromParticle(daughter);
214 }
215 output += ")";
216 }
217
218 return output;
219
220}
221
223{
224
225 std::string output;
226
227 output = getMCDecayStringFromMCParticle(p->getRelatedTo<MCParticle>());
228 // Some FSPs can have daughters, e.g. converted Photons and K-Shorts
229 if (not isFSP(p->getPDGCode())) {
230 for (auto& daughter : p->getDaughters()) {
231 output += " | " + getMCDecayStringFromParticle(daughter);
232 }
233 }
234
235 return output;
236
237}
238
240{
241
242 if (mcPMatched == nullptr)
243 return "(No match)";
244
245 // TODO Performance can be optimized, this mcPMother does not change during the construction
246 const MCParticle* mcPMother = getInitialParticle(mcPMatched);
247
248 std::string decayString = buildMCDecayString(mcPMother, mcPMatched);
249
250 if (mcPMatched->getPDG() == 10022)
251 return decayString + " (Virtual gamma match)";
252 return decayString;
253}
254
255
256std::string ParticleMCDecayStringModule::buildMCDecayString(const MCParticle* mcPMother, const MCParticle* mcPMatched)
257{
258
259 std::stringstream ss;
260 ss << " ";
261 if (mcPMother->getArrayIndex() == mcPMatched->getArrayIndex()) {
262 ss << "^";
263 }
264
265 ss << mcPMother->getPDG();
266
267 if (not isFSP(mcPMother->getPDG())) {
268 ss << " (-->";
269 for (auto daughter : mcPMother->getDaughters()) {
270 ss << buildMCDecayString(daughter, mcPMatched);
271 }
272 ss << ")";
273 }
274
275 return ss.str();
276}
277
279{
280
281 std::vector<std::string> decayStrings;
282 boost::split(decayStrings, string, boost::is_any_of("|"));
283
284 if (decayStrings.empty()) {
285 B2WARNING("ParticleMCDecayStringModule: unable to convert decay string to concise format.");
286 return;
287 }
288
289 unsigned int nParticles(decayStrings.size() - 1);
290 if (nParticles > m_identifiers.size()) {
291 B2WARNING("ParticleMCDecayStringModule: not enough identifiers have been specified to use the concise string format:"
292 << std::endl << "Number of particles in your decay mode = " << nParticles << std::endl
293 << "Available identifiers: " << m_identifiers << std::endl
294 << "Standard format will be used instead.");
295 return;
296 }
297
298 //Find positions of carets in original strings, store them, and then erase them.
299 std::string mode("");
300 std::vector<int> caretPositions;
301 for (auto& decayString : decayStrings) {
302 std::string thisString(decayString);
303 if ("" == mode) {
304 mode = thisString;
305 continue;
306 }
307
308 int caretPosition(thisString.find('^')); // -1 if no match.
309 caretPositions.push_back(caretPosition);
310 if (caretPosition > -1) {
311 decayString.erase(caretPosition, 1);
312 }
313 }
314
315 //Check if all of the decay strings are the same (except for No matches):
316 std::string theDecayString("");
317 for (auto thisString : decayStrings) {
318 if (thisString == mode) {continue;}
319
320 //last decay string does not have a space at the end, don't want this to stop a match.
321 char finalChar(thisString.back());
322 if (finalChar != ' ') {thisString = thisString + " ";}
323
324 if (" (No match) " != thisString) {
325 if ("" == theDecayString) {
326 theDecayString = thisString;
327 } else {
328 if (theDecayString != thisString) {
329 //TODO: add string format if multiple decay strings are present (e.g. pile-up events).
330 return;
331 }
332 }
333 }
334 }
335
336 std::string modifiedString(theDecayString);
337
338 //insert identifiers in positions where carets were:
339 int nStrings(caretPositions.size());
340 for (int iString(0); iString < nStrings; ++iString) {
341 std::string identifier(m_identifiers.substr(iString, 1));
342 int insertPosition(caretPositions.at(iString));
343 if (insertPosition > -1) {
344 for (int jString(0); jString < iString; ++jString) {
345 if (caretPositions.at(jString) > -1 && caretPositions.at(jString) <= caretPositions.at(iString)) {
346 ++insertPosition;
347 }
348 }
349 modifiedString.insert(insertPosition, identifier);
350 }
351 }
352
353 modifiedString = mode + "|" + modifiedString;
354
355 //add a list of the unmatched particles at the end of the string:
356 bool noMatchStringAdded(false);
357 for (int iString(0); iString < nStrings; ++iString) {
358 int insertPosition(caretPositions.at(iString));
359 if (-1 == insertPosition) {
360 if (!noMatchStringAdded) {
361 modifiedString += " | No match: ";
362 noMatchStringAdded = true;
363 }
364 modifiedString += m_identifiers.substr(iString, 1);
365 }
366 }
367
368 string = modifiedString;
369 return;
370}
In the store you can park objects that have to be accessed by various modules.
Definition: DataStore.h:51
@ c_DontWriteOut
Object/array should be NOT saved by output modules.
Definition: DataStore.h:71
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:32
std::vector< Belle2::MCParticle * > getDaughters() const
Get vector of all daughter particles, empty vector if none.
Definition: MCParticle.cc:52
int getArrayIndex() const
Get 0-based index of the particle in the corresponding MCParticle list.
Definition: MCParticle.h:244
int getPDG() const
Return PDG code of particle.
Definition: MCParticle.h:112
Base class for Modules.
Definition: Module.h:72
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
@ 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
float m_decayHash
Decay hash -> The hash of the decay string of the mother particle.
std::string getDecayStringFromParticle(const Particle *p)
get decay string of particle
virtual void initialize() override
Initialize the module.
const std::string c_ExtraInfoName
Name of the extraInfo, which is stored in each Particle.
virtual void event() override
Called for each event.
std::string m_decayString
The complete decay string.
std::string m_listName
Name of the particle list.
virtual void terminate() override
Terminate modules.
std::string m_fileName
Filename in which the hash strings are saved, if empty the strings are not saved.
std::string getMCDecayStringFromMCParticle(const MCParticle *mcPMatched)
get mc decay string from mc particle
std::hash< std::string > m_hasher
Hash function.
bool m_useConciseString
Switch to use concise format for the extended string.
StoreObjPtr< RootMergeable< TTree > > m_tree
ROOT TNtuple containing the saved hashes and strings.
std::shared_ptr< TFile > m_file
ROOT file to store the hashes and strings.
const std::string c_ExtraInfoNameExtended
Name of the extraInfo, which is stored in each Particle.
StoreObjPtr< SetMergeable< std::unordered_set< uint64_t > > > m_hashset
Mergeable unordered set containing the encountered hashes.
StoreArray< StringWrapper > m_stringWrapperArray
StoreArray of StringWrappers.
std::string m_treeName
Tree name in which the hash strings are saved.
std::string m_identifiers
Characters used to identify particles in the concise decay string format (default: alphabet).
void convertToConciseString(std::string &string)
Convert the extended string to a more concise format.
std::string buildMCDecayString(const MCParticle *mcPMother, const MCParticle *mcPMatched)
return decay string for mcPMother, highlight mcPMatched.
const MCParticle * getInitialParticle(const MCParticle *mcP)
search from mcP upwards for a particle that matches specified mother PDG codes.
float m_decayHashExtended
Extended decay hash -> The hash of the decay string of all daughter particles.
std::string getMCDecayStringFromParticle(const Particle *p)
get mc decay string from particle
StoreObjPtr< ParticleList > m_pList
input particle list
std::string getDecayString(const Particle &p)
get the decay string for p.
Class to store reconstructed particles.
Definition: Particle.h:75
void addExtraInfo(const std::string &name, double value)
Sets the user-defined data of given name to the given value.
Definition: Particle.cc:1336
static bool isOutputProcess()
Return true if the process is an output process.
Definition: ProcHandler.cc:232
static bool parallelProcessingUsed()
Returns true if multiple processes have been spawned, false in single-core mode.
Definition: ProcHandler.cc:226
void addRelationTo(const RelationsInterface< BASE > *object, float weight=1.0, const std::string &namedRelation="") const
Add a relation from this object to another object (with caching).
TO * getRelatedTo(const std::string &name="", const std::string &namedRelation="") const
Get the object to which this object has a relation.
bool isRequired(const std::string &name="")
Ensure this array/object has been registered previously.
Accessor to arrays stored in the data store.
Definition: StoreArray.h:113
bool registerRelationTo(const StoreArray< TO > &toArray, DataStore::EDurability durability=DataStore::c_Event, DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut, const std::string &namedRelation="") const
Register a relation to the given StoreArray.
Definition: StoreArray.h:140
This class is a wrapper for strings, such as MCDecayStrings, to allow them to be associated with part...
Definition: StringWrapper.h:23
void setString(const std::string &inputstring)
Set string.
Definition: StringWrapper.h:41
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:560
std::shared_ptr< TFile > getFile(std::string, bool ignoreErrors=false)
Get a file with a specific name, if is does not exist it will be created.
static RootFileCreationManager & getInstance()
Interface for the FileManager.
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
MCParticle * getMother() const
Returns a pointer to the mother particle.
Definition: MCParticle.h:600
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24