//----------Author:        Markus D. Oldenburg
//----------Last Modified: 09.03.1999
//----------Copyright:     &copy MDO Production 1996


#include "TMath.h"
#include "TRandom.h"
#include "MScatter.h"


///////////////////////////////////////////////////////////////////////////////
// MScatter                                                                  //
//                                                                           //
// Calculates the scattering of two pions. The input values are the two      //
// four momenta of the pions.                                                //
//                                                                           //
// To tell the truth - this is not a real class, because it has no data      //
// members. It is just a collection of functions which are needed to         //
// calculate a relativistic scatter process.                                 //
//                                                                           //
// On the other hand a lot of usefull functions are members of this 'class': //
//                                                                           //
//   - some random routines to obtain randomly distibuted angles             //
//   - three dimensional turns of coordinate systems of given angles         //
//   - reltivistic calculations of center of gravity, total momentum,        //
//     velocities, energies and gamma                                        //
//   - lorentz boosts                                                        //
///////////////////////////////////////////////////////////////////////////////

ClassImp(MScatter)


 MScatter::MScatter() {
  // Constructor. Creates a MScatter object. 
}


 MScatter::~MScatter() {
  // MScatter destructor. Does nothing else.
}
  

 void MScatter::Scattering(double p1[4], double p2[4])
{
  // Complete scattering. 
  // Calculates everything which is needed for one complete scattering
  // of two pions with the given four momenta p1 and p2.

  int i; 
  int parity;
  
  double P[4];
  double phi, tetha, v;
  
  Center_of_gravity(&P[0], &p1[0], &p2[0]);
  
  parity = 1;

  if (P[2] < 0.) {

    parity = -1;

    for (i = 1; i < 4; i++) {

      p1[i] = -p1[i];
      p2[i] = -p2[i];
       P[i] =  -P[i];
    }
  }
  
  Turn(&p1[0], &p2[0], phi = Get_phi(&P[0]), tetha = Get_tetha(&P[0]));
  Lorentz(v = Get_velocity(&p1[0], &p2[0]), &p1[0], &p2[0]);
  Scatter(Dice_tetha(Get_R(&p1[0], &p2[0])), Dice_phi(), &p1[0], &p2[0]);
  Lorentz(-v, &p1[0], &p2[0]);
  Turnback(&p1[0], &p2[0], phi, tetha);

  if (parity == -1) {

   for (i = 1; i < 4; i++) {

      p1[i] = -p1[i];
      p2[i] = -p2[i];
       P[i] =  -P[i];
    }
  }
  
  return;
}


 double MScatter::Get_phi(double P[4])
{
  // Returns the angle phi between a four vector and the x-axis.

  if(sqrt(P[1]*P[1] + P[2]*P[2]) == 0.) {

    return 0.;
  }

  else {

    if (P[2] > 0.) {

      return +acos(P[1] / sqrt(P[1]*P[1] + P[2]*P[2]));
    }

    if (P[2] < 0.) {

      return  -acos(P[1] / sqrt(P[1]*P[1] + P[2]*P[2]));
    }

    if (P[2] == 0.) {

      if (P[1] > 0.) {

        return +TMath::Pi()/2.;
      }

      else {

        return -TMath::Pi()/2.;
      }
    }
  }
  // The following line will never be executed.
  // I just added it for getting rid of the warning
  // control reaches end of non-void function `MScatter::Get_phi(double *)'
  return 0.;
}


 double MScatter::Get_tetha(double P[4])
{
  // Returns the angle tetha between a four vector and the x-axis.

  if(sqrt(P[1]*P[1] + P[2]*P[2] + P[3]*P[3]) == 0.) {

    return 0.;
  }

  else {

    return asin(P[3] / sqrt(P[1]*P[1] + P[2]*P[2] + P[3]*P[3]));
  }
}


 double MScatter::Get_energy(double p[4])
{
  // Calculates with the given four momentum the energy of an particle. 

  return sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2] + p[3]*p[3]);
}


 double MScatter::Get_velocity(double p1[4], double p2[4])
{  
  // Returns the velocity of the center of mass of an two particle system.
  
  return sqrt((p1[1] + p2[1]) * (p1[1] + p2[1]) +
              (p1[2] + p2[2]) * (p1[2] + p2[2]) +
              (p1[3] + p2[3]) * (p1[3] + p2[3])) /
    (Get_energy(&p1[0]) + Get_energy(&p2[0]));
}


 double MScatter::Get_R(double p1[4], double p2[4])
{
  // Calculates R (I can't remember what R is good for. 
  // May be some scattering cross section? The source of this is either
  // the particle data book or Perkins.)

  double s;
  double mr = 0.769, taur = 0.154;

  s = Get_invariant_momentum(&p1[0], &p2[0]);

  return (s - (p1[0] - p2[0]) * (p1[0] - p2[0]))/
    ((s - mr*mr) * (s - mr*mr) + mr*mr * taur*taur);
}


 double MScatter::Dice_phi(void)
{  
  // Returns a randomly distributed scatter angle phi.

  static TRandom *myRandom;
  static int FirstCall = kTRUE;

  if(FirstCall) {

    FirstCall = kFALSE;
    myRandom = new TRandom();
  }

  return myRandom->Rndm() * 2*TMath::Pi();
}


 double MScatter::Dice_tetha(double R)
{
  // Returns a randomly distributed scatter angle tetha.

  double x, y;
  double wert, height;

  static TRandom *myRandom;
  static int FirstCall = kTRUE;

  if(FirstCall) {

    FirstCall = kFALSE;
    myRandom = new TRandom();
  }


  wert=3*1.9E5*R;
  height=0.8+wert;

  do {

    x=myRandom->Rndm() * TMath::Pi();
    y=myRandom->Rndm() * height;

  } while (y > .8 + wert * cos(x)*cos(x));
  
  return x;

}


 double MScatter::Get_gamma(double v)
{
  // Calculates gamma for a given velocity.

  return 1 / (sqrt(1 - v*v));
}


 double MScatter::Get_invariant_momentum(double p1[4], double p2[4])
{
  // Returns the invariant four momentum of a two particle system.

  return ((Get_energy(&p1[0]) + Get_energy(&p2[0])) *
          (Get_energy(&p1[0]) + Get_energy(&p2[0])) -
          ((p1[1] + p2[1]) * (p1[1] + p2[1]) +
           (p1[2] + p2[2]) * (p1[2] + p2[2]) +
           (p1[3] + p2[3]) * (p1[3] + p2[3])));
}


 void MScatter::Scatter(double tetha, double phi, double p1[4], double p2[4])
{
  // Calculates the new four momenta after a scattering with the given
  // angles tetha and phi.

  double l1[4], l2[4];
  double alpha, betha;

  int j;
  
  if (sqrt(p1[1] * p1[1] + p1[2] * p1[2]) == 0.) {
    alpha = 0.;
  }

  else {
    alpha = atan(p1[3] / sqrt(p1[1]*p1[1] + p1[2]*p1[2]));
  }
  
  if (p1[1] == 0.) {
    betha = 0.;
  }

  else {
    betha = atan(p1[2] / p1[1]);
  }
  

  // Turn the coordinate system in the direction of the momenta

  p1[1] =  sqrt(p1[1]*p1[1] + p1[2]*p1[2] + p1[3]*p1[3]);
  p2[1] = -sqrt(p2[1]*p2[1] + p2[2]*p2[2] + p2[3]*p2[3]);

  p2[2] = p1[2] = 0.;
  p2[3] = p1[3] = 0.;


  // actual scattering happens here
  
  l1[1] = p1[1] * cos(tetha);
  l1[2] = p1[1] * sin(tetha) *sin(phi);
  l1[3] = p1[1] * sin(tetha) *cos(phi);

  l2[1] =p2[1] * cos(tetha);
  l2[2] =p2[1] * sin(tetha) * sin(phi);
  l2[3] =p2[1] * sin(tetha) * cos(phi);

  for (j = 1; j < 4;j++) {

    p1[j] = l1[j];
    p2[j] = l2[j];
  }

  

  // Turn the coordinate system back in the primal direction
    
  l1[1] = cos(betha) * cos(alpha) * p1[1] -
    sin(betha) * p1[2] -
    cos(betha) * sin(alpha) * p1[3];
  
  l1[2] = sin(betha) * cos(alpha) * p1[1] +
    cos(betha) * p1[2] -
    sin(betha) * sin(alpha) * p1[3];
  
  l1[3] = sin(alpha) * p1[1] +
    cos(alpha) * p1[3];

  l2[1] = cos(betha) * cos(alpha) * p2[1] -
    sin(betha) * p2[2] -
    cos(betha) * sin(alpha) * p2[3];
  
  l2[2] = sin(betha) * cos(alpha) * p2[1] +
    cos(betha) * p2[2] -
    sin(betha) * sin(alpha) * p2[3];
  
  l2[3] = sin(alpha) * p2[1] +
    cos(alpha) * p2[3];

  for (j = 1; j < 4; j++) {

    p1[j] = l1[j];
    p2[j] = l2[j];
  }  
  
  return;
}


 void MScatter::Center_of_gravity(double P[4], double p1[4], double p2[4])
{
  // Calculates the center of gravity P (the sum of two four vetors, p1 and p2).

  int i;

  for (i = 0; i < 4; i++) {
    P[i] = p1[i] + p2[i];
  }

  return;
}


 void MScatter::Turn(double p1[4], double p2[4], double phi, double tetha)
{
  // Turns the x-axis in the direction of the total momentum.

  int j;

  double pturn1[4], pturn2[4];
   
  pturn1[1] = cos(tetha) * cos(phi) * p1[1] +
    cos(tetha) * sin(phi) * p1[2] +
    sin(tetha) * p1[3];
  
  pturn1[2] = -sin(phi) * p1[1] +
    cos(phi) * p1[2];
  
  pturn1[3] = -sin(tetha) * cos(phi) * p1[1] -
    sin(tetha) * sin(phi) * p1[2] +
    cos(tetha) * p1[3];
  

  pturn2[1] = cos(tetha) * cos(phi) * p2[1] +
    cos(tetha) * sin(phi) * p2[2] +
    sin(tetha) * p2[3];
  
  pturn2[2] = -sin(phi) * p2[1] +
    cos(phi) * p2[2];
  
  pturn2[3] = -sin(tetha) * cos(phi) * p2[1] -
    sin(tetha) * sin(phi) * p2[2] +
    cos(tetha) * p2[3];
  
  for (j = 1; j < 4; j++) {

    p1[j] = pturn1[j];
    p2[j] = pturn2[j];
  }
  
  return;
} 


 void MScatter::Turnback(double p1[4], double p2[4], double phi, double tetha)
{
  // Turns the coordinate system back in the primal direction.

  int j;
  double pturn1[4], pturn2[4];
  
  pturn1[1] = cos(tetha) * cos(phi) * p1[1] -
    sin(phi) * p1[2] -
    cos(phi) * sin(tetha)*p1[3];
  
  pturn1[2] = sin(phi) * cos(tetha) * p1[1] +
    cos(phi) * p1[2] -
    sin(phi) * sin(tetha) * p1[3];
  
  pturn1[3] = sin(tetha) * p1[1] +
    cos(tetha) * p1[3];


  pturn2[1] = cos(tetha) * cos(phi) * p2[1] -
    sin(phi) * p2[2] -
    cos(phi) * sin(tetha) * p2[3];
  
  pturn2[2] = sin(phi) * cos(tetha) * p2[1]+
    cos(phi) * p2[2] -
    sin(phi) * sin(tetha) * p2[3];
  
  pturn2[3] = sin(tetha) * p2[1] +
    cos(tetha) * p2[3];

  
  for (j = 1; j < 4; j++) {

    p1[j] = pturn1[j];
    p2[j] = pturn2[j];
  }  
  
  return;
}


 void MScatter::Lorentz(double v, double p1[4], double p2[4])
{
  // Lorentz transform of two four momenta for a given velocity in 
  // the direction of the x-axis.

  p1[1] = Get_gamma(v) * (p1[1] - v * Get_energy(&p1[0]));
  p2[1] = Get_gamma(v) * (p2[1] - v * Get_energy(&p2[0]));  

  return;
}