CustomPhysics/src/Decay3Body.cc

0001 #include <cmath>
0002 #include "SimG4Core/CustomPhysics/interface/Decay3Body.h"
0003
0004 #include <CLHEP/Units/SystemOfUnits.h>
0005 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0006 #include "DataFormats/Math/interface/Vector3D.h"
0007 #include "DataFormats/Math/interface/LorentzVector.h"
0008 #include "Math/GenVector/Rotation3D.h"
0009 #include "Math/GenVector/EulerAngles.h"
0010 #include "Math/GenVector/Boost.h"
0011 #include <vector>
0012 #include <cmath>
0013
0014 #include "Randomize.hh"
0015
0016 using CLHEP::GeV;
0017
0018 Decay3Body::Decay3Body() {}
0019
0020 Decay3Body::~Decay3Body() {}
0021
0022 void Decay3Body::doDecay(const G4LorentzVector& mother,
0023                          G4LorentzVector& daughter1,
0024                          G4LorentzVector& daughter2,
0025                          G4LorentzVector& daughter3) {
0026   double m0 = mother.m();
0027   double m1 = daughter1.m();
0028   double m2 = daughter2.m();
0029   double m3 = daughter3.m();
0030   double sumM2 = m0 * m0 + m1 * m1 + m2 * m2 + m3 * m3;
0031   double tolerance = 1.0e-9;
0032
0033   math::XYZTLorentzVectorD mmm(mother.px(), mother.py(), mother.pz(), mother.e());
0034
0035   if (m0 < m1 + m2 + m3) {
0036     std::cout << "Error: Daughters too heavy!" << std::endl;
0037     std::cout << "M: " << m0 / GeV << " < m1+m2+m3: " << m1 / GeV + m2 / GeV + m3 / GeV << std::endl;
0038     return;
0039   } else {
0040     double m2_12max = sqr(m0 - m3);
0041     double m2_12min = sqr(m1 + m2);
0042     double m2_23max = sqr(m0 - m1);
0043     double m2_23min = sqr(m2 + m3);
0044
0045     double x1, x2;
0046     double m2_12 = 0.0;
0047     double m2_23 = 0.0;
0048     double E2_12, E3_12;
0049     double m2_23max_12, m2_23min_12;
0050
0051     do {
0052       // Pick values for m2_12 and m2_23 uniformly:
0053       x1 = G4UniformRand();
0054       m2_12 = m2_12min + x1 * (m2_12max - m2_12min);
0055       x2 = G4UniformRand();
0056       m2_23 = m2_23min + x2 * (m2_23max - m2_23min);
0057
0058       // From the allowed range of m2_23 (given m2_12), determine if the point is valid:
0059       // (formulae taken from PDG booklet 2004 kinematics, page 305, Eqs. 38.22a+b)
0060       E2_12 = (m2_12 - m1 * m1 + m2 * m2) / (2.0 * sqrt(m2_12));
0061       E3_12 = (m0 * m0 - m2_12 - m3 * m3) / (2.0 * sqrt(m2_12));
0062       m2_23max_12 = sqr(E2_12 + E3_12) - sqr(sqrt(sqr(E2_12) - m2 * m2) - sqrt(sqr(E3_12) - m3 * m3));
0063       m2_23min_12 = sqr(E2_12 + E3_12) - sqr(sqrt(sqr(E2_12) - m2 * m2) + sqrt(sqr(E3_12) - m3 * m3));
0064     } while ((m2_23 > m2_23max_12) || (m2_23 < m2_23min_12));
0065
0066     // Determine the value of the third invariant mass squared:
0067     double m2_13 = sumM2 - m2_12 - m2_23;
0068
0069     // Calculate the energy and size of the momentum of the three daughters:
0070     double e1 = (m0 * m0 + m1 * m1 - m2_23) / (2.0 * m0);
0071     double e2 = (m0 * m0 + m2 * m2 - m2_13) / (2.0 * m0);
0072     double e3 = (m0 * m0 + m3 * m3 - m2_12) / (2.0 * m0);
0073     double p1 = sqrt(e1 * e1 - m1 * m1);
0074     double p2 = sqrt(e2 * e2 - m2 * m2);
0075     double p3 = sqrt(e3 * e3 - m3 * m3);
0076
0077     // Calculate cosine of the relative angles between the three daughters:
0078     double cos12 = (m1 * m1 + m2 * m2 + 2.0 * e1 * e2 - m2_12) / (2.0 * p1 * p2);
0079     double cos13 = (m1 * m1 + m3 * m3 + 2.0 * e1 * e3 - m2_13) / (2.0 * p1 * p3);
0080     double cos23 = (m2 * m2 + m3 * m3 + 2.0 * e2 * e3 - m2_23) / (2.0 * p2 * p3);
0081     if (fabs(cos12) > 1.0)
0082       std::cout << "Error: Undefined angle12!" << std::endl;
0083     if (fabs(cos13) > 1.0)
0084       std::cout << "Error: Undefined angle13!" << std::endl;
0085     if (fabs(cos23) > 1.0)
0086       std::cout << "Error: Undefined angle23!" << std::endl;
0087
0088     // Find the four vectors of the particles in a chosen (i.e. simple) frame:
0089     double xi = 2.0 * pi * G4UniformRand();
0090     math::XYZVectorD q1(0.0, 0.0, p1);
0091     math::XYZVectorD q2(sin(acos(cos12)) * cos(xi) * p2, sin(acos(cos12)) * sin(xi) * p2, cos12 * p2);
0092     math::XYZVectorD q3(-sin(acos(cos13)) * cos(xi) * p3, -sin(acos(cos13)) * sin(xi) * p3, cos13 * p3);
0093
0094     // Rotate all three daughters momentum with the angles theta and phi:
0095     double theta = acos(2.0 * G4UniformRand() - 1.0);
0096     double phi = 2.0 * pi * G4UniformRand();
0097     double psi = 2.0 * pi * G4UniformRand();
0098
0099     ROOT::Math::EulerAngles ang(phi, theta, psi);
0100     ROOT::Math::Rotation3D rot(ang);
0101
0102     math::XYZVectorD q1rot = rot * q1;
0103     math::XYZVectorD q2rot = rot * q2;
0104     math::XYZVectorD q3rot = rot * q3;
0105
0106     math::XYZTLorentzVectorD daughter1_orig(q1rot.X(), q1rot.Y(), q1rot.Z(), e1);
0107     math::XYZTLorentzVectorD daughter2_orig(q2rot.X(), q2rot.Y(), q2rot.Z(), e2);
0108     math::XYZTLorentzVectorD daughter3_orig(q3rot.X(), q3rot.Y(), q3rot.Z(), e3);
0109
0110     ROOT::Math::Boost cmboost(mmm.BoostToCM());
0111
0112     // Check of total angle and momentum:
0113     if (acos(cos12) + acos(cos13) + acos(cos23) - 2.0 * pi > tolerance)
0114       std::cout << "Error: Total angle not 2pi! " << acos(cos12) + acos(cos13) + acos(cos23) - 2.0 * pi << std::endl;
0115     if (fabs(daughter1_orig.px() + daughter2_orig.px() + daughter3_orig.px()) / GeV > tolerance)
0116       std::cout << "Error: Total 3B Px not conserved! "
0117                 << (daughter1_orig.px() + daughter2_orig.px() + daughter3_orig.px()) / GeV << std::endl;
0118     if (fabs(daughter1_orig.py() + daughter2_orig.py() + daughter3_orig.py()) / GeV > tolerance)
0119       std::cout << "Error: Total 3B Py not conserved! "
0120                 << (daughter1_orig.py() + daughter2_orig.py() + daughter3_orig.py()) / GeV << std::endl;
0121     if (fabs(daughter1_orig.pz() + daughter2_orig.pz() + daughter3_orig.pz()) / GeV > tolerance)
0122       std::cout << "Error: Total 3B Pz not conserved! " << (daughter1.pz() + daughter2.pz() + daughter3.pz()) / GeV
0123                 << std::endl;
0124
0125     // Boost the daughters back to the frame of the mother:
0126
0127     math::XYZTLorentzVectorD temp1(cmboost(daughter1_orig));
0128     math::XYZTLorentzVectorD temp2(cmboost(daughter2_orig));
0129     math::XYZTLorentzVectorD temp3(cmboost(daughter3_orig));
0130
0131     daughter1.setPx(temp1.Px());
0132     daughter1.setPy(temp1.Py());
0133     daughter1.setPz(temp1.Pz());
0134     daughter1.setE(temp1.e());
0135
0136     daughter2.setPx(temp2.Px());
0137     daughter2.setPy(temp2.Py());
0138     daughter2.setPz(temp2.Pz());
0139     daughter2.setE(temp2.e());
0140
0141     daughter3.setPx(temp3.Px());
0142     daughter3.setPy(temp3.Py());
0143     daughter3.setPz(temp3.Pz());
0144     daughter3.setE(temp3.e());
0145
0146     return;
0147   }
0148 }
0149
0150 double Decay3Body::sqr(double a) { return a * a; }