Belle II Software  release-05-01-25
HelicityVariables.cc
1 /**************************************************************************
2 * BASF2 (Belle Analysis Framework 2) *
3 * Copyright(C) 2019 - Belle II Collaboration *
4 * *
5 * Author: The Belle II Collaboration *
6 * Contributors: Phillip Urquijo, Fernando Abudinen, Michael de Nuccio *
7 * *
8 * This software is provided "as is" without any warranty. *
9 **************************************************************************/
10 
11 
12 #include <analysis/variables/HelicityVariables.h>
13 #include <analysis/variables/EventVariables.h>
14 
15 #include <analysis/dataobjects/Particle.h>
16 
17 #include <analysis/utility/ReferenceFrame.h>
18 
19 #include <framework/utilities/Conversion.h>
20 
21 #include <TLorentzVector.h>
22 #include <TVector3.h>
23 #include <cmath>
24 
25 namespace Belle2 {
30  namespace Variable {
31 
32  double cosHelicityAngleMomentum(const Particle* part)
33  {
34 
35  const auto& frame = ReferenceFrame::GetCurrent();
36  TVector3 motherBoost = - frame.getMomentum(part).BoostVector();
37  TVector3 motherMomentum = frame.getMomentum(part).Vect();
38  const auto& daughters = part -> getDaughters() ;
39 
40  if (daughters.size() == 2) {
41 
42  bool isOneConversion = false;
43 
44  for (auto& idaughter : daughters) {
45  if (idaughter -> getNDaughters() == 2) {
46  if (std::abs(idaughter -> getDaughters()[0]-> getPDGCode()) == 11
47  && std::abs(idaughter -> getDaughters()[1]-> getPDGCode()) == 11) {
48  isOneConversion = true;
49  }
50  }
51  }
52 
53  if (isOneConversion) {
54  //only for pi0 decay where one gamma converts
55 
56  TLorentzVector pGamma;
57 
58  for (auto& idaughter : daughters) {
59  if (idaughter -> getNDaughters() == 2) continue;
60  else pGamma = frame.getMomentum(idaughter);
61  }
62 
63  pGamma.Boost(motherBoost);
64 
65  return std::cos(motherMomentum.Angle(pGamma.Vect()));
66 
67  } else {
68  TLorentzVector pDaughter1 = frame.getMomentum(daughters[0]);
69  TLorentzVector pDaughter2 = frame.getMomentum(daughters[1]);
70 
71  pDaughter1.Boost(motherBoost);
72  pDaughter2.Boost(motherBoost);
73 
74  TVector3 p12 = (pDaughter2 - pDaughter1).Vect();
75 
76  return std::cos(motherMomentum.Angle(p12));
77  }
78 
79  } else if (daughters.size() == 3) {
80 
81  TLorentzVector pDaughter1 = frame.getMomentum(daughters[0]);
82  TLorentzVector pDaughter2 = frame.getMomentum(daughters[1]);
83  TLorentzVector pDaughter3 = frame.getMomentum(daughters[2]);
84 
85  pDaughter1.Boost(motherBoost);
86  pDaughter2.Boost(motherBoost);
87  pDaughter3.Boost(motherBoost);
88 
89  TVector3 p12 = (pDaughter2 - pDaughter1).Vect();
90  TVector3 p13 = (pDaughter3 - pDaughter1).Vect();
91 
92  TVector3 n = p12.Cross(p13);
93 
94  return std::cos(motherMomentum.Angle(n));
95 
96  } else return std::numeric_limits<float>::quiet_NaN();
97 
98  }
99 
100  double cosHelicityAngleMomentumPi0Dalitz(const Particle* part)
101  {
102 
103  const auto& frame = ReferenceFrame::GetCurrent();
104  TVector3 motherBoost = - frame.getMomentum(part).BoostVector();
105  TVector3 motherMomentum = frame.getMomentum(part).Vect();
106  const auto& daughters = part -> getDaughters() ;
107 
108 
109  if (daughters.size() == 3) {
110 
111  TLorentzVector pGamma;
112 
113  for (auto& idaughter : daughters) {
114  if (std::abs(idaughter -> getPDGCode()) == 22) pGamma = frame.getMomentum(idaughter);
115  }
116 
117  pGamma.Boost(motherBoost);
118 
119  return std::cos(motherMomentum.Angle(pGamma.Vect()));
120 
121  } else return std::numeric_limits<float>::quiet_NaN();
122 
123  }
124 
125 
126  Manager::FunctionPtr cosHelicityAngleBeamMomentum(const std::vector<std::string>& arguments)
127  {
128  int idau = 0;
129  if (arguments.size() == 1) {
130  try {
131  idau = Belle2::convertString<int>(arguments[0]);
132  } catch (std::invalid_argument&) {
133  B2FATAL("The argument of cosHelicityAngleWrtCMSFrame must be an integer!");
134  }
135  } else {
136  B2FATAL("Wrong number of arguments for cosHelicityAngleIfCMSIsTheMother");
137  }
138  auto func = [idau](const Particle * mother) -> double {
139  const Particle* part = mother->getDaughter(idau);
140  if (!part)
141  {
142  B2FATAL("Couldn't find the " << idau << "th daughter");
143  }
144 
145  TLorentzVector beam4Vector(getBeamPx(nullptr), getBeamPy(nullptr), getBeamPz(nullptr), getBeamE(nullptr));
146  TLorentzVector part4Vector = part->get4Vector();
147  TLorentzVector mother4Vector = mother->get4Vector();
148 
149  TVector3 motherBoost = -(mother4Vector.BoostVector());
150 
151  beam4Vector.Boost(motherBoost);
152  part4Vector.Boost(motherBoost);
153 
154  return - part4Vector.Vect().Dot(beam4Vector.Vect()) / part4Vector.Vect().Mag() / beam4Vector.Vect().Mag();
155  };
156  return func;
157  }
158 
159 
160  Manager::FunctionPtr cosHelicityAngle(const std::vector<std::string>& arguments)
161  {
162  int iDau = 0;
163  int iGrandDau = 0;
164  if (arguments.size() == 2) {
165  try {
166  iDau = Belle2::convertString<int>(arguments[0]);
167  iGrandDau = Belle2::convertString<int>(arguments[1]);
168  } catch (std::invalid_argument&) {
169  B2FATAL("The two arguments of cosHelicityAngleIfRefFrameIsTheDaughter must be integers!");
170  }
171  } else {
172  B2FATAL("Wrong number of arguments for cosHelicityAngleIfRefFrameIsTheDaughter: two are needed.");
173  }
174  auto func = [iDau, iGrandDau](const Particle * mother) -> double {
175 
176  const Particle* daughter = mother->getDaughter(iDau);
177  if (!daughter)
178  B2FATAL("Couldn't find the " << iDau << "th daughter.");
179 
180  const Particle* grandDaughter = daughter->getDaughter(iGrandDau);
181  if (!grandDaughter)
182  B2FATAL("Couldn't find the " << iGrandDau << "th daughter of the " << iDau << "th daughter.");
183 
184 
185  TLorentzVector mother4Vector = mother->get4Vector();
186  TLorentzVector daughter4Vector = daughter->get4Vector();
187  TLorentzVector grandDaughter4Vector = grandDaughter->get4Vector();
188 
189  TVector3 daughterBoost = -(daughter4Vector.BoostVector());
190 
191  // We boost the momentum of the mother and of the granddaughter to the reference frame of the daughter.
192  grandDaughter4Vector.Boost(daughterBoost);
193  mother4Vector.Boost(daughterBoost);
194 
195  return - grandDaughter4Vector.Vect().Dot(mother4Vector.Vect()) / grandDaughter4Vector.Vect().Mag() / mother4Vector.Vect().Mag();
196 
197  };
198  return func;
199  }
200 
201  Manager::FunctionPtr cosAcoplanarityAngle(const std::vector<std::string>& arguments)
202  {
203  int iGrandDau1 = 0;
204  int iGrandDau2 = 0;
205  if (arguments.size() == 2) {
206  try {
207  iGrandDau1 = Belle2::convertString<int>(arguments[0]);
208  iGrandDau2 = Belle2::convertString<int>(arguments[1]);
209  } catch (std::invalid_argument&) {
210  B2FATAL("The two arguments of cosAcoplanarityAngleIfRefFrameIsTheMother must be integers!");
211  }
212  } else {
213  B2FATAL("Wrong number of arguments for cosAcoplanarityAngleIfRefFrameIsTheMother: two are needed.");
214  }
215  auto func = [iGrandDau1, iGrandDau2](const Particle * mother) -> double {
216 
217  if (mother->getNDaughters() != 2)
218  B2FATAL("cosAcoplanarityAngleIfRefFrameIsTheMother: this variable works only for two-body decays.");
219 
220  const Particle* daughter1 = mother-> getDaughter(0);
221  const Particle* daughter2 = mother-> getDaughter(1);
222 
223  const Particle* grandDaughter1 = daughter1 -> getDaughter(iGrandDau1);
224  if (!grandDaughter1)
225  B2FATAL("Couldn't find the " << iGrandDau1 << "th daughter of the first daughter.");
226 
227  const Particle* grandDaughter2 = daughter2 -> getDaughter(iGrandDau2);
228  if (!grandDaughter2)
229  B2FATAL("Couldn't find the " << iGrandDau2 << "th daughter of the second daughter.");
230 
231  TLorentzVector mother4Vector = mother->get4Vector();
232  TLorentzVector daughter4Vector1 = daughter1->get4Vector();
233  TLorentzVector daughter4Vector2 = daughter2->get4Vector();
234  TLorentzVector grandDaughter4Vector1 = grandDaughter1->get4Vector();
235  TLorentzVector grandDaughter4Vector2 = grandDaughter2->get4Vector();
236 
237  TVector3 motherBoost = -(mother4Vector.BoostVector());
238  TVector3 daughter1Boost = -(daughter4Vector1.BoostVector());
239  TVector3 daughter2Boost = -(daughter4Vector2.BoostVector());
240 
241  // Boosting daughters to reference frame of the mother
242  daughter4Vector1.Boost(motherBoost);
243  daughter4Vector2.Boost(motherBoost);
244 
245  // Boosting each granddaughter to reference frame of its mother
246  grandDaughter4Vector1.Boost(daughter1Boost);
247  grandDaughter4Vector2.Boost(daughter2Boost);
248 
249  // We calculate the normal vectors of the decay two planes
250  TVector3 normalVector1 = daughter4Vector1.Vect().Cross(grandDaughter4Vector1.Vect());
251  TVector3 normalVector2 = daughter4Vector2.Vect().Cross(grandDaughter4Vector2.Vect());
252 
253  return std::cos(normalVector1.Angle(normalVector2));
254 
255  };
256  return func;
257  }
258 
259 
260  double cosHelicityAnglePrimary(const Particle* part)
261  {
262  return part->getCosHelicity();
263  }
264 
265  Manager::FunctionPtr cosHelicityAngleDaughter(const std::vector<std::string>& arguments)
266  {
267  int iDaughter = 0;
268  int iGrandDaughter = 0;
269  if ((arguments.size() == 0) || (arguments.size() > 2)) {
270  B2FATAL("Wrong number of arguments for cosHelicityAngleDaughter: one or two are needed.");
271  }
272  try {
273  iDaughter = Belle2::convertString<int>(arguments[0]);
274  if (arguments.size() == 2) {
275  iGrandDaughter = Belle2::convertString<int>(arguments[1]);
276  }
277  } catch (std::invalid_argument&) {
278  B2FATAL("The arguments of cosHelicityAngleDaughter must be integers!");
279  }
280 
281  auto func = [iDaughter, iGrandDaughter](const Particle * part) -> double {
282  return part->getCosHelicityDaughter(iDaughter, iGrandDaughter);
283  };
284  return func;
285  }
286 
287  double acoplanarityAngle(const Particle* part)
288  {
289  return part->getAcoplanarity();
290  }
291 
292 
293  VARIABLE_GROUP("Helicity variables");
294 
295  REGISTER_VARIABLE("cosHelicityAngleMomentum",
296  cosHelicityAngleMomentum,
297  R"DOC(
298  If the given particle has two daughters: cosine of the angle between the line defined by the momentum difference of the two daughters
299  in the frame of the given particle (mother)
300  and the momentum of the given particle in the lab frame.
301 
302  If the given particle has three daughters: cosine of the angle between the normal vector of the plane defined by
303  the momenta of the three daughters in the frame of the given particle (mother)
304  and the momentum of the given particle in the lab frame.
305 
306  Otherwise, it returns 0.)DOC");
307 
308  REGISTER_VARIABLE("cosHelicityAngleMomentumPi0Dalitz",
309  cosHelicityAngleMomentumPi0Dalitz,
310  R"DOC(
311  To be used for the decay :math:`\pi^0 \to e^+ e^- \gamma`:
312  cosine of the angle between the momentum of the gamma in the frame of the given particle (mother)
313  and the momentum of the given particle in the lab frame.
314 
315  Otherwise, it returns 0.)DOC");
316 
317  REGISTER_VARIABLE("cosHelicityAngleBeamMomentum(i)", cosHelicityAngleBeamMomentum,
318  R"DOC(
319  Cosine of the helicity angle of the :math:`i`-th daughter of the particle provided,
320  assuming that the mother of the provided particle corresponds to the centre-of-mass system, whose parameters are
321  automatically loaded by the function, given the accelerator's conditions.)DOC");
322 
323  REGISTER_VARIABLE("cosHelicityAngle(i, j)", cosHelicityAngle,
324  R"DOC(
325  Cosine of the helicity angle between the momentum of the provided particle and the momentum of the selected granddaughter
326  in the reference frame of the selected daughter (:math:`\theta_1` and :math:`\theta_2` in the
327  `PDG <https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.030001>`_ 2018, p. 722).
328 
329  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.
330 
331  For example, in the decay :math:`B^0 \to \left(J/\psi \to \mu^+ \mu^-\right) \left(K^{*0} \to K^+ \pi^-\right)`,
332  if the provided particle is :math:`B^0` and the selected indices are (0, 0),
333  the variable will return the angle between the momentum of the :math:`B^0` and the momentum of the :math:`\mu^+`,
334  both momenta in the rest frame of the :math:`J/\psi`.
335 
336  This variable is needed for angular analyses of :math:`B`-meson decays into two vector particles.)DOC");
337 
338  REGISTER_VARIABLE("cosAcoplanarityAngle(i, j)", cosAcoplanarityAngle,
339  R"DOC(
340  Cosine of the acoplanarity angle (:math:`\Phi` in the `PDG Polarization Review <http://pdg.lbl.gov/2019/reviews/rpp2018-rev-b-decays-polarization.pdf>`_).
341  Given a two-body decay, the acoplanarity angle is defined as
342  the angle between the two decay planes in the reference frame of the mother.
343 
344  We calculate the acoplanarity angle as the angle between the two
345  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
346  momentum of one of the granddaughters (in the reference frame of the daughter).
347 
348  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
349  second granddaughter.
350 
351  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),
352  the variable will return the acoplanarity using the :math:`\mu^+` and the :math:`K^+` granddaughters.)DOC");
353 
354  REGISTER_VARIABLE("cosHelicityAnglePrimary", cosHelicityAnglePrimary,
355  R"DOC(
356  Cosine of the helicity angle (see``Particle::getCosHelicity``) assuming the center of mass system as mother rest frame.
357  See `PDG Polarization Review <http://pdg.lbl.gov/2019/reviews/rpp2018-rev-b-decays-polarization.pdf>`_ for the definition of the helicity angle.)DOC");
358 
359  REGISTER_VARIABLE("cosHelicityAngleDaughter(i [, j] )", cosHelicityAngleDaughter,
360  R"DOC(
361  Cosine of the helicity angle of the i-th daughter (see ``Particle::getCosHelicityDaughter``).
362  The optional second argument is the index of the granddaughter that defines the angle, default is 0.
363 
364  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,
365  the variable will return the helicity angle of the :math:`\mu^+`.
366  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}`).
367  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^-`).
368 
369  See `PDG Polarization Review <http://pdg.lbl.gov/2019/reviews/rpp2018-rev-b-decays-polarization.pdf>`_ for the definition of the helicity angle.)DOC");
370 
371  REGISTER_VARIABLE("acoplanarityAngle", acoplanarityAngle,
372  R"DOC(
373  Acoplanarity angle (see ``Particle::getAcoplanarity``) assuming a two body decay of the particle and its daughters.
374  See `PDG Polarization Review <http://pdg.lbl.gov/2019/reviews/rpp2018-rev-b-decays-polarization.pdf>`_ for the definition of the acoplanarity angle.)DOC");
375 
376  }
378 }
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::Variable::Manager::FunctionPtr
std::function< double(const Particle *)> FunctionPtr
NOTE: the python interface is documented manually in analysis/doc/Variables.rst (because we use ROOT ...
Definition: Manager.h:118
Belle2::ReferenceFrame::GetCurrent
static const ReferenceFrame & GetCurrent()
Get current rest frame.
Definition: ReferenceFrame.cc:28