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