/*  Particle.h

    Particle Class  (Sep 2008)
*/


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();
};


//------------------------------------------------------------------------------
// 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*Random::Flat()- 1.0 );
   double phi   = 2.0*M_PI *Random::Flat();

   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 );
}