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>
18 #include <analysis/utility/AnalysisConfiguration.h>
20 #include <mdst/dataobjects/MCParticle.h>
22 #include <framework/logging/Logger.h>
24 #include <TDatabasePDG.h>
26 #include <boost/variant/get.hpp>
27 #include <boost/spirit/include/qi.hpp>
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;
68 B2WARNING(
"Could not initialise mother particle " << p->m_strName);
72 }
else if (
const DecayStringDecay* d = boost::get< DecayStringDecay > (&s)) {
76 B2WARNING(
"Could not initialise mother particle " << d->m_mother.m_strName);
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=>") {
91 B2WARNING(
"Unknown arrow: " << d->m_strArrow);
97 if (d->m_daughters.empty()) {
101 int nDaughters = d->m_daughters.size();
102 for (
int iDaughter = 0; iDaughter < nDaughters; iDaughter++) {
104 m_isInitOK = daughter.init(d->m_daughters[iDaughter]);
106 B2WARNING(
"Could not initialise daughter!");
114 if ((std::find(d->m_keywords.begin(), d->m_keywords.end(),
"?nu")) != d->m_keywords.end()) {
115 m_properties |= Particle::PropertyFlags::c_IsIgnoreNeutrino;
118 if ((std::find(d->m_keywords.begin(), d->m_keywords.end(),
"?gamma")) != d->m_keywords.end()) {
119 m_properties |= Particle::PropertyFlags::c_IsIgnoreGamma;
122 if ((std::find(d->m_keywords.begin(), d->m_keywords.end(),
"...")) != d->m_keywords.end()) {
123 m_properties |= Particle::PropertyFlags::c_IsIgnoreMassive;
126 if ((std::find(d->m_keywords.begin(), d->m_keywords.end(),
"?addbrems")) != d->m_keywords.end()) {
127 m_properties |= Particle::PropertyFlags::c_IsIgnoreBrems;
144 B2WARNING(
"NULL pointer provided instead of particle.");
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();
154 B2WARNING(
"Template type not supported!");
158 int iPDGCodeCC_p = TDatabasePDG::Instance()->GetParticle(iPDGCode_p)->AntiParticle()->PdgCode();
160 if (abs(iPDGCode_d) != abs(iPDGCode_p))
return 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;
166 const std::vector<T*> daughterList = p->getDaughters();
167 int nDaughters_p = daughterList.size();
177 if (nDaughters_p == 0)
return 0;
187 vector< pair< int, int > > singlematch;
189 vector< pair< int, set<int> > > multimatch;
191 bool isAmbiguities =
false;
193 set<int> matches_global;
196 for (
int iDaughter_d = 0; iDaughter_d <
getNDaughters(); iDaughter_d++) {
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();
207 && iPDGCode_daughter_p == 22)
208 matches_global.insert(jDaughter_p);
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);
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);
227 &&
int(matches_global.size()) != nDaughters_p)
return 0;
232 bool isModified =
true;
233 for (
int iTry = 0; iTry < 20; iTry++) {
235 if (!isModified)
break;
237 for (
auto& itMulti : multimatch) {
238 for (
auto& itSingle : singlematch) {
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!");
268 if (!multimatch.empty()) isAmbiguities =
true;
269 if (isAmbiguities)
return -iCC;
287 vector<const Particle*> selparticles;
289 int motherPDG = abs(particle->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);
298 for (
int iDaughter_d = 0; iDaughter_d < nDaughters_d; ++iDaughter_d) {
300 int iDaughter_p =
m_daughters[iDaughter_d].getMatchedDaughter();
304 if (iDaughter_p < 0) iDaughter_p = iDaughter_d;
305 const Particle* daughter = particle->getDaughter(iDaughter_p);
307 B2WARNING(
"Could not find daughter!");
311 int daughterPDG = abs(daughter->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.");
319 vector<const Particle*> seldaughters =
m_daughters[iDaughter_d].getSelectionParticles(daughter);
320 selparticles.insert(selparticles.end(), seldaughters.begin(), seldaughters.end());
328 std::vector<int> decay, decaybar;
331 int pdg = daughter->getPDGCode();
332 decay.push_back(pdg);
336 std::sort(decay.begin(), decay.end());
337 std::sort(decaybar.begin(), decaybar.end());
345 vector<string> strNames;
348 vector<string> strDaughterNames = daughter.getSelectionNames();
349 int nDaughters = strDaughterNames.size();
350 for (
int iDaughter = 0; iDaughter < nDaughters; iDaughter++) {
361 strNames.insert(strNames.end(), strDaughterNames.begin(), strDaughterNames.end());
366 for (
auto itName = strNames.begin(); itName != strNames.end(); ++itName) {
367 if (count(itName, strNames.end(), *itName) == 1)
continue;
369 string strNameOld = *itName;
370 auto itOccurrence = strNames.begin();
372 while (iOccurrence <= 10) {
374 itOccurrence = find(itOccurrence, strNames.end(), strNameOld);
376 if (itOccurrence == strNames.end())
break;
378 string strNameNew = strNameOld + std::to_string(iOccurrence);
380 if (count(strNames.begin(), strNames.end(), strNameNew) == 0) {
381 *itOccurrence = strNameNew;
386 if (iOccurrence == 10) {
387 B2ERROR(
"DecayDescriptor::getSelectionNames - Something is wrong! More than 10x the same name!");
397 std::vector<std::vector<std::pair<int, std::string>>> hierarchy =
m_hierarchy;
400 std::vector<std::pair<int, std::string>> currentPath;
406 const std::vector<std::pair<int, std::string>>& currentPath)
409 for (std::size_t i = 0; i <
m_daughters.size(); i++) {
410 std::vector<std::pair<int, std::string>> newPath = currentPath;
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);
415 std::vector<std::vector<std::pair<int, std::string>>> hierarchy =
m_hierarchy;