Belle II Software  light-2205-abys
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 
18 using namespace std;
19 using namespace Belle2;
20 
21 // Register module in the framework
22 REG_MODULE(MCDecayFinder);
23 
24 MCDecayFinderModule::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("writeOut", m_writeOut,
32  "If true, the output ParticleList will be saved by RootOutput. If false, it will be ignored when writing the file.", false);
33 }
34 
36 {
38 
41 
42  B2DEBUG(10, "particle list name: " << m_listName);
43  B2DEBUG(10, "antiparticle list name: " << m_antiListName);
44 
45 
46  // Register output particle list, particle store and relation to MCParticles
47  m_mcparticles.isRequired();
48 
50  m_outputList.registerInDataStore(m_listName, flags);
51  m_particles.registerInDataStore(flags);
52  m_extraInfoMap.registerInDataStore();
53  m_particles.registerRelationTo(m_mcparticles, DataStore::c_Event, flags);
54 
56  m_antiOutputList.registerInDataStore(m_antiListName, flags);
57  }
58 }
59 
61 {
62  // Create output particle list
63  m_outputList.create();
65 
67  m_antiOutputList.create();
69  m_outputList->bindAntiParticleList(*(m_antiOutputList));
70  }
71 
72  // loop over all MCParticles
73  int nMCParticles = m_mcparticles.getEntries();
74  for (int i = 0; i < nMCParticles; i++) {
75  for (int iCC = 0; iCC < 2; iCC++) {
76  std::unique_ptr<DecayTree<MCParticle>> decay(match(m_mcparticles[i], m_decaydescriptor, iCC));
77  if (decay->getObj()) {
78  B2DEBUG(19, "Match!");
79  int iIndex = write(decay.get());
80  m_outputList->addParticle(m_particles[iIndex]);
81  }
82  }
83  }
84 }
85 
87 {
88  // Suffixes used in this method:
89  // P = Information from MCParticle
90  // D = Information from DecayDescriptor
91 
92  // Create empty DecayTree as return value
93  auto* decay = new DecayTree<MCParticle>();
94 
95  // Load PDG codes and compare,
96  int iPDGD = d->getMother()->getPDGCode();
97  int iPDGP = mcp->getPDG();
98 
99  bool isSelfConjugatedParticle = !(EvtPDLUtil::hasAntiParticle(iPDGD));
100 
101  if (!isCC && iPDGD != iPDGP) return decay;
102  else if (isCC && (iPDGD != -iPDGP && !isSelfConjugatedParticle)) return decay;
103  else if (isCC && (iPDGD != iPDGP && isSelfConjugatedParticle)) return decay;
104  B2DEBUG(19, "PDG code matched: " << iPDGP);
105 
106  // Get number of daughters in the decay descriptor.
107  // If no daughters in decay descriptor, no more checks needed.
108  int nDaughtersD = d->getNDaughters();
109  if (nDaughtersD == 0) {
110  B2DEBUG(19, "DecayDescriptor has no Daughters, everything OK!");
111  decay->setObj(const_cast<MCParticle*>(mcp));
112  return decay;
113  }
114  // Get number of daughter recursively if missing intermediate states are accepted.
115  int nDaughtersRecursiveD = nDaughtersD;
116  if (d->isIgnoreIntermediate()) {
117  nDaughtersRecursiveD = getNDaughtersRecursive(d);
118  }
119 
120  // Get daughters of MCParticle
121  vector<const MCParticle*> daughtersP;
122  if (d->isIgnoreIntermediate()) {
123  appendParticles(mcp, daughtersP);
124  } else {
125  const vector<MCParticle*>& tmpDaughtersP = mcp->getDaughters();
126  for (auto daug : tmpDaughtersP)
127  daughtersP.push_back(daug);
128  }
129 
130  // The MCParticle must have at least as many daughters as the decaydescriptor
131  int nDaughtersP = daughtersP.size();
132  if (nDaughtersRecursiveD > nDaughtersP) {
133  B2DEBUG(10, "DecayDescriptor has more daughters than MCParticle!");
134  return decay;
135  }
136 
137  // loop over all daughters of the decay descriptor
138  for (int iDD = 0; iDD < nDaughtersD; iDD++) {
139  // check if there is an unmatched particle daughter matching this decay descriptor daughter
140  bool isMatchDaughter = false;
141  auto itDP = daughtersP.begin();
142  while (itDP != daughtersP.end()) {
143  const MCParticle* daugP = *itDP;
144  DecayTree<MCParticle>* daughter = match(daugP, d->getDaughter(iDD), isCC);
145  if (!daughter->getObj()) {
146  ++itDP;
147  continue;
148  }
149  // Matching daughter found, remove it from list of unmatched particle daughters
150  decay->append(daughter);
151  isMatchDaughter = true;
152  itDP = daughtersP.erase(itDP);
153 
154  // if the matched daughter has daughters, they are also removed from the list
155  if (d->isIgnoreIntermediate() and d->getDaughter(iDD)->getNDaughters() != 0) {
156  vector<const MCParticle*> grandDaughtersP;
157  if (d->getDaughter(iDD)->isIgnoreIntermediate()) {
158  appendParticles(daugP, grandDaughtersP);
159  } else {
160  const vector<MCParticle*>& tmpGrandDaughtersP = daugP->getDaughters();
161  for (auto grandDaugP : tmpGrandDaughtersP)
162  grandDaughtersP.push_back(grandDaugP);
163  }
164 
165  for (auto grandDaugP : grandDaughtersP) {
166  auto jtDP = itDP;
167  while (jtDP != daughtersP.end()) {
168  const MCParticle* daugP_j = *jtDP;
169  // if a grand-daughter matched a daughter, remove it.
170  if (grandDaugP == daugP_j)
171  jtDP = daughtersP.erase(jtDP);
172  else
173  ++jtDP;
174  }
175  }
176  }
177 
178  break;
179  }
180  if (!isMatchDaughter) {
181  return decay;
182  }
183  }
184 
185  // Ok, it seems that everything from the DecayDescriptor could be matched.
186  // If the decay is NOT INCLUSIVE, no unmatched MCParticles should be left
187  bool isInclusive = (d->isIgnoreMassive() and d->isIgnoreNeutrino() and d->isIgnoreGamma() and d->isIgnoreRadiatedPhotons());
188  if (!isInclusive) {
189  B2DEBUG(10, "Check for left over MCParticles!\n");
190  for (auto daugP : daughtersP) {
191  if (daugP->getPDG() == Const::photon.getPDGCode()) { // gamma
192  // if gamma is FSR or produced by PHOTOS
193  if (MCMatching::isFSR(daugP) or daugP->hasStatus(MCParticle::c_IsPHOTOSPhoton)) {
194  // if the radiated photons are ignored, we can skip the particle
195  if (d->isIgnoreRadiatedPhotons()) continue;
196  }
197  // else if missing gamma is ignored, we can skip the particle
198  else if (d->isIgnoreGamma()) continue;
199  } else if ((abs(daugP->getPDG()) == 12 or abs(daugP->getPDG()) == 14 or abs(daugP->getPDG()) == 16)) { // neutrino
200  if (d->isIgnoreNeutrino()) continue;
201  } else if (MCMatching::isFSP(daugP->getPDG())) { // other final state particle
202  if (d->isIgnoreMassive()) continue;
203  } else { // intermediate
204  if (d->isIgnoreIntermediate()) continue;
205  }
206  B2DEBUG(10, "There was an additional particle left. Found " << daugP->getPDG());
207  return decay;
208  }
209  }
210  decay->setObj(const_cast<MCParticle*>(mcp));
211  B2DEBUG(19, "Match found!");
212  return decay;
213 }
214 
216 {
217  // Create new Particle in particles array
218  Particle* newParticle = m_particles.appendNew(decay->getObj());
219 
220  // set relation between the created Particle and the MCParticle
221  newParticle->addRelationTo(decay->getObj());
222  const int iIndex = m_particles.getEntries() - 1;
223 
224  // Now save also daughters of this MCParticle and set the daughter relation
225  vector< DecayTree<MCParticle>* > daughters = decay->getDaughters();
226  int nDaughers = daughters.size();
227  for (int i = 0; i < nDaughers; ++i) {
228  int iIndexDaughter = write(daughters[i]);
229  newParticle->appendDaughter(iIndexDaughter);
230  }
231  return iIndex;
232 }
233 
234 void MCDecayFinderModule::appendParticles(const MCParticle* gen, vector<const MCParticle*>& children)
235 {
236  if (MCMatching::isFSP(gen->getPDG()))
237  return; //stop at the bottom of the MC decay tree (ignore secondaries)
238 
239  const vector<MCParticle*>& genDaughters = gen->getDaughters();
240  for (auto daug : genDaughters) {
241  children.push_back(daug);
242  appendParticles(daug, children);
243  }
244 }
245 
247 {
248  const int nDaughter = d->getNDaughters();
249  if (nDaughter == 0) return nDaughter;
250 
251  int nDaughterRecursive = nDaughter;
252  for (int iDaug = 0; iDaug < nDaughter; iDaug++) {
253  const DecayDescriptor* dDaug = d->getDaughter(iDaug);
254 
255  nDaughterRecursive += getNDaughtersRecursive(dDaug);
256  }
257 
258  return nDaughterRecursive;
259 }
int getPDGCode() const
PDG code.
Definition: Const.h:354
static const ParticleType photon
photon particle
Definition: Const.h:554
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.
The DecayDescriptor stores information about a decay tree or parts of a decay tree.
bool isIgnoreRadiatedPhotons() const
Check if additional radiated photons shall be ignored.
bool init(const std::string &str)
Initialise the DecayDescriptor from given string.
const DecayDescriptorParticle * getMother() const
return mother.
bool isIgnoreNeutrino() const
Check if missing neutrinos shall be ignored.
int getNDaughters() const
return number of direct daughters.
bool isIgnoreMassive() const
Check if missing massive final state particles shall be ignored.
const DecayDescriptor * getDaughter(int i) const
return i-th daughter (0 based index).
bool isIgnoreGamma() const
Check if missing gammas shall be ignored.
bool isIgnoreIntermediate() const
Check if intermediate resonances/particles shall be ignored.
This is a helper class for the MCDecayFinderModule.
Definition: DecayTree.h:20
bool m_isSelfConjugatedParticle
Is the particle list for a self-conjugated particle.
std::string m_antiListName
Name of output anti-particle list.
int getNDaughtersRecursive(const DecayDescriptor *d)
Recursively get number of daughters of given DecayDescriptor.
virtual void initialize() override
Initialises the module.
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.
int write(DecayTree< MCParticle > *decay)
Create Particle from matched MCParticle and write to Particle list.
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.
DecayTree< MCParticle > * match(const MCParticle *mcp, const DecayDescriptor *d, bool isCC)
Search for MCParticles matching the given DecayString.
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.
StoreObjPtr< ParticleExtraInfoMap > m_extraInfoMap
object pointer to ParticleExtraInfoMaps
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:34
@ c_IsPHOTOSPhoton
bit 8: Particle is an radiative photon from PHOTOS
Definition: MCParticle.h:65
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:114
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:74
void appendDaughter(const Particle *daughter, const bool updateType=true)
Appends index of daughter to daughters index array.
Definition: Particle.cc:700
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).
T * getObj()
Return the decaying object itself, e.g.
Definition: DecayTree.h:68
std::vector< DecayTree< T > * > getDaughters()
Return list of decay daughters.
Definition: DecayTree.h:62
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:23
static bool isFSR(const Belle2::MCParticle *p)
Returns true if given MCParticle is a final state radiation (FSR) photon based on MCParticle::c_IsFSR...
Definition: MCMatching.cc:466
static bool isFSP(int pdg)
Returns true if given PDG code indicates a FSP.
Definition: MCMatching.cc:445