Belle II Software  release-05-02-19
Particle.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2012-2019 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Anze Zupanc, Marko Staric, Christian Pulvermacher, *
7  * Sam Cunliffe, Torben Ferber, Thomas Kuhr, *
8  * Umberto Tamponi *
9  * *
10  * This software is provided "as is" without any warranty. *
11  **************************************************************************/
12 
13 #include <analysis/dataobjects/Particle.h>
14 #include <analysis/dataobjects/ParticleExtraInfoMap.h>
15 
16 #include <analysis/ClusterUtility/ClusterUtils.h>
17 
18 #include <mdst/dataobjects/KLMCluster.h>
19 #include <mdst/dataobjects/MCParticle.h>
20 #include <mdst/dataobjects/PIDLikelihood.h>
21 #include <mdst/dataobjects/Track.h>
22 #include <mdst/dataobjects/TrackFitResult.h>
23 #include <mdst/dataobjects/V0.h>
24 #include <mdst/dbobjects/CollisionBoostVector.h>
25 #include <mdst/dbobjects/CollisionInvariantMass.h>
26 
27 #include <framework/datastore/StoreArray.h>
28 #include <framework/datastore/StoreObjPtr.h>
29 #include <framework/database/DBObjPtr.h>
30 #include <framework/logging/Logger.h>
31 #include <framework/utilities/HTML.h>
32 #include <framework/utilities/Conversion.h>
33 
34 #include <TClonesArray.h>
35 #include <TDatabasePDG.h>
36 #include <TMatrixFSym.h>
37 
38 #include <iostream>
39 #include <iomanip>
40 #include <stdexcept>
41 #include <queue>
42 
43 #include <boost/algorithm/string.hpp>
44 
45 using namespace Belle2;
46 
48  m_pdgCode(0), m_mass(0), m_px(0), m_py(0), m_pz(0), m_x(0), m_y(0), m_z(0),
49  m_pValue(nan("")), m_flavorType(c_Unflavored), m_particleSource(c_Undefined), m_mdstIndex(0), m_properties(0),
50  m_arrayPointer(nullptr)
51 {
53 }
54 
55 
56 Particle::Particle(const TLorentzVector& momentum, const int pdgCode) :
57  m_pdgCode(pdgCode), m_mass(0), m_px(0), m_py(0), m_pz(0), m_x(0), m_y(0), m_z(0),
58  m_pValue(-1), m_flavorType(c_Unflavored), m_particleSource(c_Undefined), m_mdstIndex(0), m_properties(0), m_arrayPointer(nullptr)
59 {
60  setFlavorType();
61  set4Vector(momentum);
63 }
64 
65 
66 Particle::Particle(const TLorentzVector& momentum,
67  const int pdgCode,
68  EFlavorType flavorType,
69  const EParticleSourceObject source,
70  const unsigned mdstIndex) :
71  m_pdgCode(pdgCode), m_mass(0), m_px(0), m_py(0), m_pz(0), m_x(0), m_y(0), m_z(0),
72  m_pValue(-1), m_flavorType(flavorType), m_particleSource(source), m_properties(0), m_arrayPointer(nullptr)
73 {
74  if (flavorType == c_Unflavored and pdgCode < 0)
75  m_pdgCode = -pdgCode;
76 
77  setMdstArrayIndex(mdstIndex);
78  set4Vector(momentum);
80 }
81 
82 
83 Particle::Particle(const TLorentzVector& momentum,
84  const int pdgCode,
85  EFlavorType flavorType,
86  const std::vector<int>& daughterIndices,
87  TClonesArray* arrayPointer) :
88  m_pdgCode(0), m_mass(0), m_px(0), m_py(0), m_pz(0), m_x(0), m_y(0), m_z(0),
89  m_pValue(-1),
90  m_daughterIndices(daughterIndices),
91  m_flavorType(c_Unflavored), m_particleSource(c_Undefined), m_mdstIndex(0),
92  m_properties(0), m_arrayPointer(arrayPointer)
93 {
94  m_pdgCode = pdgCode;
95  m_flavorType = flavorType;
96  if (flavorType == c_Unflavored and pdgCode < 0)
97  m_pdgCode = -pdgCode;
98  set4Vector(momentum);
100 
101  if (!daughterIndices.empty()) {
102  m_particleSource = c_Composite;
103  if (getArrayPointer() == nullptr) {
104  B2FATAL("Composite Particle (with daughters) was constructed outside StoreArray without specifying daughter array!");
105  }
106  for (unsigned int i = 0; i < m_daughterIndices.size(); i++) {
107  m_daughterProperties.push_back(Particle::PropertyFlags::c_Ordinary);
108  }
109  }
110 }
111 
112 Particle::Particle(const TLorentzVector& momentum,
113  const int pdgCode,
114  EFlavorType flavorType,
115  const std::vector<int>& daughterIndices,
116  const int properties,
117  TClonesArray* arrayPointer) :
118  m_pdgCode(0), m_mass(0), m_px(0), m_py(0), m_pz(0), m_x(0), m_y(0), m_z(0),
119  m_pValue(-1),
120  m_daughterIndices(daughterIndices),
121  m_flavorType(c_Unflavored), m_particleSource(c_Undefined), m_mdstIndex(0),
122  m_arrayPointer(arrayPointer)
123 {
124  m_pdgCode = pdgCode;
125  m_flavorType = flavorType;
126  if (flavorType == c_Unflavored and pdgCode < 0)
127  m_pdgCode = -pdgCode;
128  set4Vector(momentum);
130  m_properties = properties;
131 
132  if (!daughterIndices.empty()) {
133  m_particleSource = c_Composite;
134  if (getArrayPointer() == nullptr) {
135  B2FATAL("Composite Particle (with daughters) was constructed outside StoreArray without specifying daughter array!");
136  }
137  for (unsigned int i = 0; i < m_daughterIndices.size(); i++) {
138  m_daughterProperties.push_back(Particle::PropertyFlags::c_Ordinary);
139  }
140  }
141 }
142 
143 
144 Particle::Particle(const TLorentzVector& momentum,
145  const int pdgCode,
146  EFlavorType flavorType,
147  const std::vector<int>& daughterIndices,
148  const int properties,
149  const std::vector<int>& daughterProperties,
150  TClonesArray* arrayPointer) :
151  m_pdgCode(0), m_mass(0), m_px(0), m_py(0), m_pz(0), m_x(0), m_y(0), m_z(0),
152  m_pValue(-1),
153  m_daughterIndices(daughterIndices),
154  m_flavorType(c_Unflavored), m_particleSource(c_Undefined), m_mdstIndex(0),
155  m_daughterProperties(daughterProperties),
156  m_arrayPointer(arrayPointer)
157 {
158  m_pdgCode = pdgCode;
159  m_flavorType = flavorType;
160  if (flavorType == c_Unflavored and pdgCode < 0)
161  m_pdgCode = -pdgCode;
162  set4Vector(momentum);
164  m_properties = properties;
165 
166  if (!daughterIndices.empty()) {
167  m_particleSource = c_Composite;
168  if (getArrayPointer() == nullptr) {
169  B2FATAL("Composite Particle (with daughters) was constructed outside StoreArray without specifying daughter array!");
170  }
171  }
172 }
173 
174 
176  const Const::ChargedStable& chargedStable) :
177  Particle(track ? track->getArrayIndex() : 0, track ? track->getTrackFitResultWithClosestMass(chargedStable) : nullptr,
178  chargedStable)
179 {
180 }
181 
182 Particle::Particle(const int trackArrayIndex,
183  const TrackFitResult* trackFit,
184  const Const::ChargedStable& chargedStable) :
185  m_pdgCode(0), m_mass(0), m_px(0), m_py(0), m_pz(0), m_x(0), m_y(0), m_z(0),
186  m_pValue(-1), m_flavorType(c_Unflavored), m_particleSource(c_Undefined), m_mdstIndex(0), m_properties(0), m_arrayPointer(nullptr)
187 {
188  if (!trackFit) return;
189 
190  m_flavorType = c_Flavored; //tracks are charged
191  m_particleSource = c_Track;
192 
193  setMdstArrayIndex(trackArrayIndex);
194 
196  m_pdgCode = generatePDGCodeFromCharge(trackFit->getChargeSign(), chargedStable);
197 
198  // set mass
199  if (TDatabasePDG::Instance()->GetParticle(m_pdgCode) == nullptr)
200  B2FATAL("PDG=" << m_pdgCode << " ***code unknown to TDatabasePDG");
201  m_mass = TDatabasePDG::Instance()->GetParticle(m_pdgCode)->Mass() ;
202 
203  // set momentum, position and error matrix
205 }
206 
207 //FIXME: Deprecated, to be removed after release-05
208 Particle::Particle(const int trackArrayIndex,
209  const TrackFitResult* trackFit,
210  const Const::ChargedStable& chargedStable,
211  const Const::ChargedStable& chargedStableUsedForFit) :
212  m_pdgCode(0), m_mass(0), m_px(0), m_py(0), m_pz(0), m_x(0), m_y(0), m_z(0),
213  m_pValue(-1), m_flavorType(c_Unflavored), m_particleSource(c_Undefined), m_mdstIndex(0), m_properties(0), m_arrayPointer(nullptr)
214 {
215  if (!trackFit) return;
216 
217  m_flavorType = c_Flavored; //tracks are charged
218  m_particleSource = c_Track;
219 
220  setMdstArrayIndex(trackArrayIndex);
221 
222  m_pdgCodeUsedForFit = chargedStableUsedForFit.getPDGCode();
223  m_pdgCode = generatePDGCodeFromCharge(trackFit->getChargeSign(), chargedStable);
224 
225  // set mass
226  if (TDatabasePDG::Instance()->GetParticle(m_pdgCode) == nullptr)
227  B2FATAL("PDG=" << m_pdgCode << " ***code unknown to TDatabasePDG");
228  m_mass = TDatabasePDG::Instance()->GetParticle(m_pdgCode)->Mass() ;
229 
230  // set momentum, position and error matrix
232 }
233 
234 
235 Particle::Particle(const ECLCluster* eclCluster, const Const::ParticleType& type) :
236  m_pdgCode(type.getPDGCode()), m_mass(type.getMass()), m_px(0), m_py(0), m_pz(0), m_x(0), m_y(0), m_z(0),
237  m_pValue(-1), m_flavorType(c_Unflavored), m_particleSource(c_Undefined), m_mdstIndex(0), m_properties(0), m_arrayPointer(nullptr)
238 {
239  if (!eclCluster) return;
240 
241  setFlavorType();
242 
243  // returns default vertex from clusterutils (from beam parameters if available)
244  // leave it like that for the moment to make it transparent
245  ClusterUtils C;
246  const TVector3 clustervertex = C.GetIPPosition();
247  setVertex(clustervertex);
248 
249  const TLorentzVector clustermom = C.Get4MomentumFromCluster(eclCluster, clustervertex, getECLClusterEHypothesisBit());
250  m_px = clustermom.Px();
251  m_py = clustermom.Py();
252  m_pz = clustermom.Pz();
253 
254  m_particleSource = c_ECLCluster;
255  setMdstArrayIndex(eclCluster->getArrayIndex());
256 
257  // set Chi^2 probability:
258  //TODO: gamma quality can be written here
259  m_pValue = 1;
260 
261  // Get covariance matrix of IP distribution.
262  const TMatrixDSym clustervertexcovmat = C.GetIPPositionCovarianceMatrix();
263 
264  // Set error matrix.
265  TMatrixDSym clustercovmat = C.GetCovarianceMatrix7x7FromCluster(eclCluster, clustervertex, clustervertexcovmat,
267  storeErrorMatrix(clustercovmat);
268 }
269 
270 
271 Particle::Particle(const KLMCluster* klmCluster, const int pdgCode) :
272  m_pdgCode(0), m_mass(0), m_px(0), m_py(0), m_pz(0), m_x(0), m_y(0), m_z(0),
273  m_pValue(-1), m_flavorType(c_Unflavored), m_particleSource(c_Undefined), m_mdstIndex(0), m_properties(0), m_arrayPointer(nullptr)
274 {
275  if (!klmCluster) return;
276 
277  m_pdgCode = pdgCode;
278  setFlavorType();
279 
280  set4Vector(klmCluster->getMomentum());
281  setVertex(klmCluster->getPosition());
282  updateMass(m_pdgCode); // KLMCluster internally use Klong mass, overwrite here to allow neutrons
283 
284  m_particleSource = c_KLMCluster;
285  setMdstArrayIndex(klmCluster->getArrayIndex());
286 
287  // set Chi^2 probability:
288  //TODO: KL quality can be written here
289  m_pValue = -1;
290 
291  // TODO: set error matrix
293  //storeErrorMatrix(klmCluster->getErrorMatrix());
294 }
295 
296 
297 Particle::Particle(const MCParticle* mcParticle) :
298  m_pdgCode(0), m_mass(0), m_px(0), m_py(0), m_pz(0), m_x(0), m_y(0), m_z(0),
299  m_pValue(-1), m_flavorType(c_Unflavored), m_particleSource(c_Undefined), m_mdstIndex(0), m_properties(0), m_arrayPointer(nullptr)
300 {
301  if (!mcParticle) return;
302 
303  m_pdgCode = mcParticle->getPDG();
304  m_particleSource = c_MCParticle; // TODO: what about daughters if not FS particle?
305 
306  setMdstArrayIndex(mcParticle->getArrayIndex());
307 
308 
309  setFlavorType();
310 
311  // mass and momentum
312  m_mass = mcParticle->getMass();
313  m_px = mcParticle->getMomentum().Px();
314  m_py = mcParticle->getMomentum().Py();
315  m_pz = mcParticle->getMomentum().Pz();
316  // production vertex
317  // TODO: good only for FS particles, for composite we must use decay vertex
318  setVertex(mcParticle->getVertex());
319 
321 }
322 
323 
324 Particle::~Particle() = default;
325 
326 void Particle::setMdstArrayIndex(const int arrayIndex)
327 {
328  m_mdstIndex = arrayIndex;
329 
330  // set the identifier
331  if (m_particleSource == c_ECLCluster) {
332  const ECLCluster* cluster = this->getECLCluster();
333  if (cluster) {
334  const int crid = cluster->getConnectedRegionId();
335  const int clusterid = cluster->getClusterId();
336  m_identifier = 1000 * crid + clusterid;
337  } else {
338  B2ERROR("Particle is of type = ECLCluster has identifier not set and no relation to ECLCluster.\n"
339  "This has happen because old microDST is analysed with newer version of software.");
340  }
341  } else {
343  }
344 }
345 
346 
348 {
349  // Is identifier already set
350  if (m_identifier > -1)
351  return m_identifier + (m_particleSource << 24);
352 
353  // Identifier is not set.
354  int identifier = 0;
355  if (m_particleSource == c_ECLCluster) {
356  const ECLCluster* cluster = this->getECLCluster();
357  if (cluster) {
358  const int crid = cluster->getConnectedRegionId();
359  const int clusterid = cluster->getClusterId();
360  identifier = 1000 * crid + clusterid;
361  } else {
362  B2ERROR("Particle is of type = ECLCluster has identifier not set and no relation to ECLCluster.\n"
363  "This has happen because old microDST is analysed with newer version of software.");
364  }
365  } else {
366  identifier = m_mdstIndex;
367  }
368 
369  return identifier + (m_particleSource << 24);
370 }
371 
372 
373 void Particle::setMomentumVertexErrorMatrix(const TMatrixFSym& m)
374 {
375  // check if provided Error Matrix is of dimension 7x7
376  // if not, reset the error matrix and print warning
377  if (m.GetNrows() != c_DimMatrix || m.GetNcols() != c_DimMatrix) {
379  B2WARNING("Error Matrix is not 7x7 ");
380  return;
381  }
382  storeErrorMatrix(m);
383 }
384 
385 
387 {
388  TMatrixFSym m(c_DimMatrix);
389 
390  int element = 0;
391  for (int irow = 0; irow < c_DimMatrix; irow++) {
392  for (int icol = irow; icol < c_DimMatrix; icol++) {
393  m(irow, icol) = m(icol, irow) = m_errMatrix[element];
394  element++;
395  }
396  }
397  return m;
398 }
399 
400 
402 {
403  TMatrixFSym mom;
404  const TMatrixFSym& full = getMomentumVertexErrorMatrix();
405 
406  // get 4x4 (momentum) submatrix from the full error matrix
407  // momentum related elements are in [0,...,3]x[0,...,3] block
408  full.GetSub(0, 3, mom, "S");
409 
410  return mom;
411 }
412 
414 {
415  TMatrixFSym pos;
416  const TMatrixFSym& full = getMomentumVertexErrorMatrix();
417 
418  // get 3x3 (position) submatrix from the full error matrix
419  // vertex related elements are in [4,5,6]x[4,5,6] block
420  full.GetSub(4, 6, pos, "S");
421 
422  return pos;
423 }
424 
425 float Particle::getCosHelicity(const Particle* mother) const
426 {
427  // boost vector to the rest frame of the particle
428  TVector3 boost = -get4Vector().BoostVector();
429 
430  // momentum of the mother in the particle's rest frame
431  TLorentzVector pMother;
432  if (mother) {
433  pMother = mother->get4Vector();
434  } else {
435  static DBObjPtr<CollisionBoostVector> cmsBoost;
436  static DBObjPtr<CollisionInvariantMass> cmsMass;
437  pMother.SetE(cmsMass->getMass());
438  pMother.Boost(cmsBoost->getBoost());
439  }
440  pMother.Boost(boost);
441 
442  // momentum of the daughter (or normal vector) in the particle's rest frame
443  TLorentzVector pDaughter;
444  if (getNDaughters() == 2) { // two body decay
445  pDaughter = getDaughter(0)->get4Vector();
446  pDaughter.Boost(boost);
447  } else if (getNDaughters() == 3) {
448  if (getPDGCode() == Const::pi0.getPDGCode()) { // pi0 Dalitz decay
449  for (auto& daughter : getDaughters()) {
450  if (daughter->getPDGCode() == Const::photon.getPDGCode()) {
451  pDaughter = daughter->get4Vector();
452  }
453  }
454  pDaughter.Boost(boost);
455  } else { // three body decay
456  TLorentzVector pDaughter0 = getDaughter(0)->get4Vector();
457  pDaughter0.Boost(boost);
458  TLorentzVector pDaughter1 = getDaughter(1)->get4Vector();
459  pDaughter1.Boost(boost);
460  pDaughter.SetVect(pDaughter0.Vect().Cross(pDaughter1.Vect()));
461  }
462  }
463 
464  double mag2 = pMother.Vect().Mag2() * pDaughter.Vect().Mag2();
465  if (mag2 <= 0) return std::numeric_limits<float>::quiet_NaN();
466  return (-pMother.Vect()) * pDaughter.Vect() / sqrt(mag2);
467 }
468 
469 float Particle::getCosHelicityDaughter(unsigned iDaughter, unsigned iGrandDaughter) const
470 {
471  // check existence of daughter
472  if (getNDaughters() <= iDaughter) {
473  B2ERROR("No daughter of particle 'name' with index 'iDaughter' for calculation of helicity angle"
474  << LogVar("name", getName()) << LogVar("iDaughter", iDaughter));
475  return std::numeric_limits<float>::quiet_NaN();
476  }
477 
478  // boost vector to the rest frame of the daughter particle
479  const Particle* daughter = getDaughter(iDaughter);
480  TVector3 boost = -daughter->get4Vector().BoostVector();
481 
482  // momentum of the this particle in the daughter's rest frame
483  TLorentzVector pMother = get4Vector();
484  pMother.Boost(boost);
485 
486  // check existence of grand daughter
487  if (daughter->getNDaughters() <= iGrandDaughter) {
488  B2ERROR("No grand daughter of daugher 'iDaughter' of particle 'name' with index 'iGrandDaughter' for calculation of helicity angle"
489  << LogVar("name", getName()) << LogVar("iDaughter", iDaughter) << LogVar("iGrandDaughter", iGrandDaughter));
490  return std::numeric_limits<float>::quiet_NaN();
491  }
492 
493  // momentum of the grand daughter in the daughter's rest frame
494  TLorentzVector pGrandDaughter = daughter->getDaughter(iGrandDaughter)->get4Vector();
495  pGrandDaughter.Boost(boost);
496 
497  double mag2 = pMother.Vect().Mag2() * pGrandDaughter.Vect().Mag2();
498  if (mag2 <= 0) return std::numeric_limits<float>::quiet_NaN();
499  return (-pMother.Vect()) * pGrandDaughter.Vect() / sqrt(mag2);
500 }
501 
503 {
504  // check that we have a decay to two daughters and then two grand daughters each
505  if (getNDaughters() != 2) {
506  B2ERROR("Cannot calculate acoplanarity of particle 'name' because the number of daughters is not 2"
507  << LogVar("name", getName()) << LogVar("# of daughters", getNDaughters()));
508  return std::numeric_limits<float>::quiet_NaN();
509  }
510  const Particle* daughter0 = getDaughter(0);
511  const Particle* daughter1 = getDaughter(1);
512  if ((daughter0->getNDaughters() != 2) || (daughter1->getNDaughters() != 2)) {
513  B2ERROR("Cannot calculate acoplanarity of particle 'name' because the number of grand daughters is not 2"
514  << LogVar("name", getName()) << LogVar("# of grand daughters of first daughter", daughter0->getNDaughters())
515  << LogVar("# of grand daughters of second daughter", daughter1->getNDaughters()));
516  return std::numeric_limits<float>::quiet_NaN();
517  }
518 
519  // boost vector to the rest frame of the particle
520  TVector3 boost = -get4Vector().BoostVector();
521 
522  // momenta of the daughters and grand daughters in the particle's rest frame
523  TLorentzVector pDaughter0 = daughter0->get4Vector();
524  pDaughter0.Boost(boost);
525  TLorentzVector pGrandDaughter0 = daughter0->getDaughter(0)->get4Vector();
526  pGrandDaughter0.Boost(boost);
527  TLorentzVector pDaughter1 = daughter1->get4Vector();
528  pDaughter1.Boost(boost);
529  TLorentzVector pGrandDaughter1 = daughter1->getDaughter(0)->get4Vector();
530  pGrandDaughter1.Boost(boost);
531 
532  // calculate angle between normal vectors
533  TVector3 normal0 = pDaughter0.Vect().Cross(pGrandDaughter0.Vect());
534  TVector3 normal1 = -pDaughter1.Vect().Cross(pGrandDaughter1.Vect());
535  double result = normal0.Angle(normal1);
536  if (normal0.Cross(normal1) * pDaughter0.Vect() < 0) result = -result;
537 
538  return result;
539 }
540 
541 
542 /*
543 float Particle::getMassError(void) const
544 {
545  float result = 0.0;
546 
547  if(m_pValue<0)
548  return result;
549 
550  float invMass = getMass();
551 
552  TMatrixFSym covarianceMatrix = getMomentumErrorMatrix();
553  TVectorF jacobian(c_DimMomentum);
554  jacobian[0] = -1.0*getPx()/invMass;
555  jacobian[1] = -1.0*getPy()/invMass;
556  jacobian[2] = -1.0*getPz()/invMass;
557  jacobian[3] = 1.0*getEnergy()/invMass;
558 
559  result = jacobian * (covarianceMatrix * jacobian);
560 
561  covarianceMatrix.Print();
562 
563  if(result<0.0)
564  result = 0.0;
565 
566  return TMath::Sqrt(result);
567 }
568 */
569 
570 void Particle::updateMass(const int pdgCode)
571 {
572  if (TDatabasePDG::Instance()->GetParticle(pdgCode) == nullptr)
573  B2FATAL("PDG=" << pdgCode << " ***code unknown to TDatabasePDG");
574  m_mass = TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass() ;
575 }
576 
577 float Particle::getPDGMass() const
578 {
579  if (TDatabasePDG::Instance()->GetParticle(m_pdgCode) == nullptr) {
580  B2ERROR("PDG=" << m_pdgCode << " ***code unknown to TDatabasePDG");
581  return 0.0;
582  }
583  return TDatabasePDG::Instance()->GetParticle(m_pdgCode)->Mass();
584 }
585 
586 float Particle::getCharge() const
587 {
588  if (TDatabasePDG::Instance()->GetParticle(m_pdgCode) == nullptr) {
589  B2ERROR("PDG=" << m_pdgCode << " ***code unknown to TDatabasePDG");
590  return 0.0;
591  }
592  return TDatabasePDG::Instance()->GetParticle(m_pdgCode)->Charge() / 3.0;
593 }
594 
595 const Particle* Particle::getDaughter(unsigned i) const
596 {
597  if (i >= getNDaughters()) return nullptr;
598  return static_cast<Particle*>(getArrayPointer()->At(m_daughterIndices[i]));
599 }
600 
601 std::vector<Belle2::Particle*> Particle::getDaughters() const
602 {
603  const unsigned int nDaughters = getNDaughters();
604  std::vector<Particle*> daughters(nDaughters);
605 
606  const TClonesArray* array = getArrayPointer();
607  for (unsigned i = 0; i < nDaughters; i++)
608  daughters[i] = static_cast<Particle*>(array->At(m_daughterIndices[i]));
609 
610  return daughters;
611 }
612 
613 std::vector<const Belle2::Particle*> Particle::getFinalStateDaughters() const
614 {
615  std::vector<const Particle*> fspDaughters;
616  fillFSPDaughters(fspDaughters);
617 
618  return fspDaughters;
619 }
620 
622 {
623  std::vector<int> mdstIndices;
624  for (const Particle* fsp : getFinalStateDaughters()) {
625  // is this FSP daughter constructed from given MDST source
626  if (fsp->getParticleSource() == source)
627  mdstIndices.push_back(fsp->getMdstArrayIndex());
628  }
629  return mdstIndices;
630 }
631 
632 
633 void Particle::appendDaughter(const Particle* daughter, const bool updateType)
634 {
635  if (updateType) {
636  // is it a composite particle or fsr corrected?
637  m_particleSource = c_Composite;
638  }
639 
640  // add daughter index
641  m_daughterIndices.push_back(daughter->getArrayIndex());
642  m_daughterProperties.push_back(Particle::PropertyFlags::c_Ordinary);
643 }
644 
645 void Particle::removeDaughter(const Particle* daughter, const bool updateType)
646 {
647  if (getNDaughters() == 0)
648  return;
649 
650  for (unsigned i = 0; i < getNDaughters(); i++) {
651  if (m_daughterIndices[i] == daughter->getArrayIndex()) {
652  m_daughterIndices.erase(m_daughterIndices.begin() + i);
653  m_daughterProperties.erase(m_daughterProperties.begin() + i);
654  i--;
655  }
656  }
657 
658  if (getNDaughters() == 0 and updateType)
659  m_particleSource = c_Undefined;
660 }
661 
662 bool Particle::overlapsWith(const Particle* oParticle) const
663 {
664  // obtain vectors of daughter final state particles
665  std::vector<const Particle*> thisFSPs = this->getFinalStateDaughters();
666  std::vector<const Particle*> otherFSPs = oParticle->getFinalStateDaughters();
667 
668  // check if they share any of the FSPs
669  for (auto& thisFSP : thisFSPs)
670  for (auto& otherFSP : otherFSPs)
671  if (thisFSP->getMdstSource() == otherFSP->getMdstSource())
672  return true;
673 
674  return false;
675 }
676 
677 bool Particle::isCopyOf(const Particle* oParticle, bool doDetailedComparison) const
678 {
679  // the name of the game is to as quickly as possible determine
680  // that the Particles are not copies
681  if (this->getPDGCode() != oParticle->getPDGCode() and !doDetailedComparison)
682  return false;
683 
684  unsigned nDaughters = this->getNDaughters();
685  if (nDaughters != oParticle->getNDaughters())
686  return false;
687 
688  if (nDaughters) {
689  // has daughters: check if the decay chain is the same and it ends with
690  // the same FSPs
691  std::vector<int> thisDecayChain(nDaughters * 2);
692  std::vector<int> othrDecayChain(nDaughters * 2);
693 
694  for (unsigned i = 0; i < nDaughters; i++) {
695  this->getDaughter(i)->fillDecayChain(thisDecayChain);
696  oParticle->getDaughter(i)->fillDecayChain(othrDecayChain);
697  }
698 
699  // Even if the daughters have been provided in a different order in the
700  // decay string the particles should be identified as copies / duplicates.
701  // Therefore, the decay chain vectors, which contain both the pdg codes and
702  // the mdst sources, are sorted here before comparing them.
703  sort(thisDecayChain.begin(), thisDecayChain.end());
704  sort(othrDecayChain.begin(), othrDecayChain.end());
705 
706  for (unsigned i = 0; i < thisDecayChain.size(); i++)
707  if (thisDecayChain[i] != othrDecayChain[i])
708  return false;
709 
710  } else if (this->getMdstSource() != oParticle->getMdstSource() and !doDetailedComparison) {
711  // has no daughters: it's a FSP, compare MDST source and index
712  return false;
713  }
714  // Stop here if we do not want a detailed comparison
715  if (!doDetailedComparison)
716  return true;
717  //If we compare here a reconstructed Particle to a generated MCParticle
718  //it means that something went horribly wrong and we must stop
719  if ((this->getParticleSource() == EParticleSourceObject::c_MCParticle
720  and oParticle->getParticleSource() != EParticleSourceObject::c_MCParticle)
721  or (this->getParticleSource() != EParticleSourceObject::c_MCParticle
722  and oParticle->getParticleSource() == EParticleSourceObject::c_MCParticle)) {
723  B2WARNING("Something went wrong: MCParticle is being compared to a non MC Particle. Please check your script!\n"
724  " If the MCParticle <-> Particle comparison happens in the RestOfEventBuilder,\n"
725  " the Rest Of Event may contain signal side particles.");
726  return false;
727  }
728  if (this->getParticleSource() == EParticleSourceObject::c_MCParticle
729  and oParticle->getParticleSource() == EParticleSourceObject::c_MCParticle) {
730  return this->getMCParticle() == oParticle->getMCParticle();
731  }
732  if (this->getParticleSource() != oParticle->getParticleSource()) {
733  return false;
734  }
735  if (this->getMdstSource() == oParticle->getMdstSource()) {
736  return true;
737  }
738  if (this->getTrack() && oParticle->getTrack() &&
739  this->getTrack()->getArrayIndex() != oParticle->getTrack()->getArrayIndex()) {
740  return false;
741  }
742  if (this->getKLMCluster() && oParticle->getKLMCluster()
743  && this->getKLMCluster()->getArrayIndex() != oParticle->getKLMCluster()->getArrayIndex()) {
744  return false;
745  }
746 
747  // It can be a bit more complicated for ECLClusters as we might also have to ensure they are connected-region unique
748  if (this->getECLCluster() && oParticle->getECLCluster()
749  && this->getECLCluster()->getArrayIndex() != oParticle->getECLCluster()->getArrayIndex()) {
750 
751  // if either is a track then they must be different
752  if (this->getECLCluster()->isTrack() or oParticle->getECLCluster()->isTrack())
753  return false;
754 
755  // we cannot combine two particles of different hypotheses from the same
756  // connected region (as their energies overlap)
757  if (this->getECLClusterEHypothesisBit() == oParticle->getECLClusterEHypothesisBit())
758  return false;
759 
760  // in the rare case that both are neutral and the hypotheses are different,
761  // we must also check that they are from different connected regions
762  // otherwise they come from the "same" underlying ECLShower
763  if (this->getECLCluster()->getConnectedRegionId() != oParticle->getECLCluster()->getConnectedRegionId())
764  return false;
765  }
766  return true;
767 
768 }
769 
770 const Track* Particle::getTrack() const
771 {
772  if (m_particleSource == c_Track) {
773  StoreArray<Track> tracks;
774  return tracks[m_mdstIndex];
775  } else
776  return nullptr;
777 }
778 
780 {
781  // if the particle is related to a TrackFitResult then return this
782  auto* selfrelated = this->getRelatedTo<TrackFitResult>();
783  if (selfrelated)
784  return selfrelated;
785 
786  // if not get the TFR with closest mass to this particle
787  auto* selftrack = this->getTrack();
788  if (selftrack)
789  return selftrack->getTrackFitResultWithClosestMass(
790  Belle2::Const::ChargedStable(std::abs(this->getPDGCode())));
791 
792  // otherwise we're probably not a track based particle
793  return nullptr;
794 }
795 
797 {
798  if (m_particleSource == c_Track) {
799  StoreArray<Track> tracks;
800  return tracks[m_mdstIndex]->getRelated<PIDLikelihood>();
801  } else
802  return nullptr;
803 }
804 
805 const V0* Particle::getV0() const
806 {
807  if (m_particleSource == c_V0) {
808  StoreArray<V0> v0s;
809  return v0s[m_mdstIndex];
810  } else {
811  return nullptr;
812  }
813 }
814 
815 
817 {
818  if (m_particleSource == c_ECLCluster) {
819  StoreArray<ECLCluster> eclClusters;
820  return eclClusters[m_mdstIndex];
821  } else if (m_particleSource == c_Track) {
822  // a track may be matched to several clusters under different hypotheses
823  // take the most energetic of the c_nPhotons hypothesis as "the" cluster
824  StoreArray<Track> tracks;
825  const ECLCluster* bestTrackMatchedCluster = nullptr;
826  double highestEnergy = -1.0;
827  // loop over all clusters matched to this track
828  for (const ECLCluster& cluster : tracks[m_mdstIndex]->getRelationsTo<ECLCluster>()) {
829  // ignore everything except the nPhotons hypothesis
830  if (!cluster.hasHypothesis(ECLCluster::EHypothesisBit::c_nPhotons))
831  continue;
832  // check if we're most energetic thus far
833  if (cluster.getEnergy(ECLCluster::EHypothesisBit::c_nPhotons) > highestEnergy) {
834  highestEnergy = cluster.getEnergy(ECLCluster::EHypothesisBit::c_nPhotons);
835  bestTrackMatchedCluster = &cluster;
836  }
837  }
838  return bestTrackMatchedCluster;
839  } else {
840  return nullptr;
841  }
842 }
843 
845 {
846  const ECLCluster* cluster = this->getECLCluster();
847  return cluster->getEnergy(this->getECLClusterEHypothesisBit());
848 }
849 
851 {
852  if (m_particleSource == c_KLMCluster) {
853  StoreArray<KLMCluster> klmClusters;
854  return klmClusters[m_mdstIndex];
855  } else if (m_particleSource == c_Track) {
856  // If there is an associated KLMCluster, it's the closest one
857  StoreArray<Track> tracks;
858  const KLMCluster* klmCluster = tracks[m_mdstIndex]->getRelatedTo<KLMCluster>();
859  return klmCluster;
860  } else {
861  return nullptr;
862  }
863 }
864 
865 
867 {
868  if (m_particleSource == c_MCParticle) {
869  StoreArray<MCParticle> mcParticles;
870  return mcParticles[m_mdstIndex];
871  } else {
872  const MCParticle* related = this->getRelated<MCParticle>();
873  if (related)
874  return related;
875  }
876  return nullptr;
877 }
878 
879 
880 const Particle* Particle::getParticleFromGeneralizedIndexString(const std::string& generalizedIndex) const
881 {
882  // Split the generalizedIndex string in a vector of strings.
883  std::vector<std::string> generalizedIndexes;
884  boost::split(generalizedIndexes, generalizedIndex, boost::is_any_of(":"));
885 
886  if (generalizedIndexes.empty()) {
887  B2WARNING("Generalized index of daughter particle is empty. Skipping.");
888  return nullptr;
889  }
890 
891  // To explore a decay tree of unknown depth, we need a place to store
892  // both the root particle and the daughter particle at each iteration
893  const Particle* dauPart =
894  nullptr; // This will be eventually returned
895  const Particle* currentPart = this; // This is the root particle of the next iteration
896 
897  // Loop over the generalizedIndexes until you get to the particle you want
898  for (auto& indexString : generalizedIndexes) {
899  // indexString is a string. First try to convert it into an int
900  int dauIndex = 0;
901  try {
902  dauIndex = Belle2::convertString<int>(indexString);
903  } catch (std::invalid_argument&) {
904  B2WARNING("Found the string " << indexString << "instead of a daughter index.");
905  return nullptr;
906  }
907 
908  // Check that the daughter index is smaller than the number of daughters of the current root particle
909  if (dauIndex >= int(currentPart->getNDaughters()) or dauIndex < 0) {
910  B2WARNING("Daughter index " << dauIndex << " out of range");
911  B2WARNING("Trying to access non-existing particle.");
912  return nullptr;
913  } else {
914  dauPart = currentPart->getDaughter(dauIndex); // Pick the particle indicated by the generalizedIndex
915  currentPart = dauPart;
916  }
917  }
918  return dauPart;
919 }
920 
921 
922 
923 //--- private methods --------------------------------------------
924 
926 {
927  // set momenum
928  m_px = trackFit->getMomentum().Px();
929  m_py = trackFit->getMomentum().Py();
930  m_pz = trackFit->getMomentum().Pz();
931 
932  // set position at which the momentum is given (= POCA)
933  setVertex(trackFit->getPosition());
934 
935  // set Chi^2 probability
936  m_pValue = trackFit->getPValue();
937 
938  // set error matrix
939  const auto cov6 = trackFit->getCovariance6();
940  constexpr unsigned order[] = {c_X, c_Y, c_Z, c_Px, c_Py, c_Pz};
941 
942  TMatrixFSym errMatrix(c_DimMatrix);
943  for (int i = 0; i < 6; i++) {
944  for (int j = i; j < 6; j++) {
945  // although it seems to make no sense to fill all elements of the
946  // symetric matrix, it has to be (do not touch this code)
947  errMatrix(order[j], order[i]) = errMatrix(order[i], order[j]) = cov6(i, j);
948  }
949  }
950 
951  /*
952  E = sqrt(px^2 + py^2 + pz^2 + m^2) thus:
953  cov(x,E) = cov(px,x) *dE/dpx + cov(py,x) *dE/dpy + cov(pz,x) *dE/dpz
954  cov(y,E) = cov(px,y) *dE/dpx + cov(py,y) *dE/dpy + cov(pz,y) *dE/dpz
955  cov(z,E) = cov(px,z) *dE/dpx + cov(py,z) *dE/dpy + cov(pz,z) *dE/dpz
956  cov(px,E) = cov(px,px)*dE/dpx + cov(px,py)*dE/dpy + cov(px,pz)*dE/dpz
957  cov(py,E) = cov(py,px)*dE/dpx + cov(py,py)*dE/dpy + cov(py,pz)*dE/dpz
958  cov(pz,E) = cov(pz,px)*dE/dpx + cov(pz,py)*dE/dpy + cov(pz,pz)*dE/dpz
959  cov(E,E) = cov(px,px)*(dE/dpx)^2 + cov(py,py)*(dE/dpy)^2 + cov(pz,pz)*(dE/dpz)^2
960  + 2*cov(px,py)*dE/dpx*dE/dpy
961  + 2*cov(py,pz)*dE/dpy*dE/dpz
962  + 2*cov(pz,px)*dE/dpz*dE/dpx
963  dE/dpx = px/E etc.
964  */
965 
966  const float E = getEnergy();
967  const float dEdp[] = {m_px / E, m_py / E, m_pz / E};
968  constexpr unsigned compMom[] = {c_Px, c_Py, c_Pz};
969  constexpr unsigned compPos[] = {c_X, c_Y, c_Z};
970 
971  // covariances (p,E)
972  for (unsigned int i : compMom) {
973  float Cov = 0;
974  for (int k = 0; k < 3; k++) {
975  Cov += errMatrix(i, compMom[k]) * dEdp[k];
976  }
977  errMatrix(i, c_E) = Cov;
978  }
979 
980  // covariances (x,E)
981  for (unsigned int comp : compPos) {
982  float Cov = 0;
983  for (int k = 0; k < 3; k++) {
984  Cov += errMatrix(comp, compMom[k]) * dEdp[k];
985  }
986  errMatrix(c_E, comp) = Cov;
987  }
988 
989  // variance (E,E)
990  float Cov = 0;
991  for (int i = 0; i < 3; i++) {
992  Cov += errMatrix(compMom[i], compMom[i]) * dEdp[i] * dEdp[i];
993  }
994  for (int i = 0; i < 3; i++) {
995  int k = (i + 1) % 3;
996  Cov += 2 * errMatrix(compMom[i], compMom[k]) * dEdp[i] * dEdp[k];
997  }
998  errMatrix(c_E, c_E) = Cov;
999 
1000  storeErrorMatrix(errMatrix);
1001 }
1002 
1004 {
1005  for (float& i : m_errMatrix)
1006  i = 0.0;
1007 }
1008 
1009 void Particle::storeErrorMatrix(const TMatrixFSym& m)
1010 {
1011  int element = 0;
1012  for (int irow = 0; irow < c_DimMatrix; irow++) {
1013  for (int icol = irow; icol < c_DimMatrix; icol++) {
1014  m_errMatrix[element] = m(irow, icol);
1015  element++;
1016  }
1017  }
1018 }
1019 
1020 
1021 void Particle::fillFSPDaughters(std::vector<const Belle2::Particle*>& fspDaughters) const
1022 {
1023  // this is FSP
1024  if (getNDaughters() == 0) {
1025  fspDaughters.push_back(this);
1026  return;
1027  }
1028 
1029  // this is not FSP (go one level down)
1030  for (unsigned i = 0; i < getNDaughters(); i++)
1031  getDaughter(i)->fillFSPDaughters(fspDaughters);
1032 }
1033 
1034 void Particle::fillDecayChain(std::vector<int>& decayChain) const
1035 {
1036  decayChain.push_back(m_pdgCode);
1037  decayChain.push_back(getMdstSource());
1038 
1039  for (unsigned i = 0; i < getNDaughters(); i++)
1040  getDaughter(i)->fillDecayChain(decayChain);
1041 }
1042 
1043 
1045 {
1047  if (m_pdgCode < 0) return;
1048  if (m_pdgCode == 22) {m_flavorType = c_Unflavored; return;} // gamma
1049  if (m_pdgCode == 310) {m_flavorType = c_Unflavored; return;} // K_s
1050  if (m_pdgCode == 130) {m_flavorType = c_Unflavored; return;} // K_L
1051  if (m_pdgCode == 43) {m_flavorType = c_Unflavored; return;} // Xu0
1052  int nnn = m_pdgCode / 10;
1053  int q3 = nnn % 10; nnn /= 10;
1054  int q2 = nnn % 10; nnn /= 10;
1055  int q1 = nnn % 10;
1056  if (q1 == 0 && q2 == q3) m_flavorType = c_Unflavored; // unflavored meson
1057 }
1058 
1059 std::string Particle::getName() const
1060 {
1061  return TDatabasePDG::Instance()->GetParticle(m_pdgCode)->GetName();
1062 }
1063 
1064 void Particle::print() const
1065 {
1066  B2INFO(getInfo());
1067 }
1068 
1069 std::vector<std::string> Particle::getExtraInfoNames() const
1070 {
1071  std::vector<std::string> out;
1072  if (!m_extraInfo.empty()) {
1073  StoreObjPtr<ParticleExtraInfoMap> extraInfoMap;
1074  if (!extraInfoMap)
1075  B2FATAL("ParticleExtraInfoMap not available, but needed for storing extra info in Particle!");
1076  const ParticleExtraInfoMap::IndexMap& map = extraInfoMap->getMap(m_extraInfo[0]);
1077  for (auto const& ee : map) out.push_back(ee.first);
1078  }
1079  return out;
1080 }
1081 
1082 std::string Particle::getInfoHTML() const
1083 {
1084  std::stringstream stream;
1085  stream << std::setprecision(4);
1086  stream << "<b>collection</b>=" << getArrayName();
1087  stream << "<br>";
1088  stream << " <b>PDGCode</b>=" << m_pdgCode;
1089  stream << " <b>Charge</b>=" << getCharge();
1090  stream << " <b>PDGMass</b>=" << getPDGMass();
1091  stream << "<br>";
1092  stream << " <b>flavorType</b>=" << m_flavorType;
1093  stream << " <b>particleType</b>=" << m_particleSource;
1094  stream << " <b>particleTypeUsedForFit</b>=" << m_pdgCodeUsedForFit;
1095  stream << "<br>";
1096 
1097  stream << " <b>mdstIndex</b>=" << m_mdstIndex;
1098  stream << " <b>arrayIndex</b>=" << getArrayIndex();
1099  stream << " <b>identifier</b>=" << m_identifier;
1100  stream << " <b>daughterIndices</b>: ";
1101  for (int daughterIndex : m_daughterIndices) {
1102  stream << daughterIndex << ", ";
1103  }
1104  if (m_daughterIndices.empty()) stream << " (none)";
1105  stream << "<br>";
1106 
1107  if (!m_daughterIndices.empty()) {
1108  stream << " <b>daughter PDGCodes</b>: ";
1109  for (unsigned i = 0; i < m_daughterIndices.size(); i++) {
1110  const Particle* p = getDaughter(i);
1111  if (p) {stream << p->getPDGCode() << ", ";}
1112  else {stream << "?, ";}
1113  }
1114  stream << "<br>";
1115  }
1116 
1117  stream << " <b>mass</b>=" << m_mass;
1118  stream << "<br>";
1119 
1120  stream << " <b>momentum</b>=" << HTML::getString(getMomentum());
1121  stream << " <b>p</b>=" << getP();
1122  stream << "<br>";
1123 
1124  stream << " <b>position</b>=" << HTML::getString(getVertex());
1125  stream << "<br>";
1126 
1127  stream << " <b>p-value of fit</b> (if done): ";
1128  stream << m_pValue;
1129  stream << "<br>";
1130 
1131  stream << " <b>error matrix</b> (px, py, pz, E, x, y ,z):<br>";
1133 
1134  stream << " <b>extra info</b>=( ";
1135  if (!m_extraInfo.empty()) {
1136  StoreObjPtr<ParticleExtraInfoMap> extraInfoMap;
1137  if (!extraInfoMap) {
1138  B2FATAL("ParticleExtraInfoMap not available, but needed for storing extra info in Particle!");
1139  }
1140  const ParticleExtraInfoMap::IndexMap& map = extraInfoMap->getMap(m_extraInfo[0]);
1141  const unsigned int nVars = m_extraInfo.size();
1142  for (const auto& pair : map) {
1143  if (pair.second < nVars) {
1144  stream << pair.first << "=" << m_extraInfo[pair.second] << " ";
1145  }
1146  }
1147 
1148  }
1149  stream << ") " << "<br>";
1150 
1151  return stream.str();
1152 }
1153 
1154 bool Particle::hasExtraInfo(const std::string& name) const
1155 {
1156  if (m_extraInfo.empty())
1157  return false;
1158 
1159  //get index for name
1160  const auto mapID = (unsigned int)m_extraInfo[0];
1161  StoreObjPtr<ParticleExtraInfoMap> extraInfoMap;
1162  if (!extraInfoMap) {
1163  B2FATAL("ParticleExtraInfoMap not available, but needed for storing extra info in Particle!");
1164  }
1165  unsigned int index = extraInfoMap->getIndex(mapID, name);
1166  if (index == 0 or index >= m_extraInfo.size()) //actualy indices start at 1
1167  return false;
1168 
1169  return true;
1170 }
1171 
1173 {
1174  m_extraInfo.clear();
1175 }
1176 
1177 float Particle::getExtraInfo(const std::string& name) const
1178 {
1179  if (m_extraInfo.empty())
1180  throw std::runtime_error(std::string("getExtraInfo: Value '") + name + "' not found in Particle!");
1181 
1182  //get index for name
1183  const auto mapID = (unsigned int)m_extraInfo[0];
1184  StoreObjPtr<ParticleExtraInfoMap> extraInfoMap;
1185  if (!extraInfoMap) {
1186  B2FATAL("ParticleExtraInfoMap not available, but needed for storing extra info in Particle!");
1187  }
1188  unsigned int index = extraInfoMap->getIndex(mapID, name);
1189  if (index == 0 or index >= m_extraInfo.size()) //actualy indices start at 1
1190  throw std::runtime_error(std::string("getExtraInfo: Value '") + name + "' not found in Particle!");
1191 
1192  return m_extraInfo[index];
1193 
1194 }
1195 
1196 void Particle::writeExtraInfo(const std::string& name, const float value)
1197 {
1198  if (this->hasExtraInfo(name)) {
1199  this->setExtraInfo(name, value);
1200  } else {
1201  this->addExtraInfo(name, value);
1202  }
1203 }
1204 
1205 void Particle::setExtraInfo(const std::string& name, float value)
1206 {
1207  if (m_extraInfo.empty())
1208  throw std::runtime_error(std::string("setExtraInfo: Value '") + name + "' not found in Particle!");
1209 
1210  //get index for name
1211  const auto mapID = (unsigned int)m_extraInfo[0];
1212  StoreObjPtr<ParticleExtraInfoMap> extraInfoMap;
1213  if (!extraInfoMap) {
1214  B2FATAL("ParticleExtraInfoMap not available, but needed for storing extra info in Particle!");
1215  }
1216  unsigned int index = extraInfoMap->getIndex(mapID, name);
1217  if (index == 0 or index >= m_extraInfo.size()) //actualy indices start at 1
1218  throw std::runtime_error(std::string("setExtraInfo: Value '") + name + "' not found in Particle!");
1219 
1220  m_extraInfo[index] = value;
1221 
1222 }
1223 
1224 void Particle::addExtraInfo(const std::string& name, float value)
1225 {
1226  if (hasExtraInfo(name))
1227  throw std::runtime_error(std::string("addExtraInfo: Value '") + name + "' already set!");
1228 
1229  StoreObjPtr<ParticleExtraInfoMap> extraInfoMap;
1230  if (!extraInfoMap)
1231  extraInfoMap.create();
1232  if (m_extraInfo.empty()) {
1233  unsigned int mapID = extraInfoMap->getMapForNewVar(name);
1234  m_extraInfo.push_back(mapID);
1235  m_extraInfo.push_back(value);
1236  } else {
1237  unsigned int oldMapID = m_extraInfo[0];
1238  unsigned int insertIndex = m_extraInfo.size();
1239  unsigned int mapID = extraInfoMap->getMapForNewVar(name, oldMapID, insertIndex);
1240 
1241  m_extraInfo[0] = mapID; //update map
1242  m_extraInfo.push_back(value); //add value
1243  }
1244 }
1245 
1246 bool Particle::forEachDaughter(const std::function<bool(const Particle*)>& function,
1247  bool recursive, bool includeSelf) const
1248 {
1249  std::queue<const Particle*> qq;
1250  // If we include ourselves add only this, otherwise directly all children
1251  if (includeSelf) {
1252  qq.push(this);
1253  } else {
1254  for (size_t i = 0; i < getNDaughters(); ++i) qq.push(getDaughter(i));
1255  }
1256  // Now just repeat until done: take the child, run the functor, remove the
1257  // child, add all children if needed
1258  while (!qq.empty()) {
1259  const Particle* p = qq.front();
1260  if (function(p)) return true;
1261  qq.pop();
1262  // Add children if we go through all children recursively or if we look at
1263  // the current particle: we always want the direct children.
1264  if (recursive || p == this)
1265  for (size_t i = 0; i < p->getNDaughters(); ++i) qq.push(p->getDaughter(i));
1266  }
1267  return false;
1268 }
1269 
1270 int Particle::generatePDGCodeFromCharge(const int chargeSign, const Const::ChargedStable& chargedStable)
1271 {
1272  int absPDGCode = chargedStable.getPDGCode();
1273  int PDGCode = absPDGCode * chargeSign;
1274  // flip sign of PDG code for leptons: their PDG code is positive if the lepton charge is negative and vice versa
1275  if (chargedStable == Const::muon || chargedStable == Const::electron) PDGCode = -PDGCode;
1276  return PDGCode;
1277 }
Belle2::Particle::getPDGCode
int getPDGCode(void) const
Returns PDG code.
Definition: Particle.h:380
Belle2::ParticleExtraInfoMap::IndexMap
std::map< std::string, unsigned int > IndexMap
string -> index map.
Definition: ParticleExtraInfoMap.h:42
Belle2::Const::photon
static const ParticleType photon
photon particle
Definition: Const.h:547
Belle2::Particle::getCosHelicityDaughter
float getCosHelicityDaughter(unsigned iDaughter, unsigned iGrandDaughter=0) const
Returns cosine of the helicity angle of the given daughter defined by given grand daughter.
Definition: Particle.cc:469
Belle2::Particle::m_daughterProperties
std::vector< int > m_daughterProperties
daughter particle properties
Definition: Particle.h:922
Belle2::Particle::getArrayPointer
TClonesArray * getArrayPointer() const
Returns the pointer to the store array which holds the daughter particles.
Definition: Particle.h:844
Belle2::Particle::writeExtraInfo
void writeExtraInfo(const std::string &name, const float value)
Sets the user defined extraInfo.
Definition: Particle.cc:1196
Belle2::Particle::set4Vector
void set4Vector(const TLorentzVector &p4)
Sets Lorentz vector.
Definition: Particle.h:279
Belle2::Particle::getTrackFitResult
const TrackFitResult * getTrackFitResult() const
Returns the pointer to the TrackFitResult that was used to create this Particle (ParticleType == c_Tr...
Definition: Particle.cc:779
Belle2::Particle::getVertex
TVector3 getVertex() const
Returns vertex position (POCA for charged, IP for neutral FS particles)
Definition: Particle.h:529
Belle2::RelationsInterface::getInfo
std::string getInfo() const
Return a short summary of this object's contents in raw text format.
Definition: RelationsObject.h:372
Belle2::V0
Object holding information for V0s.
Definition: V0.h:40
Belle2::TrackFitResult::getMomentum
TVector3 getMomentum() const
Getter for vector of momentum at closest approach of track in r/phi projection.
Definition: TrackFitResult.h:116
Belle2::ECLCluster
ECL cluster data.
Definition: ECLCluster.h:39
Belle2::Particle::generatePDGCodeFromCharge
int generatePDGCodeFromCharge(const int chargedSign, const Const::ChargedStable &chargedStable)
Generate the PDG code with correct sign, using the charge.
Definition: Particle.cc:1270
Belle2::Const::electron
static const ChargedStable electron
electron particle
Definition: Const.h:533
Belle2::Particle::print
void print() const
Prints the contents of a Particle object to standard output.
Definition: Particle.cc:1064
Belle2::Particle::m_px
float m_px
momentum component x
Definition: Particle.h:909
Belle2::Particle::m_daughterIndices
std::vector< int > m_daughterIndices
daughter particle indices
Definition: Particle.h:917
Belle2::Particle::m_pdgCode
int m_pdgCode
PDG code.
Definition: Particle.h:906
Belle2::Particle::addExtraInfo
void addExtraInfo(const std::string &name, float value)
Sets the user-defined data of given name to the given value.
Definition: Particle.cc:1224
Belle2::TrackFitResult::getPValue
double getPValue() const
Getter for Chi2 Probability of the track fit.
Definition: TrackFitResult.h:163
Belle2::Particle::resetErrorMatrix
void resetErrorMatrix()
Resets 7x7 error matrix All elements are set to 0.0.
Definition: Particle.cc:1003
Belle2::ClusterUtils::Get4MomentumFromCluster
const TLorentzVector Get4MomentumFromCluster(const ECLCluster *cluster, ECLCluster::EHypothesisBit hypo)
Returns four momentum vector.
Definition: ClusterUtils.cc:26
Belle2::Particle::getMdstArrayIndices
std::vector< int > getMdstArrayIndices(EParticleSourceObject type) const
Returns a vector of StoreArray indices of given MDST dataobjects.
Definition: Particle.cc:621
Belle2::Particle::m_flavorType
EFlavorType m_flavorType
flavor type.
Definition: Particle.h:918
Belle2::Particle::m_extraInfo
std::vector< float > m_extraInfo
Stores associated user defined values.
Definition: Particle.h:938
Belle2::RelationsInterface::getArrayName
std::string getArrayName() const
Get name of array this object is stored in, or "" if not found.
Definition: RelationsObject.h:379
Belle2::Particle::setExtraInfo
void setExtraInfo(const std::string &name, float value)
Sets the user-defined data of given name to the given value.
Definition: Particle.cc:1205
Belle2::ECLCluster::EHypothesisBit::c_nPhotons
@ c_nPhotons
CR is split into n photons (N1)
Belle2::Particle::m_mass
float m_mass
particle (invariant) mass
Definition: Particle.h:908
Belle2::Const::ParticleType::getPDGCode
int getPDGCode() const
PDG code.
Definition: Const.h:349
Belle2::Particle::m_pValue
float m_pValue
chi^2 probability of the fit.
Definition: Particle.h:916
Belle2::Particle::m_identifier
int m_identifier
Identifier that can be used to identify whether the particle is unique or is a copy or representation...
Definition: Particle.h:931
Belle2::PIDLikelihood
Class to collect log likelihoods from TOP, ARICH, dEdx, ECL and KLM aimed for output to mdst includes...
Definition: PIDLikelihood.h:37
Belle2::Particle::fillDecayChain
void fillDecayChain(std::vector< int > &decayChain) const
Fill vector with (PDGCode, MdstSource) pairs for the entire decay chain.
Definition: Particle.cc:1034
Belle2::Particle::storeErrorMatrix
void storeErrorMatrix(const TMatrixFSym &errMatrix)
Stores 7x7 error matrix into private member m_errMatrix.
Definition: Particle.cc:1009
Belle2::Particle::getParticleFromGeneralizedIndexString
const Particle * getParticleFromGeneralizedIndexString(const std::string &generalizedIndex) const
Explores the decay tree of the particle and returns the (grand^n)daughter identified by a generalized...
Definition: Particle.cc:880
Belle2::Particle::getDaughter
const Particle * getDaughter(unsigned i) const
Returns a pointer to the i-th daughter particle.
Definition: Particle.cc:595
Belle2::MCParticle::getArrayIndex
int getArrayIndex() const
Get 0-based index of the particle in the corresponding MCParticle list.
Definition: MCParticle.h:255
Belle2::TrackFitResult
Values of the result of a track fit with a given particle hypothesis.
Definition: TrackFitResult.h:59
Belle2::Particle::getP
float getP() const
Returns momentum magnitude (same as getMomentumMagnitude but with shorter name)
Definition: Particle.h:493
Belle2::TrackFitResult::getPosition
TVector3 getPosition() const
Getter for vector of position at closest approach of track in r/phi projection.
Definition: TrackFitResult.h:109
Belle2::ClusterUtils::GetIPPosition
const TVector3 GetIPPosition()
Returns default IP position from beam parameters.
Definition: ClusterUtils.cc:182
Belle2::Particle::getMdstSource
int getMdstSource() const
Returns unique identifier of final state particle (needed in particle combiner)
Definition: Particle.cc:347
Belle2::Particle::c_Flavored
@ c_Flavored
Is either particle or antiparticle.
Definition: Particle.h:97
Belle2::DBObjPtr
Class for accessing objects in the database.
Definition: DBObjPtr.h:31
Belle2::Particle::m_particleSource
EParticleSourceObject m_particleSource
(mdst) source of particle
Definition: Particle.h:919
Belle2::Particle::getKLMCluster
const KLMCluster * getKLMCluster() const
Returns the pointer to the KLMCluster object that was used to create this Particle (ParticleType == c...
Definition: Particle.cc:850
Belle2::Particle::m_properties
int m_properties
particle property
Definition: Particle.h:921
Belle2::Particle::m_pdgCodeUsedForFit
int m_pdgCodeUsedForFit
PDG code used for the track fit.
Definition: Particle.h:907
Belle2::ECLCluster::getEnergy
double getEnergy(const EHypothesisBit &hypothesis) const
Return Energy (GeV).
Definition: ECLCluster.cc:21
Belle2::Particle::getEnergy
float getEnergy() const
Returns total energy.
Definition: Particle.h:455
Belle2::ECLCluster::isTrack
bool isTrack() const
Return true if the cluster matches with track.
Definition: ECLCluster.h:247
Belle2::Particle::removeExtraInfo
void removeExtraInfo()
Remove all stored extra info fields.
Definition: Particle.cc:1172
Belle2::Particle::getECLClusterEHypothesisBit
ECLCluster::EHypothesisBit getECLClusterEHypothesisBit() const
Returns the ECLCluster EHypothesisBit for this Particle.
Definition: Particle.h:874
Belle2::Particle::getFinalStateDaughters
std::vector< const Belle2::Particle * > getFinalStateDaughters() const
Returns a vector of pointers to Final State daughter particles.
Definition: Particle.cc:613
Belle2::Particle::getCharge
float getCharge(void) const
Returns particle charge.
Definition: Particle.cc:586
Belle2::Particle::getPIDLikelihood
const PIDLikelihood * getPIDLikelihood() const
Returns the pointer to the PIDLikelihood object that is related to the Track, which was used to creat...
Definition: Particle.cc:796
Belle2::ClusterUtils
Class to provide momentum-related information from ECLClusters.
Definition: ClusterUtils.h:44
Belle2::TrackFitResult::getParticleType
Const::ParticleType getParticleType() const
Getter for ParticleType of the mass hypothesis of the track fit.
Definition: TrackFitResult.h:154
Belle2::ClusterUtils::GetCovarianceMatrix7x7FromCluster
const TMatrixDSym GetCovarianceMatrix7x7FromCluster(const ECLCluster *cluster, ECLCluster::EHypothesisBit hypo)
Returns 7x7 covariance matrix (px, py, pz, E, x, y, z)
Definition: ClusterUtils.cc:149
Belle2::Particle::getDaughters
std::vector< Belle2::Particle * > getDaughters() const
Returns a vector of pointers to daughter particles.
Definition: Particle.cc:601
Belle2::Particle::EFlavorType
EFlavorType
describes flavor type, see getFlavorType().
Definition: Particle.h:95
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::Particle::overlapsWith
bool overlapsWith(const Particle *oParticle) const
Returns true if final state ancessors of oParticle overlap.
Definition: Particle.cc:662
Belle2::StoreObjPtr
Type-safe access to single objects in the data store.
Definition: ParticleList.h:33
Belle2::HTML::getString
std::string getString(const TMatrixFBase &matrix, int precision=2, bool color=true)
get HTML table representing a matrix.
Definition: HTML.cc:18
Belle2::MCParticle::getVertex
TVector3 getVertex() const
Return production vertex position, shorthand for getProductionVertex().
Definition: MCParticle.h:194
Belle2::Particle::getECLCluster
const ECLCluster * getECLCluster() const
Returns the pointer to the ECLCluster object that was used to create this Particle (if ParticleType =...
Definition: Particle.cc:816
Belle2::Particle::getMomentumVertexErrorMatrix
TMatrixFSym getMomentumVertexErrorMatrix() const
Returns 7x7 error matrix.
Definition: Particle.cc:386
Belle2::KLMCluster::getPosition
TVector3 getPosition() const
Get global position (TVector3 version) of the origin of KLMCluster (always return (0,...
Definition: KLMCluster.h:99
Belle2::Particle::getNDaughters
unsigned getNDaughters(void) const
Returns number of daughter particles.
Definition: Particle.h:625
Belle2::Particle::appendDaughter
void appendDaughter(const Particle *daughter, const bool updateType=true)
Appends index of daughter to daughters index array.
Definition: Particle.cc:633
Belle2::Particle::setMomentumVertexErrorMatrix
void setMomentumVertexErrorMatrix(const TMatrixFSym &errMatrix)
Sets 7x7 error matrix.
Definition: Particle.cc:373
LogVar
Class to store variables with their name which were sent to the logging service.
Definition: LogVariableStream.h:24
Belle2::MCParticle::getPDG
int getPDG() const
Return PDG code of particle.
Definition: MCParticle.h:123
Belle2::Const::pi0
static const ParticleType pi0
neutral pion particle
Definition: Const.h:548
Belle2::Particle::c_Unflavored
@ c_Unflavored
Is its own antiparticle or we don't know whether it is a particle/antiparticle.
Definition: Particle.h:96
Belle2::Particle::getExtraInfo
float getExtraInfo(const std::string &name) const
Return given value if set.
Definition: Particle.cc:1177
Belle2::MCParticle::getMass
float getMass() const
Return the particle mass in GeV.
Definition: MCParticle.h:146
Belle2::Particle::setVertex
void setVertex(const TVector3 &vertex)
Sets position (decay vertex)
Definition: Particle.h:291
Belle2::KLMCluster::getMomentum
TLorentzVector getMomentum() const
Get momentum.
Definition: KLMCluster.cc:50
Belle2::Particle::m_errMatrix
float m_errMatrix[c_SizeMatrix]
error matrix (1D representation)
Definition: Particle.h:915
Belle2::RelationsInterface::getArrayIndex
int getArrayIndex() const
Returns this object's array index (in StoreArray), or -1 if not found.
Definition: RelationsObject.h:387
Belle2::Particle::EParticleSourceObject
EParticleSourceObject
particle source enumerators
Definition: Particle.h:84
Belle2::KLMCluster
KLM cluster data.
Definition: KLMCluster.h:38
Belle2::Particle::getMomentum
TVector3 getMomentum() const
Returns momentum vector.
Definition: Particle.h:475
Belle2::Particle::isCopyOf
bool isCopyOf(const Particle *oParticle, bool doDetailedComparison=false) const
Returns true if this Particle and oParticle are copies of each other.
Definition: Particle.cc:677
Belle2::Particle::getTrack
const Track * getTrack() const
Returns the pointer to the Track object that was used to create this Particle (ParticleType == c_Trac...
Definition: Particle.cc:770
Belle2::Particle
Class to store reconstructed particles.
Definition: Particle.h:77
Belle2::Particle::getVertexErrorMatrix
TMatrixFSym getVertexErrorMatrix() const
Returns the 3x3 position error sub-matrix.
Definition: Particle.cc:413
Belle2::Particle::updateMass
void updateMass(const int pdgCode)
Updates particle mass with the mass of the particle corresponding to the given PDG.
Definition: Particle.cc:570
Belle2::Const::ParticleType
The ParticleType class for identifying different particle types.
Definition: Const.h:284
Belle2::Particle::m_pz
float m_pz
momentum component z
Definition: Particle.h:911
Belle2::Particle::removeDaughter
void removeDaughter(const Particle *daughter, const bool updateType=true)
Removes index of daughter from daughters index array.
Definition: Particle.cc:645
Belle2::Const::ChargedStable
Provides a type-safe way to pass members of the chargedStableSet set.
Definition: Const.h:465
Belle2::Particle::getCosHelicity
float getCosHelicity(const Particle *mother=nullptr) const
Returns cosine of the helicity angle The helicity angle is defined in the rest frame of the particle ...
Definition: Particle.cc:425
Belle2::MCParticle::getMomentum
TVector3 getMomentum() const
Return momentum.
Definition: MCParticle.h:209
Belle2::Particle::getPDGMass
float getPDGMass(void) const
Returns uncertaint on the invariant mass (requiers valid momentum error matrix)
Definition: Particle.cc:577
Belle2::Particle::m_mdstIndex
unsigned m_mdstIndex
0-based index of MDST store array object
Definition: Particle.h:920
Belle2::Particle::setMdstArrayIndex
void setMdstArrayIndex(const int arrayIndex)
set mdst array index
Definition: Particle.cc:326
Belle2::Const::muon
static const ChargedStable muon
muon particle
Definition: Const.h:534
Belle2::Particle::setFlavorType
void setFlavorType()
sets m_flavorType using m_pdgCode
Definition: Particle.cc:1044
Belle2::Track
Class that bundles various TrackFitResults.
Definition: Track.h:35
Belle2::MCParticle
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:43
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition: ECLMatchingPerformanceExpertModule.h:33
Belle2::Particle::m_py
float m_py
momentum component y
Definition: Particle.h:910
Belle2::Particle::Particle
Particle()
Default constructor.
Definition: Particle.cc:47
Belle2::Particle::get4Vector
TLorentzVector get4Vector() const
Returns Lorentz vector.
Definition: Particle.h:464
Belle2::Particle::getParticleSource
EParticleSourceObject getParticleSource() const
Returns particle source as defined with enum EParticleSourceObject.
Definition: Particle.h:404
Belle2::Particle::getECLClusterEnergy
double getECLClusterEnergy() const
Returns the energy of the ECLCluster for the particle.
Definition: Particle.cc:844
Belle2::Particle::getAcoplanarity
float getAcoplanarity() const
Returns acoplanarity angle defined as the angle between the decay planes of the grand daughters in th...
Definition: Particle.cc:502
Belle2::Particle::getV0
const V0 * getV0() const
Returns the pointer to the V0 object that was used to create this Particle (if ParticleType == c_V0).
Definition: Particle.cc:805
Belle2::Particle::setMomentumPositionErrorMatrix
void setMomentumPositionErrorMatrix(const TrackFitResult *trackFit)
Sets the momentum, position and error matrix for this particle (created from charged Track)
Definition: Particle.cc:925
Belle2::Particle::fillFSPDaughters
void fillFSPDaughters(std::vector< const Belle2::Particle * > &fspDaughters) const
Fill final state particle daughters into a vector.
Definition: Particle.cc:1021
Belle2::ClusterUtils::GetIPPositionCovarianceMatrix
const TMatrixDSym GetIPPositionCovarianceMatrix()
Returns default IP position covariance matrix from beam parameters.
Definition: ClusterUtils.cc:191
Belle2::TrackFitResult::getCovariance6
TMatrixDSym getCovariance6() const
Position and Momentum Covariance Matrix.
Definition: TrackFitResult.h:147
Belle2::ECLCluster::getConnectedRegionId
int getConnectedRegionId() const
Return connected region id.
Definition: ECLCluster.h:262
Belle2::Particle::getMCParticle
const MCParticle * getMCParticle() const
Returns the pointer to the MCParticle object that was used to create this Particle (ParticleType == c...
Definition: Particle.cc:866
Belle2::Particle::getName
virtual std::string getName() const
Return name of this particle.
Definition: Particle.cc:1059
Belle2::Particle::getMomentumErrorMatrix
TMatrixFSym getMomentumErrorMatrix() const
Returns the 4x4 momentum error matrix.
Definition: Particle.cc:401
Belle2::TrackFitResult::getChargeSign
short getChargeSign() const
Return track charge (1 or -1).
Definition: TrackFitResult.h:160
Belle2::Particle::getInfoHTML
virtual std::string getInfoHTML() const
Return a short summary of this object's contents in HTML format.
Definition: Particle.cc:1082
Belle2::Particle::getExtraInfoNames
std::vector< std::string > getExtraInfoNames() const
get a list of the extra info names
Definition: Particle.cc:1069
Belle2::Particle::~Particle
~Particle()
Destructor.
Belle2::Particle::hasExtraInfo
bool hasExtraInfo(const std::string &name) const
Return wether the extra info with the given name is set.
Definition: Particle.cc:1154
Belle2::Particle::forEachDaughter
bool forEachDaughter(const std::function< bool(const Particle *)> &function, bool recursive=true, bool includeSelf=true) const
Apply a function to all daughters of this particle.
Definition: Particle.cc:1246