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