Belle II Software light-2406-ragdoll
MCDecayFinderModule.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/MCDecayFinder/MCDecayFinderModule.h>
10#include <framework/gearbox/Const.h>
11#include <analysis/DecayDescriptor/ParticleListName.h>
12#include <analysis/utility/EvtPDLUtil.h>
13#include <analysis/utility/MCMatching.h>
14
15#include <string>
16#include <memory>
17
18using namespace std;
19using namespace Belle2;
20
21// Register module in the framework
22REG_MODULE(MCDecayFinder);
23
24MCDecayFinderModule::MCDecayFinderModule() : Module(), m_isSelfConjugatedParticle(false)
25{
26 //Set module properties
27 setDescription("Find decays in MCParticle list matching a given DecayString and create Particles from them.");
28 //Parameter definition
29 addParam("decayString", m_strDecay, "DecayDescriptor string.");
30 addParam("listName", m_listName, "Name of the output particle list");
31 addParam("appendAllDaughters", m_appendAllDaughters,
32 "If true, all daughters of the matched MCParticle will be added in the order defined at the MCParticle. "
33 "If false, only the daughters described in the given decayString will be appended to the output particle.",
34 false);
35 addParam("skipNonPrimaryDaughters", m_skipNonPrimaryDaughters,
36 "If true, the secondary MC daughters will be skipped to append to the output particles. Default: true",
37 true);
38 addParam("writeOut", m_writeOut,
39 "If true, the output ParticleList will be saved by RootOutput. If false, it will be ignored when writing the file.", false);
40}
41
43{
44 bool valid = m_decaydescriptor.init(m_strDecay);
45 if (!valid)
46 B2ERROR("Invalid input DecayString: " << m_strDecay);
47
50
51 B2DEBUG(10, "particle list name: " << m_listName);
52 B2DEBUG(10, "antiparticle list name: " << m_antiListName);
53
54
55 // Register output particle list, particle store and relation to MCParticles
56 m_mcparticles.isRequired();
57
59 m_outputList.registerInDataStore(m_listName, flags);
61 m_extraInfoMap.registerInDataStore();
63
65 m_antiOutputList.registerInDataStore(m_antiListName, flags);
66 }
67}
68
70{
71 if (not m_extraInfoMap)
72 m_extraInfoMap.create();
73
74 // Create output particle list
75 const int motherPDG = m_decaydescriptor.getMother()->getPDGCode();
76 m_outputList.create();
77 m_outputList->initialize(motherPDG, m_listName);
78
80 m_antiOutputList.create();
81 m_antiOutputList->initialize(-1 * motherPDG, m_antiListName);
82 m_outputList->bindAntiParticleList(*(m_antiOutputList));
83 }
84
85 // loop over all MCParticles
86 int nMCParticles = m_mcparticles.getEntries();
87 for (int i = 0; i < nMCParticles; i++) {
88
89 if (abs(m_mcparticles[i]->getPDG()) != abs(motherPDG))
90 continue;
91
92 for (int iCC = 0; iCC < 2; iCC++) {
93 int arrayIndex = -1;
94 std::unique_ptr<DecayTree<MCParticle>> decay(match(m_mcparticles[i], m_decaydescriptor, iCC, arrayIndex));
95
96 MCParticle* mcp = decay->getObj();
97 if (!mcp) continue;
98
99 B2DEBUG(19, "Match!");
100
102 // if m_appendAllDaughters is True, create a new Particle appending all daughters
103 Particle* newParticle = m_particles.appendNew(mcp);
104 newParticle->addRelationTo(mcp);
105
106 for (auto mcDaughter : mcp->getDaughters()) {
107 if (mcDaughter->hasStatus(MCParticle::c_PrimaryParticle) or not m_skipNonPrimaryDaughters
108 or abs(mcp->getPDG()) == Const::Lambda.getPDGCode()) { // Lambda's daughters are not primary but it is not FSP at mcmatching
109
110 auto partDaughter = createParticleRecursively(mcDaughter, m_skipNonPrimaryDaughters);
111 newParticle->appendDaughter(partDaughter, /* updateType = */ false); // particleSource remain c_MCParticle
112 }
113 }
114 addUniqueParticleToList(newParticle);
115
116 } else if (arrayIndex == -1) {
117 // Particle is not created when no daughter is described in decayString
118 Particle* newParticle = m_particles.appendNew(mcp);
119 newParticle->addRelationTo(mcp);
120
121 addUniqueParticleToList(newParticle);
122 } else {
123 // Particle is already created
125
126 }
127
128 }
129 }
130}
131
132DecayTree<MCParticle>* MCDecayFinderModule::match(const MCParticle* mcp, const DecayDescriptor* d, bool isCC, int& arrayIndex)
133{
134 // Suffixes used in this method:
135 // P = Information from MCParticle
136 // D = Information from DecayDescriptor
137
138 // Create empty DecayTree as return value
139 auto* decay = new DecayTree<MCParticle>();
140
141 // Load PDG codes and compare,
142 const int iPDGD = d->getMother()->getPDGCode();
143 const int iPDGP = mcp->getPDG();
144 const bool isSelfConjugatedParticle = !(EvtPDLUtil::hasAntiParticle(iPDGD));
145
146 if (!isCC && iPDGD != iPDGP) return decay;
147 else if (isCC && (iPDGD != -iPDGP && !isSelfConjugatedParticle)) return decay;
148 else if (isCC && (iPDGD != iPDGP && isSelfConjugatedParticle)) return decay;
149 B2DEBUG(19, "PDG code matched: " << iPDGP);
150
151 // Get number of daughters in the decay descriptor.
152 // If no daughters in decay descriptor, no more checks needed.
153 const int nDaughtersD = d->getNDaughters();
154 if (nDaughtersD == 0) {
155 B2DEBUG(19, "DecayDescriptor has no Daughters, everything OK!");
156 decay->setObj(const_cast<MCParticle*>(mcp));
157 return decay; // arrayIndex is not set
158 }
159
160 // Get daughters of MCParticle
161 vector<const MCParticle*> daughtersP;
162 int nDaughtersRecursiveD = 0;
163 if (d->isIgnoreIntermediate()) {
164 appendParticles(mcp, daughtersP);
165 // Get number of daughter recursively if missing intermediate states are accepted.
166 nDaughtersRecursiveD = getNDaughtersRecursive(d);
167 } else {
168 const vector<MCParticle*>& tmpDaughtersP = mcp->getDaughters();
169 for (auto daug : tmpDaughtersP)
170 daughtersP.push_back(daug);
171 nDaughtersRecursiveD = nDaughtersD;
172 }
173
174 // The MCParticle must have at least as many daughters as the decaydescriptor
175 if (nDaughtersRecursiveD > (int)daughtersP.size()) {
176 B2DEBUG(10, "DecayDescriptor has more daughters than MCParticle!");
177 return decay;
178 }
179
180 // nested daughter has to be called at first to avoid overlap
181 std::vector<std::pair<int, int>> daughtersDepthMapD; // first: -1*depth, second:iDD
182 for (int iDD = 0; iDD < nDaughtersD; iDD++) {
183 int depth = 1; // to be incremented
184 countMaxDepthOfNest(d->getDaughter(iDD), depth);
185
186 daughtersDepthMapD.push_back({-1 * depth, iDD});
187 }
188 // sorted in ascending order = the deepest nested daughter will come first
189 std::sort(daughtersDepthMapD.begin(), daughtersDepthMapD.end());
190
191 // loop over all daughters of the decay descriptor
192 for (auto pair : daughtersDepthMapD) {
193 int iDD = pair.second;
194 // check if there is an unmatched particle daughter matching this decay descriptor daughter
195 bool isMatchDaughter = false;
196 auto itDP = daughtersP.begin();
197 while (itDP != daughtersP.end()) {
198 const MCParticle* daugP = *itDP;
199
200 int tmp; // array index of a created daughter Particle (not to be used)
201 DecayTree<MCParticle>* daughter = match(daugP, d->getDaughter(iDD), isCC, tmp);
202 if (!daughter->getObj()) {
203 ++itDP;
204 delete daughter;
205 continue;
206 }
207 // Matching daughter found, remove it from list of unmatched particle daughters
208 decay->append(daughter);
209 isMatchDaughter = true;
210 itDP = daughtersP.erase(itDP);
211
212 // if the matched daughter has daughters, they are also removed from the list
213 if (d->isIgnoreIntermediate() and d->getDaughter(iDD)->getNDaughters() != 0) {
214 vector<const MCParticle*> grandDaughtersP;
215 if (d->getDaughter(iDD)->isIgnoreIntermediate()) {
216 appendParticles(daugP, grandDaughtersP);
217 } else {
218 const vector<MCParticle*>& tmpGrandDaughtersP = daugP->getDaughters();
219 for (auto grandDaugP : tmpGrandDaughtersP)
220 grandDaughtersP.push_back(grandDaugP);
221 }
222
223 for (auto grandDaugP : grandDaughtersP) {
224 auto jtDP = itDP;
225 while (jtDP != daughtersP.end()) {
226 const MCParticle* daugP_j = *jtDP;
227 // if a grand-daughter matched a daughter, remove it.
228 if (grandDaugP == daugP_j)
229 jtDP = daughtersP.erase(jtDP);
230 else
231 ++jtDP;
232 }
233 }
234 }
235
236 break;
237 }
238 if (!isMatchDaughter) {
239 return decay;
240 }
241 }
242
243 // Ok, it seems that everything from the DecayDescriptor could be matched.
244 decay->setObj(const_cast<MCParticle*>(mcp));
245
246 Particle* newParticle = buildParticleFromDecayTree(decay, d);
247 int status = MCMatching::getMCErrors(newParticle);
248 if (status != MCMatching::c_Correct) { // isSignal != 1
249 B2DEBUG(10, "isSignal is not True. There was an additional particle left.");
250 decay->setObj(nullptr);
251 return decay;
252 }
253
254 arrayIndex = newParticle->getArrayIndex();
255
256 B2DEBUG(19, "Match found!");
257 return decay;
258}
259
260void MCDecayFinderModule::appendParticles(const MCParticle* gen, vector<const MCParticle*>& children)
261{
262 if (MCMatching::isFSP(gen->getPDG())) {
263 if (gen->getPDG() != Const::Kshort.getPDGCode()) // exception for K_S0
264 return; //stop at the bottom of the MC decay tree (ignore secondaries)
265
266 // Currently the decay of "FSP" cannot be specified except for K_S0,
267 // e.g. photon-conversion: gamma -> e+ e-, decay-in-flight: K+ -> mu+ nu_mu
268 }
269
270 const vector<MCParticle*>& genDaughters = gen->getDaughters();
271 for (auto daug : genDaughters) {
272 children.push_back(daug);
273 appendParticles(daug, children);
274 }
275}
276
278{
279 const int nDaughter = d->getNDaughters();
280 if (nDaughter == 0) return nDaughter;
281
282 int nDaughterRecursive = nDaughter;
283 for (int iDaug = 0; iDaug < nDaughter; iDaug++) {
284 const DecayDescriptor* dDaug = d->getDaughter(iDaug);
285
286 nDaughterRecursive += getNDaughtersRecursive(dDaug);
287 }
288
289 return nDaughterRecursive;
290}
291
293{
294 int maxDepth = 0;
295 for (int i = 0; i < d->getNDaughters(); i++) {
296 int tmp_depth = 1;
297 countMaxDepthOfNest(d->getDaughter(i), tmp_depth);
298
299 if (tmp_depth > maxDepth)
300 maxDepth = tmp_depth;
301 }
302
303 depth += maxDepth;
304}
305
307{
308 auto newParticle = buildParticleFromDecayTree(decay, dd);
309 int status = MCMatching::getMCErrors(newParticle);
310 return (status == MCMatching::c_Correct);
311}
312
313
315{
316 const int nDaughters = dd->getNDaughters();
317 const vector<DecayTree<MCParticle>*> decayDaughters = decay->getDaughters();
318
319 // sanity check
320 if ((int)decayDaughters.size() != nDaughters) {
321 B2ERROR("MCDecayFinderModule::buildParticleFromDecayTree Inconsistency on the number daughters between DecayTree and DecayDescriptor");
322 return nullptr;
323 }
324
325 // build particle from head of DecayTree
326 MCParticle* mcp = decay->getObj();
327 Particle* newParticle = m_particles.appendNew(mcp);
328 newParticle->addRelationTo(mcp);
329
330 int property = dd->getProperty();
331 property |= dd->getMother()->getProperty();
332 newParticle->setProperty(property);
333
334 // if nDaughters is 0 but mcp is not FSP, all primary daughters are appended for mcmatching
335 if (nDaughters == 0 and not MCMatching::isFSP(mcp->getPDG())) {
336 for (auto mcDaughter : mcp->getDaughters()) {
337 if (mcDaughter->hasStatus(MCParticle::c_PrimaryParticle) or
338 (abs(mcp->getPDG()) == Const::Lambda.getPDGCode())) { // Lambda's daughters are not primary but it is not FSP at mcmatching
339
340 auto partDaughter = createParticleRecursively(mcDaughter, true); // for mcmatching non-primary should be omitted
341 newParticle->appendDaughter(partDaughter, false);
342 }
343 }
344 }
345
346
347 // Daughters of DecayTree were filled in ascending order of depth of nest.
348 std::vector<std::pair<int, int>> daughtersDepthMapD; // first: -1*depth, second:iDD
349 for (int iDD = 0; iDD < nDaughters; iDD++) {
350 auto daugD = dd->getDaughter(iDD);
351
352 int depth = 1; // to be incremented
353 countMaxDepthOfNest(daugD, depth);
354 daughtersDepthMapD.push_back({-1 * depth, iDD});
355 }
356 std::sort(daughtersDepthMapD.begin(), daughtersDepthMapD.end());
357
358 for (int iDD = 0; iDD < nDaughters; iDD++) {
359
360 int index_decayDaughter = 0;
361 for (auto pair : daughtersDepthMapD) {
362 if (pair.second == iDD) break;
363 index_decayDaughter++;
364 }
365
366 Particle* partDaughter = buildParticleFromDecayTree(decayDaughters[index_decayDaughter], dd->getDaughter(iDD));
367
368 int daughterProperty = dd->getDaughter(iDD)->getMother()->getProperty();
369 newParticle->appendDaughter(partDaughter, false, daughterProperty);
370 }
371
372 return newParticle;
373}
374
376{
377 Particle* newParticle = m_particles.appendNew(mcp);
378 newParticle->addRelationTo(mcp);
379
380 for (auto mcDaughter : mcp->getDaughters()) {
381 if (mcDaughter->hasStatus(MCParticle::c_PrimaryParticle) or not skipNonPrimaryDaughters
382 or abs(mcp->getPDG()) == Const::Lambda.getPDGCode()) { // Lambda's daughters are not primary but it is not FSP at mcmatching
383 auto partDaughter = createParticleRecursively(mcDaughter, skipNonPrimaryDaughters);
384 newParticle->appendDaughter(partDaughter, false);
385 }
386 }
387
388 return newParticle;
389}
390
392{
393 for (auto& existingParticle : *m_outputList) {
394 if (existingParticle.isCopyOf(newParticle))
395 return;
396 }
397
398 m_outputList->addParticle(newParticle);
399}
int getPDGCode() const
PDG code.
Definition: Const.h:473
static const ParticleType Lambda
Lambda particle.
Definition: Const.h:679
static const ParticleType Kshort
K^0_S particle.
Definition: Const.h:677
EStoreFlags
Flags describing behaviours of objects etc.
Definition: DataStore.h:69
@ c_WriteOut
Object/array should be saved by output modules.
Definition: DataStore.h:70
@ c_DontWriteOut
Object/array should be NOT saved by output modules.
Definition: DataStore.h:71
@ c_Event
Different object in each event, all objects/arrays are invalidated after event() function has been ca...
Definition: DataStore.h:59
int getPDGCode() const
Return PDG code.
int getProperty() const
return property of the particle.
The DecayDescriptor stores information about a decay tree or parts of a decay tree.
bool init(const std::string &str)
Initialise the DecayDescriptor from given string.
const DecayDescriptor * getDaughter(int i) const
return i-th daughter (0 based index).
int getNDaughters() const
return number of direct daughters.
int getProperty() const
return property of the particle.
bool isIgnoreIntermediate() const
Check if intermediate resonances/particles shall be ignored.
const DecayDescriptorParticle * getMother() const
return mother.
This is a helper class for the MCDecayFinderModule.
Definition: DecayTree.h:20
Particle * createParticleRecursively(const MCParticle *mcp, bool skipNonPrimaryDaughters)
Create a Particle from the given MCParticle appending all daughters of the MCParticle.
bool m_isSelfConjugatedParticle
Is the particle list for a self-conjugated particle.
std::string m_antiListName
Name of output anti-particle list.
DecayTree< MCParticle > * match(const MCParticle *mcp, const DecayDescriptor *d, bool isCC, int &arrayIndex)
Search for MCParticles matching the given DecayString.
Particle * buildParticleFromDecayTree(const DecayTree< MCParticle > *decay, const DecayDescriptor *dd)
Build a Particle from the given DecayTree and DecayDescriptor.
int getNDaughtersRecursive(const DecayDescriptor *d)
Recursively get number of daughters of given DecayDescriptor.
virtual void initialize() override
Initialises the module.
bool performMCMatching(const DecayTree< MCParticle > *decay, const DecayDescriptor *dd)
Perform the MCMatching on the Particle built from the given DecayTree and DecayDescriptor.
virtual void event() override
Method called for each event.
StoreArray< MCParticle > m_mcparticles
StoreArray of MCParticles.
std::string m_listName
Name of output particle list.
StoreArray< Particle > m_particles
StoreArray of Particles.
StoreObjPtr< ParticleList > m_outputList
output particle list
std::string m_strDecay
Decay string to build the decay descriptor.
void appendParticles(const MCParticle *gen, std::vector< const MCParticle * > &children)
Recursively gather all MC daughters of gen.
void addUniqueParticleToList(Particle *newParticle)
Add a new Particle into the output ParticleList if it does not exist already.
bool m_appendAllDaughters
if true, all daughters are appended
bool m_skipNonPrimaryDaughters
toggle skip of secondary MC daughters
StoreObjPtr< ParticleList > m_antiOutputList
output anti-particle list
DecayDescriptor m_decaydescriptor
Decay descriptor of decays to look for.
bool m_writeOut
toggle output particle list btw.
void countMaxDepthOfNest(const DecayDescriptor *d, int &depth)
Count the max depth of nest from the given DecayDescriptor.
StoreObjPtr< ParticleExtraInfoMap > m_extraInfoMap
object pointer to ParticleExtraInfoMaps
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:32
@ c_PrimaryParticle
bit 0: Particle is primary particle.
Definition: MCParticle.h:47
std::vector< Belle2::MCParticle * > getDaughters() const
Get vector of all daughter particles, empty vector if none.
Definition: MCParticle.cc:52
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
Class to store reconstructed particles.
Definition: Particle.h:75
void setProperty(const int properties)
sets m_properties
Definition: Particle.h:374
void appendDaughter(const Particle *daughter, const bool updateType=true, const int daughterProperty=c_Ordinary)
Appends index of daughter to daughters index array.
Definition: Particle.cc:676
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).
int getArrayIndex() const
Returns this object's array index (in StoreArray), or -1 if not found.
bool registerInDataStore(DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut)
Register the object/array in the DataStore.
T * appendNew()
Construct a new T object at the end of the array.
Definition: StoreArray.h:246
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
T * getObj() const
Return the decaying object itself, e.g.
Definition: DecayTree.h:65
std::vector< DecayTree< T > * > getDaughters() const
Return list of decay daughters.
Definition: DecayTree.h:58
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
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
bool hasAntiParticle(int pdgCode)
Checks if the particle with given pdg code has an anti-particle or not.
Definition: EvtPDLUtil.cc:12
std::string antiParticleListName(const std::string &listName)
Returns name of anti-particle-list corresponding to listName.
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24
STL namespace.
@ c_Correct
This Particle and all its daughters are perfectly reconstructed.
Definition: MCMatching.h:34
static int getMCErrors(const Belle2::Particle *particle, const Belle2::MCParticle *mcParticle=nullptr)
Returns quality indicator of the match as a bit pattern where the individual bits indicate the the ty...
Definition: MCMatching.cc:282
static bool isFSP(int pdg)
Returns true if given PDG code indicates a FSP.
Definition: MCMatching.cc:448