CMSSW Reference Manual

Decay3Body::Decay3Body ( )

Definition at line 16 of file Decay3Body.cc.

00016                        {
00017 }

Decay3Body::~Decay3Body ( )

Definition at line 19 of file Decay3Body.cc.

00019                         {
00020 }

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

Definition at line 23 of file Decay3Body.cc.

References funct::cos(), GenMuonPlsPt100GeV_cfg::cout, e1, e2, e3, lat::endl(), m1, p1, p2, p3, phi, pi, rot, funct::sin(), sqr(), funct::sqrt(), and theta.

00026                                                             {
00027 
00028   double m0 = mother.m();
00029   double m1 = daughter1.m();
00030   double m2 = daughter2.m();
00031   double m3 = daughter3.m();
00032   double sumM2 = m0*m0 + m1*m1 + m2*m2 + m3*m3;
00033   double tolerance = 1.0e-9;
00034 
00035   math::XYZTLorentzVectorD mmm(mother.px(),mother.py(),mother.pz(),mother.e());
00036 
00037   if (m0 < m1+m2+m3) {
00038         std::cout << "Error: Daughters too heavy!" << std::endl;
00039     std::cout << "M: " << m0/GeV <<
00040        " < m1+m2+m3: " << m1/GeV + m2/GeV + m3/GeV << std::endl;
00041     return;
00042   } else {
00043     double m2_12max = sqr(m0-m3);
00044     double m2_12min = sqr(m1+m2);
00045     double m2_23max = sqr(m0-m1);
00046     double m2_23min = sqr(m2+m3);
00047 
00048     double x1,x2;
00049     double m2_12 = 0.0;
00050     double m2_23 = 0.0;
00051     double E2_12,E3_12;
00052         double m2_23max_12,m2_23min_12;
00053 
00054     do {
00055 // Pick values for m2_12 and m2_23 uniformly:
00056       x1 = G4UniformRand();
00057       m2_12 = m2_12min + x1*(m2_12max-m2_12min);
00058       x2 = G4UniformRand();
00059       m2_23 = m2_23min + x2*(m2_23max-m2_23min);
00060 
00061 // From the allowed range of m2_23 (given m2_12), determine if the point is valid:
00062 // (formulae taken from PDG booklet 2004 kinematics, page 305, Eqs. 38.22a+b)
00063       E2_12 = (m2_12 - m1*m1 + m2*m2)/(2.0*sqrt(m2_12));
00064       E3_12 = (m0*m0 - m2_12 - m3*m3)/(2.0*sqrt(m2_12));
00065       m2_23max_12 = sqr(E2_12+E3_12)-sqr(sqrt(sqr(E2_12)-m2*m2)-sqrt(sqr(E3_12)-m3*m3));
00066       m2_23min_12 = sqr(E2_12+E3_12)-sqr(sqrt(sqr(E2_12)-m2*m2)+sqrt(sqr(E3_12)-m3*m3));
00067     } while ((m2_23 > m2_23max_12) || (m2_23 < m2_23min_12));
00068 
00069 // Determine the value of the third invariant mass squared:
00070     double m2_13 = sumM2 - m2_12 - m2_23;
00071 
00072 // Calculate the energy and size of the momentum of the three daughters:  
00073     double e1 = (m0*m0 + m1*m1 - m2_23)/(2.0*m0);
00074     double e2 = (m0*m0 + m2*m2 - m2_13)/(2.0*m0);
00075     double e3 = (m0*m0 + m3*m3 - m2_12)/(2.0*m0);
00076     double p1 = sqrt(e1*e1 - m1*m1);
00077     double p2 = sqrt(e2*e2 - m2*m2);
00078     double p3 = sqrt(e3*e3 - m3*m3);
00079 
00080 // Calculate cosine of the relative angles between the three daughters:
00081     double cos12 = (m1*m1 + m2*m2 + 2.0*e1*e2 - m2_12)/(2.0*p1*p2);
00082     double cos13 = (m1*m1 + m3*m3 + 2.0*e1*e3 - m2_13)/(2.0*p1*p3);
00083     double cos23 = (m2*m2 + m3*m3 + 2.0*e2*e3 - m2_23)/(2.0*p2*p3);
00084     if (fabs(cos12) > 1.0) std::cout << "Error: Undefined angle12!" << std::endl;
00085     if (fabs(cos13) > 1.0) std::cout << "Error: Undefined angle13!" << std::endl;
00086     if (fabs(cos23) > 1.0) std::cout << "Error: Undefined angle23!" << std::endl;
00087 
00088 // Find the four vectors of the particles in a chosen (i.e. simple) frame:
00089     double xi    = 2.0 * pi * G4UniformRand();
00090     math::XYZVectorD q1(0.0,0.0,p1);
00091     math::XYZVectorD q2( sin(acos(cos12))*cos(xi)*p2, sin(acos(cos12))*sin(xi)*p2,cos12*p2);
00092     math::XYZVectorD q3(-sin(acos(cos13))*cos(xi)*p3,-sin(acos(cos13))*sin(xi)*p3,cos13*p3);
00093 
00094 // Rotate all three daughters momentum with the angles theta and phi:
00095     double theta = acos(2.0 * G4UniformRand() - 1.0);
00096     double phi   = 2.0 * pi * G4UniformRand();
00097     double psi   = 2.0 * pi * G4UniformRand();
00098     
00099     ROOT::Math::EulerAngles ang(phi,theta,psi);
00100     ROOT::Math::Rotation3D rot(ang);
00101 
00102     math::XYZVectorD q1rot = rot*q1;
00103     math::XYZVectorD q2rot = rot*q2;
00104     math::XYZVectorD q3rot = rot*q3;
00105 
00106     math::XYZTLorentzVectorD daughter1_orig(q1rot.X(),q1rot.Y(),q1rot.Z(),e1);
00107     math::XYZTLorentzVectorD daughter2_orig(q2rot.X(),q2rot.Y(),q2rot.Z(),e2);
00108     math::XYZTLorentzVectorD daughter3_orig(q3rot.X(),q3rot.Y(),q3rot.Z(),e3);
00109 
00110     ROOT::Math::Boost cmboost(mmm.BoostToCM());
00111 
00112     // Check of total angle and momentum:
00113     if (acos(cos12)+acos(cos13)+acos(cos23)-2.0*pi > tolerance)
00114       std::cout << "Error: Total angle not 2pi! " <<
00115         acos(cos12)+acos(cos13)+acos(cos23)-2.0*pi << std::endl;
00116     if (fabs(daughter1_orig.px()+daughter2_orig.px()+daughter3_orig.px())/GeV > tolerance)
00117       std::cout << "Error: Total 3B Px not conserved! " << 
00118         (daughter1_orig.px()+daughter2_orig.px()+daughter3_orig.px())/GeV << std::endl;
00119     if (fabs(daughter1_orig.py()+daughter2_orig.py()+daughter3_orig.py())/GeV > tolerance)
00120       std::cout << "Error: Total 3B Py not conserved! " << 
00121         (daughter1_orig.py()+daughter2_orig.py()+daughter3_orig.py())/GeV << std::endl;
00122     if (fabs(daughter1_orig.pz()+daughter2_orig.pz()+daughter3_orig.pz())/GeV > tolerance)
00123       std::cout << "Error: Total 3B Pz not conserved! " << 
00124         (daughter1.pz()+daughter2.pz()+daughter3.pz())/GeV << std::endl;
00125     
00126     // Boost the daughters back to the frame of the mother:
00127 
00128     math::XYZTLorentzVectorD temp1(cmboost(daughter1_orig));
00129     math::XYZTLorentzVectorD temp2(cmboost(daughter2_orig));
00130     math::XYZTLorentzVectorD temp3(cmboost(daughter3_orig));
00131     
00132     daughter1.setPx(temp1.Px());
00133     daughter1.setPy(temp1.Py());
00134     daughter1.setPz(temp1.Pz());
00135     daughter1.setE(temp1.e());
00136 
00137     daughter2.setPx(temp2.Px());
00138     daughter2.setPy(temp2.Py());
00139     daughter2.setPz(temp2.Pz());
00140     daughter2.setE(temp2.e());
00141 
00142     daughter3.setPx(temp3.Px());
00143     daughter3.setPy(temp3.Py());
00144     daughter3.setPz(temp3.Pz());
00145     daughter3.setE(temp3.e());
00146 
00147     return;
00148   }
00149 }

double Decay3Body::sqr ( double a ) [inline, private]

Definition at line 152 of file Decay3Body.cc.

Referenced by doDecay().

00152                                {
00153   return a*a;
00154 }


Public Member Functions
	Decay3Body ()
void	doDecay (const G4LorentzVector &mother, G4LorentzVector &daughter1, G4LorentzVector &daughter2, G4LorentzVector &daughter3)
	~Decay3Body ()
Private Member Functions
double	sqr (double a)

Decay3Body Class Reference

Public Member Functions

Private Member Functions

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation