/*  twoBody.cpp

    This is an example of two consecutive two-body decays. 
    First, a charged rho meson decays to a charged and a neutral pion.
    Then,  the neutral pion decays to a pair of photons.
    The decays are as follows:

          rho+-  -->  pi+-  +   pi0
                                 |
                                 +---> gamma1 + gamma2

    The two-body decays are implemented by a twoBodyDecay() routine defined 
    in Particle class.

    Program prints components of 4-momentums (px, py, pz, E) of the parent and 
    decay products in the lab frame in GeV. An example output is given below:

    P4( rho+-  ): ( 0.387106,  -0.876443,  -0.286352,  1.265403)
    P4( pi+-   ): ( 0.487285,  -0.544223,  -0.438312,  0.863263)
    P4( pi0    ): (-0.100180,  -0.332220,   0.151961,  0.402140)
    P4( gamma1 ): ( 0.000397,   0.009283,   0.020730,  0.022717)
    P4( gamma2 ): (-0.100577,  -0.341503,   0.131231,  0.379423)

    Aug 2008
*/

#include <iostream>
#include <cmath>
#include <cstdlib>
using namespace std;

//------------------------------------------------------------------------------
// Particle Class
//
class Particle{

  public:
    Particle();
    Particle(double);
    double   px, py, pz, E, m, pAbs, pt, p[4]; 
    void     p4(double, double, double, double);
    void     boost(Particle);
    void     twoBodyDecay(Particle&, Particle&);
    void     print();
};
//------------------------------------------------------------------------------
// Returns a random real number between [0,1]
//
double Rndm(){
   return ( double(rand())/RAND_MAX );
}
//------------------------------------------------------------------------------

int main(void)
{
  const int evMax = 10;

  // Define the particle names and set (nominal) rest masses in GeV/c2.
  // Data is taken fom PDG. <durpdg.dur.ac.uk>  
  Particle  rho(0.77540000); // rho meson (peak mass only)
  Particle  pi0(0.13497660); // neutral pion
  Particle  pic(0.13957018); // charged pion
  Particle  gamma1, gamma2;  // photons

  for(int event = 1; event <= evMax; event++){

     cout << "Event #: " << event << endl;

     //***
     //*** First decay: rho+- --> pi0 + pi+- (rho is at rest)
     //***
     rho.twoBodyDecay(pi0, pic);

     // Isotropic angles to give a random direction to the rho+-
     double theta = acos( 2.0*Rndm() - 1.0 );
     double phi   = 2.0 *M_PI*Rndm();

     // rho 4-momentum components in lab frame
     double P  = 1.0;  // GeV/c
     double e  = sqrt(P*P + rho.m*rho.m);
     double pX = P*sin(theta)*cos(phi);
     double pY = P*sin(theta)*sin(phi);
     double pZ = P*cos(theta);

     rho.p4( pX, pY, pZ, e );

     // boost the pions
     pi0.boost( rho );
     pic.boost( rho );

     //***
     //*** Second decay: pi0 --> gamma1 + gamma2 (pi0 rest frame)
     //***
     pi0.twoBodyDecay(gamma1, gamma2);

     // boost the gammas
     gamma1.boost( pi0 );
     gamma2.boost( pi0 );
    
     cout << "P4( rho+-  ): "; rho.print(); 
     cout << "P4( pi+-   ): "; pic.print(); 
     cout << "P4( pi0    ): "; pi0.print();
     cout << "P4( gamma1 ): "; gamma1.print();
     cout << "P4( gamma2 ): "; gamma2.print();
     
     cout << "________________________________________________________" << endl;
   }

} // main


//*****************************************************************************
//                                                                             *
//    MEMBERS functions of the Particle Class                                  * 
//                                                                             *
//*****************************************************************************

//
//*** Default constructor ------------------------------------------------------
//
Particle::Particle(){
  px = py = pz = E = m = pAbs = pt = 0.0;
  p[0] = p[1] = p[2] = p[3] = 0.0;
}
Particle::Particle(double mass){
  m = mass;
  px = py = pz = E = pAbs = pt = 0.0;
  p[0] = p[1] = p[2] = p[3] = 0.0;
}
//
//*** Sets components of 4-momentum vector -------------------------------------
//
void Particle::p4(double momx, double momy, double momz, double energy){
  // components of 4-momenta
  px = p[0] = momx;
  py = p[1] = momy;
  pz = p[2] = momz;
  E  = p[3] = energy;
  // transverse momentum and the magnitude of the space momentum
  pt      = sqrt(momx*momx + momy*momy);
  pAbs    = sqrt(momx*momx + momy*momy + momz*momz);
}
//
//*** Prints 4-vector ----------------------------------------------------------
//
void Particle::print(){
  cout << "(" << p[0] <<",\t" << p[1] <<",\t"<< p[2] <<",\t"<< p[3] << ")" << endl;
}
//
//*** Evaluates 4-vector of decay product in the rest frame of parent. ---------
//
void Particle::twoBodyDecay(Particle& prod1, Particle& prod2) {

   double m1 = prod1.m;
   double m2 = prod2.m;

   // check if particle m can decay
   if( m < m1+m2 ){
     cout << "Particle::twoBodyDecay  parent mass is less than sum of products." 
          << endl;
      prod1.p4 ( 0., 0., 0., m1);
      prod2.p4 ( 0., 0., 0., m2);
      return;
   }

   // CM energies and momentum
   double e1 = (m*m + m1*m1 - m2*m2) / (2.0*m);
   double e2 = (m*m - m1*m1 + m2*m2) / (2.0*m);
   double P  = sqrt(e1*e1 - m1*m1);

   // Isotropic random angles
   double theta = acos( 2.0*Rndm() - 1.0 );
   double phi   = 2.0*M_PI *Rndm();

   double pX = P*sin(theta)*cos(phi);
   double pY = P*sin(theta)*sin(phi);
   double pZ = P*cos(theta);

   // set the 4-momenta
   prod1.p4(  pX,  pY,  pZ, e1 );
   prod2.p4( -pX, -pY, -pZ, e2 );
}
//
//*** Lorentz Boost  -----------------------------------------------------------
//
void Particle::boost(Particle parent){

  // beta and gamma values
  double betax = parent.px / parent.E;
  double betay = parent.py / parent.E;
  double betaz = parent.pz / parent.E;
  double beta2 = betax*betax + betay*betay + betaz*betaz;
  double gamma = 1.0/sqrt(1.0-beta2);
  double dot   = betax*px + betay*py + betaz*pz;
  double prod  = gamma*( gamma*dot/(1.0+gamma) + E );

  double pX = px + betax*prod;
  double pY = py + betay*prod;
  double pZ = pz + betaz*prod;
  double e  = gamma*(E + dot);

  p4( pX, pY, pZ, e );
}