13 #include <analysis/dataobjects/Particle.h>
14 #include <analysis/dataobjects/ParticleExtraInfoMap.h>
16 #include <analysis/ClusterUtility/ClusterUtils.h>
18 #include <mdst/dataobjects/KLMCluster.h>
19 #include <mdst/dataobjects/MCParticle.h>
20 #include <mdst/dataobjects/PIDLikelihood.h>
21 #include <mdst/dataobjects/Track.h>
22 #include <mdst/dataobjects/TrackFitResult.h>
23 #include <mdst/dataobjects/V0.h>
24 #include <mdst/dbobjects/CollisionBoostVector.h>
25 #include <mdst/dbobjects/CollisionInvariantMass.h>
27 #include <framework/datastore/StoreArray.h>
28 #include <framework/datastore/StoreObjPtr.h>
29 #include <framework/database/DBObjPtr.h>
30 #include <framework/logging/Logger.h>
31 #include <framework/utilities/HTML.h>
32 #include <framework/utilities/Conversion.h>
34 #include <TClonesArray.h>
35 #include <TDatabasePDG.h>
36 #include <TMatrixFSym.h>
43 #include <boost/algorithm/string.hpp>
48 m_pdgCode(0), m_mass(0), m_px(0), m_py(0), m_pz(0), m_x(0), m_y(0), m_z(0),
49 m_pValue(nan(
"")), m_flavorType(c_Unflavored), m_particleSource(c_Undefined), m_mdstIndex(0), m_properties(0),
50 m_arrayPointer(nullptr)
57 m_pdgCode(pdgCode), m_mass(0), m_px(0), m_py(0), m_pz(0), m_x(0), m_y(0), m_z(0),
58 m_pValue(-1), m_flavorType(c_Unflavored), m_particleSource(c_Undefined), m_mdstIndex(0), m_properties(0), m_arrayPointer(nullptr)
70 const unsigned mdstIndex) :
71 m_pdgCode(pdgCode), m_mass(0), m_px(0), m_py(0), m_pz(0), m_x(0), m_y(0), m_z(0),
72 m_pValue(-1), m_flavorType(flavorType), m_particleSource(source), m_properties(0), m_arrayPointer(nullptr)
86 const std::vector<int>& daughterIndices,
87 TClonesArray* arrayPointer) :
88 m_pdgCode(0), m_mass(0), m_px(0), m_py(0), m_pz(0), m_x(0), m_y(0), m_z(0),
90 m_daughterIndices(daughterIndices),
91 m_flavorType(c_Unflavored), m_particleSource(c_Undefined), m_mdstIndex(0),
92 m_properties(0), m_arrayPointer(arrayPointer)
101 if (!daughterIndices.empty()) {
104 B2FATAL(
"Composite Particle (with daughters) was constructed outside StoreArray without specifying daughter array!");
115 const std::vector<int>& daughterIndices,
116 const int properties,
117 TClonesArray* arrayPointer) :
118 m_pdgCode(0), m_mass(0), m_px(0), m_py(0), m_pz(0), m_x(0), m_y(0), m_z(0),
120 m_daughterIndices(daughterIndices),
121 m_flavorType(c_Unflavored), m_particleSource(c_Undefined), m_mdstIndex(0),
122 m_arrayPointer(arrayPointer)
132 if (!daughterIndices.empty()) {
135 B2FATAL(
"Composite Particle (with daughters) was constructed outside StoreArray without specifying daughter array!");
147 const std::vector<int>& daughterIndices,
148 const int properties,
149 const std::vector<int>& daughterProperties,
150 TClonesArray* arrayPointer) :
151 m_pdgCode(0), m_mass(0), m_px(0), m_py(0), m_pz(0), m_x(0), m_y(0), m_z(0),
153 m_daughterIndices(daughterIndices),
154 m_flavorType(c_Unflavored), m_particleSource(c_Undefined), m_mdstIndex(0),
155 m_daughterProperties(daughterProperties),
156 m_arrayPointer(arrayPointer)
166 if (!daughterIndices.empty()) {
169 B2FATAL(
"Composite Particle (with daughters) was constructed outside StoreArray without specifying daughter array!");
177 Particle(track ? track->getArrayIndex() : 0, track ? track->getTrackFitResultWithClosestMass(chargedStable) : nullptr,
185 m_pdgCode(0), m_mass(0), m_px(0), m_py(0), m_pz(0), m_x(0), m_y(0), m_z(0),
186 m_pValue(-1), m_flavorType(c_Unflavored), m_particleSource(c_Undefined), m_mdstIndex(0), m_properties(0), m_arrayPointer(nullptr)
188 if (!trackFit)
return;
199 if (TDatabasePDG::Instance()->GetParticle(
m_pdgCode) ==
nullptr)
200 B2FATAL(
"PDG=" <<
m_pdgCode <<
" ***code unknown to TDatabasePDG");
212 m_pdgCode(0), m_mass(0), m_px(0), m_py(0), m_pz(0), m_x(0), m_y(0), m_z(0),
213 m_pValue(-1), m_flavorType(c_Unflavored), m_particleSource(c_Undefined), m_mdstIndex(0), m_properties(0), m_arrayPointer(nullptr)
215 if (!trackFit)
return;
226 if (TDatabasePDG::Instance()->GetParticle(
m_pdgCode) ==
nullptr)
227 B2FATAL(
"PDG=" <<
m_pdgCode <<
" ***code unknown to TDatabasePDG");
236 m_pdgCode(type.getPDGCode()), m_mass(type.getMass()), m_px(0), m_py(0), m_pz(0), m_x(0), m_y(0), m_z(0),
237 m_pValue(-1), m_flavorType(c_Unflavored), m_particleSource(c_Undefined), m_mdstIndex(0), m_properties(0), m_arrayPointer(nullptr)
239 if (!eclCluster)
return;
250 m_px = clustermom.Px();
251 m_py = clustermom.Py();
252 m_pz = clustermom.Pz();
272 m_pdgCode(0), m_mass(0), m_px(0), m_py(0), m_pz(0), m_x(0), m_y(0), m_z(0),
273 m_pValue(-1), m_flavorType(c_Unflavored), m_particleSource(c_Undefined), m_mdstIndex(0), m_properties(0), m_arrayPointer(nullptr)
275 if (!klmCluster)
return;
298 m_pdgCode(0), m_mass(0), m_px(0), m_py(0), m_pz(0), m_x(0), m_y(0), m_z(0),
299 m_pValue(-1), m_flavorType(c_Unflavored), m_particleSource(c_Undefined), m_mdstIndex(0), m_properties(0), m_arrayPointer(nullptr)
301 if (!mcParticle)
return;
334 const int crid = cluster->getConnectedRegionId();
335 const int clusterid = cluster->getClusterId();
338 B2ERROR(
"Particle is of type = ECLCluster has identifier not set and no relation to ECLCluster.\n"
339 "This has happen because old microDST is analysed with newer version of software.");
358 const int crid = cluster->getConnectedRegionId();
359 const int clusterid = cluster->getClusterId();
360 identifier = 1000 * crid + clusterid;
362 B2ERROR(
"Particle is of type = ECLCluster has identifier not set and no relation to ECLCluster.\n"
363 "This has happen because old microDST is analysed with newer version of software.");
377 if (m.GetNrows() != c_DimMatrix || m.GetNcols() != c_DimMatrix) {
379 B2WARNING(
"Error Matrix is not 7x7 ");
388 TMatrixFSym m(c_DimMatrix);
391 for (
int irow = 0; irow < c_DimMatrix; irow++) {
392 for (
int icol = irow; icol < c_DimMatrix; icol++) {
393 m(irow, icol) = m(icol, irow) =
m_errMatrix[element];
408 full.GetSub(0, 3, mom,
"S");
420 full.GetSub(4, 6, pos,
"S");
431 TLorentzVector pMother;
433 pMother = mother->get4Vector();
437 pMother.SetE(cmsMass->getMass());
438 pMother.Boost(cmsBoost->getBoost());
440 pMother.Boost(boost);
443 TLorentzVector pDaughter;
446 pDaughter.Boost(boost);
451 pDaughter = daughter->get4Vector();
454 pDaughter.Boost(boost);
457 pDaughter0.Boost(boost);
459 pDaughter1.Boost(boost);
460 pDaughter.SetVect(pDaughter0.Vect().Cross(pDaughter1.Vect()));
464 double mag2 = pMother.Vect().Mag2() * pDaughter.Vect().Mag2();
465 if (mag2 <= 0)
return std::numeric_limits<float>::quiet_NaN();
466 return (-pMother.Vect()) * pDaughter.Vect() / sqrt(mag2);
473 B2ERROR(
"No daughter of particle 'name' with index 'iDaughter' for calculation of helicity angle"
475 return std::numeric_limits<float>::quiet_NaN();
480 TVector3 boost = -daughter->get4Vector().BoostVector();
484 pMother.Boost(boost);
487 if (daughter->getNDaughters() <= iGrandDaughter) {
488 B2ERROR(
"No grand daughter of daugher 'iDaughter' of particle 'name' with index 'iGrandDaughter' for calculation of helicity angle"
490 return std::numeric_limits<float>::quiet_NaN();
494 TLorentzVector pGrandDaughter = daughter->getDaughter(iGrandDaughter)->get4Vector();
495 pGrandDaughter.Boost(boost);
497 double mag2 = pMother.Vect().Mag2() * pGrandDaughter.Vect().Mag2();
498 if (mag2 <= 0)
return std::numeric_limits<float>::quiet_NaN();
499 return (-pMother.Vect()) * pGrandDaughter.Vect() / sqrt(mag2);
506 B2ERROR(
"Cannot calculate acoplanarity of particle 'name' because the number of daughters is not 2"
508 return std::numeric_limits<float>::quiet_NaN();
513 B2ERROR(
"Cannot calculate acoplanarity of particle 'name' because the number of grand daughters is not 2"
516 return std::numeric_limits<float>::quiet_NaN();
523 TLorentzVector pDaughter0 = daughter0->
get4Vector();
524 pDaughter0.Boost(boost);
526 pGrandDaughter0.Boost(boost);
527 TLorentzVector pDaughter1 = daughter1->
get4Vector();
528 pDaughter1.Boost(boost);
530 pGrandDaughter1.Boost(boost);
533 TVector3 normal0 = pDaughter0.Vect().Cross(pGrandDaughter0.Vect());
534 TVector3 normal1 = -pDaughter1.Vect().Cross(pGrandDaughter1.Vect());
535 double result = normal0.Angle(normal1);
536 if (normal0.Cross(normal1) * pDaughter0.Vect() < 0) result = -result;
572 if (TDatabasePDG::Instance()->GetParticle(pdgCode) ==
nullptr)
573 B2FATAL(
"PDG=" << pdgCode <<
" ***code unknown to TDatabasePDG");
574 m_mass = TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass() ;
579 if (TDatabasePDG::Instance()->GetParticle(
m_pdgCode) ==
nullptr) {
580 B2ERROR(
"PDG=" <<
m_pdgCode <<
" ***code unknown to TDatabasePDG");
583 return TDatabasePDG::Instance()->GetParticle(
m_pdgCode)->Mass();
588 if (TDatabasePDG::Instance()->GetParticle(
m_pdgCode) ==
nullptr) {
589 B2ERROR(
"PDG=" <<
m_pdgCode <<
" ***code unknown to TDatabasePDG");
592 return TDatabasePDG::Instance()->GetParticle(
m_pdgCode)->Charge() / 3.0;
604 std::vector<Particle*> daughters(nDaughters);
607 for (
unsigned i = 0; i < nDaughters; i++)
615 std::vector<const Particle*> fspDaughters;
623 std::vector<int> mdstIndices;
626 if (fsp->getParticleSource() == source)
627 mdstIndices.push_back(fsp->getMdstArrayIndex());
669 for (
auto& thisFSP : thisFSPs)
670 for (
auto& otherFSP : otherFSPs)
671 if (thisFSP->getMdstSource() == otherFSP->getMdstSource())
691 std::vector<int> thisDecayChain(nDaughters * 2);
692 std::vector<int> othrDecayChain(nDaughters * 2);
694 for (
unsigned i = 0; i < nDaughters; i++) {
703 sort(thisDecayChain.begin(), thisDecayChain.end());
704 sort(othrDecayChain.begin(), othrDecayChain.end());
706 for (
unsigned i = 0; i < thisDecayChain.size(); i++)
707 if (thisDecayChain[i] != othrDecayChain[i])
715 if (!doDetailedComparison)
723 B2WARNING(
"Something went wrong: MCParticle is being compared to a non MC Particle. Please check your script!\n"
724 " If the MCParticle <-> Particle comparison happens in the RestOfEventBuilder,\n"
725 " the Rest Of Event may contain signal side particles.");
782 auto* selfrelated = this->getRelatedTo<TrackFitResult>();
789 return selftrack->getTrackFitResultWithClosestMass(
825 const ECLCluster* bestTrackMatchedCluster =
nullptr;
826 double highestEnergy = -1.0;
835 bestTrackMatchedCluster = &cluster;
838 return bestTrackMatchedCluster;
872 const MCParticle* related = this->getRelated<MCParticle>();
883 std::vector<std::string> generalizedIndexes;
884 boost::split(generalizedIndexes, generalizedIndex, boost::is_any_of(
":"));
886 if (generalizedIndexes.empty()) {
887 B2WARNING(
"Generalized index of daughter particle is empty. Skipping.");
898 for (
auto& indexString : generalizedIndexes) {
902 dauIndex = Belle2::convertString<int>(indexString);
903 }
catch (std::invalid_argument&) {
904 B2WARNING(
"Found the string " << indexString <<
"instead of a daughter index.");
909 if (dauIndex >=
int(currentPart->
getNDaughters()) or dauIndex < 0) {
910 B2WARNING(
"Daughter index " << dauIndex <<
" out of range");
911 B2WARNING(
"Trying to access non-existing particle.");
915 currentPart = dauPart;
940 constexpr
unsigned order[] = {c_X, c_Y, c_Z, c_Px, c_Py, c_Pz};
942 TMatrixFSym errMatrix(c_DimMatrix);
943 for (
int i = 0; i < 6; i++) {
944 for (
int j = i; j < 6; j++) {
947 errMatrix(order[j], order[i]) = errMatrix(order[i], order[j]) = cov6(i, j);
968 constexpr
unsigned compMom[] = {c_Px, c_Py, c_Pz};
969 constexpr
unsigned compPos[] = {c_X, c_Y, c_Z};
972 for (
unsigned int i : compMom) {
974 for (
int k = 0; k < 3; k++) {
975 Cov += errMatrix(i, compMom[k]) * dEdp[k];
977 errMatrix(i, c_E) = Cov;
981 for (
unsigned int comp : compPos) {
983 for (
int k = 0; k < 3; k++) {
984 Cov += errMatrix(comp, compMom[k]) * dEdp[k];
986 errMatrix(c_E, comp) = Cov;
991 for (
int i = 0; i < 3; i++) {
992 Cov += errMatrix(compMom[i], compMom[i]) * dEdp[i] * dEdp[i];
994 for (
int i = 0; i < 3; i++) {
996 Cov += 2 * errMatrix(compMom[i], compMom[k]) * dEdp[i] * dEdp[k];
998 errMatrix(c_E, c_E) = Cov;
1012 for (
int irow = 0; irow < c_DimMatrix; irow++) {
1013 for (
int icol = irow; icol < c_DimMatrix; icol++) {
1025 fspDaughters.push_back(
this);
1053 int q3 = nnn % 10; nnn /= 10;
1054 int q2 = nnn % 10; nnn /= 10;
1061 return TDatabasePDG::Instance()->GetParticle(
m_pdgCode)->GetName();
1071 std::vector<std::string> out;
1075 B2FATAL(
"ParticleExtraInfoMap not available, but needed for storing extra info in Particle!");
1077 for (
auto const& ee : map) out.push_back(ee.first);
1084 std::stringstream stream;
1085 stream << std::setprecision(4);
1088 stream <<
" <b>PDGCode</b>=" <<
m_pdgCode;
1089 stream <<
" <b>Charge</b>=" <<
getCharge();
1090 stream <<
" <b>PDGMass</b>=" <<
getPDGMass();
1100 stream <<
" <b>daughterIndices</b>: ";
1102 stream << daughterIndex <<
", ";
1108 stream <<
" <b>daughter PDGCodes</b>: ";
1111 if (p) {stream << p->getPDGCode() <<
", ";}
1112 else {stream <<
"?, ";}
1117 stream <<
" <b>mass</b>=" <<
m_mass;
1121 stream <<
" <b>p</b>=" <<
getP();
1127 stream <<
" <b>p-value of fit</b> (if done): ";
1131 stream <<
" <b>error matrix</b> (px, py, pz, E, x, y ,z):<br>";
1134 stream <<
" <b>extra info</b>=( ";
1137 if (!extraInfoMap) {
1138 B2FATAL(
"ParticleExtraInfoMap not available, but needed for storing extra info in Particle!");
1142 for (
const auto& pair : map) {
1143 if (pair.second < nVars) {
1144 stream << pair.first <<
"=" <<
m_extraInfo[pair.second] <<
" ";
1149 stream <<
") " <<
"<br>";
1151 return stream.str();
1162 if (!extraInfoMap) {
1163 B2FATAL(
"ParticleExtraInfoMap not available, but needed for storing extra info in Particle!");
1165 unsigned int index = extraInfoMap->getIndex(mapID, name);
1180 throw std::runtime_error(std::string(
"getExtraInfo: Value '") + name +
"' not found in Particle!");
1185 if (!extraInfoMap) {
1186 B2FATAL(
"ParticleExtraInfoMap not available, but needed for storing extra info in Particle!");
1188 unsigned int index = extraInfoMap->getIndex(mapID, name);
1190 throw std::runtime_error(std::string(
"getExtraInfo: Value '") + name +
"' not found in Particle!");
1208 throw std::runtime_error(std::string(
"setExtraInfo: Value '") + name +
"' not found in Particle!");
1213 if (!extraInfoMap) {
1214 B2FATAL(
"ParticleExtraInfoMap not available, but needed for storing extra info in Particle!");
1216 unsigned int index = extraInfoMap->getIndex(mapID, name);
1218 throw std::runtime_error(std::string(
"setExtraInfo: Value '") + name +
"' not found in Particle!");
1227 throw std::runtime_error(std::string(
"addExtraInfo: Value '") + name +
"' already set!");
1231 extraInfoMap.create();
1233 unsigned int mapID = extraInfoMap->getMapForNewVar(name);
1239 unsigned int mapID = extraInfoMap->getMapForNewVar(name, oldMapID, insertIndex);
1247 bool recursive,
bool includeSelf)
const
1249 std::queue<const Particle*> qq;
1258 while (!qq.empty()) {
1260 if (
function(p))
return true;
1264 if (recursive || p ==
this)
1265 for (
size_t i = 0; i < p->getNDaughters(); ++i) qq.push(p->getDaughter(i));
1273 int PDGCode = absPDGCode * chargeSign;