CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
List of all members | Public Member Functions | Private Member Functions
Decay3Body Class Reference

#include <Decay3Body.h>

Public Member Functions

 Decay3Body ()
 
void doDecay (const G4LorentzVector &mother, G4LorentzVector &daughter1, G4LorentzVector &daughter2, G4LorentzVector &daughter3)
 
 ~Decay3Body ()
 

Private Member Functions

double sqr (double a)
 

Detailed Description

Definition at line 6 of file Decay3Body.h.

Constructor & Destructor Documentation

Decay3Body::Decay3Body ( )

Definition at line 16 of file Decay3Body.cc.

16 {}
Decay3Body::~Decay3Body ( )

Definition at line 18 of file Decay3Body.cc.

18 {}

Member Function Documentation

void Decay3Body::doDecay ( const G4LorentzVector &  mother,
G4LorentzVector &  daughter1,
G4LorentzVector &  daughter2,
G4LorentzVector &  daughter3 
)

Definition at line 20 of file Decay3Body.cc.

References funct::cos(), gather_cfg::cout, GeV, callgraph::m2, fireworks::p1, fireworks::p2, phi, pi, makeMuonMisalignmentScenario::rot, funct::sin(), sqr(), mathSSE::sqrt(), theta(), and tolerance.

23  {
24  double m0 = mother.m();
25  double m1 = daughter1.m();
26  double m2 = daughter2.m();
27  double m3 = daughter3.m();
28  double sumM2 = m0 * m0 + m1 * m1 + m2 * m2 + m3 * m3;
29  double tolerance = 1.0e-9;
30 
31  math::XYZTLorentzVectorD mmm(mother.px(), mother.py(), mother.pz(), mother.e());
32 
33  if (m0 < m1 + m2 + m3) {
34  std::cout << "Error: Daughters too heavy!" << std::endl;
35  std::cout << "M: " << m0 / GeV << " < m1+m2+m3: " << m1 / GeV + m2 / GeV + m3 / GeV << std::endl;
36  return;
37  } else {
38  double m2_12max = sqr(m0 - m3);
39  double m2_12min = sqr(m1 + m2);
40  double m2_23max = sqr(m0 - m1);
41  double m2_23min = sqr(m2 + m3);
42 
43  double x1, x2;
44  double m2_12 = 0.0;
45  double m2_23 = 0.0;
46  double E2_12, E3_12;
47  double m2_23max_12, m2_23min_12;
48 
49  do {
50  // Pick values for m2_12 and m2_23 uniformly:
51  x1 = G4UniformRand();
52  m2_12 = m2_12min + x1 * (m2_12max - m2_12min);
53  x2 = G4UniformRand();
54  m2_23 = m2_23min + x2 * (m2_23max - m2_23min);
55 
56  // From the allowed range of m2_23 (given m2_12), determine if the point is valid:
57  // (formulae taken from PDG booklet 2004 kinematics, page 305, Eqs. 38.22a+b)
58  E2_12 = (m2_12 - m1 * m1 + m2 * m2) / (2.0 * sqrt(m2_12));
59  E3_12 = (m0 * m0 - m2_12 - m3 * m3) / (2.0 * sqrt(m2_12));
60  m2_23max_12 = sqr(E2_12 + E3_12) - sqr(sqrt(sqr(E2_12) - m2 * m2) - sqrt(sqr(E3_12) - m3 * m3));
61  m2_23min_12 = sqr(E2_12 + E3_12) - sqr(sqrt(sqr(E2_12) - m2 * m2) + sqrt(sqr(E3_12) - m3 * m3));
62  } while ((m2_23 > m2_23max_12) || (m2_23 < m2_23min_12));
63 
64  // Determine the value of the third invariant mass squared:
65  double m2_13 = sumM2 - m2_12 - m2_23;
66 
67  // Calculate the energy and size of the momentum of the three daughters:
68  double e1 = (m0 * m0 + m1 * m1 - m2_23) / (2.0 * m0);
69  double e2 = (m0 * m0 + m2 * m2 - m2_13) / (2.0 * m0);
70  double e3 = (m0 * m0 + m3 * m3 - m2_12) / (2.0 * m0);
71  double p1 = sqrt(e1 * e1 - m1 * m1);
72  double p2 = sqrt(e2 * e2 - m2 * m2);
73  double p3 = sqrt(e3 * e3 - m3 * m3);
74 
75  // Calculate cosine of the relative angles between the three daughters:
76  double cos12 = (m1 * m1 + m2 * m2 + 2.0 * e1 * e2 - m2_12) / (2.0 * p1 * p2);
77  double cos13 = (m1 * m1 + m3 * m3 + 2.0 * e1 * e3 - m2_13) / (2.0 * p1 * p3);
78  double cos23 = (m2 * m2 + m3 * m3 + 2.0 * e2 * e3 - m2_23) / (2.0 * p2 * p3);
79  if (fabs(cos12) > 1.0)
80  std::cout << "Error: Undefined angle12!" << std::endl;
81  if (fabs(cos13) > 1.0)
82  std::cout << "Error: Undefined angle13!" << std::endl;
83  if (fabs(cos23) > 1.0)
84  std::cout << "Error: Undefined angle23!" << std::endl;
85 
86  // Find the four vectors of the particles in a chosen (i.e. simple) frame:
87  double xi = 2.0 * pi * G4UniformRand();
88  math::XYZVectorD q1(0.0, 0.0, p1);
89  math::XYZVectorD q2(sin(acos(cos12)) * cos(xi) * p2, sin(acos(cos12)) * sin(xi) * p2, cos12 * p2);
90  math::XYZVectorD q3(-sin(acos(cos13)) * cos(xi) * p3, -sin(acos(cos13)) * sin(xi) * p3, cos13 * p3);
91 
92  // Rotate all three daughters momentum with the angles theta and phi:
93  double theta = acos(2.0 * G4UniformRand() - 1.0);
94  double phi = 2.0 * pi * G4UniformRand();
95  double psi = 2.0 * pi * G4UniformRand();
96 
97  ROOT::Math::EulerAngles ang(phi, theta, psi);
98  ROOT::Math::Rotation3D rot(ang);
99 
100  math::XYZVectorD q1rot = rot * q1;
101  math::XYZVectorD q2rot = rot * q2;
102  math::XYZVectorD q3rot = rot * q3;
103 
104  math::XYZTLorentzVectorD daughter1_orig(q1rot.X(), q1rot.Y(), q1rot.Z(), e1);
105  math::XYZTLorentzVectorD daughter2_orig(q2rot.X(), q2rot.Y(), q2rot.Z(), e2);
106  math::XYZTLorentzVectorD daughter3_orig(q3rot.X(), q3rot.Y(), q3rot.Z(), e3);
107 
108  ROOT::Math::Boost cmboost(mmm.BoostToCM());
109 
110  // Check of total angle and momentum:
111  if (acos(cos12) + acos(cos13) + acos(cos23) - 2.0 * pi > tolerance)
112  std::cout << "Error: Total angle not 2pi! " << acos(cos12) + acos(cos13) + acos(cos23) - 2.0 * pi << std::endl;
113  if (fabs(daughter1_orig.px() + daughter2_orig.px() + daughter3_orig.px()) / GeV > tolerance)
114  std::cout << "Error: Total 3B Px not conserved! "
115  << (daughter1_orig.px() + daughter2_orig.px() + daughter3_orig.px()) / GeV << std::endl;
116  if (fabs(daughter1_orig.py() + daughter2_orig.py() + daughter3_orig.py()) / GeV > tolerance)
117  std::cout << "Error: Total 3B Py not conserved! "
118  << (daughter1_orig.py() + daughter2_orig.py() + daughter3_orig.py()) / GeV << std::endl;
119  if (fabs(daughter1_orig.pz() + daughter2_orig.pz() + daughter3_orig.pz()) / GeV > tolerance)
120  std::cout << "Error: Total 3B Pz not conserved! " << (daughter1.pz() + daughter2.pz() + daughter3.pz()) / GeV
121  << std::endl;
122 
123  // Boost the daughters back to the frame of the mother:
124 
125  math::XYZTLorentzVectorD temp1(cmboost(daughter1_orig));
126  math::XYZTLorentzVectorD temp2(cmboost(daughter2_orig));
127  math::XYZTLorentzVectorD temp3(cmboost(daughter3_orig));
128 
129  daughter1.setPx(temp1.Px());
130  daughter1.setPy(temp1.Py());
131  daughter1.setPz(temp1.Pz());
132  daughter1.setE(temp1.e());
133 
134  daughter2.setPx(temp2.Px());
135  daughter2.setPy(temp2.Py());
136  daughter2.setPz(temp2.Pz());
137  daughter2.setE(temp2.e());
138 
139  daughter3.setPx(temp3.Px());
140  daughter3.setPy(temp3.Py());
141  daughter3.setPz(temp3.Pz());
142  daughter3.setE(temp3.e());
143 
144  return;
145  }
146 }
GeV
const double GeV
Definition: MathUtil.h:16
math::XYZTLorentzVectorD
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > XYZTLorentzVectorD
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:14
fireworks::p2
const TString p2
Definition: fwPaths.cc:13
tolerance
const double tolerance
Definition: HGCalGeomParameters.cc:29
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Decay3Body::sqr
double sqr(double a)
Definition: Decay3Body.cc:148
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:150
psi
std::map< std::string, int, std::less< std::string > > psi
Definition: CountProcessesAction.h:15
callgraph.m2
tuple m2
Definition: callgraph.py:57
pi
const Double_t pi
Definition: trackSplitPlot.h:36
math::XYZVectorD
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double > > XYZVectorD
spatial vector with cartesian internal representation
Definition: Vector3D.h:8
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
fireworks::p1
const TString p1
Definition: fwPaths.cc:12
align::EulerAngles
AlgebraicVector EulerAngles
Definition: Definitions.h:34
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:144
double Decay3Body::sqr ( double  a)
inlineprivate

Definition at line 148 of file Decay3Body.cc.

References a.

Referenced by doDecay().

148 { return a * a; }
a
double a
Definition: hdecay.h:119
Generated for CMSSW Reference Manual by doxygen 1.8.5