11 #include <analysis/variables/ParameterVariables.h>
12 #include <analysis/VariableManager/Manager.h>
13 #include <analysis/dataobjects/Particle.h>
14 #include <analysis/utility/PCmsLabTransform.h>
15 #include <analysis/utility/ReferenceFrame.h>
17 #include <framework/logging/Logger.h>
18 #include <framework/datastore/StoreArray.h>
20 #include <mdst/dataobjects/MCParticle.h>
22 #include <mdst/dataobjects/Track.h>
23 #include <mdst/dataobjects/TrackFitResult.h>
25 #include <TLorentzVector.h>
39 bool almostContains(
const std::vector<double>& vector,
int value)
41 for (
const auto& item : vector)
42 if (std::abs(value - item) < 1e-3)
47 double RandomChoice(
const Particle*,
const std::vector<double>& choices)
49 int r = std::rand() % choices.size() + 1;
50 auto it = choices.begin();
55 double NumberOfMCParticlesInEvent(
const Particle*,
const std::vector<double>& pdgs)
57 StoreArray<MCParticle> mcParticles;
59 for (
int i = 0; i < mcParticles.getEntries(); ++i) {
66 double isAncestorOf(
const Particle* part,
const std::vector<double>& daughterIDs)
69 return std::numeric_limits<float>::quiet_NaN();
72 const MCParticle* mcpart = part->getRelatedTo<MCParticle>();
73 if (mcpart ==
nullptr)
76 if (daughterIDs.empty())
77 B2FATAL(
"Wrong number of arguments for parameter function isAncestorOf. At least one needed!");
80 const Particle* curParticle = part;
81 double isAncestor = 0.0;
83 for (
unsigned int i = 0; i < daughterIDs.size(); i++) {
84 int nCurDaughters = curParticle->getNDaughters();
85 if (nCurDaughters == 0)
86 B2FATAL(
"Assumed mother of particle at argument " << i <<
" has no daughters!");
87 if (daughterIDs[i] >= nCurDaughters)
88 B2FATAL(
"Assumed mother of particle at argument " << i <<
" has only " << nCurDaughters
89 <<
" daughters, but daughter at position " << daughterIDs[i] <<
" expected!");
90 const Particle* curDaughter = curParticle->getDaughter(daughterIDs[i]);
91 if (curDaughter ==
nullptr)
92 return std::numeric_limits<float>::quiet_NaN();
93 curParticle = curDaughter;
97 const MCParticle* finalMCDaughter = curParticle->
getRelatedTo<MCParticle>();
98 if (finalMCDaughter ==
nullptr)
102 const MCParticle* curMCParticle = finalMCDaughter;
104 while (curMCParticle !=
nullptr) {
105 const MCParticle* curMCMother = curMCParticle->getMother();
106 if (curMCMother ==
nullptr)
109 if (curMCMother->getArrayIndex() == mcpart->getArrayIndex()) {
113 curMCParticle = curMCMother;
121 double hasAncestor(
const Particle* part,
const std::vector<double>& args)
124 return std::numeric_limits<float>::quiet_NaN();
127 const MCParticle* mcpart = part->getRelatedTo<MCParticle>();
128 if (mcpart ==
nullptr)
131 int m_PDG, m_sign = 0;
134 B2FATAL(
"Wrong number of arguments for variable hasAncestor!");
135 else if (args.size() == 1) {
137 B2FATAL(
"PDG code in variable hasAncestor is 0!");
140 }
else if (args.size() == 2) {
141 if (args[0] == 0 or (args[1] != 0 and args[1] != 1))
142 B2FATAL(
"PDG code in variable hasAncestor is 0 or second argument is not 0 or 1!");
148 B2FATAL(
"Too many arguments for variable hasAncestor!");
151 unsigned int nLevels = 0;
153 const MCParticle* curMCParticle = mcpart;
155 while (curMCParticle !=
nullptr) {
156 const MCParticle* curMCMother = curMCParticle->getMother();
157 if (curMCMother ==
nullptr)
160 int pdg = curMCMother->getPDG();
169 curMCParticle = curMCMother;
176 double daughterInvariantMass(
const Particle* particle,
const std::vector<double>& daughter_indexes)
179 return std::numeric_limits<float>::quiet_NaN();
182 const auto& daughters = particle->getDaughters();
183 int nDaughters =
static_cast<int>(daughters.size());
185 for (
auto& double_daughter : daughter_indexes) {
186 long daughter = std::lround(double_daughter);
187 if (daughter >= nDaughters)
188 return std::numeric_limits<float>::quiet_NaN();
190 sum += daughters[daughter]->get4Vector();
196 double daughterMCInvariantMass(
const Particle* particle,
const std::vector<double>& daughter_indexes)
199 return std::numeric_limits<float>::quiet_NaN();
202 const auto& daughters = particle->getDaughters();
203 int nDaughters =
static_cast<int>(daughters.size());
205 for (
auto& double_daughter : daughter_indexes) {
206 long daughter = std::lround(double_daughter);
207 if (daughter >= nDaughters)
208 return std::numeric_limits<float>::quiet_NaN();
210 const MCParticle* mcdaughter = daughters[daughter]->
getRelated<MCParticle>();
212 return std::numeric_limits<float>::quiet_NaN();
214 sum += mcdaughter->get4Vector();
221 double massDifference(
const Particle* particle,
const std::vector<double>& daughters)
224 return std::numeric_limits<float>::quiet_NaN();
226 long daughter = std::lround(daughters[0]);
227 if (daughter >=
static_cast<int>(particle->getNDaughters()))
228 return std::numeric_limits<float>::quiet_NaN();
230 double motherMass = particle->getMass();
231 double daughterMass = particle->getDaughter(daughter)->getMass();
233 return motherMass - daughterMass;
236 double massDifferenceError(
const Particle* particle,
const std::vector<double>& daughters)
239 return std::numeric_limits<float>::quiet_NaN();
241 long daughter = std::lround(daughters[0]);
242 if (daughter >=
static_cast<int>(particle->getNDaughters()))
243 return std::numeric_limits<float>::quiet_NaN();
247 TLorentzVector thisDaughterMomentum = particle->getDaughter(daughter)->get4Vector();
249 TMatrixFSym thisDaughterCovM(Particle::c_DimMomentum);
250 thisDaughterCovM = particle->getDaughter(daughter)->getMomentumErrorMatrix();
251 TMatrixFSym othrDaughterCovM(Particle::c_DimMomentum);
253 for (
int j = 0; j < int(particle->getNDaughters()); ++j) {
257 othrDaughterCovM += particle->getDaughter(j)->getMomentumErrorMatrix();
260 TMatrixFSym covarianceMatrix(2 * Particle::c_DimMomentum);
261 covarianceMatrix.SetSub(0, thisDaughterCovM);
262 covarianceMatrix.SetSub(4, othrDaughterCovM);
264 double motherMass = particle->getMass();
265 double daughterMass = particle->getDaughter(daughter)->getMass();
267 TVectorF jacobian(2 * Particle::c_DimMomentum);
268 jacobian[0] = thisDaughterMomentum.Px() / daughterMass - particle->getPx() / motherMass;
269 jacobian[1] = thisDaughterMomentum.Py() / daughterMass - particle->getPy() / motherMass;
270 jacobian[2] = thisDaughterMomentum.Pz() / daughterMass - particle->getPz() / motherMass;
271 jacobian[3] = particle->getEnergy() / motherMass - thisDaughterMomentum.E() / daughterMass;
272 jacobian[4] = -1.0 * particle->getPx() / motherMass;
273 jacobian[5] = -1.0 * particle->getPy() / motherMass;
274 jacobian[6] = -1.0 * particle->getPz() / motherMass;
275 jacobian[7] = 1.0 * particle->getEnergy() / motherMass;
277 result = jacobian * (covarianceMatrix * jacobian);
282 return TMath::Sqrt(result);
285 double massDifferenceSignificance(
const Particle* particle,
const std::vector<double>& daughters)
288 return std::numeric_limits<float>::quiet_NaN();
290 long daughter = std::lround(daughters[0]);
291 if (daughter >=
static_cast<int>(particle->getNDaughters()))
292 return std::numeric_limits<float>::quiet_NaN();
294 double massDiff = massDifference(particle, daughters);
295 double massDiffErr = massDifferenceError(particle, daughters);
297 double massDiffNominal = particle->getPDGMass() - particle->getDaughter(daughter)->getPDGMass();
299 return (massDiff - massDiffNominal) / massDiffErr;
303 double particleDecayAngle(
const Particle* particle,
const std::vector<double>& daughters)
306 return std::numeric_limits<float>::quiet_NaN();
309 TLorentzVector m = - T.getBeamFourMomentum();
311 TLorentzVector motherMomentum = particle->get4Vector();
312 TVector3 motherBoost = -(motherMomentum.BoostVector());
314 long daughter = std::lround(daughters[0]);
315 if (daughter >=
static_cast<int>(particle->getNDaughters()))
316 return std::numeric_limits<float>::quiet_NaN();
318 TLorentzVector daugMomentum = particle->getDaughter(daughter)->get4Vector();
319 daugMomentum.Boost(motherBoost);
321 m.Boost(motherBoost);
323 return daugMomentum.Angle(m.Vect());
326 double pointingAngle(
const Particle* particle,
const std::vector<double>& daughters)
329 return std::numeric_limits<float>::quiet_NaN();
331 long daughter = std::lround(daughters[0]);
332 if (daughter >=
static_cast<int>(particle->getNDaughters()))
333 return std::numeric_limits<float>::quiet_NaN();
335 if (particle->getDaughter(daughter)->getNDaughters() < 2)
336 return std::numeric_limits<float>::quiet_NaN();
338 TVector3 productionVertex = particle->getVertex();
339 TVector3 decayVertex = particle->getDaughter(daughter)->getVertex();
341 TVector3 vertexDiffVector = decayVertex - productionVertex;
344 TVector3 daughterMomentumVector = frame.getMomentum(particle->getDaughter(daughter)).Vect();
346 return daughterMomentumVector.Angle(vertexDiffVector);
349 double azimuthalAngleInDecayPlane(
const Particle* particle,
const std::vector<double>& daughters)
352 return std::numeric_limits<float>::quiet_NaN();
354 int nDaughters =
static_cast<int>(particle->getNDaughters());
356 long daughter1 = std::lround(daughters[0]);
357 long daughter2 = std::lround(daughters[1]);
358 if (daughter1 >= nDaughters || daughter2 >= nDaughters)
359 return std::numeric_limits<float>::quiet_NaN();
362 TLorentzVector m = T.getBeamFourMomentum();
363 TLorentzVector p = particle->get4Vector();
364 TLorentzVector d1 = particle->getDaughter(daughter1)->get4Vector();
365 TLorentzVector d2 = particle->getDaughter(daughter2)->get4Vector();
368 l.SetX(p.Py() * (d1.Pz() * d2.E() - d1.E() * d2.Pz()) + p.Pz() * (d1.E() * d2.Py() - d1.Py() * d2.E())
369 + p.E() * (d1.Py() * d2.Pz() - d1.Pz() * d2.Py()));
370 l.SetY(p.Px() * (d1.E() * d2.Pz() - d1.Pz() * d2.E()) + p.Pz() * (d1.Px() * d2.E() - d1.E() * d2.Px())
371 + p.E() * (d1.Pz() * d2.Px() - d1.Px() * d2.Pz()));
372 l.SetZ(p.Px() * (d1.Py() * d2.E() - d1.E() * d2.Py()) + p.Py() * (d1.E() * d2.Px() - d1.Px() * d2.E())
373 + p.E() * (d1.Px() * d2.Py() - d1.Py() * d2.Px()));
374 l.SetE(-(p.Px() * (d1.Pz() * d2.Py() - d1.Py() * d2.Pz()) + p.Py() * (d1.Px() * d2.Pz() - d1.Pz() * d2.Px())
375 + p.Pz() * (d1.Py() * d2.Px() - d1.Px() * d2.Py())));
377 double m_times_p = m * p;
378 double m_times_l = m * l;
379 double m_times_d1 = m * d1;
380 double l_times_d1 = l * d1;
381 double d1_times_p = d1 * p;
382 double m_abs = TMath::Sqrt(pow(m_times_p / p.M(), 2) - m.M2());
383 double d1_abs = TMath::Sqrt(pow(d1_times_p / p.M(), 2) - d1.M2());
384 double cos_phi = -m_times_l / (m_abs * TMath::Sqrt(-l.M2()));
385 double m_parallel_abs = m_abs * TMath::Sqrt(1 - cos_phi * cos_phi);
386 double m_parallel_times_d1 = m_times_p * d1_times_p / p.M2() + m_times_l * l_times_d1 / l.M2() - m_times_d1;
388 return TMath::ACos(-m_parallel_times_d1 / (m_parallel_abs * d1_abs));
391 double v0DaughterD0(
const Particle* particle,
const std::vector<double>& daughterID)
394 return std::numeric_limits<float>::quiet_NaN();
396 TVector3 v0Vertex = particle->getVertex();
398 const Particle* daug = particle->getDaughter(daughterID[0]);
400 const Track* track = daug->getTrack();
401 if (!track)
return std::numeric_limits<float>::quiet_NaN();
403 const TrackFitResult* trackFit = track->getTrackFitResultWithClosestMass(Const::ChargedStable(abs(daug->getPDGCode())));
404 if (!trackFit)
return std::numeric_limits<float>::quiet_NaN();
406 UncertainHelix helix = trackFit->getUncertainHelix();
407 helix.passiveMoveBy(v0Vertex);
409 return helix.getD0();
412 double v0DaughterD0Diff(
const Particle* particle)
414 return v0DaughterD0(particle, {0}) - v0DaughterD0(particle, {1});
417 double v0DaughterZ0(
const Particle* particle,
const std::vector<double>& daughterID)
420 return std::numeric_limits<float>::quiet_NaN();
422 TVector3 v0Vertex = particle->getVertex();
424 const Particle* daug = particle->getDaughter(daughterID[0]);
426 const Track* track = daug->getTrack();
427 if (!track)
return std::numeric_limits<float>::quiet_NaN();
429 const TrackFitResult* trackFit = track->getTrackFitResultWithClosestMass(Const::ChargedStable(abs(daug->getPDGCode())));
430 if (!trackFit)
return std::numeric_limits<float>::quiet_NaN();
432 UncertainHelix helix = trackFit->getUncertainHelix();
433 helix.passiveMoveBy(v0Vertex);
435 return helix.getZ0();
438 double v0DaughterZ0Diff(
const Particle* particle)
440 return v0DaughterZ0(particle, {0}) - v0DaughterZ0(particle, {1});
443 double Constant(
const Particle*,
const std::vector<double>& constant)
450 VARIABLE_GROUP(
"ParameterFunctions");
451 REGISTER_VARIABLE(
"NumberOfMCParticlesInEvent(pdgcode)", NumberOfMCParticlesInEvent , R
"DOC(
452 Returns number of MC Particles (including anti-particles) with the given pdgcode in the event.
454 Used in the FEI to determine to calculate reconstruction efficiencies.
456 The variable is event-based and does not need a valid particle pointer as input.)DOC");
457 REGISTER_VARIABLE("isAncestorOf(i, j, ...)", isAncestorOf, R
"DOC(
458 Returns a positive integer if daughter at position particle->daughter(i)->daughter(j)... is an ancestor of the related MC particle, 0 otherwise.
460 Positive integer represents the number of steps needed to get from final MC daughter to ancestor.
461 If any particle or MCparticle is a nullptr, NaN is returned. If MC relations of any particle doesn't exist, -1.0 is returned.)DOC");
462 REGISTER_VARIABLE("hasAncestor(PDG, abs)", hasAncestor, R
"DOC(
464 Returns a positive integer if an ancestor with the given PDG code is found, 0 otherwise.
466 The integer is the level where the ancestor was found, 1: first mother, 2: grandmother, etc.
468 Second argument is optional, 1 means that the sign of the PDG code is taken into account, default is 0.
470 If there is no MC relations found, -1 is returned. In case of nullptr particle, NaN is returned.)DOC");
471 REGISTER_VARIABLE("daughterInvariantMass(i, j, ...)", daughterInvariantMass , R
"DOC(
472 Returns invariant mass of the given daughter particles. E.g.:
474 * daughterInvariantMass(0, 1) returns the invariant mass of the first and second daughter.
475 * daughterInvariantMass(0, 1, 2) returns the invariant mass of the first, second and third daughter.
477 Useful to identify intermediate resonances in a decay, which weren't reconstructed explicitly.
479 Returns NaN if particle is nullptr or if the given daughter-index is out of bound (>= amount of daughters).)DOC");
480 REGISTER_VARIABLE("daughterMCInvariantMass(i, j, ...)", daughterMCInvariantMass ,
481 "Returns true invariant mass of the given daughter particles, same behaviour as daughterInvariantMass variable.");
482 REGISTER_VARIABLE(
"decayAngle(i)", particleDecayAngle,
483 "Angle in the mother's rest frame between the reverted CMS momentum vector and the direction of the i-th daughter");
484 REGISTER_VARIABLE(
"pointingAngle(i)", pointingAngle, R
"DOC(
485 Angle between i-th daughter's momentum vector and vector connecting production and decay vertex of i-th daughter.
486 This makes only sense if the i-th daughter has itself daughter particles and therefore a properly defined vertex.)DOC");
487 REGISTER_VARIABLE("azimuthalAngleInDecayPlane(i, j)", azimuthalAngleInDecayPlane, R
"DOC(
488 Azimuthal angle of i-th daughter in decay plane towards projection of particle momentum into decay plane.
490 First we define the following symbols:
492 * P: four-momentum vector of decaying particle in whose decay plane the azimuthal angle is measured
493 * M: "mother" of p, however not necessarily the direct mother but any higher state, here the CMS itself is chosen
494 * D1: daughter for which the azimuthal angle is supposed to be calculated
495 * D2: another daughter needed to span the decay plane
496 * L: normal to the decay plane (four-component vector)
498 L can be defined via the following relation:
500 .. math:: L^{\sigma} = \delta^{\sigma\nu} \epsilon_{\mu\nu\alpha\beta} P^{\mu}D1^{\alpha}D2^{\beta}
502 The azimuthal angle is given by
504 .. math:: \phi \equiv \cos^{-1} \left(\frac{-\vec{M_{\parallel}} \cdot \vec{D1}}{|\vec{M_{\parallel}}| \cdot |\vec{D1}|}\right)
506 For a frame independent formulation the three component vectors need to be written via invariant four-momentum vectors.
510 -\vec{M_{\parallel}} \cdot \vec{D1} &= \biggl[M - \frac{(M \cdot L)L}{L^2}\biggr] \cdot D1 - \frac{(M \cdot P)(D1 \cdot P)}{m^2_P}\\
511 |\vec{M_{\parallel}}| &= |\vec{M}| \sqrt{1 - \cos^2 \psi}\\
512 |\vec{M}| &= \sqrt{\frac{(M \cdot P)^2}{m^2_P} - m^2_M}\\
513 \cos \psi &= \frac{\vec{M} \cdot \vec{L}}{|\vec{M}| \cdot |\vec{L}|} = \frac{-M \cdot L}{|\vec{M}| \cdot \sqrt{-L^2}}\\
514 |\vec{D1}| &= \sqrt{\frac{(D1 \cdot P)^2}{m^2_P} - m^2_{D1}}
518 REGISTER_VARIABLE("massDifference(i)", massDifference,
"Difference in invariant masses of this particle and its i-th daughter");
519 REGISTER_VARIABLE(
"massDifferenceError(i)", massDifferenceError,
520 "Estimated uncertainty on difference in invariant masses of this particle and its i-th daughter");
521 REGISTER_VARIABLE(
"massDifferenceSignificance(i)", massDifferenceSignificance,
522 "Signed significance of the deviation from the nominal mass difference of this particle and its i-th daughter [(massDiff - NOMINAL_MASS_DIFF)/ErrMassDiff]");
524 REGISTER_VARIABLE(
"V0d0(id)", v0DaughterD0,
525 "Return the d0 impact parameter of a V0's daughter with daughterID index with the V0 vertex point as a pivot for the track.");
526 REGISTER_VARIABLE(
"V0Deltad0", v0DaughterD0Diff,
527 "Return the difference between d0 impact parameters of V0's daughters with the V0 vertex point as a pivot for the track.");
528 REGISTER_VARIABLE(
"V0z0(id)", v0DaughterZ0,
529 "Return the z0 impact parameter of a V0's daughter with daughterID index with the V0 vertex point as a pivot for the track.");
530 REGISTER_VARIABLE(
"V0Deltaz0", v0DaughterZ0Diff,
531 "Return the difference between z0 impact parameters of V0's daughters with the V0 vertex point as a pivot for the track.");
533 REGISTER_VARIABLE(
"constant(float i)", Constant, R
"DOC(
536 Useful for debugging purposes and in conjunction with the formula meta-variable.)DOC");
538 REGISTER_VARIABLE("randomChoice(i, j, ...)", RandomChoice, R
"DOC(
539 Returns random element of given numbers.
541 Useful for testing purposes.)DOC");