Belle II Software development
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 {
27
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),
113 {
114
115 DecayDescriptor decaydescriptor;
116 bool valid = decaydescriptor.init(decayString);
117 if (!valid) {
118 B2ERROR("Invalid input DecayString: " << decayString);
119 return;
120 }
121
122 m_properties = decaydescriptor.getProperty();
123
124 // Mother particle
125 const DecayDescriptorParticle* mother = decaydescriptor.getMother();
126 m_pdgCode = mother->getPDGCode();
127 m_properties |= mother->getProperty();
128
129 // Daughters
130 m_numberOfLists = decaydescriptor.getNDaughters();
131 for (unsigned int i = 0; i < m_numberOfLists; ++i) {
132 // Get the mother of the subdecaystring of the ith daughter
133 // eg. "B -> [D -> K pi] [tau -> pi pi pi]". The 0th daughter is the
134 // *decaystring* D -> K pi whose mother is the D.
135 const DecayDescriptorParticle* daughter = decaydescriptor.getDaughter(i)->getMother();
136 StoreObjPtr<ParticleList> list(daughter->getFullName());
137 m_plists.push_back(list);
138
139 int daughterProperty = daughter->getProperty();
140 m_daughterProperties.push_back(daughterProperty);
141 }
142
143 m_cut = Variable::Cut::compile(cutParameter);
144
145 m_isSelfConjugated = decaydescriptor.isSelfConjugated();
146
147 // check if input lists can contain copies
148 // remember (list_i, list_j) pairs, where list_i and list_j can contain copies
149 m_inputListsCollide = false;
150 m_collidingLists.clear();
151 for (unsigned int i = 0; i < m_numberOfLists; ++i) {
152 const DecayDescriptorParticle* daughter_i = decaydescriptor.getDaughter(i)->getMother();
153 for (unsigned int j = i + 1; j < m_numberOfLists; ++j) {
154 const DecayDescriptorParticle* daughter_j = decaydescriptor.getDaughter(j)->getMother();
155
156 if (abs(daughter_i->getPDGCode()) != abs(daughter_j->getPDGCode()))
157 continue;
158
159 if (daughter_i->getLabel() != daughter_j->getLabel()) {
160 m_inputListsCollide = true;
161 m_collidingLists.emplace_back(i, j);
162 }
163 }
164 }
165
166 }
167
168
169 ParticleGenerator::ParticleGenerator(const DecayDescriptor& decaydescriptor, const std::string& cutParameter) : m_iParticleType(0),
172 {
173 bool valid = decaydescriptor.isInitOK();
174 if (!valid)
175 B2ERROR("Given decaydescriptor failed to initialized");
176
177 m_properties = decaydescriptor.getProperty();
178
179 // Mother particle
180 const DecayDescriptorParticle* mother = decaydescriptor.getMother();
181 m_pdgCode = mother->getPDGCode();
182 m_properties |= mother->getProperty();
183
184 // Daughters
185 m_numberOfLists = decaydescriptor.getNDaughters();
186 for (unsigned int i = 0; i < m_numberOfLists; ++i) {
187 // Get the mother of the subdecaystring of the ith daughter
188 // eg. "B -> [D -> K pi] [tau -> pi pi pi]". The 0th daughter is the
189 // *decaystring* D -> K pi whose mother is the D.
190 const DecayDescriptorParticle* daughter = decaydescriptor.getDaughter(i)->getMother();
191 StoreObjPtr<ParticleList> list(daughter->getFullName());
192 m_plists.push_back(list);
193
194 int daughterProperty = daughter->getProperty();
195 m_daughterProperties.push_back(daughterProperty);
196 }
197
198 m_cut = Variable::Cut::compile(cutParameter);
199
200 m_isSelfConjugated = decaydescriptor.isSelfConjugated();
201
202 // check if input lists can contain copies
203 // remember (list_i, list_j) pairs, where list_i and list_j can contain copies
204 m_inputListsCollide = false;
205 m_collidingLists.clear();
206 for (unsigned int i = 0; i < m_numberOfLists; ++i) {
207 const DecayDescriptorParticle* daughter_i = decaydescriptor.getDaughter(i)->getMother();
208 for (unsigned int j = i + 1; j < m_numberOfLists; ++j) {
209 const DecayDescriptorParticle* daughter_j = decaydescriptor.getDaughter(j)->getMother();
210
211 if (abs(daughter_i->getPDGCode()) != abs(daughter_j->getPDGCode()))
212 continue;
213
214 if (daughter_i->getLabel() != daughter_j->getLabel()) {
215 m_inputListsCollide = true;
216 m_collidingLists.emplace_back(i, j);
217 }
218 }
219 }
220
221 }
222
224 {
225 m_iParticleType = 0;
226
227 m_particleIndexGenerator.init(std::vector<unsigned int> {}); // ParticleIndexGenerator will be initialised on first call
228 m_listIndexGenerator.init(m_numberOfLists); // ListIndexGenerator must be initialised here!
229 m_usedCombinations.clear();
232
235 }
236
238 {
239 m_indicesToUniqueIDs.clear();
240
241 unsigned inputParticlesCount = 0;
242 for (unsigned int i = 0; i < m_numberOfLists; ++i)
243 inputParticlesCount += m_plists[i]->getListSize();
244
245 m_indicesToUniqueIDs.reserve(inputParticlesCount);
246
247 int uniqueID = 1;
248
249 for (auto& collidingList : m_collidingLists) {
250 StoreObjPtr<ParticleList> listA = m_plists[collidingList.first];
251 StoreObjPtr<ParticleList> listB = m_plists[collidingList.second];
252
253 bool sameSign = (listA->getPDGCode() == listB->getPDGCode());
254
255 // if sameSign == true then
256 // 1. compare FS to FS particles in lists A and B
257 // 2. compare anti-FS to anti-FS particles in lists A and B
258 // else
259 // 1. compare FS to anti-FS particles in lists A and B
260 // 2. compare anti-FS to FS particles in lists A and B
261 // and in either case compare
262 // 3. compare SC to SC particles in lists A and B
263 if (sameSign) {
264 fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_FlavorSpecificParticle, false),
265 listB->getList(ParticleList::c_FlavorSpecificParticle, false), uniqueID);
266 fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_FlavorSpecificParticle, true),
267 listB->getList(ParticleList::c_FlavorSpecificParticle, true), uniqueID);
268 } else {
269 fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_FlavorSpecificParticle, false),
270 listB->getList(ParticleList::c_FlavorSpecificParticle, true), uniqueID);
271 fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_FlavorSpecificParticle, true),
272 listB->getList(ParticleList::c_FlavorSpecificParticle, false), uniqueID);
273 }
274
275 fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_SelfConjugatedParticle),
276 listB->getList(ParticleList::c_SelfConjugatedParticle), uniqueID);
277 }
278
279 // assign unique ids to others as well
280 for (unsigned i = 0; i < m_numberOfLists; i++) {
282
283 fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_FlavorSpecificParticle, true), uniqueID);
284 fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_FlavorSpecificParticle, false), uniqueID);
285 fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_SelfConjugatedParticle), uniqueID);
286 }
287 }
288
289 void ParticleGenerator::fillIndicesToUniqueIDMap(const std::vector<int>& listA, const std::vector<int>& listB, int& uniqueID)
290 {
291 const Particle* A, *B;
292 bool copies = false;
293 for (int i : listA) {
294 bool aIsAlreadyIn = m_indicesToUniqueIDs.count(i) ? true : false;
295
296 if (not aIsAlreadyIn)
297 m_indicesToUniqueIDs[ i ] = uniqueID++;
298
299 for (int j : listB) {
300 bool bIsAlreadyIn = m_indicesToUniqueIDs.count(j) ? true : false;
301
302 if (bIsAlreadyIn)
303 continue;
304
305 // are these two particles copies
306 A = m_particleArray[ i ];
307 B = m_particleArray[ j ];
308 copies = B->isCopyOf(A);
309
310 if (copies)
312 }
313 }
314 }
315
316
317 void ParticleGenerator::fillIndicesToUniqueIDMap(const std::vector<int>& listA, int& uniqueID)
318 {
319 for (int i : listA) {
320 bool aIsAlreadyIn = m_indicesToUniqueIDs.count(i) ? true : false;
321
322 if (not aIsAlreadyIn)
323 m_indicesToUniqueIDs[ i ] = uniqueID++;
324 }
325 }
326
327
328 bool ParticleGenerator::loadNext(bool loadAntiParticle)
329 {
330
331 bool loadedNext = false;
338 while (true) {
339 if (m_iParticleType == 0) {
340 loadedNext = loadNextParticle(false);
341 } else if (m_iParticleType == 1 and loadAntiParticle) {
342 loadedNext = loadNextParticle(true);
343 } else if (m_iParticleType == 2) {
344 loadedNext = loadNextSelfConjugatedParticle();
345 } else {
346 return false;
347 }
348
349 if (loadedNext) return true;
350 else ++m_iParticleType;
351
352 if (m_iParticleType == 2) {
353 std::vector<unsigned int> sizes(m_numberOfLists);
354 for (unsigned int i = 0; i < m_numberOfLists; ++i) {
355 sizes[i] = m_plists[i]->getList(ParticleList::c_SelfConjugatedParticle, false).size();
356 }
357 m_particleIndexGenerator.init(sizes);
358 } else {
360 }
361
362 }
363 }
364
365 bool ParticleGenerator::loadNextParticle(bool useAntiParticle)
366 {
367 while (true) {
368
369 // Load next index combination if available
370 if (m_particleIndexGenerator.loadNext()) {
371
372 const auto& m_types = m_listIndexGenerator.getCurrentIndices();
373 const auto& indices = m_particleIndexGenerator.getCurrentIndices();
374
375 for (unsigned int i = 0; i < m_numberOfLists; i++) {
376 const auto& list = m_plists[i]->getList(m_types[i], m_types[i] == ParticleList::c_FlavorSpecificParticle ? useAntiParticle : false);
377 m_indices[i] = list[ indices[i] ];
379 }
380
381 if (not currentCombinationHasDifferentSources()) continue;
382
384 if (!m_cut->check(&m_current_particle))
385 continue;
386
387 if (not currentCombinationIsUnique()) continue;
388
389 if (not currentCombinationIsECLCRUnique()) continue;
390
391 return true;
392 }
393
394 // Load next list combination if available and reset indexCombiner
395 if (m_listIndexGenerator.loadNext()) {
396 const auto& m_types = m_listIndexGenerator.getCurrentIndices();
397 std::vector<unsigned int> sizes(m_numberOfLists);
398 for (unsigned int i = 0; i < m_numberOfLists; ++i) {
399 sizes[i] = m_plists[i]->getList(m_types[i], m_types[i] == ParticleList::c_FlavorSpecificParticle ? useAntiParticle : false).size();
400 }
401 m_particleIndexGenerator.init(sizes);
402 continue;
403 }
404 return false;
405 }
406
407 }
408
410 {
411 while (true) {
412
413 // Load next index combination if available
414 if (m_particleIndexGenerator.loadNext()) {
415
416 const auto& indices = m_particleIndexGenerator.getCurrentIndices();
417
418 for (unsigned int i = 0; i < m_numberOfLists; i++) {
419 m_indices[i] = m_plists[i]->getList(ParticleList::c_SelfConjugatedParticle, false) [ indices[i] ];
421 }
422
423 if (not currentCombinationHasDifferentSources()) continue;
424
426 if (!m_cut->check(&m_current_particle))
427 continue;
428
429 if (not currentCombinationIsUnique()) continue;
430
431 if (not currentCombinationIsECLCRUnique()) continue;
432
433 return true;
434 }
435
436 return false;
437 }
438
439 }
440
442 {
443 double px = 0;
444 double py = 0;
445 double pz = 0;
446 double E = 0;
447 for (const Particle* d : m_particles) {
448 px += d->getPx();
449 py += d->getPy();
450 pz += d->getPz();
451 E += d->getEnergy();
452 }
453 const ROOT::Math::PxPyPzEVector vec(px, py, pz, E);
454
455 switch (m_iParticleType) {
458 m_particleArray.getPtr());
461 m_particleArray.getPtr());
464 m_particleArray.getPtr());
465 default: B2FATAL("You called getCurrentParticle although loadNext should have returned false!");
466 }
467
468 return Particle(); // This should never happen
469 }
470
475
477 {
478 std::vector<Particle*> stack = m_particles;
479 static std::vector<int> sources; // stack for particle sources
480 sources.clear();
481
482 // recursively check all daughters and daughters of daughters
483 while (!stack.empty()) {
484 Particle* p = stack.back();
485 stack.pop_back();
486 const std::vector<int>& daughters = p->getDaughterIndices();
487
488 if (daughters.empty()) {
489 int source = p->getMdstSource();
490 for (int i : sources) {
491 if (source == i) return false;
492 }
493 sources.push_back(source);
494 } else {
495 for (int daughter : daughters) stack.push_back(m_particleArray[daughter]);
496 }
497 }
498 return true;
499 }
500
501
503 {
504 std::set<int> indexSet;
505 if (not m_inputListsCollide)
506 indexSet.insert(m_indices.begin(), m_indices.end());
507 else
508 for (unsigned int i = 0; i < m_numberOfLists; i++)
509 indexSet.insert(m_indicesToUniqueIDs.at(m_indices[i]));
510
511 bool elementInserted = m_usedCombinations.insert(indexSet).second;
512 return elementInserted;
513 }
514
515 bool ParticleGenerator::inputListsCollide(const std::pair<unsigned, unsigned>& pair) const
516 {
517 for (const auto& collidingList : m_collidingLists)
518 if (pair == collidingList)
519 return true;
520
521 return false;
522 }
523
525 {
526 unsigned nECLSource = 0;
527 std::vector<Particle*> stack = m_particles;
528 static std::vector<int> connectedregions;
529 static std::vector<ECLCluster::EHypothesisBit> hypotheses;
530 connectedregions.clear();
531 hypotheses.clear();
532
533 // recursively check all daughters and daughters of daughters
534 while (!stack.empty()) {
535 Particle* p = stack.back();
536 stack.pop_back();
537 const std::vector<int>& daughters = p->getDaughterIndices();
538
539 if (daughters.empty()) {
540 // Only test if the particle was created from an ECLCluster at source.
541 // This CAN CHANGE if we change the cluster <--> track matching,
542 // (currently we match nPhotons clusters to all track type particles,
543 // i.e. electrons, pions, kaons etc. We might gain by matching
544 // neutralHadron hypothesis clusters to kaons, for example).
545 //
546 // In the above case one would need to extend the check to ALL
547 // particles with an associated ECLCluster. Then replace the following
548 // active line with two lines...
549 //
550 // auto cluster = p->getECLCluster();
551 // if (cluster) { // then do stuff
552 if (p->getParticleSource() == Particle::EParticleSourceObject::c_ECLCluster) {
553 nECLSource++;
554 auto* cluster = p->getECLCluster();
555 connectedregions.push_back(cluster->getConnectedRegionId());
556 hypotheses.push_back(p->getECLClusterEHypothesisBit());
557 }
558 } else {
559 for (int daughter : daughters) stack.push_back(m_particleArray[daughter]);
560 }
561 }
562
563 // less than two particles from an ECL source is fine
564 // (unless cluster <--> track matching changes)
565 if (nECLSource < 2) return true;
566
567 // yes this is a nested for loop but it's fast, we promise
568 for (unsigned icr = 0; icr < connectedregions.size(); ++icr)
569 for (unsigned jcr = icr + 1; jcr < connectedregions.size(); ++jcr)
570 if (connectedregions[icr] == connectedregions[jcr])
571 if (hypotheses[icr] != hypotheses[jcr]) return false;
572
573 return true;
574 }
575
577 {
578 if (not m_inputListsCollide)
579 return 0;
580
581 if (not m_indicesToUniqueIDs.count(index))
582 return -1;
583
584 return m_indicesToUniqueIDs.at(index);
585 }
586
587}
R E
internal precision of FFTW codelets
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)
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
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
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:76
@ c_Unflavored
Is its own antiparticle or we don't know whether it is a particle/antiparticle.
Definition Particle.h:97
@ c_Flavored
Is either particle or antiparticle.
Definition Particle.h:98
Type-safe access to single objects in the data store.
Definition StoreObjPtr.h:96
STL class.
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.