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 <mdst/dataobjects/Track.h>
25#include <mdst/dataobjects/TrackFitResult.h>
26
27#include <Math/Boost.h>
28#include <Math/Vector4D.h>
29#include <TVectorF.h>
30
31#include <cmath>
32
33
34namespace Belle2 {
39 namespace Variable {
40
41 bool almostContains(const std::vector<double>& vector, int value)
42 {
43 for (const auto& item : vector)
44 if (std::abs(value - item) < 1e-3)
45 return true;
46 return false;
47 }
48
49 double RandomChoice(const Particle*, const std::vector<double>& choices)
50 {
51 int r = std::rand() % choices.size() + 1;
52 auto it = choices.begin();
53 std::advance(it, r);
54 return *it;
55 }
56
57 int NumberOfMCParticlesInEvent(const Particle*, const std::vector<double>& pdgs)
58 {
59 StoreArray<MCParticle> mcParticles;
60 int counter = 0;
61 for (int i = 0; i < mcParticles.getEntries(); ++i) {
62 if (mcParticles[i]->getStatus(MCParticle::c_PrimaryParticle) and almostContains(pdgs, std::abs(mcParticles[i]->getPDG())))
63 counter++;
64 }
65 return counter;
66 }
67
68 double isAncestorOf(const Particle* part, const std::vector<double>& daughterIDs)
69 {
70 if (part == nullptr)
71 return Const::doubleNaN;
72
73 // If particle has no MC relation, MC chain doesn't exist
74 const MCParticle* mcpart = part->getMCParticle();
75 if (mcpart == nullptr)
76 return Const::doubleNaN;
77
78 if (daughterIDs.empty())
79 B2FATAL("Wrong number of arguments for parameter function isAncestorOf. At least one needed!");
80
81 // Get to the daughter of interest
82 const Particle* curParticle = part;
83 double isAncestor = 0.0;
84
85 for (unsigned int i = 0; i < daughterIDs.size(); i++) {
86 int nCurDaughters = curParticle->getNDaughters();
87 if (nCurDaughters == 0)
88 B2FATAL("Assumed mother of particle at argument " << i << " has no daughters!");
89 if (daughterIDs[i] >= nCurDaughters)
90 B2FATAL("Assumed mother of particle at argument " << i << " has only " << nCurDaughters
91 << " daughters, but daughter at position " << daughterIDs[i] << " expected!");
92 const Particle* curDaughter = curParticle->getDaughter(daughterIDs[i]);
93 if (curDaughter == nullptr)
94 return Const::doubleNaN;
95 curParticle = curDaughter;
96 }
97
98 // Daughter obtained, get MC particle of daughter
99 const MCParticle* finalMCDaughter = curParticle->getMCParticle();
100 if (finalMCDaughter == nullptr)
101 return Const::doubleNaN;
102
103 // Go up the MC chain, check for ancestor
104 const MCParticle* curMCParticle = finalMCDaughter;
105
106 while (curMCParticle != nullptr) {
107 const MCParticle* curMCMother = curMCParticle->getMother();
108 if (curMCMother == nullptr)
109 return 0.0;
110 else {
111 if (curMCMother->getArrayIndex() == mcpart->getArrayIndex()) {
112 isAncestor++;
113 break;
114 } else {
115 curMCParticle = curMCMother;
116 isAncestor++;
117 }
118 }
119 }
120 return isAncestor;
121 }
122
123 double hasAncestor(const Particle* part, const std::vector<double>& args)
124 {
125 if (part == nullptr)
126 return Const::doubleNaN;
127
128 // If particle has no MC relation, MC chain doesn't exist
129 const MCParticle* mcpart = part->getMCParticle();
130 if (mcpart == nullptr)
131 return Const::doubleNaN;
132
133 int m_PDG, m_sign = 0;
134
135 if (args.empty())
136 B2FATAL("Wrong number of arguments for variable hasAncestor!");
137 else if (args.size() == 1) {
138 if (args[0] == 0)
139 B2FATAL("PDG code in variable hasAncestor is 0!");
140 else
141 m_PDG = args[0];
142 } else if (args.size() == 2) {
143 if (args[0] == 0 or (args[1] != 0 and args[1] != 1))
144 B2FATAL("PDG code in variable hasAncestor is 0 or second argument is not 0 or 1!");
145 else {
146 m_PDG = args[0];
147 m_sign = args[1];
148 }
149 } else {
150 B2FATAL("Too many arguments for variable hasAncestor!");
151 }
152
153 unsigned int nLevels = 0;
154
155 const MCParticle* curMCParticle = mcpart;
156
157 while (curMCParticle != nullptr) {
158 const MCParticle* curMCMother = curMCParticle->getMother();
159 if (curMCMother == nullptr)
160 return 0;
161
162 int pdg = curMCMother->getPDG();
163 if (m_sign == 0)
164 pdg = abs(pdg);
165
166 if (pdg == m_PDG) {
167 ++nLevels;
168 break;
169 } else {
170 ++nLevels;
171 curMCParticle = curMCMother;
172 }
173 }
174 return nLevels;
175 }
176
177
178 double daughterInvariantMass(const Particle* particle, const std::vector<double>& daughter_indexes)
179 {
180 if (!particle)
181 return Const::doubleNaN;
182
183 ROOT::Math::PxPyPzEVector sum;
184 const auto& daughters = particle->getDaughters();
185 int nDaughters = static_cast<int>(daughters.size());
186
187 for (auto& double_daughter : daughter_indexes) {
188 long daughter = std::lround(double_daughter);
189 if (daughter >= nDaughters)
190 return Const::doubleNaN;
191
192 sum += daughters[daughter]->get4Vector();
193 }
194
195 return sum.M();
196 }
197
198 double daughterMCInvariantMass(const Particle* particle, const std::vector<double>& daughter_indexes)
199 {
200 if (!particle)
201 return Const::doubleNaN;
202
203 ROOT::Math::PxPyPzEVector sum;
204 const auto& daughters = particle->getDaughters();
205 int nDaughters = static_cast<int>(daughters.size());
206
207 for (auto& double_daughter : daughter_indexes) {
208 long daughter = std::lround(double_daughter);
209 if (daughter >= nDaughters)
210 return Const::doubleNaN;
211
212 const MCParticle* mcdaughter = daughters[daughter]->getMCParticle();
213 if (!mcdaughter)
214 return Const::doubleNaN;
215
216 sum += mcdaughter->get4Vector();
217 }
218
219 return sum.M();
220 }
221
222
223 double massDifference(const Particle* particle, const std::vector<double>& daughters)
224 {
225 if (!particle)
226 return Const::doubleNaN;
227
228 long daughter = std::lround(daughters[0]);
229 if (daughter >= static_cast<int>(particle->getNDaughters()))
230 return Const::doubleNaN;
231
232 double motherMass = particle->getMass();
233 double daughterMass = particle->getDaughter(daughter)->getMass();
234
235 return motherMass - daughterMass;
236 }
237
238 double massDifferenceError(const Particle* particle, const std::vector<double>& daughters)
239 {
240 if (!particle)
241 return Const::doubleNaN;
242
243 long daughter = std::lround(daughters[0]);
244 if (daughter >= static_cast<int>(particle->getNDaughters()))
245 return Const::doubleNaN;
246
247 float result = 0.0;
248
249 ROOT::Math::PxPyPzEVector thisDaughterMomentum = particle->getDaughter(daughter)->get4Vector();
250
251 TMatrixFSym thisDaughterCovM(Particle::c_DimMomentum);
252 thisDaughterCovM = particle->getDaughter(daughter)->getMomentumErrorMatrix();
253 TMatrixFSym othrDaughterCovM(Particle::c_DimMomentum);
254
255 for (int j = 0; j < int(particle->getNDaughters()); ++j) {
256 if (j == daughter)
257 continue;
258
259 othrDaughterCovM += particle->getDaughter(j)->getMomentumErrorMatrix();
260 }
261
262 TMatrixFSym covarianceMatrix(2 * Particle::c_DimMomentum);
263 covarianceMatrix.SetSub(0, thisDaughterCovM);
264 covarianceMatrix.SetSub(4, othrDaughterCovM);
265
266 double motherMass = particle->getMass();
267 double daughterMass = particle->getDaughter(daughter)->getMass();
268
269 TVectorF jacobian(2 * Particle::c_DimMomentum);
270 jacobian[0] = thisDaughterMomentum.Px() / daughterMass - particle->getPx() / motherMass;
271 jacobian[1] = thisDaughterMomentum.Py() / daughterMass - particle->getPy() / motherMass;
272 jacobian[2] = thisDaughterMomentum.Pz() / daughterMass - particle->getPz() / motherMass;
273 jacobian[3] = particle->getEnergy() / motherMass - thisDaughterMomentum.E() / daughterMass;
274 jacobian[4] = -1.0 * particle->getPx() / motherMass;
275 jacobian[5] = -1.0 * particle->getPy() / motherMass;
276 jacobian[6] = -1.0 * particle->getPz() / motherMass;
277 jacobian[7] = 1.0 * particle->getEnergy() / motherMass;
278
279 result = jacobian * (covarianceMatrix * jacobian);
280
281 if (result < 0.0)
282 result = 0.0;
283
284 return TMath::Sqrt(result);
285 }
286
287 double massDifferenceSignificance(const Particle* particle, const std::vector<double>& daughters)
288 {
289 if (!particle)
290 return Const::doubleNaN;
291
292 long daughter = std::lround(daughters[0]);
293 if (daughter >= static_cast<int>(particle->getNDaughters()))
294 return Const::doubleNaN;
295
296 double massDiff = massDifference(particle, daughters);
297 double massDiffErr = massDifferenceError(particle, daughters);
298
299 double massDiffNominal = particle->getPDGMass() - particle->getDaughter(daughter)->getPDGMass();
300
301 return (massDiff - massDiffNominal) / massDiffErr;
302 }
303
304 // Decay Kinematics -------------------------------------------------------
305 double particleDecayAngle(const Particle* particle, const std::vector<double>& daughters)
306 {
307 if (!particle)
308 return Const::doubleNaN;
309
310 PCmsLabTransform T;
311 ROOT::Math::PxPyPzEVector m = - T.getBeamFourMomentum();
312
313 ROOT::Math::PxPyPzEVector motherMomentum = particle->get4Vector();
314 B2Vector3D motherBoost = motherMomentum.BoostToCM();
315
316 long daughter = std::lround(daughters[0]);
317 if (daughter >= static_cast<int>(particle->getNDaughters()))
318 return Const::doubleNaN;
319
320 ROOT::Math::PxPyPzEVector daugMomentum = particle->getDaughter(daughter)->get4Vector();
321 daugMomentum = ROOT::Math::Boost(motherBoost) * daugMomentum;
322
323 m = ROOT::Math::Boost(motherBoost) * m;
324
325 return B2Vector3D(daugMomentum.Vect()).Angle(B2Vector3D(m.Vect()));
326 }
327
328 double pointingAngle(const Particle* particle, const std::vector<double>& daughters)
329 {
330 if (!particle)
331 return Const::doubleNaN;
332
333 long daughter = std::lround(daughters[0]);
334 if (daughter >= static_cast<int>(particle->getNDaughters()))
335 return Const::doubleNaN;
336
337 if (particle->getDaughter(daughter)->getNDaughters() < 2)
338 return Const::doubleNaN;
339
340 B2Vector3D productionVertex = particle->getVertex();
341 B2Vector3D decayVertex = particle->getDaughter(daughter)->getVertex();
342
343 B2Vector3D vertexDiffVector = decayVertex - productionVertex;
344
345 const auto& frame = ReferenceFrame::GetCurrent();
346 B2Vector3D daughterMomentumVector = frame.getMomentum(particle->getDaughter(daughter)).Vect();
347
348 return daughterMomentumVector.Angle(vertexDiffVector);
349 }
350
351 double azimuthalAngleInDecayPlane(const Particle* particle, const std::vector<double>& daughters)
352 {
353 if (!particle)
354 return Const::doubleNaN;
355
356 int nDaughters = static_cast<int>(particle->getNDaughters());
357
358 long daughter1 = std::lround(daughters[0]);
359 long daughter2 = std::lround(daughters[1]);
360 if (daughter1 >= nDaughters || daughter2 >= nDaughters)
361 return Const::doubleNaN;
362
363 PCmsLabTransform T;
364 ROOT::Math::PxPyPzEVector m = T.getBeamFourMomentum();
365 ROOT::Math::PxPyPzEVector p = particle->get4Vector();
366 ROOT::Math::PxPyPzEVector d1 = particle->getDaughter(daughter1)->get4Vector();
367 ROOT::Math::PxPyPzEVector d2 = particle->getDaughter(daughter2)->get4Vector();
368
369 ROOT::Math::PxPyPzEVector l;
370 l.SetPx(p.Py() * (d1.Pz() * d2.E() - d1.E() * d2.Pz()) + p.Pz() * (d1.E() * d2.Py() - d1.Py() * d2.E())
371 + p.E() * (d1.Py() * d2.Pz() - d1.Pz() * d2.Py()));
372 l.SetPy(p.Px() * (d1.E() * d2.Pz() - d1.Pz() * d2.E()) + p.Pz() * (d1.Px() * d2.E() - d1.E() * d2.Px())
373 + p.E() * (d1.Pz() * d2.Px() - d1.Px() * d2.Pz()));
374 l.SetPz(p.Px() * (d1.Py() * d2.E() - d1.E() * d2.Py()) + p.Py() * (d1.E() * d2.Px() - d1.Px() * d2.E())
375 + p.E() * (d1.Px() * d2.Py() - d1.Py() * d2.Px()));
376 l.SetE(-(p.Px() * (d1.Pz() * d2.Py() - d1.Py() * d2.Pz()) + p.Py() * (d1.Px() * d2.Pz() - d1.Pz() * d2.Px())
377 + p.Pz() * (d1.Py() * d2.Px() - d1.Px() * d2.Py())));
378
379 double m_times_p = m.Dot(p);
380 double m_times_l = m.Dot(l);
381 double m_times_d1 = m.Dot(d1);
382 double l_times_d1 = l.Dot(d1);
383 double d1_times_p = d1.Dot(p);
384 double m_abs = TMath::Sqrt(pow(m_times_p / p.M(), 2) - m.M2());
385 double d1_abs = TMath::Sqrt(pow(d1_times_p / p.M(), 2) - d1.M2());
386 double cos_phi = -m_times_l / (m_abs * TMath::Sqrt(-l.M2()));
387 double m_parallel_abs = m_abs * TMath::Sqrt(1 - cos_phi * cos_phi);
388 double m_parallel_times_d1 = m_times_p * d1_times_p / p.M2() + m_times_l * l_times_d1 / l.M2() - m_times_d1;
389
390 return TMath::ACos(-m_parallel_times_d1 / (m_parallel_abs * d1_abs));
391 }
392
393 double Constant(const Particle*, const std::vector<double>& constant)
394 {
395 return constant[0];
396 }
397
398
399
400 VARIABLE_GROUP("ParameterFunctions");
401 REGISTER_VARIABLE("NumberOfMCParticlesInEvent(pdgcode)", NumberOfMCParticlesInEvent, R"DOC(
402 [Eventbased] Returns number of MC Particles (including anti-particles) with the given pdgcode in the event.
403
404 Used in the FEI to determine to calculate reconstruction efficiencies.
405
406 The variable is event-based and does not need a valid particle pointer as input.)DOC");
407 REGISTER_VARIABLE("isAncestorOf(i, j, ...)", isAncestorOf, R"DOC(
408 Returns a positive integer if daughter at position particle->daughter(i)->daughter(j)... is an ancestor of the related MC particle, 0 otherwise.
409
410 Positive integer represents the number of steps needed to get from final MC daughter to ancestor.
411 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");
412 REGISTER_VARIABLE("hasAncestor(PDG, abs)", hasAncestor, R"DOC(
413
414 Returns a positive integer if an ancestor with the given PDG code is found, 0 otherwise.
415
416 The integer is the level where the ancestor was found, 1: first mother, 2: grandmother, etc.
417
418 Second argument is optional, 1 means that the sign of the PDG code is taken into account, default is 0.
419
420 If there is no MC relations found, -1 is returned. In case of nullptr particle, NaN is returned.)DOC");
421 REGISTER_VARIABLE("daughterInvariantMass(i, j, ...)", daughterInvariantMass, R"DOC(
422Returns invariant mass of the given daughter particles. E.g.:
423
424* daughterInvariantMass(0, 1) returns the invariant mass of the first and second daughter.
425* daughterInvariantMass(0, 1, 2) returns the invariant mass of the first, second and third daughter.
426
427Useful to identify intermediate resonances in a decay, which weren't reconstructed explicitly.
428
429Returns NaN if particle is nullptr or if the given daughter-index is out of bound (>= number of daughters).
430
431)DOC", "GeV/:math:`\\text{c}^2`");
432 MAKE_DEPRECATED("daughterInvariantMass", false, "light-2203-zeus", R"DOC(
433 The variable `daughterInvM` provides exactly the same functionality.)DOC");
434 REGISTER_VARIABLE("daughterMCInvariantMass(i, j, ...)", daughterMCInvariantMass, R"DOC(
435Returns true invariant mass of the given daughter particles, same behaviour as daughterInvariantMass variable.
436
437)DOC", "GeV/:math:`\\text{c}^2`");
438 REGISTER_VARIABLE("decayAngle(i)", particleDecayAngle, R"DOC(
439Angle in the mother's rest frame between the reverted CMS momentum vector and the direction of the i-th daughter
440
441)DOC", "rad");
442 REGISTER_VARIABLE("pointingAngle(i)", pointingAngle, R"DOC(
443Angle between i-th daughter's momentum vector and vector connecting production and decay vertex of i-th daughter.
444This makes only sense if the i-th daughter has itself daughter particles and therefore a properly defined vertex.
445
446)DOC", "rad");
447 REGISTER_VARIABLE("azimuthalAngleInDecayPlane(i, j)", azimuthalAngleInDecayPlane, R"DOC(
448Azimuthal angle of i-th daughter in decay plane towards projection of particle momentum into decay plane.
449
450First we define the following symbols:
451
452* P: four-momentum vector of decaying particle in whose decay plane the azimuthal angle is measured
453* M: "mother" of p, however not necessarily the direct mother but any higher state, here the CMS itself is chosen
454* D1: daughter for which the azimuthal angle is supposed to be calculated
455* D2: another daughter needed to span the decay plane
456* L: normal to the decay plane (four-component vector)
457
458L can be defined via the following relation:
459
460.. math:: L^{\sigma} = \delta^{\sigma\nu} \epsilon_{\mu\nu\alpha\beta} P^{\mu}D1^{\alpha}D2^{\beta}
461
462The azimuthal angle is given by
463
464.. math:: \phi \equiv \cos^{-1} \left(\frac{-\vec{M_{\parallel}} \cdot \vec{D1}}{|\vec{M_{\parallel}}| \cdot |\vec{D1}|}\right)
465
466For a frame independent formulation the three component vectors need to be written via invariant four-momentum vectors.
467
468.. math::
469
470 -\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}\\
471 |\vec{M_{\parallel}}| &= |\vec{M}| \sqrt{1 - \cos^2 \psi}\\
472 |\vec{M}| &= \sqrt{\frac{(M \cdot P)^2}{m^2_P} - m^2_M}\\
473 \cos \psi &= \frac{\vec{M} \cdot \vec{L}}{|\vec{M}| \cdot |\vec{L}|} = \frac{-M \cdot L}{|\vec{M}| \cdot \sqrt{-L^2}}\\
474 |\vec{D1}| &= \sqrt{\frac{(D1 \cdot P)^2}{m^2_P} - m^2_{D1}}
475
476)DOC", "rad");
477
478 REGISTER_VARIABLE("massDifference(i)", massDifference, "Difference in invariant masses of this particle and its i-th daughter\n\n",
479 "GeV/:math:`\\text{c}^2`");
480 REGISTER_VARIABLE("massDifferenceError(i)", massDifferenceError,
481 "Estimated uncertainty on difference in invariant masses of this particle and its i-th daughter\n\n", "GeV/:math:`\\text{c}^2`");
482 REGISTER_VARIABLE("massDifferenceSignificance(i)", massDifferenceSignificance,
483 "Signed significance of the deviation from the nominal mass difference of this particle and its i-th daughter [(massDiff - NOMINAL_MASS_DIFF)/ErrMassDiff]");
484
485 REGISTER_VARIABLE("constant(float i)", Constant, R"DOC(
486 Returns i.
487
488 Useful for debugging purposes and in conjunction with the formula meta-variable.)DOC");
489
490 REGISTER_VARIABLE("randomChoice(i, j, ...)", RandomChoice, R"DOC(
491 Returns random element of given numbers.
492
493 Useful for testing purposes.)DOC");
494
495
496 }
498}
DataType Angle(const B2Vector3< DataType > &q) const
The angle w.r.t.
Definition: B2Vector3.h:302
static const double doubleNaN
quiet_NaN
Definition: Const.h:703
@ c_PrimaryParticle
bit 0: Particle is primary particle.
Definition: MCParticle.h:47
const MCParticle * getMCParticle() const
Returns the pointer to the MCParticle object that was used to create this Particle (ParticleType == c...
Definition: Particle.cc:942
static const ReferenceFrame & GetCurrent()
Get current rest frame.
#define MAKE_DEPRECATED(name, make_fatal, version, description)
Registers a variable as deprecated.
Definition: Manager.h:443
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition: B2Vector3.h:516
Abstract base class for different kinds of events.