10#include <analysis/variables/HelicityVariables.h>
12#include <analysis/variables/EventVariables.h>
14#include <analysis/dataobjects/Particle.h>
16#include <analysis/utility/ReferenceFrame.h>
17#include <analysis/VariableManager/Manager.h>
19#include <framework/gearbox/Const.h>
21#include <Math/Boost.h>
22#include <Math/Vector4D.h>
23#include <Math/VectorUtil.h>
24using namespace ROOT::Math;
34 double cosHelicityAngleMomentum(
const Particle* part)
38 XYZVector motherBoost = frame.getMomentum(part).BoostToCM();
39 PxPyPzEVector motherMomentum = frame.getMomentum(part);
40 const auto& daughters = part -> getDaughters() ;
42 if (daughters.size() == 2) {
45 bool isOneConversion =
false;
47 for (
auto& idaughter : daughters) {
49 if (idaughter -> getPDGCode() !=
Const::photon.getPDGCode()) {
50 isOneConversion =
false;
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()) {
57 isOneConversion =
true;
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.");
70 for (
auto& idaughter : daughters) {
71 if (idaughter -> getNDaughters() == 2)
continue;
72 else pGamma = frame.getMomentum(idaughter);
75 pGamma = Boost(motherBoost) * pGamma;
77 return VectorUtil::CosTheta(motherMomentum, pGamma);
80 PxPyPzEVector pDaughter1 = frame.getMomentum(daughters[0]);
81 PxPyPzEVector pDaughter2 = frame.getMomentum(daughters[1]);
83 pDaughter1 = Boost(motherBoost) * pDaughter1;
84 pDaughter2 = Boost(motherBoost) * pDaughter2;
86 PxPyPzEVector p12 = pDaughter2 - pDaughter1;
88 return VectorUtil::CosTheta(motherMomentum, p12);
91 }
else if (daughters.size() == 3) {
93 PxPyPzEVector pDaughter1 = frame.getMomentum(daughters[0]);
94 PxPyPzEVector pDaughter2 = frame.getMomentum(daughters[1]);
95 PxPyPzEVector pDaughter3 = frame.getMomentum(daughters[2]);
97 pDaughter1 = Boost(motherBoost) * pDaughter1;
98 pDaughter2 = Boost(motherBoost) * pDaughter2;
99 pDaughter3 = Boost(motherBoost) * pDaughter3;
101 XYZVector p12 = (pDaughter2 - pDaughter1).Vect();
102 XYZVector p13 = (pDaughter3 - pDaughter1).Vect();
104 XYZVector n = p12.Cross(p13);
106 return VectorUtil::CosTheta(motherMomentum, n);
112 double cosHelicityAngleMomentumPi0Dalitz(
const Particle* part)
116 XYZVector motherBoost = frame.getMomentum(part).BoostToCM();
117 PxPyPzEVector motherMomentum = frame.getMomentum(part);
118 const auto& daughters = part -> getDaughters() ;
121 if (daughters.size() == 3) {
123 PxPyPzEVector pGamma;
125 for (
auto& idaughter : daughters) {
126 if (std::abs(idaughter -> getPDGCode()) ==
Const::photon.getPDGCode()) {
127 pGamma = frame.getMomentum(idaughter);
131 pGamma = Boost(motherBoost) * pGamma;
133 return VectorUtil::CosTheta(motherMomentum, pGamma);
135 }
else if (daughters.size() == 2) {
137 PxPyPzEVector pGamma;
140 if (daughters[0] -> getPDGCode() !=
Const::photon.getPDGCode() or
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()) {
147 pGamma = frame.getMomentum(daughters[1]);
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()) {
154 pGamma = frame.getMomentum(daughters[0]);
162 pGamma = Boost(motherBoost) * pGamma;
164 return VectorUtil::CosTheta(motherMomentum, pGamma);
171 double cosHelicityAngleBeamMomentum(
const Particle* mother,
const std::vector<double>& index)
173 if (index.size() != 1) {
174 B2FATAL(
"Wrong number of arguments for cosHelicityAngleIfCMSIsTheMother");
177 int idau = std::lround(index[0]);
179 const Particle* part = mother->getDaughter(idau);
181 B2FATAL(
"Couldn't find the " << idau <<
"th daughter");
184 PxPyPzEVector beam4Vector(getBeamPx(
nullptr), getBeamPy(
nullptr), getBeamPz(
nullptr), getBeamE(
nullptr));
185 PxPyPzEVector part4Vector = part->get4Vector();
186 PxPyPzEVector mother4Vector = mother->get4Vector();
188 XYZVector motherBoost = mother4Vector.BoostToCM();
190 beam4Vector = Boost(motherBoost) * beam4Vector;
191 part4Vector = Boost(motherBoost) * part4Vector;
193 return - VectorUtil::CosTheta(part4Vector, beam4Vector);
197 double cosHelicityAngle(
const Particle* mother,
const std::vector<double>& indices)
199 if (indices.size() != 2) {
200 B2FATAL(
"Wrong number of arguments for cosHelicityAngleIfRefFrameIsTheDaughter: two are needed.");
203 int iDau = std::lround(indices[0]);
204 int iGrandDau = std::lround(indices[1]);
206 const Particle* daughter = mother->getDaughter(iDau);
208 B2FATAL(
"Couldn't find the " << iDau <<
"th daughter.");
210 const Particle* grandDaughter = daughter->getDaughter(iGrandDau);
212 B2FATAL(
"Couldn't find the " << iGrandDau <<
"th daughter of the " << iDau <<
"th daughter.");
214 PxPyPzEVector mother4Vector = mother->get4Vector();
215 PxPyPzEVector daughter4Vector = daughter->get4Vector();
216 PxPyPzEVector grandDaughter4Vector = grandDaughter->get4Vector();
218 XYZVector daughterBoost = daughter4Vector.BoostToCM();
221 grandDaughter4Vector = Boost(daughterBoost) * grandDaughter4Vector;
222 mother4Vector = Boost(daughterBoost) * mother4Vector;
224 return - VectorUtil::CosTheta(grandDaughter4Vector, mother4Vector);
227 double cosAcoplanarityAngle(
const Particle* mother,
const std::vector<double>& granddaughters)
229 if (granddaughters.size() != 2) {
230 B2FATAL(
"Wrong number of arguments for cosAcoplanarityAngleIfRefFrameIsTheMother: two are needed.");
233 if (mother->getNDaughters() != 2)
234 B2FATAL(
"cosAcoplanarityAngleIfRefFrameIsTheMother: this variable works only for two-body decays.");
236 int iGrandDau1 = std::lround(granddaughters[0]);
237 int iGrandDau2 = std::lround(granddaughters[1]);
239 const Particle* daughter1 = mother->getDaughter(0);
240 const Particle* daughter2 = mother->getDaughter(1);
242 const Particle* grandDaughter1 = daughter1->getDaughter(iGrandDau1);
244 B2FATAL(
"Couldn't find the " << iGrandDau1 <<
"th daughter of the first daughter.");
246 const Particle* grandDaughter2 = daughter2->getDaughter(iGrandDau2);
248 B2FATAL(
"Couldn't find the " << iGrandDau2 <<
"th daughter of the second daughter.");
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();
256 XYZVector motherBoost = mother4Vector.BoostToCM();
257 XYZVector daughter1Boost = daughter4Vector1.BoostToCM();
258 XYZVector daughter2Boost = daughter4Vector2.BoostToCM();
261 daughter4Vector1 = Boost(motherBoost) * daughter4Vector1;
262 daughter4Vector2 = Boost(motherBoost) * daughter4Vector2;
265 grandDaughter4Vector1 = Boost(daughter1Boost) * grandDaughter4Vector1;
266 grandDaughter4Vector2 = Boost(daughter2Boost) * grandDaughter4Vector2;
269 XYZVector normalVector1 = daughter4Vector1.Vect().Cross(grandDaughter4Vector1.Vect());
270 XYZVector normalVector2 = daughter4Vector2.Vect().Cross(grandDaughter4Vector2.Vect());
272 return VectorUtil::CosTheta(normalVector1, normalVector2);
275 double cosHelicityAnglePrimary(
const Particle* part)
277 return part->getCosHelicity();
280 double cosHelicityAngleDaughter(
const Particle* part,
const std::vector<double>& indices)
282 if ((indices.size() == 0) || (indices.size() > 2)) {
283 B2FATAL(
"Wrong number of arguments for cosHelicityAngleDaughter: one or two are needed.");
286 int iDaughter = std::lround(indices[0]);
287 int iGrandDaughter = 0;
288 if (indices.size() == 2) {
289 iGrandDaughter = std::lround(indices[1]);
292 return part->getCosHelicityDaughter(iDaughter, iGrandDaughter);
295 double acoplanarityAngle(
const Particle* part)
297 return part->getAcoplanarity();
301 double cosHelicityAngleForQuasiTwoBodyDecay(
const Particle* mother,
const std::vector<double>& indices)
303 if (indices.size() != 2) {
304 B2FATAL(
"Wrong number of arguments for cosHelicityAngleForQuasiTwoBodyDecay: two are needed.");
307 if (mother->getNDaughters() != 3)
310 int iDau = std::lround(indices[0]);
311 int jDau = std::lround(indices[1]);
313 const Particle* iDaughter = mother->getDaughter(iDau);
317 const Particle* jDaughter = mother->getDaughter(jDau);
321 PxPyPzEVector mother4Vector = mother->get4Vector();
322 PxPyPzEVector iDaughter4Vector = iDaughter->get4Vector();
323 PxPyPzEVector jDaughter4Vector = jDaughter->get4Vector();
325 PxPyPzEVector resonance4Vector = iDaughter4Vector + jDaughter4Vector;
326 XYZVector resonanceBoost = resonance4Vector.BoostToCM();
328 iDaughter4Vector = Boost(resonanceBoost) * iDaughter4Vector;
329 mother4Vector = Boost(resonanceBoost) * mother4Vector;
331 return - VectorUtil::CosTheta(iDaughter4Vector, mother4Vector);
336 if (arguments.size() != 3) {
337 B2FATAL(
"Wrong number of arguments for momentaTripleProduct: three (particles) are needed.");
341 auto func = [arguments](
const Particle * mother) ->
double {
342 auto iDau = arguments[0];
343 auto jDau = arguments[1];
344 auto kDau = arguments[2];
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.");
353 PxPyPzEVector mother4Vector = mother->get4Vector();
354 PxPyPzEVector iDaughter4Vector = iDaughter->get4Vector();
355 PxPyPzEVector jDaughter4Vector = jDaughter->get4Vector();
356 PxPyPzEVector kDaughter4Vector = kDaughter->get4Vector();
358 XYZVector motherBoost = mother4Vector.BoostToCM();
361 iDaughter4Vector = Boost(motherBoost) * iDaughter4Vector;
362 jDaughter4Vector = Boost(motherBoost) * jDaughter4Vector;
363 kDaughter4Vector = Boost(motherBoost) * kDaughter4Vector;
366 XYZVector jkDaughterCrossProduct = jDaughter4Vector.Vect().Cross(kDaughter4Vector.Vect());
368 return iDaughter4Vector.Vect().Dot(jkDaughterCrossProduct) ;
374 VARIABLE_GROUP(
"Helicity variables");
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.
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.
385 Otherwise, it returns 0.)DOC");
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.
392 One can call the variable for the decay :math:`\pi^0 \to \gamma \gamma, \gamma \to e^+ e^-` as well.
394 Otherwise, it returns 0.)DOC");
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");
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).
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.
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`.
413 This variable is needed for angular analyses of :math:`B`-meson decays into two vector particles.)DOC");
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.
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).
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.
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");
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");
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.
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^-`).
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");
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.
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).
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.
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");
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);
int getPDGCode() const
PDG code.
static const ParticleType pi0
neutral pion particle
static const double doubleNaN
quiet_NaN
static const ParticleType photon
photon particle
static const ChargedStable electron
electron particle
static const ReferenceFrame & GetCurrent()
Get current rest frame.
std::function< VarVariant(const Particle *)> FunctionPtr
functions stored take a const Particle* and return VarVariant.
Abstract base class for different kinds of events.