Belle II Software development
HelicityVariables.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/HelicityVariables.h>
11
12#include <analysis/variables/EventVariables.h>
13
14#include <analysis/dataobjects/Particle.h>
15
16#include <analysis/utility/ReferenceFrame.h>
17#include <analysis/VariableManager/Manager.h>
18
19#include <framework/gearbox/Const.h>
20
21#include <Math/Boost.h>
22#include <Math/Vector4D.h>
23#include <Math/VectorUtil.h>
24using namespace ROOT::Math;
25#include <cmath>
26
27namespace Belle2 {
32 namespace Variable {
33
34 double cosHelicityAngleMomentum(const Particle* part)
35 {
36
37 const auto& frame = ReferenceFrame::GetCurrent();
38 XYZVector motherBoost = frame.getMomentum(part).BoostToCM();
39 PxPyPzEVector motherMomentum = frame.getMomentum(part);
40 const auto& daughters = part -> getDaughters() ;
41
42 if (daughters.size() == 2) {
43
44 // Only for pi0 -> gamma gamma, gamma -> e+ e-
45 bool isOneConversion = false;
46 if (part->getPDGCode() == Const::pi0.getPDGCode()) {
47 for (auto& idaughter : daughters) {
48 // both daughter must be gamma
49 if (idaughter -> getPDGCode() != Const::photon.getPDGCode()) {
50 isOneConversion = false;
51 break;
52 }
53 // check if one of gammas has two daughters
54 if (idaughter -> getNDaughters() == 2) {
55 if (std::abs(idaughter -> getDaughters()[0]-> getPDGCode()) == Const::electron.getPDGCode()
56 && std::abs(idaughter -> getDaughters()[1]-> getPDGCode()) == Const::electron.getPDGCode()) { // e+ e-
57 isOneConversion = true;
58 }
59 }
60 }
61 }
62
63 if (isOneConversion) {
64 B2WARNING("cosHelicityAngleMomentum: Special treatment for pi0->gamma gamma, gamma -> e+ e-, is called. "
65 "This treatment is going to be deprecated and we recommend using ``cosHelicityAngleMomentumPi0Dalitz`` "
66 "If you find this message in another case, it must be a bug. Please report it to the software mailing list.");
67
68 //only for pi0 decay where one gamma converts
69 PxPyPzEVector pGamma;
70 for (auto& idaughter : daughters) {
71 if (idaughter -> getNDaughters() == 2) continue;
72 else pGamma = frame.getMomentum(idaughter);
73 }
74
75 pGamma = Boost(motherBoost) * pGamma;
76
77 return VectorUtil::CosTheta(motherMomentum, pGamma);
78
79 } else {
80 PxPyPzEVector pDaughter1 = frame.getMomentum(daughters[0]);
81 PxPyPzEVector pDaughter2 = frame.getMomentum(daughters[1]);
82
83 pDaughter1 = Boost(motherBoost) * pDaughter1;
84 pDaughter2 = Boost(motherBoost) * pDaughter2;
85
86 PxPyPzEVector p12 = pDaughter2 - pDaughter1;
87
88 return VectorUtil::CosTheta(motherMomentum, p12);
89 }
90
91 } else if (daughters.size() == 3) {
92
93 PxPyPzEVector pDaughter1 = frame.getMomentum(daughters[0]);
94 PxPyPzEVector pDaughter2 = frame.getMomentum(daughters[1]);
95 PxPyPzEVector pDaughter3 = frame.getMomentum(daughters[2]);
96
97 pDaughter1 = Boost(motherBoost) * pDaughter1;
98 pDaughter2 = Boost(motherBoost) * pDaughter2;
99 pDaughter3 = Boost(motherBoost) * pDaughter3;
100
101 XYZVector p12 = (pDaughter2 - pDaughter1).Vect();
102 XYZVector p13 = (pDaughter3 - pDaughter1).Vect();
103
104 XYZVector n = p12.Cross(p13);
105
106 return VectorUtil::CosTheta(motherMomentum, n);
107
108 } else return Const::doubleNaN;
109
110 }
111
112 double cosHelicityAngleMomentumPi0Dalitz(const Particle* part)
113 {
114
115 const auto& frame = ReferenceFrame::GetCurrent();
116 XYZVector motherBoost = frame.getMomentum(part).BoostToCM();
117 PxPyPzEVector motherMomentum = frame.getMomentum(part);
118 const auto& daughters = part -> getDaughters() ;
119
120
121 if (daughters.size() == 3) {
122
123 PxPyPzEVector pGamma;
124
125 for (auto& idaughter : daughters) {
126 if (std::abs(idaughter -> getPDGCode()) == Const::photon.getPDGCode()) {
127 pGamma = frame.getMomentum(idaughter);
128 break;
129 }
130 }
131 pGamma = Boost(motherBoost) * pGamma;
132
133 return VectorUtil::CosTheta(motherMomentum, pGamma);
134
135 } else if (daughters.size() == 2) { // only for pi0 -> gamma gamma, gamma -> e+ e-
136
137 PxPyPzEVector pGamma;
138
139 // both daughters must be gamma
140 if (daughters[0] -> getPDGCode() != Const::photon.getPDGCode() or
141 daughters[1] -> getPDGCode() != Const::photon.getPDGCode())
142 return Const::doubleNaN;
143
144 if (daughters[0] -> getNDaughters() == 2 and daughters[1] -> getNDaughters() == 0) {
145 if (std::abs(daughters[0] -> getDaughters()[0]-> getPDGCode()) == Const::electron.getPDGCode()
146 && std::abs(daughters[0] -> getDaughters()[1]-> getPDGCode()) == Const::electron.getPDGCode()) { // e+ e-
147 pGamma = frame.getMomentum(daughters[1]);
148 } else {
149 return Const::doubleNaN;
150 }
151 } else if (daughters[0] -> getNDaughters() == 0 and daughters[1] -> getNDaughters() == 2) {
152 if (std::abs(daughters[1] -> getDaughters()[0]-> getPDGCode()) == Const::electron.getPDGCode()
153 && std::abs(daughters[1] -> getDaughters()[1]-> getPDGCode()) == Const::electron.getPDGCode()) { // e+ e-
154 pGamma = frame.getMomentum(daughters[0]);
155 } else {
156 return Const::doubleNaN;
157 }
158 } else {
159 return Const::doubleNaN;
160 }
161
162 pGamma = Boost(motherBoost) * pGamma;
163
164 return VectorUtil::CosTheta(motherMomentum, pGamma);
165
166 } else return Const::doubleNaN;
167
168 }
169
170
171 double cosHelicityAngleBeamMomentum(const Particle* mother, const std::vector<double>& index)
172 {
173 if (index.size() != 1) {
174 B2FATAL("Wrong number of arguments for cosHelicityAngleIfCMSIsTheMother");
175 }
176
177 int idau = std::lround(index[0]);
178
179 const Particle* part = mother->getDaughter(idau);
180 if (!part) {
181 B2FATAL("Couldn't find the " << idau << "th daughter");
182 }
183
184 PxPyPzEVector beam4Vector(getBeamPx(nullptr), getBeamPy(nullptr), getBeamPz(nullptr), getBeamE(nullptr));
185 PxPyPzEVector part4Vector = part->get4Vector();
186 PxPyPzEVector mother4Vector = mother->get4Vector();
187
188 XYZVector motherBoost = mother4Vector.BoostToCM();
189
190 beam4Vector = Boost(motherBoost) * beam4Vector;
191 part4Vector = Boost(motherBoost) * part4Vector;
192
193 return - VectorUtil::CosTheta(part4Vector, beam4Vector);
194 }
195
196
197 double cosHelicityAngle(const Particle* mother, const std::vector<double>& indices)
198 {
199 if (indices.size() != 2) {
200 B2FATAL("Wrong number of arguments for cosHelicityAngleIfRefFrameIsTheDaughter: two are needed.");
201 }
202
203 int iDau = std::lround(indices[0]);
204 int iGrandDau = std::lround(indices[1]);
205
206 const Particle* daughter = mother->getDaughter(iDau);
207 if (!daughter)
208 B2FATAL("Couldn't find the " << iDau << "th daughter.");
209
210 const Particle* grandDaughter = daughter->getDaughter(iGrandDau);
211 if (!grandDaughter)
212 B2FATAL("Couldn't find the " << iGrandDau << "th daughter of the " << iDau << "th daughter.");
213
214 PxPyPzEVector mother4Vector = mother->get4Vector();
215 PxPyPzEVector daughter4Vector = daughter->get4Vector();
216 PxPyPzEVector grandDaughter4Vector = grandDaughter->get4Vector();
217
218 XYZVector daughterBoost = daughter4Vector.BoostToCM();
219
220 // We boost the momentum of the mother and of the granddaughter to the reference frame of the daughter.
221 grandDaughter4Vector = Boost(daughterBoost) * grandDaughter4Vector;
222 mother4Vector = Boost(daughterBoost) * mother4Vector;
223
224 return - VectorUtil::CosTheta(grandDaughter4Vector, mother4Vector);
225 }
226
227 double cosAcoplanarityAngle(const Particle* mother, const std::vector<double>& granddaughters)
228 {
229 if (granddaughters.size() != 2) {
230 B2FATAL("Wrong number of arguments for cosAcoplanarityAngleIfRefFrameIsTheMother: two are needed.");
231 }
232
233 if (mother->getNDaughters() != 2)
234 B2FATAL("cosAcoplanarityAngleIfRefFrameIsTheMother: this variable works only for two-body decays.");
235
236 int iGrandDau1 = std::lround(granddaughters[0]);
237 int iGrandDau2 = std::lround(granddaughters[1]);
238
239 const Particle* daughter1 = mother->getDaughter(0);
240 const Particle* daughter2 = mother->getDaughter(1);
241
242 const Particle* grandDaughter1 = daughter1->getDaughter(iGrandDau1);
243 if (!grandDaughter1)
244 B2FATAL("Couldn't find the " << iGrandDau1 << "th daughter of the first daughter.");
245
246 const Particle* grandDaughter2 = daughter2->getDaughter(iGrandDau2);
247 if (!grandDaughter2)
248 B2FATAL("Couldn't find the " << iGrandDau2 << "th daughter of the second daughter.");
249
250 PxPyPzEVector mother4Vector = mother->get4Vector();
251 PxPyPzEVector daughter4Vector1 = daughter1->get4Vector();
252 PxPyPzEVector daughter4Vector2 = daughter2->get4Vector();
253 PxPyPzEVector grandDaughter4Vector1 = grandDaughter1->get4Vector();
254 PxPyPzEVector grandDaughter4Vector2 = grandDaughter2->get4Vector();
255
256 XYZVector motherBoost = mother4Vector.BoostToCM();
257 XYZVector daughter1Boost = daughter4Vector1.BoostToCM();
258 XYZVector daughter2Boost = daughter4Vector2.BoostToCM();
259
260 // Boosting daughters to reference frame of the mother
261 daughter4Vector1 = Boost(motherBoost) * daughter4Vector1;
262 daughter4Vector2 = Boost(motherBoost) * daughter4Vector2;
263
264 // Boosting each granddaughter to reference frame of its mother
265 grandDaughter4Vector1 = Boost(daughter1Boost) * grandDaughter4Vector1;
266 grandDaughter4Vector2 = Boost(daughter2Boost) * grandDaughter4Vector2;
267
268 // We calculate the normal vectors of the decay two planes
269 XYZVector normalVector1 = daughter4Vector1.Vect().Cross(grandDaughter4Vector1.Vect());
270 XYZVector normalVector2 = daughter4Vector2.Vect().Cross(grandDaughter4Vector2.Vect());
271
272 return VectorUtil::CosTheta(normalVector1, normalVector2);
273 }
274
275 double cosHelicityAnglePrimary(const Particle* part)
276 {
277 return part->getCosHelicity();
278 }
279
280 double cosHelicityAngleDaughter(const Particle* part, const std::vector<double>& indices)
281 {
282 if ((indices.size() == 0) || (indices.size() > 2)) {
283 B2FATAL("Wrong number of arguments for cosHelicityAngleDaughter: one or two are needed.");
284 }
285
286 int iDaughter = std::lround(indices[0]);
287 int iGrandDaughter = 0;
288 if (indices.size() == 2) {
289 iGrandDaughter = std::lround(indices[1]);
290 }
291
292 return part->getCosHelicityDaughter(iDaughter, iGrandDaughter);
293 }
294
295 double acoplanarityAngle(const Particle* part)
296 {
297 return part->getAcoplanarity();
298 }
299
300
301 double cosHelicityAngleForQuasiTwoBodyDecay(const Particle* mother, const std::vector<double>& indices)
302 {
303 if (indices.size() != 2) {
304 B2FATAL("Wrong number of arguments for cosHelicityAngleForQuasiTwoBodyDecay: two are needed.");
305 }
306
307 if (mother->getNDaughters() != 3)
308 return Const::doubleNaN;
309
310 int iDau = std::lround(indices[0]);
311 int jDau = std::lround(indices[1]);
312
313 const Particle* iDaughter = mother->getDaughter(iDau);
314 if (!iDaughter)
315 return Const::doubleNaN;
316
317 const Particle* jDaughter = mother->getDaughter(jDau);
318 if (!jDaughter)
319 return Const::doubleNaN;
320
321 PxPyPzEVector mother4Vector = mother->get4Vector();
322 PxPyPzEVector iDaughter4Vector = iDaughter->get4Vector();
323 PxPyPzEVector jDaughter4Vector = jDaughter->get4Vector();
324
325 PxPyPzEVector resonance4Vector = iDaughter4Vector + jDaughter4Vector;
326 XYZVector resonanceBoost = resonance4Vector.BoostToCM();
327
328 iDaughter4Vector = Boost(resonanceBoost) * iDaughter4Vector;
329 mother4Vector = Boost(resonanceBoost) * mother4Vector;
330
331 return - VectorUtil::CosTheta(iDaughter4Vector, mother4Vector);
332 }
333
334 Manager::FunctionPtr momentaTripleProduct(const std::vector<std::string>& arguments)
335 {
336 if (arguments.size() != 3) {
337 B2FATAL("Wrong number of arguments for momentaTripleProduct: three (particles) are needed.");
338 }
339
340 // wrap with func and return it
341 auto func = [arguments](const Particle * mother) -> double {
342 auto iDau = arguments[0];
343 auto jDau = arguments[1];
344 auto kDau = arguments[2];
345
346 const Particle* iDaughter = mother->getParticleFromGeneralizedIndexString(iDau);
347 if (!iDaughter) B2FATAL("Couldn't find the " << iDau << "th daughter.");
348 const Particle* jDaughter = mother->getParticleFromGeneralizedIndexString(jDau);
349 if (!jDaughter) B2FATAL("Couldn't find the " << jDau << "th daughter.");
350 const Particle* kDaughter = mother->getParticleFromGeneralizedIndexString(kDau);
351 if (!kDaughter) B2FATAL("Couldn't find the " << kDau << "th daughter.");
352
353 PxPyPzEVector mother4Vector = mother->get4Vector();
354 PxPyPzEVector iDaughter4Vector = iDaughter->get4Vector();
355 PxPyPzEVector jDaughter4Vector = jDaughter->get4Vector();
356 PxPyPzEVector kDaughter4Vector = kDaughter->get4Vector();
357
358 XYZVector motherBoost = mother4Vector.BoostToCM();
359
360 // We boost the momenta of offspring to the reference frame of the mother.
361 iDaughter4Vector = Boost(motherBoost) * iDaughter4Vector;
362 jDaughter4Vector = Boost(motherBoost) * jDaughter4Vector;
363 kDaughter4Vector = Boost(motherBoost) * kDaughter4Vector;
364
365 // cross product: p_j x p_k
366 XYZVector jkDaughterCrossProduct = jDaughter4Vector.Vect().Cross(kDaughter4Vector.Vect());
367 // triple product: p_i * (p_j x p_k)
368 return iDaughter4Vector.Vect().Dot(jkDaughterCrossProduct) ;
369 };
370 return func;
371 }
372
373
374 VARIABLE_GROUP("Helicity variables");
375
376 REGISTER_VARIABLE("cosHelicityAngleMomentum", cosHelicityAngleMomentum, R"DOC(
377 If the given particle has two daughters: cosine of the angle between the line defined by the momentum difference of the two daughters
378 in the frame of the given particle (mother)
379 and the momentum of the given particle in the lab frame.
380
381 If the given particle has three daughters: cosine of the angle between the normal vector of the plane defined by
382 the momenta of the three daughters in the frame of the given particle (mother)
383 and the momentum of the given particle in the lab frame.
384
385 Otherwise, it returns 0.)DOC");
386
387 REGISTER_VARIABLE("cosHelicityAngleMomentumPi0Dalitz", cosHelicityAngleMomentumPi0Dalitz, R"DOC(
388 To be used for the decay :math:`\pi^0 \to e^+ e^- \gamma`:
389 cosine of the angle between the momentum of the gamma in the frame of the given particle (mother)
390 and the momentum of the given particle in the lab frame.
391
392 One can call the variable for the decay :math:`\pi^0 \to \gamma \gamma, \gamma \to e^+ e^-` as well.
393
394 Otherwise, it returns 0.)DOC");
395
396 REGISTER_VARIABLE("cosHelicityAngleBeamMomentum(i)", cosHelicityAngleBeamMomentum, R"DOC(
397 Cosine of the helicity angle of the :math:`i`-th daughter of the particle provided,
398 assuming that the mother of the provided particle corresponds to the centre-of-mass system, whose parameters are
399 automatically loaded by the function, given the accelerator's conditions.)DOC");
400
401 REGISTER_VARIABLE("cosHelicityAngle(i, j)", cosHelicityAngle, R"DOC(
402 Cosine of the helicity angle between the momentum of the selected granddaughter and the direction opposite to the momentum of the provided particle
403 in the reference frame of the selected daughter (:math:`\theta_1` and :math:`\theta_2` in the
404 `PDG <https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.030001>`_ 2018, p. 722).
405
406 This variable needs two integer arguments: the first one, ``i``, is the index of the daughter and the second one, ``j`` is the index of the granddaughter.
407
408 For example, in the decay :math:`B^0 \to \left(J/\psi \to \mu^+ \mu^-\right) \left(K^{*0} \to K^+ \pi^-\right)`,
409 if the provided particle is :math:`B^0` and the selected indices are (0, 0),
410 the variable will return the angle between the momentum of the :math:`\mu^+` and the direction opposite to the momentum of
411 the :math:`B^0`, both momenta in the rest frame of the :math:`J/\psi`.
412
413 This variable is needed for angular analyses of :math:`B`-meson decays into two vector particles.)DOC");
414
415 REGISTER_VARIABLE("cosAcoplanarityAngle(i, j)", cosAcoplanarityAngle, R"DOC(
416 Cosine of the acoplanarity angle (:math:`\Phi` in the `PDG Polarization Review <http://pdg.lbl.gov/2019/reviews/rpp2018-rev-b-decays-polarization.pdf>`_).
417 Given a two-body decay, the acoplanarity angle is defined as
418 the angle between the two decay planes in the reference frame of the mother.
419
420 We calculate the acoplanarity angle as the angle between the two
421 normal vectors of the decay planes. Each normal vector is the cross product of the momentum of one daughter (in the frame of the mother) and the
422 momentum of one of the granddaughters (in the reference frame of the daughter).
423
424 This variable needs two integer arguments: the first one, ``i`` is the index of the first granddaughter, and the second one, ``j`` the index of the
425 second granddaughter.
426
427 For example, in the decay :math:`B^0 \to \left(J/\psi \to \mu^+ \mu^-\right) \left(K^{*0} \to K^+ \pi^-\right)`, if the provided particle is :math:`B^0` and the selected indices are (0, 0),
428 the variable will return the acoplanarity using the :math:`\mu^+` and the :math:`K^+` granddaughters.)DOC");
429
430 REGISTER_VARIABLE("cosHelicityAnglePrimary", cosHelicityAnglePrimary, R"DOC(
431 Cosine of the helicity angle (see``Particle::getCosHelicity``) assuming the center of mass system as mother rest frame.
432 See `PDG Polarization Review <http://pdg.lbl.gov/2019/reviews/rpp2018-rev-b-decays-polarization.pdf>`_ for the definition of the helicity angle.)DOC");
433
434 REGISTER_VARIABLE("cosHelicityAngleDaughter(i [, j] )", cosHelicityAngleDaughter, R"DOC(
435 Cosine of the helicity angle of the i-th daughter (see ``Particle::getCosHelicityDaughter``).
436 The optional second argument is the index of the granddaughter that defines the angle, default is 0.
437
438 For example, in the decay: :math:`B^0 \to \left(J/\psi \to \mu^+ \mu^-\right) \left(K^{*0} \to K^+ \pi^-\right)`, if the provided particle is :math:`B^0` and the selected index is 0,
439 the variable will return the helicity angle of the :math:`\mu^+`.
440 If the selected index is 1 the variable will return the helicity angle of the :math:`K^+` (defined via the rest frame of the :math:`K^{*0}`).
441 In rare cases if one wanted the helicity angle of the second granddaughter, indices 1,1 would return the helicity angle of the :math:`\pi^-`).
442
443 See `PDG Polarization Review <http://pdg.lbl.gov/2019/reviews/rpp2018-rev-b-decays-polarization.pdf>`_ for the definition of the helicity angle.)DOC");
444
445 REGISTER_VARIABLE("acoplanarityAngle", acoplanarityAngle, R"DOC(
446Acoplanarity angle (see ``Particle::getAcoplanarity``) assuming a two body decay of the particle and its daughters.
447See `PDG Polarization Review <http://pdg.lbl.gov/2019/reviews/rpp2018-rev-b-decays-polarization.pdf>`_ for the definition of the acoplanarity angle.
448
449)DOC", "rad");
450
451 REGISTER_VARIABLE("cosHelicityAngleForQuasiTwoBodyDecay(i, j)", cosHelicityAngleForQuasiTwoBodyDecay, R"DOC(
452 Cosine of the helicity angle between the momentum of the provided particle and the momentum of the first selected
453 daughter (i-th) in the reference frame of the sum of two selected daughters (i-th + j-th).
454
455 The variable is supposed to be used for the analysis of a quasi-two-body decay. The number of daughters of the given
456 particle must be three. Otherwise, the variable returns NaN.
457
458 For example, in the decay :math:`\bar{B}^0 \to D^+ K^- K^{*0}`, if the provided particle is :math:`\bar{B}^0` and
459 the selected indices are (1, 2), the variable will return the angle between the momentum of the :math:`\bar{B}^0`
460 and the momentum of the :math:`K^-`, both momenta in the rest frame of the :math:`K^- K^{*0}`.)DOC");
461
462 REGISTER_METAVARIABLE("momentaTripleProduct(i,j,k)", momentaTripleProduct, R"DOC(
463a triple-product of three momenta of offspring in the mother rest frame: :math:`C_T=\vec{p}_i\cdot(\vec{p}_j\times\vec{p}_k)`. For examples,
464In a four-body decay M->D1D2D3D4, momentaTripleProduct(0,1,2) returns CT using the momenta of D1D2D3 particles.
465In other decays involving secondary decay, e.g. for M->(R->D1D2)D3D4, momentaTripleProduct(0:0,1,2) returns C_T using momenta of D1D3D4 particles.
466)DOC", Manager::VariableDataType::c_double);
467
468 }
470}
int getPDGCode() const
PDG code.
Definition Const.h:473
static const ParticleType pi0
neutral pion particle
Definition Const.h:674
static const double doubleNaN
quiet_NaN
Definition Const.h:703
static const ParticleType photon
photon particle
Definition Const.h:673
static const ChargedStable electron
electron particle
Definition Const.h:659
static const ReferenceFrame & GetCurrent()
Get current rest frame.
std::function< VarVariant(const Particle *)> FunctionPtr
functions stored take a const Particle* and return VarVariant.
Definition Manager.h:112
Abstract base class for different kinds of events.