10#include <analysis/variables/ParameterVariables.h>
13#include <analysis/VariableManager/Manager.h>
15#include <analysis/dataobjects/Particle.h>
16#include <analysis/utility/PCmsLabTransform.h>
17#include <analysis/utility/ReferenceFrame.h>
19#include <framework/logging/Logger.h>
20#include <framework/datastore/StoreArray.h>
22#include <mdst/dataobjects/MCParticle.h>
24#include <Math/Boost.h>
25#include <Math/Vector4D.h>
26#include <Math/VectorUtil.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 int 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)
72 const MCParticle* mcpart = part->getMCParticle();
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)
93 curParticle = curDaughter;
97 const MCParticle* finalMCDaughter = curParticle->getMCParticle();
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)
127 const MCParticle* mcpart = part->getMCParticle();
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)
181 ROOT::Math::PxPyPzEVector sum;
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)
190 sum += daughters[daughter]->get4Vector();
196 double daughterMCInvariantMass(
const Particle* particle,
const std::vector<double>& daughter_indexes)
201 ROOT::Math::PxPyPzEVector sum;
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)
210 const MCParticle* mcdaughter = daughters[daughter]->getMCParticle();
214 sum += mcdaughter->get4Vector();
221 double massDifference(
const Particle* particle,
const std::vector<double>& daughters)
226 long daughter = std::lround(daughters[0]);
227 if (daughter >=
static_cast<int>(particle->getNDaughters()))
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)
241 long daughter = std::lround(daughters[0]);
242 if (daughter >=
static_cast<int>(particle->getNDaughters()))
247 ROOT::Math::PxPyPzEVector 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)
290 long daughter = std::lround(daughters[0]);
291 if (daughter >=
static_cast<int>(particle->getNDaughters()))
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)
309 ROOT::Math::PxPyPzEVector m = - T.getBeamFourMomentum();
311 ROOT::Math::PxPyPzEVector motherMomentum = particle->get4Vector();
312 ROOT::Math::XYZVector motherBoost = motherMomentum.BoostToCM();
314 long daughter = std::lround(daughters[0]);
315 if (daughter >=
static_cast<int>(particle->getNDaughters()))
318 ROOT::Math::PxPyPzEVector daugMomentum = particle->getDaughter(daughter)->get4Vector();
319 daugMomentum = ROOT::Math::Boost(motherBoost) * daugMomentum;
321 m = ROOT::Math::Boost(motherBoost) * m;
323 return ROOT::Math::VectorUtil::Angle(daugMomentum, m);
326 double pointingAngle(
const Particle* particle,
const std::vector<double>& daughters)
331 long daughter = std::lround(daughters[0]);
332 if (daughter >=
static_cast<int>(particle->getNDaughters()))
335 if (particle->getDaughter(daughter)->getNDaughters() < 2)
338 ROOT::Math::XYZVector productionVertex = particle->getVertex();
339 ROOT::Math::XYZVector decayVertex = particle->getDaughter(daughter)->getVertex();
341 ROOT::Math::XYZVector vertexDiffVector = decayVertex - productionVertex;
344 ROOT::Math::PxPyPzEVector daughterMomentumVector = frame.getMomentum(particle->getDaughter(daughter));
346 return ROOT::Math::VectorUtil::Angle(daughterMomentumVector, vertexDiffVector);
349 double azimuthalAngleInDecayPlane(
const Particle* particle,
const std::vector<double>& daughters)
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)
362 ROOT::Math::PxPyPzEVector m = T.getBeamFourMomentum();
363 ROOT::Math::PxPyPzEVector p = particle->get4Vector();
364 ROOT::Math::PxPyPzEVector d1 = particle->getDaughter(daughter1)->get4Vector();
365 ROOT::Math::PxPyPzEVector d2 = particle->getDaughter(daughter2)->get4Vector();
367 ROOT::Math::PxPyPzEVector l;
368 l.SetPx(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.SetPy(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.SetPz(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.Dot(p);
378 double m_times_l = m.Dot(l);
379 double m_times_d1 = m.Dot(d1);
380 double l_times_d1 = l.Dot(d1);
381 double d1_times_p = d1.Dot(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 Constant(
const Particle*,
const std::vector<double>& constant)
398 VARIABLE_GROUP(
"ParameterFunctions");
399 REGISTER_VARIABLE(
"NumberOfMCParticlesInEvent(pdgcode)", NumberOfMCParticlesInEvent, R
"DOC(
400 [Eventbased] Returns number of MC Particles (including anti-particles) with the given pdgcode in the event.
402 Used in the FEI to determine to calculate reconstruction efficiencies.
404 The variable is event-based and does not need a valid particle pointer as input.)DOC");
405 REGISTER_VARIABLE("isAncestorOf(i, j, ...)", isAncestorOf, R
"DOC(
406 Returns a positive integer if daughter at position particle->daughter(i)->daughter(j)... is an ancestor of the related MC particle, 0 otherwise.
408 Positive integer represents the number of steps needed to get from final MC daughter to ancestor.
409 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");
410 REGISTER_VARIABLE("hasAncestor(PDG, abs)", hasAncestor, R
"DOC(
412 Returns a positive integer if an ancestor with the given PDG code is found, 0 otherwise.
414 The integer is the level where the ancestor was found, 1: first mother, 2: grandmother, etc.
416 Second argument is optional, 1 means that the sign of the PDG code is taken into account, default is 0.
418 If there is no MC relations found, -1 is returned. In case of nullptr particle, NaN is returned.)DOC");
419 REGISTER_VARIABLE("daughterInvariantMass(i, j, ...)", daughterInvariantMass, R
"DOC(
420Returns invariant mass of the given daughter particles. E.g.:
422* daughterInvariantMass(0, 1) returns the invariant mass of the first and second daughter.
423* daughterInvariantMass(0, 1, 2) returns the invariant mass of the first, second and third daughter.
425Useful to identify intermediate resonances in a decay, which weren't reconstructed explicitly.
427Returns NaN if particle is nullptr or if the given daughter-index is out of bound (>= number of daughters).
429)DOC", "GeV/:math:`\\text{c}^2`");
430 MAKE_DEPRECATED(
"daughterInvariantMass",
true,
"light-2203-zeus", R
"DOC(
431 The variable `daughterInvM` provides exactly the same functionality.)DOC");
432 REGISTER_VARIABLE("daughterMCInvariantMass(i, j, ...)", daughterMCInvariantMass, R
"DOC(
433Returns true invariant mass of the given daughter particles, same behaviour as daughterInvariantMass variable.
435)DOC", "GeV/:math:`\\text{c}^2`");
436 REGISTER_VARIABLE(
"decayAngle(i)", particleDecayAngle, R
"DOC(
437Angle in the mother's rest frame between the reverted CMS momentum vector and the direction of the i-th daughter
440 REGISTER_VARIABLE(
"pointingAngle(i)", pointingAngle, R
"DOC(
441Angle between i-th daughter's momentum vector and vector connecting production and decay vertex of i-th daughter.
442This makes only sense if the i-th daughter has itself daughter particles and therefore a properly defined vertex.
445 REGISTER_VARIABLE(
"azimuthalAngleInDecayPlane(i, j)", azimuthalAngleInDecayPlane, R
"DOC(
446Azimuthal angle of i-th daughter in decay plane towards projection of particle momentum into decay plane.
448First we define the following symbols:
450* P: four-momentum vector of decaying particle in whose decay plane the azimuthal angle is measured
451* M: "mother" of p, however not necessarily the direct mother but any higher state, here the CMS itself is chosen
452* D1: daughter for which the azimuthal angle is supposed to be calculated
453* D2: another daughter needed to span the decay plane
454* L: normal to the decay plane (four-component vector)
456L can be defined via the following relation:
458.. math:: L^{\sigma} = \delta^{\sigma\nu} \epsilon_{\mu\nu\alpha\beta} P^{\mu}D1^{\alpha}D2^{\beta}
460The azimuthal angle is given by
462.. math:: \phi \equiv \cos^{-1} \left(\frac{-\vec{M_{\parallel}} \cdot \vec{D1}}{|\vec{M_{\parallel}}| \cdot |\vec{D1}|}\right)
464For a frame independent formulation the three component vectors need to be written via invariant four-momentum vectors.
468 -\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}\\
469 |\vec{M_{\parallel}}| &= |\vec{M}| \sqrt{1 - \cos^2 \psi}\\
470 |\vec{M}| &= \sqrt{\frac{(M \cdot P)^2}{m^2_P} - m^2_M}\\
471 \cos \psi &= \frac{\vec{M} \cdot \vec{L}}{|\vec{M}| \cdot |\vec{L}|} = \frac{-M \cdot L}{|\vec{M}| \cdot \sqrt{-L^2}}\\
472 |\vec{D1}| &= \sqrt{\frac{(D1 \cdot P)^2}{m^2_P} - m^2_{D1}}
476 REGISTER_VARIABLE(
"massDifference(i)", massDifference,
"Difference in invariant masses of this particle and its i-th daughter\n\n",
477 "GeV/:math:`\\text{c}^2`");
478 REGISTER_VARIABLE(
"massDifferenceError(i)", massDifferenceError,
479 "Estimated uncertainty on difference in invariant masses of this particle and its i-th daughter\n\n",
"GeV/:math:`\\text{c}^2`");
480 REGISTER_VARIABLE(
"massDifferenceSignificance(i)", massDifferenceSignificance,
481 "Signed significance of the deviation from the nominal mass difference of this particle and its i-th daughter [(massDiff - NOMINAL_MASS_DIFF)/ErrMassDiff]");
483 REGISTER_VARIABLE(
"constant(float i)", Constant, R
"DOC(
486 Useful for debugging purposes and in conjunction with the formula meta-variable.)DOC");
488 REGISTER_VARIABLE("randomChoice(i, j, ...)", RandomChoice, R
"DOC(
489 Returns random element of given numbers.
491 Useful for testing purposes.)DOC");
static const double doubleNaN
quiet_NaN
@ c_PrimaryParticle
bit 0: Particle is primary particle.
static const ReferenceFrame & GetCurrent()
Get current rest frame.
#define MAKE_DEPRECATED(name, make_fatal, version, description)
Registers a variable as deprecated.
Abstract base class for different kinds of events.