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