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