Belle II Software light-2406-ragdoll
FourVector.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * Forked from https://github.com/iLCSoft/MarlinKinfit *
6 * *
7 * Further information about the fit engine and the user interface *
8 * provided in MarlinKinfit can be found at *
9 * https://www.desy.de/~blist/kinfit/doc/html/ *
10 * and in the LCNotes LC-TOOL-2009-001 and LC-TOOL-2009-004 available *
11 * from http://www-flc.desy.de/lcnotes/ *
12 * *
13 * See git log for contributors and copyright holders. *
14 * This file is licensed under LGPL-3.0, see LICENSE.md. *
15 **************************************************************************/
16
17#ifdef MARLIN_USE_ROOT
18
19#include "analysis/OrcaKInFit/FourVector.h"
20
21#include <cmath>
22
23#undef NDEBUG
24#include <cassert>
25
26#include <TRandom3.h>
27
28namespace Belle2 {
33 namespace OrcaKinFit {
34
35 static TRandom* rnd = 0;
36
37 FourVector& FourVector::boost(const FourVector& P)
38 {
39 // See CERNLIB U101 for a description
40
41 double pP = -(p * P.p);
42 double e = getE();
43 double M = P.getM();
44
45 E = (e * P.getE() - pP) / M;
46 p = p - ((pP / (P.getE() + M) - e) / M) * P.p;
47
48 return *this;
49 }
50
51 void FourVector::decayto(FourVector& d1, FourVector& d2) const
52 {
53 // Let this particle decay isotropically into 4-vectors d1 and d2;
54 // d1 and d2 must have definite mass at beginning
55 using std::abs;
56 using std::sqrt;
57 using std::pow;
58
59 double M2 = getM2();
60 double M = getM();;
61 double m1 = d1.getM();
62 double m2 = d2.getM();
63
64 double randoms[2];
65 if (rnd == 0) rnd = new TRandom3();
66 rnd->RndmArray(2, randoms);
67
68
69 assert(m1 + m2 <= M);
70
71 double pstar = 0.5 * sqrt(abs((M2 - pow(m1 + m2, 2)) * (M2 - pow(m1 - m2, 2)))) / M;
72 double phistar = 2 * M_PI * randoms[0];
73 double costhetastar = 2 * randoms[1] - 1;
74 double sinthetastar = sqrt(abs(1 - costhetastar * costhetastar));
75 double E1 = sqrt(m1 * m1 + pstar * pstar);
76 double E2 = sqrt(m2 * m2 + pstar * pstar);
77
78
79 d1 = FourVector(E1, pstar * sinthetastar * cos(phistar),
80 pstar * sinthetastar * sin(phistar),
81 pstar * costhetastar);
82 d2 = FourVector(E2, -pstar * sinthetastar * cos(phistar),
83 -pstar * sinthetastar * sin(phistar),
84 -pstar * costhetastar);
85
86 d1.boost(*this);
87 d2.boost(*this);
88
89 }
90
91 }// end OrcaKinFit namespace
93} // end Belle2 namespace
94
95#endif // MARLIN_USE_ROOT
FourVector()
Default constructor.
Definition: FourVector.h:111
double E
The energy / 0 component.
Definition: FourVector.h:107
double getM2() const
Returns the mass squared / magnitude squared.
Definition: FourVector.h:138
double getM() const
Returns the mass / magnitude.
Definition: FourVector.h:139
ThreeVector p
The momentum three vector.
Definition: FourVector.h:108
double getE() const
Returns the energy / 0 component.
Definition: FourVector.h:127
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24