12 #include <analysis/variables/Variables.h>
13 #include <analysis/VariableManager/Manager.h>
14 #include <analysis/utility/PCmsLabTransform.h>
15 #include <analysis/utility/ReferenceFrame.h>
17 #include <analysis/utility/MCMatching.h>
18 #include <analysis/ClusterUtility/ClusterUtils.h>
21 #include <framework/datastore/StoreArray.h>
22 #include <framework/datastore/StoreObjPtr.h>
25 #include <analysis/dataobjects/Particle.h>
26 #include <analysis/dataobjects/EventExtraInfo.h>
27 #include <analysis/dataobjects/EventShapeContainer.h>
29 #include <mdst/dataobjects/MCParticle.h>
30 #include <mdst/dataobjects/Track.h>
31 #include <mdst/dataobjects/ECLCluster.h>
33 #include <mdst/dbobjects/BeamSpot.h>
36 #include <framework/logging/Logger.h>
38 #include <TLorentzVector.h>
57 double particleP(
const Particle* part)
59 const auto& frame = ReferenceFrame::GetCurrent();
60 return frame.getMomentum(part).P();
64 double particleE(
const Particle* part)
66 const auto& frame = ReferenceFrame::GetCurrent();
67 return frame.getMomentum(part).E();
70 double particleClusterEUncertainty(
const Particle* part)
72 const ECLCluster* cluster = part->getECLCluster();
73 const auto EPhiThetaCov = cluster->getCovarianceMatrix3x3();
75 return std::sqrt(EPhiThetaCov[0][0]);
77 return std::numeric_limits<double>::quiet_NaN();
80 double particlePx(
const Particle* part)
82 const auto& frame = ReferenceFrame::GetCurrent();
83 return frame.getMomentum(part).Px();
86 double particlePy(
const Particle* part)
88 const auto& frame = ReferenceFrame::GetCurrent();
89 return frame.getMomentum(part).Py();
92 double particlePz(
const Particle* part)
94 const auto& frame = ReferenceFrame::GetCurrent();
95 return frame.getMomentum(part).Pz();
98 double particlePt(
const Particle* part)
100 const auto& frame = ReferenceFrame::GetCurrent();
101 return frame.getMomentum(part).Pt();
104 double covMatrixElement(
const Particle* part,
const std::vector<double>& element)
106 int elementI = int(std::lround(element[0]));
107 int elementJ = int(std::lround(element[1]));
109 if (elementI < 0 || elementI > 6) {
110 B2WARNING(
"Requested particle's momentumVertex covariance matrix element is out of boundaries [0 - 6]: i = " << elementI);
111 return std::numeric_limits<double>::quiet_NaN();
113 if (elementJ < 0 || elementJ > 6) {
114 B2WARNING(
"Requested particle's momentumVertex covariance matrix element is out of boundaries [0 - 6]: j = " << elementJ);
115 return std::numeric_limits<double>::quiet_NaN();
118 return part->getMomentumVertexErrorMatrix()(elementI, elementJ);
121 double particleEUncertainty(
const Particle* part)
123 const auto& frame = ReferenceFrame::GetCurrent();
125 double errorSquared = frame.getMomentumErrorMatrix(part)(3, 3);
127 if (errorSquared >= 0.0)
128 return sqrt(errorSquared);
130 return std::numeric_limits<double>::quiet_NaN();
133 double particlePErr(
const Particle* part)
135 TMatrixD jacobianRot(3, 3);
138 double cosPhi = cos(particlePhi(part));
139 double sinPhi = sin(particlePhi(part));
140 double cosTheta = particleCosTheta(part);
141 double sinTheta = sin(acos(cosTheta));
142 double p = particleP(part);
144 jacobianRot(0, 0) = sinTheta * cosPhi;
145 jacobianRot(0, 1) = sinTheta * sinPhi;
146 jacobianRot(1, 0) = cosTheta * cosPhi / p;
147 jacobianRot(1, 1) = cosTheta * sinPhi / p;
148 jacobianRot(0, 2) = cosTheta;
149 jacobianRot(2, 0) = -sinPhi / sinTheta / p;
150 jacobianRot(1, 2) = -sinTheta / p;
151 jacobianRot(2, 1) = cosPhi / sinTheta / p;
153 const auto& frame = ReferenceFrame::GetCurrent();
155 double errorSquared = frame.getMomentumErrorMatrix(part).GetSub(0, 2, 0, 2,
" ").Similarity(jacobianRot)(0, 0);
157 if (errorSquared >= 0.0)
158 return sqrt(errorSquared);
160 return std::numeric_limits<double>::quiet_NaN();
163 double particlePxErr(
const Particle* part)
165 const auto& frame = ReferenceFrame::GetCurrent();
167 double errorSquared = frame.getMomentumErrorMatrix(part)(0, 0);
169 if (errorSquared >= 0.0)
170 return sqrt(errorSquared);
172 return std::numeric_limits<double>::quiet_NaN();
175 double particlePyErr(
const Particle* part)
177 const auto& frame = ReferenceFrame::GetCurrent();
178 double errorSquared = frame.getMomentumErrorMatrix(part)(1, 1);
180 if (errorSquared >= 0.0)
181 return sqrt(errorSquared);
183 return std::numeric_limits<double>::quiet_NaN();
186 double particlePzErr(
const Particle* part)
188 const auto& frame = ReferenceFrame::GetCurrent();
189 double errorSquared = frame.getMomentumErrorMatrix(part)(2, 2);
191 if (errorSquared >= 0.0)
192 return sqrt(errorSquared);
194 return std::numeric_limits<double>::quiet_NaN();
197 double particlePtErr(
const Particle* part)
199 TMatrixD jacobianRot(3, 3);
202 double px = particlePx(part);
203 double py = particlePy(part);
204 double pt = particlePt(part);
206 jacobianRot(0, 0) = px / pt;
207 jacobianRot(0, 1) = py / pt;
208 jacobianRot(1, 0) = -py / (pt * pt);
209 jacobianRot(1, 1) = px / (pt * pt);
210 jacobianRot(2, 2) = 1;
212 const auto& frame = ReferenceFrame::GetCurrent();
213 double errorSquared = frame.getMomentumErrorMatrix(part).GetSub(0, 2, 0, 2,
" ").Similarity(jacobianRot)(0, 0);
215 if (errorSquared >= 0.0)
216 return sqrt(errorSquared);
218 return std::numeric_limits<double>::quiet_NaN();
222 double momentumDeviationChi2(
const Particle* part)
224 double result = std::numeric_limits<double>::quiet_NaN();
227 if (part->getPValue() < 0.0)
231 const MCParticle* mcp = part->getRelated<MCParticle>();
236 result += TMath::Power(part->getPx() - mcp->getMomentum().Px(), 2.0) / part->getMomentumVertexErrorMatrix()(0, 0);
237 result += TMath::Power(part->getPy() - mcp->getMomentum().Py(), 2.0) / part->getMomentumVertexErrorMatrix()(1, 1);
238 result += TMath::Power(part->getPz() - mcp->getMomentum().Pz(), 2.0) / part->getMomentumVertexErrorMatrix()(2, 2);
243 double particleTheta(
const Particle* part)
245 const auto& frame = ReferenceFrame::GetCurrent();
246 return acos(frame.getMomentum(part).CosTheta());
249 double particleThetaErr(
const Particle* part)
251 TMatrixD jacobianRot(3, 3);
254 double cosPhi = cos(particlePhi(part));
255 double sinPhi = sin(particlePhi(part));
256 double cosTheta = particleCosTheta(part);
257 double sinTheta = sin(acos(cosTheta));
258 double p = particleP(part);
260 jacobianRot(0, 0) = sinTheta * cosPhi;
261 jacobianRot(0, 1) = sinTheta * sinPhi;
262 jacobianRot(1, 0) = cosTheta * cosPhi / p;
263 jacobianRot(1, 1) = cosTheta * sinPhi / p;
264 jacobianRot(0, 2) = cosTheta;
265 jacobianRot(2, 0) = -sinPhi / sinTheta / p;
266 jacobianRot(1, 2) = -sinTheta / p;
267 jacobianRot(2, 1) = cosPhi / sinTheta / p;
269 const auto& frame = ReferenceFrame::GetCurrent();
270 double errorSquared = frame.getMomentumErrorMatrix(part).GetSub(0, 2, 0, 2,
" ").Similarity(jacobianRot)(1, 1);
272 if (errorSquared >= 0.0)
273 return sqrt(errorSquared);
275 return std::numeric_limits<double>::quiet_NaN();
278 double particleCosTheta(
const Particle* part)
280 const auto& frame = ReferenceFrame::GetCurrent();
281 return frame.getMomentum(part).CosTheta();
284 double particleCosThetaErr(
const Particle* part)
286 return fabs(particleThetaErr(part) * sin(particleTheta(part)));
289 double particlePhi(
const Particle* part)
291 const auto& frame = ReferenceFrame::GetCurrent();
292 return frame.getMomentum(part).Phi();
295 double particlePhiErr(
const Particle* part)
297 TMatrixD jacobianRot(3, 3);
300 double px = particlePx(part);
301 double py = particlePy(part);
302 double pt = particlePt(part);
304 jacobianRot(0, 0) = px / pt;
305 jacobianRot(0, 1) = py / pt;
306 jacobianRot(1, 0) = -py / (pt * pt);
307 jacobianRot(1, 1) = px / (pt * pt);
308 jacobianRot(2, 2) = 1;
310 const auto& frame = ReferenceFrame::GetCurrent();
311 double errorSquared = frame.getMomentumErrorMatrix(part).GetSub(0, 2, 0, 2,
" ").Similarity(jacobianRot)(1, 1);
313 if (errorSquared >= 0.0)
314 return sqrt(errorSquared);
316 return std::numeric_limits<double>::quiet_NaN();
320 double particleXp(
const Particle* part)
323 TLorentzVector p4 = part -> get4Vector();
324 TLorentzVector p4CMS = T.rotateLabToCms() * p4;
325 float s = T.getCMSEnergy();
326 float M = part->getMass();
327 return p4CMS.P() / TMath::Sqrt(s * s / 4 - M * M);
330 double particlePDGCode(
const Particle* part)
332 return part->getPDGCode();
335 double cosAngleBetweenMomentumAndVertexVectorInXYPlane(
const Particle* part)
337 static DBObjPtr<BeamSpot> beamSpotDB;
338 double px = part->getMomentum().Px();
339 double py = part->getMomentum().Py();
341 double xV = part->getVertex().X();
342 double yV = part->getVertex().Y();
344 double xIP = (beamSpotDB->getIPPosition()).X();
345 double yIP = (beamSpotDB->getIPPosition()).Y();
350 double cosangle = (px * x + py * y) / (sqrt(px * px + py * py) * sqrt(x * x + y * y));
354 double cosAngleBetweenMomentumAndVertexVector(
const Particle* part)
356 static DBObjPtr<BeamSpot> beamSpotDB;
357 return std::cos((part->getVertex() - beamSpotDB->getIPPosition()).Angle(part->getMomentum()));
360 double cosThetaBetweenParticleAndNominalB(
const Particle* part)
363 int particlePDG = abs(part->getPDGCode());
364 if (particlePDG != 511 and particlePDG != 521)
365 B2FATAL(
"The Variables cosThetaBetweenParticleAndNominalB is only meant to be used on B mesons!");
369 double e_Beam = 1.0579400E+1 / 2.0;
370 double m_B = part->getPDGMass();
372 double p_B = std::sqrt(e_Beam * e_Beam - m_B * m_B);
374 TLorentzVector p = T.rotateLabToCms() * part->get4Vector();
377 double p_d = p.Rho();
379 double theta_BY = (2 * e_Beam * e_d - m_B * m_B - m_d * m_d)
384 double cosToThrustOfEvent(
const Particle* part)
386 StoreObjPtr<EventShapeContainer> evtShape;
388 B2WARNING(
"Cannot find thrust of event information, did you forget to load the event shape calculation?");
389 return std::numeric_limits<float>::quiet_NaN();
392 TVector3 th = evtShape->getThrustAxis();
393 TVector3 particleMomentum = (T.rotateLabToCms() * part -> get4Vector()).Vect();
394 return std::cos(th.Angle(particleMomentum));
397 double ImpactXY(
const Particle* particle)
399 static DBObjPtr<BeamSpot> beamSpotDB;
401 double px = particle->getPx();
402 double py = particle->getPy();
404 if (py == py && px == px) {
406 double x = particle->getX() - (beamSpotDB->getIPPosition()).X();
407 double y = particle->getY() - (beamSpotDB->getIPPosition()).Y();
409 double pt = sqrt(px * px + py * py);
414 double a = -0.2998 * 1.5 * particle->getCharge();
415 double T = TMath::Sqrt(pt * pt - 2 * a * (x * py - y * px) + a * a * (x * x + y * y));
417 return TMath::Abs((-2 * (x * py - y * px) + a * (x * x + y * y)) / (T + pt));
418 }
else return std::numeric_limits<double>::quiet_NaN();
421 double ArmenterosLongitudinalMomentumAsymmetry(
const Particle* part)
423 const auto& frame = ReferenceFrame::GetCurrent();
424 int nDaughters = part -> getNDaughters();
426 B2FATAL(
"You are trying to use an Armenteros variable. The mother particle is required to have exactly two daughters");
428 const auto& daughters = part -> getDaughters();
429 TVector3 motherMomentum = frame.getMomentum(part).Vect();
430 TVector3 daughter1Momentum = frame.getMomentum(daughters[0]).Vect();
431 TVector3 daughter2Momentum = frame.getMomentum(daughters[1]).Vect();
433 int daughter1Charge = daughters[0] -> getCharge();
434 int daughter2Charge = daughters[1] -> getCharge();
435 double daughter1Ql = daughter1Momentum.Dot(motherMomentum) / motherMomentum.Mag();
436 double daughter2Ql = daughter2Momentum.Dot(motherMomentum) / motherMomentum.Mag();
439 if (daughter2Charge > daughter1Charge)
440 Arm_alpha = (daughter2Ql - daughter1Ql) / (daughter2Ql + daughter1Ql);
442 Arm_alpha = (daughter1Ql - daughter2Ql) / (daughter1Ql + daughter2Ql);
447 double ArmenterosDaughter1Qt(
const Particle* part)
449 const auto& frame = ReferenceFrame::GetCurrent();
450 int nDaughters = part -> getNDaughters();
452 B2FATAL(
"You are trying to use an Armenteros variable. The mother particle is required to have exactly two daughters.");
454 const auto& daughters = part -> getDaughters();
455 TVector3 motherMomentum = frame.getMomentum(part).Vect();
456 TVector3 daughter1Momentum = frame.getMomentum(daughters[0]).Vect();
457 double qt = daughter1Momentum.Perp(motherMomentum);
462 double ArmenterosDaughter2Qt(
const Particle* part)
464 const auto& frame = ReferenceFrame::GetCurrent();
465 int nDaughters = part -> getNDaughters();
467 B2FATAL(
"You are trying to use an Armenteros variable. The mother particle is required to have exactly two daughters.");
469 const auto& daughters = part -> getDaughters();
470 TVector3 motherMomentum = frame.getMomentum(part).Vect();
471 TVector3 daughter2Momentum = frame.getMomentum(daughters[1]).Vect();
472 double qt = daughter2Momentum.Perp(motherMomentum);
480 double particleMass(
const Particle* part)
482 return part->getMass();
485 double particleDMass(
const Particle* part)
487 return part->getMass() - part->getPDGMass();
490 double particleInvariantMassFromDaughters(
const Particle* part)
492 const std::vector<Particle*> daughters = part->getDaughters();
493 if (daughters.size() > 0) {
495 for (
auto daughter : daughters)
496 sum += daughter->get4Vector();
500 return part->getMass();
504 double particleInvariantMassLambda(
const Particle* part)
506 const std::vector<Particle*> daughters = part->getDaughters();
507 if (daughters.size() == 2) {
510 TLorentzVector dtsum;
511 double mpi = Const::pionMass;
512 double mpr = Const::protonMass;
513 dt1 = daughters[0]->get4Vector();
514 dt2 = daughters[1]->get4Vector();
515 double E1 = hypot(mpi, dt1.P());
516 double E2 = hypot(mpr, dt2.P());
518 return sqrt((E1 + E2) * (E1 + E2) - dtsum.P() * dtsum.P());
521 return part->getMass();
525 double particleInvariantMassError(
const Particle* part)
527 float invMass = part->getMass();
528 TMatrixFSym covarianceMatrix = part->getMomentumErrorMatrix();
530 TVectorF jacobian(Particle::c_DimMomentum);
531 jacobian[0] = -1.0 * part->getPx() / invMass;
532 jacobian[1] = -1.0 * part->getPy() / invMass;
533 jacobian[2] = -1.0 * part->getPz() / invMass;
534 jacobian[3] = 1.0 * part->getEnergy() / invMass;
536 double result = jacobian * (covarianceMatrix * jacobian);
539 return std::numeric_limits<double>::quiet_NaN();
541 return TMath::Sqrt(result);
544 double particleInvariantMassSignificance(
const Particle* part)
546 return particleDMass(part) / particleInvariantMassError(part);
549 double particleMassSquared(
const Particle* part)
551 TLorentzVector p4 = part->get4Vector();
555 double b2bTheta(
const Particle* part)
558 TLorentzVector pcms = T.rotateLabToCms() * part->get4Vector();
559 TLorentzVector b2bcms(-pcms.Px(), -pcms.Py(), -pcms.Pz(), pcms.E());
560 TLorentzVector b2blab = T.rotateCmsToLab() * b2bcms;
561 return b2blab.Theta();
564 double b2bPhi(
const Particle* part)
567 TLorentzVector pcms = T.rotateLabToCms() * part->get4Vector();
568 TLorentzVector b2bcms(-pcms.Px(), -pcms.Py(), -pcms.Pz(), pcms.E());
569 TLorentzVector b2blab = T.rotateCmsToLab() * b2bcms;
573 double b2bClusterTheta(
const Particle* part)
576 const ECLCluster* cluster = part->getECLCluster();
577 const ECLCluster::EHypothesisBit clusterHypothesis = part->getECLClusterEHypothesisBit();
578 if (!cluster)
return std::numeric_limits<float>::quiet_NaN();
582 TLorentzVector p4Cluster = clutls.Get4MomentumFromCluster(cluster, clusterHypothesis);
586 TLorentzVector pcms = T.rotateLabToCms() * p4Cluster;
587 TLorentzVector b2bcms(-pcms.Px(), -pcms.Py(), -pcms.Pz(), pcms.E());
588 TLorentzVector b2blab = T.rotateCmsToLab() * b2bcms;
589 return b2blab.Theta();
592 double b2bClusterPhi(
const Particle* part)
595 const ECLCluster* cluster = part->getECLCluster();
596 const ECLCluster::EHypothesisBit clusterHypothesis = part->getECLClusterEHypothesisBit();
597 if (!cluster)
return std::numeric_limits<float>::quiet_NaN();
601 TLorentzVector p4Cluster = clutls.Get4MomentumFromCluster(cluster, clusterHypothesis);
605 TLorentzVector pcms = T.rotateLabToCms() * p4Cluster;
606 TLorentzVector b2bcms(-pcms.Px(), -pcms.Py(), -pcms.Pz(), pcms.E());
607 TLorentzVector b2blab = T.rotateCmsToLab() * b2bcms;
614 double particleQ(
const Particle* part)
616 float m = part->getMass();
617 for (
unsigned i = 0; i < part->getNDaughters(); i++) {
618 const Particle* child = part->getDaughter(i);
620 m -= child->getMass();
625 double particleDQ(
const Particle* part)
627 float m = part->getMass() - part->getPDGMass();
628 for (
unsigned i = 0; i < part->getNDaughters(); i++) {
629 const Particle* child = part->getDaughter(i);
631 m -= (child->getMass() - child->getPDGMass());
638 double particleMbc(
const Particle* part)
641 TLorentzVector vec = T.rotateLabToCms() * part->get4Vector();
642 double E = T.getCMSEnergy() / 2;
643 double m2 = E * E - vec.Vect().Mag2();
644 double mbc = m2 >= 0 ? sqrt(m2) : std::numeric_limits<double>::quiet_NaN();
648 double particleDeltaE(
const Particle* part)
651 TLorentzVector vec = T.rotateLabToCms() * part->get4Vector();
652 return vec.E() - T.getCMSEnergy() / 2;
658 void printParticleInternal(
const Particle* p,
int depth)
661 for (
int i = 0; i < depth; i++) {
664 s << p->getPDGCode();
665 const MCParticle* mcp = p->getRelated<MCParticle>();
667 unsigned int flags = MCMatching::getMCErrors(p, mcp);
668 s <<
" -> MC: " << mcp->getPDG() <<
", mcErrors: " << flags <<
" ("
669 << MCMatching::explainFlags(flags) <<
")";
670 s <<
", mc-index " << mcp->getIndex();
672 s <<
" (no MC match)";
674 s <<
", mdst-source " << p->getMdstSource();
676 for (
const auto* daughter : p->getDaughters()) {
677 printParticleInternal(daughter, depth + 1);
681 double printParticle(
const Particle* p)
683 printParticleInternal(p, 0);
688 double particleMCMomentumTransfer2(
const Particle* part)
691 const MCParticle* mcB = part->getRelated<MCParticle>();
694 return std::numeric_limits<double>::quiet_NaN();
696 TLorentzVector pB = mcB->get4Vector();
698 std::vector<MCParticle*> mcDaug = mcB->getDaughters();
701 return std::numeric_limits<double>::quiet_NaN();
707 for (
auto mcTemp : mcDaug) {
708 if (abs(mcTemp->getPDG()) <= 16)
711 pX += mcTemp->get4Vector();
714 TLorentzVector q = pB - pX;
720 double recoilPx(
const Particle* particle)
725 TLorentzVector pIN = T.getBeamFourMomentum();
728 const auto& frame = ReferenceFrame::GetCurrent();
729 return frame.getMomentum(pIN - particle->get4Vector()).Px();
732 double recoilPy(
const Particle* particle)
737 TLorentzVector pIN = T.getBeamFourMomentum();
740 const auto& frame = ReferenceFrame::GetCurrent();
741 return frame.getMomentum(pIN - particle->get4Vector()).Py();
744 double recoilPz(
const Particle* particle)
749 TLorentzVector pIN = T.getBeamFourMomentum();
752 const auto& frame = ReferenceFrame::GetCurrent();
753 return frame.getMomentum(pIN - particle->get4Vector()).Pz();
757 double recoilMomentum(
const Particle* particle)
762 TLorentzVector pIN = T.getBeamFourMomentum();
765 const auto& frame = ReferenceFrame::GetCurrent();
766 return frame.getMomentum(pIN - particle->get4Vector()).P();
769 double recoilMomentumTheta(
const Particle* particle)
774 TLorentzVector pIN = T.getBeamFourMomentum();
777 const auto& frame = ReferenceFrame::GetCurrent();
778 return frame.getMomentum(pIN - particle->get4Vector()).Theta();
781 double recoilMomentumPhi(
const Particle* particle)
786 TLorentzVector pIN = T.getBeamFourMomentum();
789 const auto& frame = ReferenceFrame::GetCurrent();
790 return frame.getMomentum(pIN - particle->get4Vector()).Phi();
793 double recoilEnergy(
const Particle* particle)
798 TLorentzVector pIN = T.getBeamFourMomentum();
801 const auto& frame = ReferenceFrame::GetCurrent();
802 return frame.getMomentum(pIN - particle->get4Vector()).E();
805 double recoilMass(
const Particle* particle)
810 TLorentzVector pIN = T.getBeamFourMomentum();
813 const auto& frame = ReferenceFrame::GetCurrent();
814 return frame.getMomentum(pIN - particle->get4Vector()).M();
817 double recoilMassSquared(
const Particle* particle)
822 TLorentzVector pIN = T.getBeamFourMomentum();
825 const auto& frame = ReferenceFrame::GetCurrent();
826 return frame.getMomentum(pIN - particle->get4Vector()).M2();
829 double m2RecoilSignalSide(
const Particle* part)
832 double beamEnergy = T.getCMSEnergy() / 2.;
833 TLorentzVector tagVec = T.rotateLabToCms()
834 * part->getDaughter(0)->get4Vector();
835 TLorentzVector sigVec = T.rotateLabToCms()
836 * part->getDaughter(1)->get4Vector();
837 tagVec.SetE(-beamEnergy);
838 return (-tagVec - sigVec).M2();
841 double recoilMCDecayType(
const Particle* particle)
843 auto* mcp = particle->getRelatedTo<MCParticle>();
846 return std::numeric_limits<double>::quiet_NaN();
848 MCParticle* mcMother = mcp->getMother();
851 return std::numeric_limits<double>::quiet_NaN();
853 std::vector<MCParticle*> daughters = mcMother->getDaughters();
855 if (daughters.size() != 2)
856 return std::numeric_limits<double>::quiet_NaN();
858 MCParticle* recoilMC =
nullptr;
859 if (daughters[0]->getArrayIndex() == mcp->getArrayIndex())
860 recoilMC = daughters[1];
862 recoilMC = daughters[0];
864 if (!recoilMC->hasStatus(MCParticle::c_PrimaryParticle))
865 return std::numeric_limits<double>::quiet_NaN();
868 checkMCParticleDecay(recoilMC, decayType,
false);
871 checkMCParticleDecay(recoilMC, decayType,
true);
876 void checkMCParticleDecay(MCParticle* mcp,
int& decayType,
bool recursive)
878 int nHadronicParticles = 0;
879 int nPrimaryParticleDaughters = 0;
880 std::vector<MCParticle*> daughters = mcp->getDaughters();
883 for (
auto& daughter : daughters) {
884 if (!daughter->hasStatus(MCParticle::c_PrimaryParticle))
887 nPrimaryParticleDaughters++;
888 if (abs(daughter->getPDG()) > 22)
889 nHadronicParticles++;
892 if (nPrimaryParticleDaughters > 1) {
893 for (
auto& daughter : daughters) {
894 if (!daughter->hasStatus(MCParticle::c_PrimaryParticle))
897 if (abs(daughter->getPDG()) == 12 or abs(daughter->getPDG()) == 14 or abs(daughter->getPDG()) == 16) {
899 if (nHadronicParticles == 0) {
913 checkMCParticleDecay(daughter, decayType, recursive);
918 double nRemainingTracksInEvent(
const Particle* particle)
921 StoreArray<Track> tracks;
922 int event_tracks = tracks.getEntries();
925 const auto& daughters = particle->getFinalStateDaughters();
926 for (
const auto& daughter : daughters) {
927 int pdg = abs(daughter->getPDGCode());
928 if (pdg == 11 or pdg == 13 or pdg == 211 or pdg == 321 or pdg == 2212)
931 return event_tracks - par_tracks;
934 double trackMatchType(
const Particle* particle)
937 double result = std::numeric_limits<double>::quiet_NaN();
939 const ECLCluster* cluster = particle->getECLCluster();
943 if (cluster->isTrack()) {
952 double False(
const Particle*)
957 double True(
const Particle*)
962 double infinity(
const Particle*)
964 double inf = std::numeric_limits<double>::infinity();
968 double random(
const Particle*)
970 return gRandom->Uniform(0, 1);
973 double eventRandom(
const Particle*)
975 std::string key =
"__eventRandom";
976 StoreObjPtr<EventExtraInfo> eventExtraInfo;
977 if (not eventExtraInfo.isValid())
978 eventExtraInfo.create();
979 if (eventExtraInfo->hasExtraInfo(key)) {
980 return eventExtraInfo->getExtraInfo(key);
982 double value = gRandom->Uniform(0, 1);
983 eventExtraInfo->addExtraInfo(key, value);
988 VARIABLE_GROUP(
"Kinematics");
989 REGISTER_VARIABLE(
"p", particleP,
"momentum magnitude");
990 REGISTER_VARIABLE(
"E", particleE,
"energy");
992 REGISTER_VARIABLE(
"E_uncertainty", particleEUncertainty,
"energy uncertainty (sqrt(sigma2))");
993 REGISTER_VARIABLE(
"ECLClusterE_uncertainty", particleClusterEUncertainty,
994 "energy uncertainty as given by the underlying ECL cluster.");
995 REGISTER_VARIABLE(
"px", particlePx,
"momentum component x");
996 REGISTER_VARIABLE(
"py", particlePy,
"momentum component y");
997 REGISTER_VARIABLE(
"pz", particlePz,
"momentum component z");
998 REGISTER_VARIABLE(
"pt", particlePt,
"transverse momentum");
999 REGISTER_VARIABLE(
"xp", particleXp,
1000 "scaled momentum: the momentum of the particle in the CMS as a fraction of its maximum available momentum in the collision");
1001 REGISTER_VARIABLE(
"pErr", particlePErr,
"error of momentum magnitude");
1002 REGISTER_VARIABLE(
"pxErr", particlePxErr,
"error of momentum component x");
1003 REGISTER_VARIABLE(
"pyErr", particlePyErr,
"error of momentum component y");
1004 REGISTER_VARIABLE(
"pzErr", particlePzErr,
"error of momentum component z");
1005 REGISTER_VARIABLE(
"ptErr", particlePtErr,
"error of transverse momentum");
1006 REGISTER_VARIABLE(
"momVertCovM(i,j)", covMatrixElement,
1007 "returns the (i,j)-th element of the MomentumVertex Covariance Matrix (7x7).\n"
1008 "Order of elements in the covariance matrix is: px, py, pz, E, x, y, z.");
1009 REGISTER_VARIABLE(
"momDevChi2", momentumDeviationChi2,
1010 "momentum deviation chi^2 value calculated as"
1011 "chi^2 = sum_i (p_i - mc(p_i))^2/sigma(p_i)^2, where sum runs over i = px, py, pz and"
1012 "mc(p_i) is the mc truth value and sigma(p_i) is the estimated error of i-th component of momentum vector")
1013 REGISTER_VARIABLE("theta", particleTheta, "polar angle in radians");
1014 REGISTER_VARIABLE("thetaErr", particleThetaErr, "error of polar angle in radians");
1015 REGISTER_VARIABLE("cosTheta", particleCosTheta, "momentum cosine of polar angle");
1016 REGISTER_VARIABLE("cosThetaErr", particleCosThetaErr, "error of momentum cosine of polar angle");
1017 REGISTER_VARIABLE("phi", particlePhi, "momentum azimuthal angle in radians");
1018 REGISTER_VARIABLE("phiErr", particlePhiErr, "error of momentum azimuthal angle in radians");
1019 REGISTER_VARIABLE("PDG", particlePDGCode, "PDG code");
1021 REGISTER_VARIABLE("cosAngleBetweenMomentumAndVertexVectorInXYPlane",
1022 cosAngleBetweenMomentumAndVertexVectorInXYPlane,
1023 "cosine of the angle between momentum and vertex vector (vector connecting ip and fitted vertex) of this particle in xy-plane");
1024 REGISTER_VARIABLE("cosAngleBetweenMomentumAndVertexVector",
1025 cosAngleBetweenMomentumAndVertexVector,
1026 "cosine of the angle between momentum and vertex vector (vector connecting ip and fitted vertex) of this particle");
1027 REGISTER_VARIABLE("cosThetaBetweenParticleAndNominalB",
1028 cosThetaBetweenParticleAndNominalB,
1029 "cosine of the angle in CMS between momentum the particle and a nominal B particle. It is somewhere between -1 and 1 if only a massless particle like a neutrino is missing in the reconstruction.");
1030 REGISTER_VARIABLE("cosToThrustOfEvent", cosToThrustOfEvent,
1031 "Returns the cosine of the angle between the particle and the thrust axis of the event, as calculate by the EventShapeCalculator module. buildEventShape() must be run before calling this variable")
1033 REGISTER_VARIABLE("ImpactXY" , ImpactXY , "The impact parameter of the given particle in the xy plane");
1035 REGISTER_VARIABLE("M", particleMass, R"DOC(
1036 The particle's mass.
1038 Note that this is context-dependent variable and can take different values depending on the situation. This should be the "best" value possible with the information provided.
1040 - If this particle is track- or cluster- based, then this is the value of the mass hypothesis.
1041 - If this particle is an MC particle then this is the mass of that particle.
1042 - If this particle is composite, then *initially* this takes the value of the invariant mass of the daughters.
1043 - If this particle is composite and a *mass or vertex fit* has been performed then this may be updated by the fit.
1045 * You will see a difference between this mass and the :b2:var:`InvM`.
1047 REGISTER_VARIABLE("dM", particleDMass, "mass minus nominal mass");
1048 REGISTER_VARIABLE("Q", particleQ, "released energy in decay");
1049 REGISTER_VARIABLE("dQ", particleDQ,
1050 "released energy in decay minus nominal one");
1051 REGISTER_VARIABLE("Mbc", particleMbc, "beam constrained mass");
1052 REGISTER_VARIABLE("deltaE", particleDeltaE, "energy difference");
1053 REGISTER_VARIABLE("M2", particleMassSquared,
1054 "The particle's mass squared.");
1056 REGISTER_VARIABLE("InvM", particleInvariantMassFromDaughters,
1057 "invariant mass (determined from particle's daughter 4-momentum vectors). If this particle has no daughters, defaults to :b2:var:`M`.");
1058 REGISTER_VARIABLE("InvMLambda", particleInvariantMassLambda,
1059 "Invariant mass (determined from particle's daughter 4-momentum vectors), assuming the first daughter is a pion and the second daughter is a proton.\n"
1060 "If the particle has not 2 daughters, it returns just the mass value.");
1062 REGISTER_VARIABLE("ErrM", particleInvariantMassError,
1063 "uncertainty of invariant mass");
1064 REGISTER_VARIABLE("SigM", particleInvariantMassSignificance,
1065 "
signed deviation of particle's invariant mass from its nominal mass in units of the uncertainty on the invariant mass (dM/ErrM)");
1067 REGISTER_VARIABLE("pxRecoil", recoilPx,
1068 "component x of 3-momentum recoiling against given Particle");
1069 REGISTER_VARIABLE("pyRecoil", recoilPy,
1070 "component y of 3-momentum recoiling against given Particle");
1071 REGISTER_VARIABLE("pzRecoil", recoilPz,
1072 "component z of 3-momentum recoiling against given Particle");
1074 REGISTER_VARIABLE("pRecoil", recoilMomentum,
1075 "magnitude of 3 - momentum recoiling against given Particle");
1076 REGISTER_VARIABLE("pRecoilTheta", recoilMomentumTheta,
1077 "Polar angle of a particle's missing momentum");
1078 REGISTER_VARIABLE("pRecoilPhi", recoilMomentumPhi,
1079 "Azimutal angle of a particle's missing momentum");
1080 REGISTER_VARIABLE("eRecoil", recoilEnergy,
1081 "energy recoiling against given Particle");
1082 REGISTER_VARIABLE("mRecoil", recoilMass,
1083 "Invariant mass of the system recoiling against given Particle");
1084 REGISTER_VARIABLE("m2Recoil", recoilMassSquared,
1085 "invariant mass squared of the system recoiling against given Particle");
1086 REGISTER_VARIABLE("m2RecoilSignalSide", m2RecoilSignalSide,
1087 "Squared recoil mass of the signal side which is calculated in the CMS frame under the assumption that the signal and tag side are produced back to back and the tag side energy equals the beam energy. The variable must be applied to the Upsilon and the tag side must be the first, the signal side the second daughter ");
1089 REGISTER_VARIABLE("b2bTheta", b2bTheta,
1090 "Polar angle in the lab system that is back-to-back to the particle in the CMS. Useful for low multiplicity studies.")
1091 REGISTER_VARIABLE("b2bPhi", b2bPhi,
1092 "Azimuthal angle in the lab system that is back-to-back to the particle in the CMS. Useful for low multiplicity studies.")
1093 REGISTER_VARIABLE("b2bClusterTheta", b2bClusterTheta,
1094 "Polar angle in the lab system that is back-to-back to the particle's associated ECLCluster in the CMS. Returns NAN if no cluster is found. Useful for low multiplicity studies.")
1095 REGISTER_VARIABLE("b2bClusterPhi", b2bClusterPhi,
1096 "Azimuthal angle in the lab system that is back-to-back to the particle's associated ECLCluster in the CMS. Returns NAN if no cluster is found. Useful for low multiplicity studies.")
1097 REGISTER_VARIABLE("ArmenterosLongitudinalMomentumAsymmetry", ArmenterosLongitudinalMomentumAsymmetry,
1098 "Longitudinal momentum asymmetry of V0's daughters.\n"
1099 "The mother (V0) is required to have exactly two daughters");
1100 REGISTER_VARIABLE("ArmenterosDaughter1Qt", ArmenterosDaughter1Qt,
1101 "Transverse momentum of the first daughter with respect to the V0 mother.\n"
1102 "The mother is required to have exactly two daughters");
1103 REGISTER_VARIABLE("ArmenterosDaughter2Qt", ArmenterosDaughter2Qt,
1104 "Transverse momentum of the second daughter with respect to the V0 mother.\n"
1105 "The mother is required to have exactly two daughters");
1107 VARIABLE_GROUP("Miscellaneous");
1108 REGISTER_VARIABLE("nRemainingTracksInEvent", nRemainingTracksInEvent,
1109 "Number of tracks in the event - Number of tracks( = charged FSPs) of particle.");
1110 REGISTER_VARIABLE("trackMatchType", trackMatchType,
1111 "-1 particle has no ECL cluster, 0 particle has no associated track, 1 there is a matched track"
1112 "called connected - region(CR) track match");
1114 REGISTER_VARIABLE("decayTypeRecoil", recoilMCDecayType,
1115 "type of the particle decay(no related mcparticle = -1, hadronic = 0, direct leptonic = 1, direct semileptonic = 2,"
1116 "lower level leptonic = 3.");
1118 REGISTER_VARIABLE("printParticle", printParticle,
1119 "For debugging, print Particle and daughter PDG codes, plus MC match. Returns 0.");
1120 REGISTER_VARIABLE("mcMomTransfer2", particleMCMomentumTransfer2,
1121 "Return the true momentum transfer to lepton pair in a B(semi -) leptonic B meson decay.");
1122 REGISTER_VARIABLE("False", False,
1123 "returns always 0, used for testing and debugging.");
1124 REGISTER_VARIABLE("True", True,
1125 "returns always 1, used for testing and debugging.");
1126 REGISTER_VARIABLE("infinity", infinity,
1127 "returns std::numeric_limits<
double>::infinity()");
1128 REGISTER_VARIABLE("random", random,
1129 "return a random number between 0 and 1 for each candidate. Can be used, e.g. for picking a random"
1130 "candidate in the best candidate selection.");
1131 REGISTER_VARIABLE("eventRandom", eventRandom,
1132 "[Eventbased] Returns a random number between 0 and 1 for this event. Can be used, e.g. for applying an event prescale.");