//----------Author:        Markus D. Oldenburg
//----------Last Modified: 23.04.1999
//----------Copyright:     &copy MDO Production 1999

#include "MFormulary.h"
#include "TMath.h"


///////////////////////////////////////////////////////////////////////////////
// MFormulary                                                                //
//                                                                           //
// The MFormulary class improves the TMath class by adding some often used   // 
// functions, constants, and calulations, especially for physicists,         //
// i.e. for me :)                                                            //
///////////////////////////////////////////////////////////////////////////////

ClassImp(MFormulary)


 MFormulary::MFormulary() : TMath() {
  // Constructor. Creates a MFormulary object.
} 


 MFormulary::~MFormulary() {
  // MFormulary destructor.
}


// mathematics

 Double_t MFormulary::Mt(const Double_t *m, const Double_t *pt) 
{
  // Calculates the transverse mass of a particle with a given mass and 
  // transverse momentum.

  return Sqrt(*m * *m + *pt * *pt);
}


 Double_t MFormulary::Pt(const Double_t *p) 
{
  // Calculates the transverse momentum of a particle with a given momentum.

  return Sqrt(p[0]*p[0] + p[1]*p[1]);
}


 Double_t MFormulary::Y(const Double_t *p, const Char_t t)
{
  // Calculates the rapidity y. In the second argument you have to specify 
  // the nature of the first argument: 
  //                        v for an ordinary velocity
  //                        b for betha, the velocity in units of c
  //                     or p for a four-momentum
  
  if (t == 'v') {
    return 1./2. * Log((1. + *p/C()) / (1. - *p/C()));
  }
  
  else if (t == 'b') {
    return 1./2. * Log((1. + *p) / (1. - *p));
  }
  
  else if (t == 'p') { 
    return 1./2. * Log((p[0] + p[3]) / (p[0] - p[3]));
  }
  
  else {
    cout << "Please specify the nature of the input velocity:" << endl;
    cout << "     v for an ordinary velocity," << endl;
    cout << "     b for betha, the velocity in units of c," <<endl;
    cout << "  or p for a four-momentum." << endl;
    
    exit(-1);
  }
}


 Double_t MFormulary::Y(const Double_t *eta, const Double_t *pt, const Double_t *m)
{
  // Calculates the rapidity of a particle with a given pseudorapidity eta,
  // transverse momentum pt, and mass m.

  return 1./2. * Log(Sqrt(*pt * *pt * CosH(*eta) * CosH(*eta) + *m * *m) + *pt * SinH(*eta) /
		     Sqrt(*pt * *pt * CosH(*eta) * CosH(*eta) + *m * *m) - *pt * SinH(*eta));
}


 Double_t MFormulary::Eta(const Double_t *p, const Char_t t)
{
  // Calculates the pseudorapidity eta. In the second argument you have 
  // to specify the nature of the first argument: 
  //                               a for the angle of the particle relative
  //                                 to the beam axis
  //                            or p for the (three dimensional) momentum
  
  if (t == 'a') {
    return -Log(Tan(*p/2.));
  }
  
  else if (t == 'p') {
    return 1./2. * Log((Abs(p, 3) + p[2]) / (Abs(p, 3) - p[2]));
  }
  
  else {
    cout << "Please specify the nature of the first argument:" << endl;
    cout << "     a for the angle of the particle relative to the beam axis," << endl;
    cout << "  or p for the (three dimensional) momentum." << endl;
    
    exit(-1);
  }
}


 Double_t MFormulary::Eta(const Double_t *y, const Double_t *mt, const Double_t *m)
{
  // Calculates the pseudorapidity of a particle with a given rapidity y,
  // transverse mass mt, and mass m.

  return 1./2. * Log(Sqrt(*mt * *mt * CosH(*y) * CosH(*y) - *m * *m) + *mt * SinH(*y) /
		     Sqrt(*mt * *mt * CosH(*y) * CosH(*y) - *m * *m) - *mt * SinH(*y));
}


 Double_t MFormulary::Energy(const Double_t 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_t MFormulary::Gamma(const Double_t *v)
{
  // Calculates gamma for a given velocity (v has to be in units of c!).

  return 1. / (Sqrt(1. - Power(*v, 2)));
}


 Double_t MFormulary::RelativisticSum(const Double_t *v1, const Double_t *v2)
{
  // Calculates the sum of two velocities with respect to the special
  // relativity. The function will return a warning if at least one of the 
  // input velocities is greater than c, but regardless of that it will 
  // return the (mathematically correct, but non-physical) result.

  if (*v1 > C() || *v2 > C()) {
    cout << "WARNING: The velocities should be less or equal c = 299792458 m/s." << endl;
  }

  // Attention: I've found a wierd feature of ROOT! *a * *b is interpreted
  // as *a ** b which means pow(*a, b). Because b is an address its value 
  // is quite high, which gives "Inf" as result in most cases!
  // The workaround is easy: just place the second address in brackets *a * (*b).
 
  return (*v1 + *v2) / (1. + ((*v1 * (*v2)) / (C() * C())));
}


 Double_t MFormulary::InvDist_squared(const Double_t x1[4], const Double_t x2[4])
{
  // This function returns the relativistic invariant distance squared of 
  // two four vectors.

  return C()*C()*Power(x2[0]-x1[0], 2) - Dist_squared(&x1[1], &x2[1], 3);
}


 Double_t MFormulary::InvDist(const Double_t x1[4], const Double_t x2[4])
{
  // This function returns the relativistic invariant distance of two four vectors.
  // It gives a warning, if the result inherits a factor i, but it will still 
  // return a value which is the result divided by i.

  Double_t abst_squared = InvDist_squared(x1, x2);

  if (abst_squared >= 0.) {
    return Sqrt(abst_squared);
  }
  
  else {
    cout << "WARNING: The true result is the returned result times i!" << endl;
    return Sqrt(-abst_squared);
  }
}


 Double_t MFormulary::RelativisticScalarProd(const Double_t x1[4], const Double_t x2[4])
{ 
  // This function returns the scalar product of two four vectors in Minkowsky metrik.

  return Power(x2[0]-x1[0], 2) - Dist_squared(&x1[1], &x2[1], 3); 
}


 void MFormulary::Lorentz(const Double_t *v, Double_t p[4])
{
  // Lorentz transform for a given velocity in the direction of the x-axis.

  p[1] = Gamma(v) * (p[1] - *v * Energy(&p[0]));

  return;
}


 void MFormulary::LorentzBoost(const Double_t x[4], Double_t x_prime[4], const Double_t v[3])
{
  // Calculates a Lorentz boost with the given velocity v of the input four 
  // vector x. The result is stored in x_prime.

  Double_t v_squared = Square(v, 3);
  Double_t v_absolut_divc = Sqrt(v_squared)/C();
  Double_t gamma = Gamma(&v_absolut_divc);

  x_prime[0] = gamma*x[0] - (v[1]*x[1] + v[2]*x[2] + v[3]*x[3]) / (C()*C());
  x_prime[1] = gamma*v[1]*x[0] + (1.+(gamma-1.)*v[1]*v[1]/v_squared)*x[1];
  x_prime[2] = gamma*v[2]*x[0] + (1.+(gamma-1.)*v[2]*v[2]/v_squared)*x[2];
  x_prime[3] = gamma*v[3]*x[0] + (1.+(gamma-1.)*v[3]*v[3]/v_squared)*x[3];

  return;
}


 Double_t MFormulary::S(const Double_t *p_a, const Double_t *p_b)
{ 
  // Returns the square of the center of mass energy of two colliding particles
  // with the four momenta p_a and p_b. This is the Mandelstam variable s.

  Double_t s[4];
  
  Sum(p_a, p_b, s, 4);
  return RelativisticScalarProd(s, s);
}				


 Double_t MFormulary::T(const Double_t p_a[4], const Double_t p_c[4])
{
  // Calculates the Mandelstam variable t.

  Double_t t[4];
  
  Diff(p_a, p_c, t, 4);
  return RelativisticScalarProd(t, t);			
}

 Double_t MFormulary::U(const Double_t p_a[4], const Double_t p_d[4])								
{
  // This is the Mandelstam variable u.

  return T(p_a, p_d);
}


// experimental

 Double_t MFormulary::Sigma_tot_pp(Double_t *p)
{
  // This function calculates the total cross section of a proton-proton-reaction.
  // p specifies the momentum of the incident proton in the laboratory frame.
  // The result is in units of mb, the momentum has to be in GeV.
  // The formula applies only for high energies!

  return 48. + 0.522 * Log(*p) * Log(*p) - 4.51 * Log(*p);
}


 Double_t MFormulary::Sigma_el_pp(Double_t *p)
{
  // This function calculates the elastic cross section of a proton-proton-reaction.
  // p specifies the momentum of the incident proton in the laboratory frame.
  // The result is in units of mb, the momentum has to be in GeV.
  // The formula applies only for high energies!

  return 11.9 + 26.9 * Power(*p, -1.21) + 0.169 * Log(*p) * Log(*p) - 1.85 * Log(*p);
}


 Double_t MFormulary::N_charged(Double_t *s) 
{
  // This function calculates the charged multiplicity of a proton-proton collision
  // at high energies. s specifies the center of mass energy squared (!) in GeV^2.

  return 0.88 + 0.44 * Log(*s) + 0.118 * Log(*s) * Log(*s);
}



// true mathematiScs

 void MFormulary::Pq(const Double_t *p, const Double_t *q, Double_t x[2]) 
{
  // Calculates the x(y=0)-values of a parabola of the form
  // x^2 + p * x + q = 0.

  Int_t i;

  Double_t minus_p_halbe = -(*p)/2.;
  Double_t wurzel = Sqrt(-(*q) + Power((*p)/2., 2));

  for (i = 0; i < 2; i++) {
    x[i] = minus_p_halbe + Power(-1., i) * wurzel;
  }
  
  return;
}


 Double_t MFormulary::Angle(const Double_t *x1, const Double_t *x2, Int_t dim = 3)
{
  // Returns the angle between to vectors (in dim dimensions).

  return TMath::ACos(ScalarProd(x1, x2, dim)/(Abs(x1, 3) * Abs(x2, dim)));
}


 Double_t MFormulary::Dist_squared(const Double_t *p1, const Double_t *p2, Int_t dim = 3)
{
  // Calculates the distance squared of two points in euclidian geometry.

  Int_t i;
  Double_t minus[dim];

  for (i=0; i<dim; i++) {
    minus[i] = p2[i] - p1[i];
  }

  return Square(minus, dim);
}


 Double_t MFormulary::Dist(const Double_t *p1,  const Double_t *p2, Int_t dim = 3)
{
  // Calculates the distance of two points in euclidian geometry.
  
  return Sqrt(Dist_squared(p1, p2, dim));
}


 void MFormulary::VectorProd(const Double_t *a, const Double_t *b, Double_t *c)
{
  // Calculates the (euclidian) vector product of two vectors in three dimensions.

  c[1] = a[2]*b[3] - a[3]*b[2];
  c[2] = a[3]*b[1] - a[1]*b[3];
  c[3] = a[1]*b[2] - a[2]*b[1];
  
  return;
}


 Double_t MFormulary::ScalarProd(const Double_t *a, const Double_t *b, Int_t dim = 3)
{
  // Calculates the (euclidian) scalar product of two vectors

  Int_t i;
  Double_t result = 0.;

  for (i=0; i<dim; i++) {
    result += a[i] * b[i];
  }

  return result;
}


 Double_t MFormulary::Square(const Double_t *p, Int_t dim = 3)
{
  // Calculates the scalar product of the given vector with itself.

  return ScalarProd(p, p, dim);
}


 void MFormulary::Sum(const Double_t *p1, const Double_t *p2, Double_t *P, Int_t n = 3)
{
  // Calculates the sum of the two n-dimensional vectors and stores the 
  // result in the vector P.

  Int_t i;

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


 void MFormulary::Diff(const Double_t *p1, const Double_t *p2, Double_t *P, Int_t n = 3)
{
  // Calculates the difference of the two n-dimensional vectors and stores the 
  // result in the vector P.

  Int_t i;

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


 Double_t MFormulary::Abs(const Double_t *p, Int_t dim = 3)
{
  // Calculates the length of a vector in ordinary euclidian geometry.
    
  return Sqrt(Square(p, dim));
}