Belle II Software light-2406-ragdoll
ParticleCombiner.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/ParticleCombiner/ParticleCombiner.h>
10
11#include <analysis/DecayDescriptor/DecayDescriptor.h>
12#include <analysis/DecayDescriptor/DecayDescriptorParticle.h>
13
14#include <framework/logging/Logger.h>
15
16#include <mdst/dataobjects/ECLCluster.h>
17
18#include <Math/Vector4D.h>
19
20#include <algorithm>
21
22namespace Belle2 {
28 void ParticleIndexGenerator::init(const std::vector<unsigned>& _sizes)
29 {
30
31 m_numberOfLists = _sizes.size();
32 sizes = _sizes;
33
35
37 for (unsigned int i = 0; i < m_numberOfLists; ++i) {
38 indices[i] = 0;
39 }
40
42 for (unsigned int i = 0; i < m_numberOfLists; ++i)
44
45 if (m_numberOfLists == 0) {
47 }
48
49 }
50
52 {
53
55 return false;
56 }
57
58 if (m_iCombination > 0) {
59
60 // TF SC this does not yet account for double counting so will produce:
61 // { 000, 100, 200, ... 010, 110, .... } even if the first and second
62 // place are the same particle list
63 for (unsigned int i = 0; i < m_numberOfLists; i++) {
64 indices[i]++;
65 if (indices[i] < sizes[i]) break;
66 indices[i] = 0;
67 }
68
69 }
70
72 return true;
73 }
74
75 const std::vector<unsigned int>& ParticleIndexGenerator::getCurrentIndices() const
76 {
77 return indices;
78 }
79
80 void ListIndexGenerator::init(unsigned int _numberOfLists)
81 {
82 m_numberOfLists = _numberOfLists;
84
87
88 }
89
91 {
92
94 return false;
95
96 for (unsigned int i = 0; i < m_numberOfLists; ++i) {
97 bool useSelfConjugatedDaughterList = (m_iCombination & (1 << i));
98 m_types[i] = useSelfConjugatedDaughterList ? ParticleList::c_SelfConjugatedParticle : ParticleList::c_FlavorSpecificParticle;
99 }
100
102 return true;
103 }
104
105 const std::vector<ParticleList::EParticleType>& ListIndexGenerator::getCurrentIndices() const
106 {
107 return m_types;
108 }
109
110 ParticleGenerator::ParticleGenerator(const std::string& decayString, const std::string& cutParameter) : m_iParticleType(0),
111 m_listIndexGenerator(),
112 m_particleIndexGenerator()
113 {
114
115 DecayDescriptor decaydescriptor;
116 bool valid = decaydescriptor.init(decayString);
117 if (!valid)
118 B2ERROR("Invalid input DecayString: " << decayString);
119
120 m_properties = decaydescriptor.getProperty();
121
122 // Mother particle
123 const DecayDescriptorParticle* mother = decaydescriptor.getMother();
124 m_pdgCode = mother->getPDGCode();
125 m_properties |= mother->getProperty();
126
127 // Daughters
128 m_numberOfLists = decaydescriptor.getNDaughters();
129 for (unsigned int i = 0; i < m_numberOfLists; ++i) {
130 // Get the mother of the subdecaystring of the ith daughter
131 // eg. "B -> [D -> K pi] [tau -> pi pi pi]". The 0th daughter is the
132 // *decaystring* D -> K pi whose mother is the D.
133 const DecayDescriptorParticle* daughter = decaydescriptor.getDaughter(i)->getMother();
134 StoreObjPtr<ParticleList> list(daughter->getFullName());
135 m_plists.push_back(list);
136
137 int daughterProperty = daughter->getProperty();
138 m_daughterProperties.push_back(daughterProperty);
139 }
140
141 m_cut = Variable::Cut::compile(cutParameter);
142
143 m_isSelfConjugated = decaydescriptor.isSelfConjugated();
144
145 // check if input lists can contain copies
146 // remember (list_i, list_j) pairs, where list_i and list_j can contain copies
147 m_inputListsCollide = false;
148 m_collidingLists.clear();
149 for (unsigned int i = 0; i < m_numberOfLists; ++i) {
150 const DecayDescriptorParticle* daughter_i = decaydescriptor.getDaughter(i)->getMother();
151 for (unsigned int j = i + 1; j < m_numberOfLists; ++j) {
152 const DecayDescriptorParticle* daughter_j = decaydescriptor.getDaughter(j)->getMother();
153
154 if (abs(daughter_i->getPDGCode()) != abs(daughter_j->getPDGCode()))
155 continue;
156
157 if (daughter_i->getLabel() != daughter_j->getLabel()) {
158 m_inputListsCollide = true;
159 m_collidingLists.emplace_back(i, j);
160 }
161 }
162 }
163
164 }
165
166
167 ParticleGenerator::ParticleGenerator(const DecayDescriptor& decaydescriptor, const std::string& cutParameter) : m_iParticleType(0),
168 m_listIndexGenerator(),
169 m_particleIndexGenerator()
170 {
171 bool valid = decaydescriptor.isInitOK();
172 if (!valid)
173 B2ERROR("Given decaydescriptor failed to initialized");
174
175 m_properties = decaydescriptor.getProperty();
176
177 // Mother particle
178 const DecayDescriptorParticle* mother = decaydescriptor.getMother();
179 m_pdgCode = mother->getPDGCode();
180 m_properties |= mother->getProperty();
181
182 // Daughters
183 m_numberOfLists = decaydescriptor.getNDaughters();
184 for (unsigned int i = 0; i < m_numberOfLists; ++i) {
185 // Get the mother of the subdecaystring of the ith daughter
186 // eg. "B -> [D -> K pi] [tau -> pi pi pi]". The 0th daughter is the
187 // *decaystring* D -> K pi whose mother is the D.
188 const DecayDescriptorParticle* daughter = decaydescriptor.getDaughter(i)->getMother();
189 StoreObjPtr<ParticleList> list(daughter->getFullName());
190 m_plists.push_back(list);
191
192 int daughterProperty = daughter->getProperty();
193 m_daughterProperties.push_back(daughterProperty);
194 }
195
196 m_cut = Variable::Cut::compile(cutParameter);
197
198 m_isSelfConjugated = decaydescriptor.isSelfConjugated();
199
200 // check if input lists can contain copies
201 // remember (list_i, list_j) pairs, where list_i and list_j can contain copies
202 m_inputListsCollide = false;
203 m_collidingLists.clear();
204 for (unsigned int i = 0; i < m_numberOfLists; ++i) {
205 const DecayDescriptorParticle* daughter_i = decaydescriptor.getDaughter(i)->getMother();
206 for (unsigned int j = i + 1; j < m_numberOfLists; ++j) {
207 const DecayDescriptorParticle* daughter_j = decaydescriptor.getDaughter(j)->getMother();
208
209 if (abs(daughter_i->getPDGCode()) != abs(daughter_j->getPDGCode()))
210 continue;
211
212 if (daughter_i->getLabel() != daughter_j->getLabel()) {
213 m_inputListsCollide = true;
214 m_collidingLists.emplace_back(i, j);
215 }
216 }
217 }
218
219 }
220
222 {
223 m_iParticleType = 0;
224
225 m_particleIndexGenerator.init(std::vector<unsigned int> {}); // ParticleIndexGenerator will be initialised on first call
226 m_listIndexGenerator.init(m_numberOfLists); // ListIndexGenerator must be initialised here!
227 m_usedCombinations.clear();
230
233 }
234
236 {
237 m_indicesToUniqueIDs.clear();
238
239 unsigned inputParticlesCount = 0;
240 for (unsigned int i = 0; i < m_numberOfLists; ++i)
241 inputParticlesCount += m_plists[i]->getListSize();
242
243 m_indicesToUniqueIDs.reserve(inputParticlesCount);
244
245 int uniqueID = 1;
246
247 for (auto& collidingList : m_collidingLists) {
248 StoreObjPtr<ParticleList> listA = m_plists[collidingList.first];
249 StoreObjPtr<ParticleList> listB = m_plists[collidingList.second];
250
251 bool sameSign = (listA->getPDGCode() == listB->getPDGCode());
252
253 // if sameSign == true then
254 // 1. compare FS to FS particles in lists A and B
255 // 2. compare anti-FS to anti-FS particles in lists A and B
256 // else
257 // 1. compare FS to anti-FS particles in lists A and B
258 // 2. compare anti-FS to FS particles in lists A and B
259 // and in either case compare
260 // 3. compare SC to SC particles in lists A and B
261 if (sameSign) {
262 fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_FlavorSpecificParticle, false),
263 listB->getList(ParticleList::c_FlavorSpecificParticle, false), uniqueID);
264 fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_FlavorSpecificParticle, true),
265 listB->getList(ParticleList::c_FlavorSpecificParticle, true), uniqueID);
266 } else {
267 fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_FlavorSpecificParticle, false),
268 listB->getList(ParticleList::c_FlavorSpecificParticle, true), uniqueID);
269 fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_FlavorSpecificParticle, true),
270 listB->getList(ParticleList::c_FlavorSpecificParticle, false), uniqueID);
271 }
272
273 fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_SelfConjugatedParticle),
274 listB->getList(ParticleList::c_SelfConjugatedParticle), uniqueID);
275 }
276
277 // assign unique ids to others as well
278 for (unsigned i = 0; i < m_numberOfLists; i++) {
280
281 fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_FlavorSpecificParticle, true), uniqueID);
282 fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_FlavorSpecificParticle, false), uniqueID);
283 fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_SelfConjugatedParticle), uniqueID);
284 }
285 }
286
287 void ParticleGenerator::fillIndicesToUniqueIDMap(const std::vector<int>& listA, const std::vector<int>& listB, int& uniqueID)
288 {
289 const Particle* A, *B;
290 bool copies = false;
291 for (int i : listA) {
292 bool aIsAlreadyIn = m_indicesToUniqueIDs.count(i) ? true : false;
293
294 if (not aIsAlreadyIn)
295 m_indicesToUniqueIDs[ i ] = uniqueID++;
296
297 for (int j : listB) {
298 bool bIsAlreadyIn = m_indicesToUniqueIDs.count(j) ? true : false;
299
300 if (bIsAlreadyIn)
301 continue;
302
303 // are these two particles copies
304 A = m_particleArray[ i ];
305 B = m_particleArray[ j ];
306 copies = B->isCopyOf(A);
307
308 if (copies)
310 }
311 }
312 }
313
314
315 void ParticleGenerator::fillIndicesToUniqueIDMap(const std::vector<int>& listA, int& uniqueID)
316 {
317 for (int i : listA) {
318 bool aIsAlreadyIn = m_indicesToUniqueIDs.count(i) ? true : false;
319
320 if (not aIsAlreadyIn)
321 m_indicesToUniqueIDs[ i ] = uniqueID++;
322 }
323 }
324
325
326 bool ParticleGenerator::loadNext(bool loadAntiParticle)
327 {
328
329 bool loadedNext = false;
336 while (true) {
337 if (m_iParticleType == 0) {
338 loadedNext = loadNextParticle(false);
339 } else if (m_iParticleType == 1 and loadAntiParticle) {
340 loadedNext = loadNextParticle(true);
341 } else if (m_iParticleType == 2) {
342 loadedNext = loadNextSelfConjugatedParticle();
343 } else {
344 return false;
345 }
346
347 if (loadedNext) return true;
348 else ++m_iParticleType;
349
350 if (m_iParticleType == 2) {
351 std::vector<unsigned int> sizes(m_numberOfLists);
352 for (unsigned int i = 0; i < m_numberOfLists; ++i) {
353 sizes[i] = m_plists[i]->getList(ParticleList::c_SelfConjugatedParticle, false).size();
354 }
356 } else {
358 }
359
360 }
361 }
362
363 bool ParticleGenerator::loadNextParticle(bool useAntiParticle)
364 {
365 while (true) {
366
367 // Load next index combination if available
369
370 const auto& m_types = m_listIndexGenerator.getCurrentIndices();
371 const auto& indices = m_particleIndexGenerator.getCurrentIndices();
372
373 for (unsigned int i = 0; i < m_numberOfLists; i++) {
374 const auto& list = m_plists[i]->getList(m_types[i], m_types[i] == ParticleList::c_FlavorSpecificParticle ? useAntiParticle : false);
375 m_indices[i] = list[ indices[i] ];
377 }
378
379 if (not currentCombinationHasDifferentSources()) continue;
380
382 if (!m_cut->check(&m_current_particle))
383 continue;
384
385 if (not currentCombinationIsUnique()) continue;
386
387 if (not currentCombinationIsECLCRUnique()) continue;
388
389 return true;
390 }
391
392 // Load next list combination if available and reset indexCombiner
394 const auto& m_types = m_listIndexGenerator.getCurrentIndices();
395 std::vector<unsigned int> sizes(m_numberOfLists);
396 for (unsigned int i = 0; i < m_numberOfLists; ++i) {
397 sizes[i] = m_plists[i]->getList(m_types[i], m_types[i] == ParticleList::c_FlavorSpecificParticle ? useAntiParticle : false).size();
398 }
400 continue;
401 }
402 return false;
403 }
404
405 }
406
408 {
409 while (true) {
410
411 // Load next index combination if available
413
414 const auto& indices = m_particleIndexGenerator.getCurrentIndices();
415
416 for (unsigned int i = 0; i < m_numberOfLists; i++) {
417 m_indices[i] = m_plists[i]->getList(ParticleList::c_SelfConjugatedParticle, false) [ indices[i] ];
419 }
420
421 if (not currentCombinationHasDifferentSources()) continue;
422
424 if (!m_cut->check(&m_current_particle))
425 continue;
426
427 if (not currentCombinationIsUnique()) continue;
428
429 if (not currentCombinationIsECLCRUnique()) continue;
430
431 return true;
432 }
433
434 return false;
435 }
436
437 }
438
440 {
441 double px = 0;
442 double py = 0;
443 double pz = 0;
444 double E = 0;
445 for (const Particle* d : m_particles) {
446 px += d->getPx();
447 py += d->getPy();
448 pz += d->getPz();
449 E += d->getEnergy();
450 }
451 const ROOT::Math::PxPyPzEVector vec(px, py, pz, E);
452
453 switch (m_iParticleType) {
463 default: B2FATAL("You called getCurrentParticle although loadNext should have returned false!");
464 }
465
466 return Particle(); // This should never happen
467 }
468
470 {
471 return m_current_particle;
472 }
473
475 {
476 std::vector<Particle*> stack = m_particles;
477 static std::vector<int> sources; // stack for particle sources
478 sources.clear();
479
480 // recursively check all daughters and daughters of daughters
481 while (!stack.empty()) {
482 Particle* p = stack.back();
483 stack.pop_back();
484 const std::vector<int>& daughters = p->getDaughterIndices();
485
486 if (daughters.empty()) {
487 int source = p->getMdstSource();
488 for (int i : sources) {
489 if (source == i) return false;
490 }
491 sources.push_back(source);
492 } else {
493 for (int daughter : daughters) stack.push_back(m_particleArray[daughter]);
494 }
495 }
496 return true;
497 }
498
499
501 {
502 std::set<int> indexSet;
503 if (not m_inputListsCollide)
504 indexSet.insert(m_indices.begin(), m_indices.end());
505 else
506 for (unsigned int i = 0; i < m_numberOfLists; i++)
507 indexSet.insert(m_indicesToUniqueIDs.at(m_indices[i]));
508
509 bool elementInserted = m_usedCombinations.insert(indexSet).second;
510 return elementInserted;
511 }
512
513 bool ParticleGenerator::inputListsCollide(const std::pair<unsigned, unsigned>& pair) const
514 {
515 for (const auto& collidingList : m_collidingLists)
516 if (pair == collidingList)
517 return true;
518
519 return false;
520 }
521
523 {
524 unsigned nECLSource = 0;
525 std::vector<Particle*> stack = m_particles;
526 static std::vector<int> connectedregions;
527 static std::vector<ECLCluster::EHypothesisBit> hypotheses;
528 connectedregions.clear();
529 hypotheses.clear();
530
531 // recursively check all daughters and daughters of daughters
532 while (!stack.empty()) {
533 Particle* p = stack.back();
534 stack.pop_back();
535 const std::vector<int>& daughters = p->getDaughterIndices();
536
537 if (daughters.empty()) {
538 // Only test if the particle was created from an ECLCluster at source.
539 // This CAN CHANGE if we change the cluster <--> track matching,
540 // (currently we match nPhotons clusters to all track type particles,
541 // i.e. electrons, pions, kaons etc. We might gain by matching
542 // neutralHadron hypothesis clusters to kaons, for example).
543 //
544 // In the above case one would need to extend the check to ALL
545 // particles with an associated ECLCluster. Then replace the following
546 // active line with two lines...
547 //
548 // auto cluster = p->getECLCluster();
549 // if (cluster) { // then do stuff
550 if (p->getParticleSource() == Particle::EParticleSourceObject::c_ECLCluster) {
551 nECLSource++;
552 auto* cluster = p->getECLCluster();
553 connectedregions.push_back(cluster->getConnectedRegionId());
554 hypotheses.push_back(p->getECLClusterEHypothesisBit());
555 }
556 } else {
557 for (int daughter : daughters) stack.push_back(m_particleArray[daughter]);
558 }
559 }
560
561 // less than two particles from an ECL source is fine
562 // (unless cluster <--> track matching changes)
563 if (nECLSource < 2) return true;
564
565 // yes this is a nested for loop but it's fast, we promise
566 for (unsigned icr = 0; icr < connectedregions.size(); ++icr)
567 for (unsigned jcr = icr + 1; jcr < connectedregions.size(); ++jcr)
568 if (connectedregions[icr] == connectedregions[jcr])
569 if (hypotheses[icr] != hypotheses[jcr]) return false;
570
571 return true;
572 }
573
575 {
576 if (not m_inputListsCollide)
577 return 0;
578
579 if (not m_indicesToUniqueIDs.count(index))
580 return -1;
581
582 return m_indicesToUniqueIDs.at(index);
583 }
585}
Represents a particle in the DecayDescriptor.
int getPDGCode() const
Return PDG code.
int getProperty() const
return property of the particle.
std::string getLabel() const
The label of this particle, "default" returned, when no label set.
The DecayDescriptor stores information about a decay tree or parts of a decay tree.
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).
bool isInitOK() const
Check if the object initialized correctly.
int getNDaughters() const
return number of direct daughters.
int getProperty() const
return property of the particle.
const DecayDescriptorParticle * getMother() const
return mother.
static std::unique_ptr< GeneralCut > compile(const std::string &cut)
Creates an instance of a cut and returns a unique_ptr to it, if you need a copy-able object instead y...
Definition: GeneralCut.h:84
unsigned int m_iCombination
The current position of the combination.
unsigned int m_numberOfLists
Number of lists which are combined.
std::vector< ParticleList::EParticleType > m_types
The current types of sublist of the ParticleLists for this combination.
unsigned int m_nCombinations
The total amount of combinations.
ListIndexGenerator m_listIndexGenerator
listIndexGenerator makes the combinations of the types of sublists of the ParticleLists
unsigned int m_iParticleType
The type of particle which is currently generated.
bool inputListsCollide() const
True if input lists collide (can contain copies of particles in the input lists).
std::unordered_set< std::set< int > > m_usedCombinations
already used combinations (as sets of indices or unique IDs).
Particle m_current_particle
The current Particle object generated by this combiner.
unsigned int m_numberOfLists
Number of lists which are combined.
std::vector< Particle * > m_particles
Pointers to the particle objects of the current combination.
std::unordered_map< int, int > m_indicesToUniqueIDs
map of store array indices of input Particles to their unique IDs.
std::unique_ptr< Variable::Cut > m_cut
cut object which performs the cuts
const StoreArray< Particle > m_particleArray
Global list of particles.
std::vector< int > m_daughterProperties
Daughter's particle properties.
ParticleIndexGenerator m_particleIndexGenerator
particleIndexGenerator makes the combinations of indices stored in the sublists of the ParticleLists
std::vector< StoreObjPtr< ParticleList > > m_plists
particle lists
int m_pdgCode
PDG Code of the particle which is combined.
int m_properties
Particle property.
bool m_inputListsCollide
True if the daughter lists can contain copies of Particles.
std::vector< int > m_indices
Indices stored in the ParticleLists of the current combination.
std::vector< std::pair< unsigned, unsigned > > m_collidingLists
pairs of lists that can contain copies.
bool m_isSelfConjugated
True if the combined particle is self-conjugated.
unsigned int m_iCombination
The current position of the combination.
unsigned int m_numberOfLists
Number of lists which are combined.
std::vector< unsigned int > indices
The indices of the current loaded combination.
std::vector< unsigned int > sizes
The sizes of the particle lists which are combined.
unsigned int m_nCombinations
The total amount of combinations.
Class to store reconstructed particles.
Definition: Particle.h:75
double getPx() const
Returns x component of momentum.
Definition: Particle.h:587
bool isCopyOf(const Particle *oParticle, bool doDetailedComparison=false) const
Returns true if this Particle and oParticle are copies of each other.
Definition: Particle.cc:752
@ c_Unflavored
Is its own antiparticle or we don't know whether it is a particle/antiparticle.
Definition: Particle.h:95
@ c_Flavored
Is either particle or antiparticle.
Definition: Particle.h:96
TClonesArray * getPtr() const
Raw access to the underlying TClonesArray.
Definition: StoreArray.h:311
Type-safe access to single objects in the data store.
Definition: StoreObjPtr.h:96
void init()
Initialises the generator to produce the given type of sublist.
void fillIndicesToUniqueIDMap(const std::vector< int > &listA, const std::vector< int > &listB, int &uniqueID)
Assigns unique IDs to all particles in list A, which do not have the unique ID already assigned.
int getUniqueID(int index) const
Returns the unique ID assigned to Particle with given index from the IndicesToUniqueID map.
void init(const std::vector< unsigned int > &_sizes)
Initialises the generator to produce combinations with the given sizes of each particle list.
bool currentCombinationHasDifferentSources()
Check that all FS particles of a combination differ.
bool currentCombinationIsECLCRUnique()
Check that: if the current combination has at least two particles from an ECL source,...
const std::vector< unsigned int > & getCurrentIndices() const
Returns theindices of the current loaded combination.
bool loadNextSelfConjugatedParticle()
Loads the next combination.
bool currentCombinationIsUnique()
Check that the combination is unique.
bool loadNext(bool loadAntiParticle=true)
Loads the next combination.
void initIndicesToUniqueIDMap()
In the case input daughter particle lists collide (two or more lists contain copies of Particles) the...
const std::vector< ParticleList::EParticleType > & getCurrentIndices() const
Returns the type of the sublist of the current loaded combination.
Particle getCurrentParticle() const
Returns the particle.
ParticleGenerator(const std::string &decayString, const std::string &cutParameter="")
Initialises the generator to produce the given type of sublist.
void init(unsigned int _numberOfLists)
Initialises the generator to produce the given type of sublist.
Particle createCurrentParticle() const
Create current particle object.
bool loadNextParticle(bool useAntiParticle)
Loads the next combination.
bool loadNext()
Loads the next combination.
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24