Belle II Software development
ParameterVariables.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// Own header.
10#include <analysis/variables/ParameterVariables.h>
11
12// include VariableManager
13#include <analysis/VariableManager/Manager.h>
14
15#include <analysis/dataobjects/Particle.h>
16#include <analysis/utility/PCmsLabTransform.h>
17#include <analysis/utility/ReferenceFrame.h>
18
19#include <framework/logging/Logger.h>
20#include <framework/datastore/StoreArray.h>
21
22#include <mdst/dataobjects/MCParticle.h>
23
24#include <Math/Boost.h>
25#include <Math/Vector4D.h>
26#include <Math/VectorUtil.h>
27#include <TVectorF.h>
28
29#include <cmath>
30
31
32namespace Belle2 {
37 namespace Variable {
38
39 bool almostContains(const std::vector<double>& vector, int value)
40 {
41 for (const auto& item : vector)
42 if (std::abs(value - item) < 1e-3)
43 return true;
44 return false;
45 }
46
47 double RandomChoice(const Particle*, const std::vector<double>& choices)
48 {
49 int r = std::rand() % choices.size() + 1;
50 auto it = choices.begin();
51 std::advance(it, r);
52 return *it;
53 }
54
55 int NumberOfMCParticlesInEvent(const Particle*, const std::vector<double>& pdgs)
56 {
57 StoreArray<MCParticle> mcParticles;
58 int counter = 0;
59 for (int i = 0; i < mcParticles.getEntries(); ++i) {
60 if (mcParticles[i]->getStatus(MCParticle::c_PrimaryParticle) and almostContains(pdgs, std::abs(mcParticles[i]->getPDG())))
61 counter++;
62 }
63 return counter;
64 }
65
66 double isAncestorOf(const Particle* part, const std::vector<double>& daughterIDs)
67 {
68 if (part == nullptr)
69 return Const::doubleNaN;
70
71 // If particle has no MC relation, MC chain doesn't exist
72 const MCParticle* mcpart = part->getMCParticle();
73 if (mcpart == nullptr)
74 return Const::doubleNaN;
75
76 if (daughterIDs.empty())
77 B2FATAL("Wrong number of arguments for parameter function isAncestorOf. At least one needed!");
78
79 // Get to the daughter of interest
80 const Particle* curParticle = part;
81 double isAncestor = 0.0;
82
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 Const::doubleNaN;
93 curParticle = curDaughter;
94 }
95
96 // Daughter obtained, get MC particle of daughter
97 const MCParticle* finalMCDaughter = curParticle->getMCParticle();
98 if (finalMCDaughter == nullptr)
99 return Const::doubleNaN;
100
101 // Go up the MC chain, check for ancestor
102 const MCParticle* curMCParticle = finalMCDaughter;
103
104 while (curMCParticle != nullptr) {
105 const MCParticle* curMCMother = curMCParticle->getMother();
106 if (curMCMother == nullptr)
107 return 0.0;
108 else {
109 if (curMCMother->getArrayIndex() == mcpart->getArrayIndex()) {
110 isAncestor++;
111 break;
112 } else {
113 curMCParticle = curMCMother;
114 isAncestor++;
115 }
116 }
117 }
118 return isAncestor;
119 }
120
121 double hasAncestor(const Particle* part, const std::vector<double>& args)
122 {
123 if (part == nullptr)
124 return Const::doubleNaN;
125
126 // If particle has no MC relation, MC chain doesn't exist
127 const MCParticle* mcpart = part->getMCParticle();
128 if (mcpart == nullptr)
129 return Const::doubleNaN;
130
131 int m_PDG, m_sign = 0;
132
133 if (args.empty())
134 B2FATAL("Wrong number of arguments for variable hasAncestor!");
135 else if (args.size() == 1) {
136 if (args[0] == 0)
137 B2FATAL("PDG code in variable hasAncestor is 0!");
138 else
139 m_PDG = args[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!");
143 else {
144 m_PDG = args[0];
145 m_sign = args[1];
146 }
147 } else {
148 B2FATAL("Too many arguments for variable hasAncestor!");
149 }
150
151 unsigned int nLevels = 0;
152
153 const MCParticle* curMCParticle = mcpart;
154
155 while (curMCParticle != nullptr) {
156 const MCParticle* curMCMother = curMCParticle->getMother();
157 if (curMCMother == nullptr)
158 return 0;
159
160 int pdg = curMCMother->getPDG();
161 if (m_sign == 0)
162 pdg = abs(pdg);
163
164 if (pdg == m_PDG) {
165 ++nLevels;
166 break;
167 } else {
168 ++nLevels;
169 curMCParticle = curMCMother;
170 }
171 }
172 return nLevels;
173 }
174
175
176 double daughterInvariantMass(const Particle* particle, const std::vector<double>& daughter_indexes)
177 {
178 if (!particle)
179 return Const::doubleNaN;
180
181 ROOT::Math::PxPyPzEVector sum;
182 const auto& daughters = particle->getDaughters();
183 int nDaughters = static_cast<int>(daughters.size());
184
185 for (auto& double_daughter : daughter_indexes) {
186 long daughter = std::lround(double_daughter);
187 if (daughter >= nDaughters)
188 return Const::doubleNaN;
189
190 sum += daughters[daughter]->get4Vector();
191 }
192
193 return sum.M();
194 }
195
196 double daughterMCInvariantMass(const Particle* particle, const std::vector<double>& daughter_indexes)
197 {
198 if (!particle)
199 return Const::doubleNaN;
200
201 ROOT::Math::PxPyPzEVector sum;
202 const auto& daughters = particle->getDaughters();
203 int nDaughters = static_cast<int>(daughters.size());
204
205 for (auto& double_daughter : daughter_indexes) {
206 long daughter = std::lround(double_daughter);
207 if (daughter >= nDaughters)
208 return Const::doubleNaN;
209
210 const MCParticle* mcdaughter = daughters[daughter]->getMCParticle();
211 if (!mcdaughter)
212 return Const::doubleNaN;
213
214 sum += mcdaughter->get4Vector();
215 }
216
217 return sum.M();
218 }
219
220
221 double massDifference(const Particle* particle, const std::vector<double>& daughters)
222 {
223 if (!particle)
224 return Const::doubleNaN;
225
226 long daughter = std::lround(daughters[0]);
227 if (daughter >= static_cast<int>(particle->getNDaughters()))
228 return Const::doubleNaN;
229
230 double motherMass = particle->getMass();
231 double daughterMass = particle->getDaughter(daughter)->getMass();
232
233 return motherMass - daughterMass;
234 }
235
236 double massDifferenceError(const Particle* particle, const std::vector<double>& daughters)
237 {
238 if (!particle)
239 return Const::doubleNaN;
240
241 long daughter = std::lround(daughters[0]);
242 if (daughter >= static_cast<int>(particle->getNDaughters()))
243 return Const::doubleNaN;
244
245 float result = 0.0;
246
247 ROOT::Math::PxPyPzEVector thisDaughterMomentum = particle->getDaughter(daughter)->get4Vector();
248
249 TMatrixFSym thisDaughterCovM(Particle::c_DimMomentum);
250 thisDaughterCovM = particle->getDaughter(daughter)->getMomentumErrorMatrix();
251 TMatrixFSym othrDaughterCovM(Particle::c_DimMomentum);
252
253 for (int j = 0; j < int(particle->getNDaughters()); ++j) {
254 if (j == daughter)
255 continue;
256
257 othrDaughterCovM += particle->getDaughter(j)->getMomentumErrorMatrix();
258 }
259
260 TMatrixFSym covarianceMatrix(2 * Particle::c_DimMomentum);
261 covarianceMatrix.SetSub(0, thisDaughterCovM);
262 covarianceMatrix.SetSub(4, othrDaughterCovM);
263
264 double motherMass = particle->getMass();
265 double daughterMass = particle->getDaughter(daughter)->getMass();
266
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;
276
277 result = jacobian * (covarianceMatrix * jacobian);
278
279 if (result < 0.0)
280 result = 0.0;
281
282 return TMath::Sqrt(result);
283 }
284
285 double massDifferenceSignificance(const Particle* particle, const std::vector<double>& daughters)
286 {
287 if (!particle)
288 return Const::doubleNaN;
289
290 long daughter = std::lround(daughters[0]);
291 if (daughter >= static_cast<int>(particle->getNDaughters()))
292 return Const::doubleNaN;
293
294 double massDiff = massDifference(particle, daughters);
295 double massDiffErr = massDifferenceError(particle, daughters);
296
297 double massDiffNominal = particle->getPDGMass() - particle->getDaughter(daughter)->getPDGMass();
298
299 return (massDiff - massDiffNominal) / massDiffErr;
300 }
301
302 // Decay Kinematics -------------------------------------------------------
303 double particleDecayAngle(const Particle* particle, const std::vector<double>& daughters)
304 {
305 if (!particle)
306 return Const::doubleNaN;
307
308 PCmsLabTransform T;
309 ROOT::Math::PxPyPzEVector m = - T.getBeamFourMomentum();
310
311 ROOT::Math::PxPyPzEVector motherMomentum = particle->get4Vector();
312 ROOT::Math::XYZVector motherBoost = motherMomentum.BoostToCM();
313
314 long daughter = std::lround(daughters[0]);
315 if (daughter >= static_cast<int>(particle->getNDaughters()))
316 return Const::doubleNaN;
317
318 ROOT::Math::PxPyPzEVector daugMomentum = particle->getDaughter(daughter)->get4Vector();
319 daugMomentum = ROOT::Math::Boost(motherBoost) * daugMomentum;
320
321 m = ROOT::Math::Boost(motherBoost) * m;
322
323 return ROOT::Math::VectorUtil::Angle(daugMomentum, m);
324 }
325
326 double pointingAngle(const Particle* particle, const std::vector<double>& daughters)
327 {
328 if (!particle)
329 return Const::doubleNaN;
330
331 long daughter = std::lround(daughters[0]);
332 if (daughter >= static_cast<int>(particle->getNDaughters()))
333 return Const::doubleNaN;
334
335 if (particle->getDaughter(daughter)->getNDaughters() < 2)
336 return Const::doubleNaN;
337
338 ROOT::Math::XYZVector productionVertex = particle->getVertex();
339 ROOT::Math::XYZVector decayVertex = particle->getDaughter(daughter)->getVertex();
340
341 ROOT::Math::XYZVector vertexDiffVector = decayVertex - productionVertex;
342
343 const auto& frame = ReferenceFrame::GetCurrent();
344 ROOT::Math::PxPyPzEVector daughterMomentumVector = frame.getMomentum(particle->getDaughter(daughter));
345
346 return ROOT::Math::VectorUtil::Angle(daughterMomentumVector, vertexDiffVector);
347 }
348
349 double azimuthalAngleInDecayPlane(const Particle* particle, const std::vector<double>& daughters)
350 {
351 if (!particle)
352 return Const::doubleNaN;
353
354 int nDaughters = static_cast<int>(particle->getNDaughters());
355
356 long daughter1 = std::lround(daughters[0]);
357 long daughter2 = std::lround(daughters[1]);
358 if (daughter1 >= nDaughters || daughter2 >= nDaughters)
359 return Const::doubleNaN;
360
361 PCmsLabTransform T;
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();
366
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())));
376
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;
387
388 return TMath::ACos(-m_parallel_times_d1 / (m_parallel_abs * d1_abs));
389 }
390
391 double Constant(const Particle*, const std::vector<double>& constant)
392 {
393 return constant[0];
394 }
395
396
397
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.
401
402 Used in the FEI to determine to calculate reconstruction efficiencies.
403
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.
407
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(
411
412 Returns a positive integer if an ancestor with the given PDG code is found, 0 otherwise.
413
414 The integer is the level where the ancestor was found, 1: first mother, 2: grandmother, etc.
415
416 Second argument is optional, 1 means that the sign of the PDG code is taken into account, default is 0.
417
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.:
421
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.
424
425Useful to identify intermediate resonances in a decay, which weren't reconstructed explicitly.
426
427Returns NaN if particle is nullptr or if the given daughter-index is out of bound (>= number of daughters).
428
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.
434
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
438
439)DOC", "rad");
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.
443
444)DOC", "rad");
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.
447
448First we define the following symbols:
449
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)
455
456L can be defined via the following relation:
457
458.. math:: L^{\sigma} = \delta^{\sigma\nu} \epsilon_{\mu\nu\alpha\beta} P^{\mu}D1^{\alpha}D2^{\beta}
459
460The azimuthal angle is given by
461
462.. math:: \phi \equiv \cos^{-1} \left(\frac{-\vec{M_{\parallel}} \cdot \vec{D1}}{|\vec{M_{\parallel}}| \cdot |\vec{D1}|}\right)
463
464For a frame independent formulation the three component vectors need to be written via invariant four-momentum vectors.
465
466.. math::
467
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}}
473
474)DOC", "rad");
475
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]");
482
483 REGISTER_VARIABLE("constant(float i)", Constant, R"DOC(
484 Returns i.
485
486 Useful for debugging purposes and in conjunction with the formula meta-variable.)DOC");
487
488 REGISTER_VARIABLE("randomChoice(i, j, ...)", RandomChoice, R"DOC(
489 Returns random element of given numbers.
490
491 Useful for testing purposes.)DOC");
492
493
494 }
496}
static const double doubleNaN
quiet_NaN
Definition Const.h:703
@ c_PrimaryParticle
bit 0: Particle is primary particle.
Definition MCParticle.h:47
static const ReferenceFrame & GetCurrent()
Get current rest frame.
#define MAKE_DEPRECATED(name, make_fatal, version, description)
Registers a variable as deprecated.
Definition Manager.h:448
Abstract base class for different kinds of events.