Belle II Software light-2406-ragdoll
DecayDescriptor.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/DecayDescriptor/DecayDescriptor.h>
10#include <analysis/DecayDescriptor/DecayString.h>
11#include <analysis/DecayDescriptor/DecayStringDecay.h>
12#include <analysis/DecayDescriptor/DecayStringGrammar.h>
13#include <analysis/utility/EvtPDLUtil.h>
14#include <analysis/dataobjects/Particle.h>
15
16#include <analysis/utility/AnalysisConfiguration.h>
17
18#include <mdst/dataobjects/MCParticle.h>
19
20#include <framework/gearbox/Const.h>
21#include <framework/logging/Logger.h>
22
23#include <TDatabasePDG.h>
24
25#include <boost/variant/get.hpp>
26#include <boost/spirit/include/qi.hpp>
27#include <algorithm>
28#include <set>
29#include <utility>
30
31using namespace Belle2;
32using namespace std;
33
35
37 m_mother(),
38 m_iDaughter_p(-1),
39 m_daughters(),
40 m_properties(0),
41 m_isNULL(false),
42 m_isInitOK(false)
43{
44}
45
46bool DecayDescriptor::init(const std::string& str)
47{
48 // The decay string grammar
51 std::string::const_iterator iter = str.begin();
52 std::string::const_iterator end = str.end();
53 bool r = phrase_parse(iter, end, g, boost::spirit::unicode::space, s);
54 if (!r || iter != end) return false;
55 return init(s);
56}
57
59{
60 // The DecayString is a hybrid, it can be
61 // a) DecayStringParticleList
62 // b) DecayStringDecay
63
64 if (const DecayStringParticle* p = boost::get< DecayStringParticle >(&s)) {
66 if (!m_isInitOK) {
67 B2WARNING("Could not initialise mother particle " << p->m_strName);
68 return false;
69 }
70 return true;
71 } else if (const DecayStringDecay* d = boost::get< DecayStringDecay > (&s)) {
72 // Initialise list of mother particles
73 m_isInitOK = m_mother.init(d->m_mother);
74 if (!m_isInitOK) {
75 B2WARNING("Could not initialise mother particle " << d->m_mother.m_strName);
76 return false;
77 }
78
79 // Identify arrow type
80 if (d->m_strArrow == "->") {
83 } else if (d->m_strArrow == "=norad=>") {
85 } else if (d->m_strArrow == "=direct=>") {
87 } else if (d->m_strArrow == "=exact=>") {
88 // do nothing
89 } else {
90 B2WARNING("Unknown arrow: " << d->m_strArrow);
91 m_isInitOK = false;
92 return false;
93 }
94
95 // Initialise list of daughters
96 if (d->m_daughters.empty()) {
97 m_isInitOK = false;
98 return false;
99 }
100 int nDaughters = d->m_daughters.size();
101 for (int iDaughter = 0; iDaughter < nDaughters; iDaughter++) {
102 DecayDescriptor daughter;
103 m_isInitOK = daughter.init(d->m_daughters[iDaughter]);
104 if (!m_isInitOK) {
105 B2WARNING("Could not initialise daughter!");
106 return false;
107 }
108 m_daughters.push_back(daughter);
109 }
110
111 // Initialise list of keywords
112 // For neutrino
113 if ((std::find(d->m_keywords.begin(), d->m_keywords.end(), "?nu")) != d->m_keywords.end()) {
115 }
116 // For gamma
117 if ((std::find(d->m_keywords.begin(), d->m_keywords.end(), "?gamma")) != d->m_keywords.end()) {
119 }
120 // For massive FSP
121 if ((std::find(d->m_keywords.begin(), d->m_keywords.end(), "...")) != d->m_keywords.end()) {
123 }
124 // For brems photons
125 if ((std::find(d->m_keywords.begin(), d->m_keywords.end(), "?addbrems")) != d->m_keywords.end()) {
127 }
128
129 return true;
130 }
131 m_isInitOK = false;
132 return false;
133}
134
135template <class T>
136int DecayDescriptor::match(const T* p, int iDaughter_p)
137{
138 // this DecayDescriptor was not matched or
139 // it is not the daughter of another DecayDescriptor
140 m_iDaughter_p = -1;
141
142 if (!p) {
143 B2WARNING("NULL pointer provided instead of particle.");
144 return 0;
145 }
146
147 int iPDGCode_p = 0;
148 if (const auto* part_test = dynamic_cast<const Particle*>(p))
149 iPDGCode_p = part_test->getPDGCode();
150 else if (const auto* mc_test = dynamic_cast<const MCParticle*>(p))
151 iPDGCode_p = mc_test->getPDG();
152 else {
153 B2WARNING("Template type not supported!");
154 return 0;
155 }
156
157 int iPDGCodeCC_p = TDatabasePDG::Instance()->GetParticle(iPDGCode_p)->AntiParticle()->PdgCode();
158 int iPDGCode_d = m_mother.getPDGCode();
159 if (abs(iPDGCode_d) != abs(iPDGCode_p)) return 0;
160 int iCC = 0;
161 if (iPDGCode_p == iPDGCodeCC_p) iCC = 3;
162 else if (iPDGCode_d == iPDGCode_p) iCC = 1;
163 else if (iPDGCode_d == iPDGCodeCC_p) iCC = 2;
164
165 const std::vector<T*> daughterList = p->getDaughters();
166 int nDaughters_p = daughterList.size();
167
168 // 1st case: the descriptor has no daughters => nothing to check
169 if (getNDaughters() == 0) {
170 m_iDaughter_p = iDaughter_p;
171 return iCC;
172 }
173
174 // 2nd case: the descriptor has daughters, but not the particle
175 // => that is not allowed!
176 if (nDaughters_p == 0) return 0;
177
178 // 3rd case: the descriptor and the particle have daughters
179 // There are two cases that can happen when matching the
180 // DecayDescriptor daughters to the particle daughters:
181 // 1. The match is unambiguous -> no problem
182 // 2. Multiple particle daughters match the same DecayDescriptor daughter
183 // -> in the latter case the ambiguity is resolved later
184
185 // 1. DecayDescriptor -> Particle relation for the cases where only one particle matches
186 vector< pair< int, int > > singlematch;
187 // 2. DecayDescriptor -> Particle relation for the cases where multiple particles match
188 vector< pair< int, set<int> > > multimatch;
189 // Are there ambiguities in the match?
190 bool isAmbiguities = false;
191 // The particle daughters that have been matched
192 set<int> matches_global;
193
194 // check if the daughters match
195 for (int iDaughter_d = 0; iDaughter_d < getNDaughters(); iDaughter_d++) {
196 set<int> matches;
197 for (int jDaughter_p = 0; jDaughter_p < nDaughters_p; jDaughter_p++) {
198 const T* daughter = daughterList[jDaughter_p];
199 int iPDGCode_daughter_p = 0;
200 if (const auto* part_test = dynamic_cast<const Particle*>(daughter))
201 iPDGCode_daughter_p = part_test->getPDGCode();
202 else if (const auto* mc_test = dynamic_cast<const MCParticle*>(daughter))
203 iPDGCode_daughter_p = mc_test->getPDG();
204
205 if (iDaughter_d == 0 && (this->isIgnoreRadiatedPhotons() or this->isIgnoreGamma() or this->isIgnoreBrems())
206 && iPDGCode_daughter_p == Const::photon.getPDGCode())
207 matches_global.insert(jDaughter_p);
208
209 int iMatchResult = m_daughters[iDaughter_d].match(daughter, jDaughter_p);
210 if (iMatchResult < 0) isAmbiguities = true;
211 if (abs(iMatchResult) == 2 && iCC == 1) continue;
212 if (abs(iMatchResult) == 1 && iCC == 2) continue;
213 if (abs(iMatchResult) == 2 && iCC == 3) continue;
214 matches.insert(jDaughter_p);
215 matches_global.insert(jDaughter_p);
216 }
217 if (matches.empty()) return 0;
218 if (matches.size() == 1) {
219 int jDaughter_p = *(matches.begin());
220 singlematch.emplace_back(iDaughter_d, jDaughter_p);
221 } else multimatch.emplace_back(iDaughter_d, matches);
222 }
223
224 // Now, all daughters of the particles should be matched to at least one DecayDescriptor daughter
225 if (!(this->isIgnoreIntermediate() or this->isIgnoreMassive() or this->isIgnoreNeutrino())
226 && int(matches_global.size()) != nDaughters_p) return 0;
227
228 // In case that there are DecayDescriptor daughters with multiple matches, try to solve the problem
229 // by removing the daughter candidates which are already used in other unambiguous relations.
230 // This is done iteratively. We limit the maximum number of attempts to 20 to avoid an infinite loop.
231 bool isModified = true;
232 for (int iTry = 0; iTry < 20; iTry++) {
233 if (int(singlematch.size()) == getNDaughters()) break;
234 if (!isModified) break;
235 isModified = false;
236 for (auto& itMulti : multimatch) {
237 for (auto& itSingle : singlematch) {
238 // try to remove particle from the multimatch list
239 if (itMulti.second.erase(itSingle.second)) {
240 B2FATAL("Trying to execute part of the code with known bug, which is not fixed yet! Send email to anze.zupanc@ijs.si with notification that this happens!");
241 /*
242 This part of the code is commented, because of the following error:
243 Iterator 'itMulti' used after element has been erased.
244
245 // if multimatch list contains only one particle candidate, move the entry to the singlematch list
246 if (itMulti->second.size() == 1) {
247 int iDaughter_d = itMulti->first;
248 int iDaughter_p = *(itMulti->second.begin());
249 singlematch.push_back(make_pair(iDaughter_d, iDaughter_p));
250 multimatch.erase(itMulti);
251 // call match function again to set the correct daughter
252 if (!isAmbiguities) {
253 const T* daughter = daughterList[iDaughter_p];
254 if (!daughter) continue;
255 m_daughters[iDaughter_d].match(daughter, iDaughter_p);
256 }
257 --itMulti;
258 isModified = true;
259 break;
260 }
261 */
262 }
263 }
264 }
265 }
266
267 if (!multimatch.empty()) isAmbiguities = true;
268 if (isAmbiguities) return -iCC;
269 else {
270 m_iDaughter_p = iDaughter_p;
271 return iCC;
272 }
273 return 0;
274}
275
277{
278 m_iDaughter_p = -1;
279 int nDaughters = m_daughters.size();
280 for (int iDaughter = 0; iDaughter < nDaughters; iDaughter++) m_daughters[iDaughter].resetMatch();
281}
282
283vector<const Particle*> DecayDescriptor::getSelectionParticles(const Particle* particle)
284{
285 // Create vector for output
286 vector<const Particle*> selparticles;
287 if (m_mother.isSelected()) {
288 int motherPDG = abs(particle->getPDGCode());
289 int decayDescriptorMotherPDG = abs(m_mother.getPDGCode());
290 if (motherPDG != decayDescriptorMotherPDG)
291 B2ERROR("The PDG code of the mother particle (" << motherPDG <<
292 ") does not match the PDG code of the DecayDescriptor mother (" << decayDescriptorMotherPDG <<
293 ")! Check the order of the decay string is the same you expect in the reconstructed Particles.");
294 selparticles.push_back(particle);
295 }
296 int nDaughters_d = getNDaughters();
297 for (int iDaughter_d = 0; iDaughter_d < nDaughters_d; ++iDaughter_d) {
298 // retrieve the particle daughter ID from this DecayDescriptor daughter
299 int iDaughter_p = m_daughters[iDaughter_d].getMatchedDaughter();
300 // If the particle daughter ID is below one, the match function was not called before
301 // or the match was ambiguous. In this case try to use the daughter ID of the DecayDescriptor.
302 // This corresponds to using the particle order in the decay string.
303 if (iDaughter_p < 0) iDaughter_p = iDaughter_d;
304 const Particle* daughter = particle->getDaughter(iDaughter_p);
305 if (!daughter) {
306 B2WARNING("Could not find daughter!");
307 continue;
308 }
309 // check if the daughter has the correct PDG code
310 int daughterPDG = abs(daughter->getPDGCode());
311 int decayDescriptorDaughterPDG = abs(m_daughters[iDaughter_d].getMother()->getPDGCode());
312 if (daughterPDG != decayDescriptorDaughterPDG) {
313 B2ERROR("The PDG code of the particle daughter (" << daughterPDG <<
314 ") does not match the PDG code of the DecayDescriptor daughter (" << decayDescriptorDaughterPDG <<
315 ")! Check the order of the decay string is the same you expect in the reconstructed Particles.");
316 break;
317 }
318 vector<const Particle*> seldaughters = m_daughters[iDaughter_d].getSelectionParticles(daughter);
319 selparticles.insert(selparticles.end(), seldaughters.begin(), seldaughters.end());
320 }
321 return selparticles;
322}
323
325{
326
327 std::vector<int> decay, decaybar;
328 for (int i = 0; i < getNDaughters(); ++i) {
329 const DecayDescriptorParticle* daughter = getDaughter(i)->getMother();
330 int pdg = daughter->getPDGCode();
331 decay.push_back(pdg);
332 decaybar.push_back(Belle2::EvtPDLUtil::hasAntiParticle(pdg) ? -pdg : pdg);
333 }
334
335 std::sort(decay.begin(), decay.end());
336 std::sort(decaybar.begin(), decaybar.end());
337
338 return (not Belle2::EvtPDLUtil::hasAntiParticle(getMother()->getPDGCode())) || (decay == decaybar);
339
340}
341
343{
344 vector<string> strNames;
345 if (m_mother.isSelected()) strNames.push_back(m_mother.getNameSimple());
346 for (auto& daughter : m_daughters) {
347 vector<string> strDaughterNames = daughter.getSelectionNames();
348 int nDaughters = strDaughterNames.size();
349 for (int iDaughter = 0; iDaughter < nDaughters; iDaughter++) {
350 // Checking variable naming scheme from AnalysisConfiguratin
351 // For example, effect of possible schemes for PX variable
352 // of pi0 from D in decay B->(D->pi0 pi) pi0:
353 // default: B_D_pi0_PX
354 // semidefault: D_pi0_PX
355 // laconic: pi01_PX
356 if (AnalysisConfiguration::instance()->getTupleStyle() == "laconic") continue;
357 if ((AnalysisConfiguration::instance()->getTupleStyle() == "semilaconic") && (iDaughter == nDaughters)) continue;
358 strDaughterNames[iDaughter] = m_mother.getNameSimple() + "_" + strDaughterNames[iDaughter];
359 }
360 strNames.insert(strNames.end(), strDaughterNames.begin(), strDaughterNames.end());
361 }
362
363 // search for multiple occurrence of the same name and then distinguish by attaching a number
364
365 for (auto itName = strNames.begin(); itName != strNames.end(); ++itName) {
366 if (count(itName, strNames.end(), *itName) == 1) continue;
367 // multiple occurrence found!
368 string strNameOld = *itName;
369 auto itOccurrence = strNames.begin();
370 int iOccurrence = 0;
371 while (iOccurrence <= 10) {
372 // find next occurrence of the identical particle name defined in DecayDescriptor
373 itOccurrence = find(itOccurrence, strNames.end(), strNameOld);
374 // stop, if nothing found
375 if (itOccurrence == strNames.end()) break;
376 // create new particle name by attaching a number
377 string strNameNew = strNameOld + std::to_string(iOccurrence);
378 // check if the new particle name exists already, if not, then it is OK to use it
379 if (count(strNames.begin(), strNames.end(), strNameNew) == 0) {
380 *itOccurrence = strNameNew;
381 ++itOccurrence;
382 }
383 iOccurrence++;
384 }
385 if (iOccurrence == 10) {
386 B2ERROR("DecayDescriptor::getSelectionNames - Something is wrong! More than 10x the same name!");
387 break;
388 }
389 }
390 return strNames;
391}
392
394{
395 vector<int> listPDG;
396 if (m_mother.isSelected()) listPDG.push_back(m_mother.getPDGCode());
397 for (auto& daughter : m_daughters) {
398 vector<int> listPDGDaughters = daughter.getSelectionPDGCodes();
399 listPDG.insert(listPDG.end(), listPDGDaughters.begin(), listPDGDaughters.end());
400 }
401 return listPDG;
402}
403
404
405std::vector<std::vector<std::pair<int, std::string>>> DecayDescriptor::getHierarchyOfSelected()
406{
407 if (not m_hierarchy.empty()) {
408 std::vector<std::vector<std::pair<int, std::string>>> hierarchy = m_hierarchy;
409 return hierarchy;
410 }
411 std::vector<std::pair<int, std::string>> currentPath;
412 currentPath.emplace_back(0, m_mother.getNameSimple());
413 return getHierarchyOfSelected(currentPath);
414}
415
416std::vector<std::vector<std::pair<int, std::string>>> DecayDescriptor::getHierarchyOfSelected(
417 const std::vector<std::pair<int, std::string>>& currentPath)
418{
419 if (m_mother.isSelected()) m_hierarchy.push_back(currentPath);
420 for (std::size_t i = 0; i < m_daughters.size(); i++) {
421 std::vector<std::pair<int, std::string>> newPath = currentPath;
422 newPath.emplace_back(i, m_daughters[i].getMother()->getNameSimple());
423 std::vector<std::vector<std::pair<int, std::string>>> foundPathes = m_daughters[i].getHierarchyOfSelected(newPath);
424 for (auto& path : foundPathes) m_hierarchy.push_back(path);
425 }
426 std::vector<std::vector<std::pair<int, std::string>>> hierarchy = m_hierarchy;
427 return hierarchy;
428}
static AnalysisConfiguration * instance()
Returns a pointer to the singleton instance.
int getPDGCode() const
PDG code.
Definition: Const.h:473
static const ParticleType photon
photon particle
Definition: Const.h:673
Represents a particle in the DecayDescriptor.
int getPDGCode() const
Return PDG code.
std::string getNameSimple() const
Return the name from getName() without + - * or anti-.
bool init(const DecayStringParticle &p)
initialise member variables from std::string member variables contained in a DecayStringParticle stru...
bool isSelected() const
Is the particle selected in the decay string?
The DecayDescriptor stores information about a decay tree or parts of a decay tree.
DecayDescriptor()
Default ctor.
bool isIgnoreBrems() const
Check if added Brems gammas shall be ignored.
bool isIgnoreRadiatedPhotons() const
Check if additional radiated photons shall be ignored.
bool init(const std::string &str)
Initialise the DecayDescriptor from given string.
bool isSelfConjugated() const
Is the decay or the particle self conjugated.
const DecayDescriptor * getDaughter(int i) const
return i-th daughter (0 based index).
DecayDescriptorParticle m_mother
Mother of the decay ('left side').
bool m_isInitOK
Is this object initialized correctly?
bool isIgnoreNeutrino() const
Check if missing neutrinos shall be ignored.
std::vector< int > getSelectionPDGCodes()
Return list of PDG codes of selected particles.
int getNDaughters() const
return number of direct daughters.
void resetMatch()
Reset results from previous call of the match() function.
bool isIgnoreMassive() const
Check if missing massive final state particles shall be ignored.
std::vector< std::vector< std::pair< int, std::string > > > getHierarchyOfSelected()
Function to get hierarchy of selected particles and their names (for python use)
bool isIgnoreGamma() const
Check if missing gammas shall be ignored.
static const DecayDescriptor & s_NULL
Singleton object representing NULL.
std::vector< std::string > getSelectionNames()
Return list of human readable names of selected particles.
std::vector< DecayDescriptor > m_daughters
Direct daughters of the decaying particle.
std::vector< std::vector< std::pair< int, std::string > > > m_hierarchy
Collection of hierarchy paths of selected particles.
std::vector< const Particle * > getSelectionParticles(const Particle *particle)
Get a vector of pointers with selected daughters in the decay tree.
int m_iDaughter_p
ID of the Daughter Particle* matched to this DecayDescriptor.
int match(const T *p, int iDaughter_p)
Internally called by match(Particle*) and match(MCParticle*) function.
bool isIgnoreIntermediate() const
Check if intermediate resonances/particles shall be ignored.
const DecayDescriptorParticle * getMother() const
return mother.
int m_properties
Particle property.
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:32
Class to store reconstructed particles.
Definition: Particle.h:75
int getPDGCode(void) const
Returns PDG code.
Definition: Particle.h:454
@ c_IsIgnoreNeutrino
Is the particle MC matched with the ignore missing neutrino flag set?
Definition: Particle.h:122
@ c_IsIgnoreRadiatedPhotons
Is the particle MC matched with the ignore radiated photon flag set?
Definition: Particle.h:119
@ c_IsIgnoreGamma
Is the particle MC matched with the ignore missing gamma flag set?
Definition: Particle.h:123
@ c_IsIgnoreBrems
Is the particle MC matched with the ignore added Brems gamma flag set?
Definition: Particle.h:124
@ c_IsIgnoreIntermediate
Is the particle MC matched with the ignore intermediate resonances flag set?
Definition: Particle.h:120
@ c_IsIgnoreMassive
Is the particle MC matched with the ignore missing massive particle flag set?
Definition: Particle.h:121
const Particle * getDaughter(unsigned i) const
Returns a pointer to the i-th daughter particle.
Definition: Particle.cc:631
boost::variant< boost::recursive_wrapper< DecayStringDecay >, DecayStringParticle > DecayString
The DecayStringElement can be either a DecayStringDecay or a vector of mother particles.
Definition: DecayString.h:23
bool hasAntiParticle(int pdgCode)
Checks if the particle with given pdg code has an anti-particle or not.
Definition: EvtPDLUtil.cc:12
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24
STL namespace.
Holds the information of a decay.
This class describes the grammar and the syntax elements of decay strings.
Holds the information of a particle in the decay string.