Belle II Software light-2406-ragdoll
MCMatcherParticlesModule.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/MCMatcherParticles/MCMatcherParticlesModule.h>
10
11// utility
12#include <analysis/utility/MCMatching.h>
13#include <analysis/utility/AnalysisConfiguration.h>
14
15// map
16#include <unordered_map>
17#include <algorithm>
18typedef std::unordered_map<unsigned int, unsigned int> CounterMap;
19
20using namespace std;
21using namespace Belle2;
22
23//-----------------------------------------------------------------
24// Register module
25//-----------------------------------------------------------------
26
27REG_MODULE(MCMatcherParticles);
28
29//-----------------------------------------------------------------
30// Implementation
31//-----------------------------------------------------------------
32
34{
35 setDescription("Performs MC matching (sets relation Particle->MCParticle) for all particles\n"
36 "(and its (grand)^N-daughter particles) in the ParticleList. The relation can\n"
37 "be used in conjunction with MCMatching::MCErrorFlags flags, e.g. using the\n"
38 "isSignal or mcPDG & mcErrors variables.\n"
39 "\n"
40 "In addition to the usual mc matching algorithm the module can run also loose mc\n"
41 "matching. The difference between loose and normal mc matching algorithm is that\n"
42 "the loose algorithm will find the common mother of the majority of daughter\n"
43 "particles while the normal algorithm finds the common mother of all daughters.\n"
44 "The results of loose mc matching algorithm are stored to the following extraInfo\n"
45 "items:\n\n"
46 "- looseMCMotherPDG: PDG code of most common mother\n"
47 "- looseMCMotherIndex: 1-based StoreArray<MCParticle> index of most common mother\n"
48 "- looseMCWrongDaughterN: number of daughters that don't originate from the most"
49 " common mother\n"
50 "- looseMCWrongDaughterPDG: PDG code of the daughter that doesn't originate from\n"
51 " the most common mother (only if looseMCWrongDaughterN = 1)\n"
52 "- looseMCWrongDaughterBiB: 1 if the wrong daughter is Beam Induced Background\n"
53 " Particle");
54
56
57 addParam("listName", m_listName, "Name of the input ParticleList.");
58 addParam("looseMCMatching", m_looseMatching, "Perform loose mc matching", false);
59}
60
61
63{
64 // check that there are MCParticles: shout if not
65 if (!m_mcparticles.isValid()) {
66 B2WARNING("No MCParticles array found!"
67 << " This is obvously fine if you're analysing real data,"
68 << " but you have added the MCMatcher module to your path,"
69 << " did you mean to do this?");
70 return;
71 }
72
73 // if we have MCParticles then continue with the initialisation
76 m_plist.isRequired(m_listName);
77
78 bool legacyAlgorithm = AnalysisConfiguration::instance()->useLegacyMCMatching();
79 B2INFO("MCMatcher module will search for Particle -> MCParticle associations for the ParticleList " << m_listName << ".");
80 if (legacyAlgorithm)
81 B2INFO(" - The MCMatcher will use legacy algorithm suitable for analysis of Belle MC.");
82 else
83 B2INFO(" - The MCMatcher will use default algorithm suitable for analysis of Belle II MC.");
84}
85
86
88{
89 // if no MCParticles then silently skip
90 if (!m_mcparticles.isValid())
91 return;
92 if (!m_plist) {
93 B2ERROR("ParticleList " << m_listName << " not found");
94 return;
95 }
96
97 const unsigned int n = m_plist->getListSize();
98 for (unsigned i = 0; i < n; i++) {
99 const Particle* part = m_plist->getParticle(i);
100
102
103 if (m_looseMatching)
104 setLooseMCMatch(part);
105 }
106}
107
109{
110 if (particle->hasExtraInfo("looseMCMotherPDG")) // nothing to do; already set
111 return;
112
113 // get all FS daughters
114 vector<const Belle2::Particle*> fsDaughters = particle->getFinalStateDaughters();
115
116 // map for counting how many times given mcparticle is mother of daughters
117 CounterMap motherCount;
118
119 for (auto daughter : fsDaughters) {
120 const MCParticle* mcDaughter = daughter->getRelatedTo<MCParticle>();
121 if (!mcDaughter)
122 continue;
123
124 vector<int> genMothers;
125 MCMatching::fillGenMothers(mcDaughter, genMothers);
126
127 for (auto motherIndex : genMothers) {
128 // exclude ROOT particles: Upsilon(nS), virtual photon
129 int motherPDG = m_mcparticles[motherIndex - 1]->getPDG();
130 if ((motherPDG == 553) ||
131 (motherPDG == 100553) ||
132 (motherPDG == 200553) ||
133 (motherPDG == 300553) ||
134 (motherPDG == 9000553) ||
135 (motherPDG == 9010553) ||
136 (motherPDG == 10022))
137 continue;
138
139 motherCount[motherIndex]++;
140 }
141 }
142
143 // find first most common mother
144 auto commonMother = std::max_element
145 (
146 std::begin(motherCount), std::end(motherCount),
147 [](std::pair <unsigned int, unsigned int> p1, std::pair <unsigned int, unsigned int> p2) {
148 bool returnValue = false;
149 if (p1.second < p2.second)
150 returnValue = true;
151 else if (p1.second == p2.second)
152 returnValue = p2.first > p1.first;
153
154 return returnValue;
155 }
156 );
157
158 // No common mother found, all daughters have no associated MC Particle
159 if (commonMother == std::end(motherCount)) {
160 Particle* thisParticle = m_particles[particle->getArrayIndex()];
161 thisParticle->addExtraInfo("looseMCMotherPDG", -1);
162 thisParticle->addExtraInfo("looseMCMotherIndex", -1);
163 thisParticle->addExtraInfo("looseMCWrongDaughterN", -1);
164 thisParticle->addExtraInfo("looseMCWrongDaughterPDG", -1);
165 thisParticle->addExtraInfo("looseMCWrongDaughterBiB", -1);
166 return;
167 }
168
169 const MCParticle* mcMother = m_mcparticles[commonMother->first - 1];
170
171 Particle* thisParticle = m_particles[particle->getArrayIndex()];
172 thisParticle->addExtraInfo("looseMCMotherPDG", mcMother->getPDG());
173 thisParticle->addExtraInfo("looseMCMotherIndex", mcMother->getArrayIndex());
174 thisParticle->addExtraInfo("looseMCWrongDaughterN", fsDaughters.size() - commonMother->second);
175
176 // find out what kind of particle was wrongly added
177 // only for the case where there is only one such particle
178 // This is the most interesting case. If two or more
179 // particles are wrongly added, then a candidate looks
180 // like background
181 int wrongParticlePDG = 0; // PDG code of the wrongly matched particle
182 int wrongParticleBiB = 0; // true (false) if particle doesn't (does) have MCParticle relation
183 if (fsDaughters.size() - commonMother->second == 1) {
184 for (auto daughter : fsDaughters) {
185 const MCParticle* mcDaughter = daughter->getRelatedTo<MCParticle>();
186 if (!mcDaughter) {
187 wrongParticlePDG = daughter->getPDGCode();
188 wrongParticleBiB = 1;
189 }
190
191 vector<int> genMothers;
192 MCMatching::fillGenMothers(mcDaughter, genMothers);
193
194 // check if current daughter descends from common mother
195 if (find(genMothers.begin(), genMothers.end(), commonMother->first) != genMothers.end())
196 continue;
197
198 // daughter is not a child of common mother
199 wrongParticlePDG = daughter->getPDGCode();
200 }
201 }
202
203 thisParticle->addExtraInfo("looseMCWrongDaughterPDG", wrongParticlePDG);
204 thisParticle->addExtraInfo("looseMCWrongDaughterBiB", wrongParticleBiB);
205
206}
void useLegacyMCMatching(const bool flag)
Determines whether to use the legacy MCMatching algorithm (true) or not (false).
static AnalysisConfiguration * instance()
Returns a pointer to the singleton instance.
virtual void initialize() override
Initialize the Module.
virtual void event() override
Event processor.
StoreArray< MCParticle > m_mcparticles
the array of MCParticles.
std::string m_listName
steering variable: name of the input ParticleList
StoreArray< Particle > m_particles
the array of Particles.
bool m_looseMatching
perform loose mc matching
void setLooseMCMatch(const Particle *particle)
Finds common mother of the majority of daughters.
StoreObjPtr< ParticleList > m_plist
the input ParticleList.
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:32
Base class for Modules.
Definition: Module.h:72
void setDescription(const std::string &description)
Sets the description of the module.
Definition: Module.cc:214
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
Definition: Module.cc:208
@ c_ParallelProcessingCertified
This module can be run in parallel processing mode safely (All I/O must be done through the data stor...
Definition: Module.h:80
Class to store reconstructed particles.
Definition: Particle.h:75
std::vector< const Belle2::Particle * > getFinalStateDaughters() const
Returns a vector of pointers to Final State daughter particles.
Definition: Particle.cc:649
bool hasExtraInfo(const std::string &name) const
Return whether the extra info with the given name is set.
Definition: Particle.cc:1266
void addExtraInfo(const std::string &name, double value)
Sets the user-defined data of given name to the given value.
Definition: Particle.cc:1336
int getArrayIndex() const
Returns this object's array index (in StoreArray), or -1 if not found.
bool isRequired(const std::string &name="")
Ensure this array/object has been registered previously.
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
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
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24
STL namespace.
static bool setMCTruth(const Belle2::Particle *particle)
This is the main function of MC matching algorithm.
Definition: MCMatching.cc:86
static void fillGenMothers(const Belle2::MCParticle *mcP, std::vector< int > &genMCPMothers)
Fills vector with array (1-based) indices of all generator ancestors of given MCParticle.
Definition: MCMatching.cc:61