// requires ROOT: http://root.cern.ch/
//
//    usage:            compile using ROOT by: ".L library_kin.C+" (note the plus sign)
//
//      By:             Jaideep Singh
//
//      Date:           02/18/07
//
//      descp:          Kinematic Variable Calculator
//
//                      To use this library, place the following line in your macro to load the lib:
//
//                              gSystem->Load("library_kin_C.so");
//
//                      Make sure that it has been compiled as described in the usage above.
//                      Make sure that the path is correct in the "Load" command
//
//                      converts between mass, energy, momentum, W, nu, Q^2, and beta
//                      also calculates the minimum and maximum values of nu, Q^2, W, and final energy
// 
//      reference:      (unfinished technote on kinematics)   

// head files needed to compile macro
#include "TMath.h"
#include "Riostream.h"

// relativistic momentum
Double_t kin_p(Double_t E , Double_t M)
{ return TMath::Sqrt(E*E-M*M); }

// relativistic speed wrt c
Double_t kin_beta(Double_t p , Double_t E)
{ return p/E; }

// relativistic energy
Double_t kin_E(Double_t p , Double_t M)
{ return TMath::Sqrt(p*p+M*M); }

// invariant mass
Double_t kin_M(Double_t E , Double_t p)
{ return TMath::Sqrt(E*E-p*p); }

// Es = incident electron
// Ep = scattered electron
// exact formula for Q^2 for electron scattering
Double_t kin_Q2(Double_t Es, Double_t Ep, Double_t theta, Double_t m)
{
        Double_t ps = kin_p(Es,m);
        Double_t pp = kin_p(Ep,m);
        return 2.0*(Es*Ep-ps*pp*TMath::Cos(theta)-m*m);
}

// energy lost by incident electron
Double_t kin_nu(Double_t Es , Double_t Ep)
{ return Es-Ep; }

// energy of scattered electron
Double_t kin_Ep(Double_t Es , Double_t nu)
{ return Es-nu; }

// invariant mass of the product
Double_t kin_W(Double_t M , Double_t nu , Double_t Q2)
{ return TMath::Sqrt(M*M + 2.0*nu*M - Q2); }

// energy lost by incident electron from W
Double_t kin_W2nu(Double_t M , Double_t m , Double_t W , Double_t Es , Double_t theta)
{ 
        Double_t ps = kin_p(Es,m);
	Double_t wm = W*W-M*M;
	Double_t s = TMath::Sin(theta);
	Double_t A = M*M+m*m+2.0*Es*M+ps*ps*s*s;
	Double_t B = ((M+Es)*wm+2.0*ps*ps*(M+Es*s*s));
	Double_t C = wm*wm/4+ps*ps*wm+ps*ps*ps*ps*s*s;
	Double_t D = B*B-4.0*A*C;
	Double_t nu = -Es;
	if (D < 0.0)
	{	cerr << "(kin_W2nu) this quadratic has complex roots, returning garbage:" << nu << endl;	}
	else if (D == 0.0)
	{	nu = B/2.0/A;	}
	else
	{
		D = TMath::Sqrt(D);
		nu = 2.0*C/(D+B);
	}
	return nu;
}

// energy loss from the incident electron from Q^2
Double_t kin_Q22nu(Double_t Q2 , Double_t Es , Double_t theta , Double_t m)
{ 
	Double_t m2 = m*m;
	Double_t Es2 = Es*Es;
        Double_t s = TMath::Sin(theta);
        Double_t c = TMath::Cos(theta);
	Double_t ps = kin_p(Es,m);
	Double_t a = 2.0*m2/Q2;
	Double_t b = 2.0*(Es2*s*s+m2*c*c)/Q2;
	Double_t one = (1.0+a)*(1.0+a)-a*b;
	Double_t Ep = Es;
	Double_t nu = -Es;
	if (one < 0.0) 
	{	cerr << "(kin_Q22nu) square root of negative number, returning garbage: " << nu << endl;	}
	else
	{ 
		Ep *= (1.0+a+ps/Es*c*TMath::Sqrt(one))/b;
		nu = kin_nu(Es,Ep);
	}
	return nu;
}

// minimum energy lost by incident electron, elastic scattering with target mass recoil
Double_t kin_numin(Double_t M , Double_t m , Double_t Es , Double_t theta)
{	
        Double_t ps = kin_p(Es,m);
	Double_t c = TMath::Cos(theta);
	Double_t s = TMath::Sin(theta);
	Double_t s2 = TMath::Sin(theta/2.0);
	Double_t one = TMath::Sqrt(1.0-m*m*s*s/M/M);
	Double_t top = (1.0+c)*2.0*ps*ps/M/Es*s2*s2;
	Double_t bot = 1.0 + c*one + (1.0+c)*2*Es/M*s2*s2;
	return Es*top/bot;
}


// maximum energy lost by incident electron, assuming it is not annihilated
Double_t kin_numax(Double_t Es , Double_t m)
{ return Es-m; }

// minimum Q^2 depends only on the incident energy and incident particle mass
Double_t kin_Q2min(Double_t Es, Double_t m)
{
        return 2.0*m*(Es-m);
}

// maximum Q^2 for a given beam energy and angle occurs for elastic scattering
Double_t kin_Q2max(Double_t Es, Double_t M , Double_t theta, Double_t m)
{
	Double_t numin = kin_numin(M , m , Es , theta);
	return 2.0*M*numin;
}

// minimum invariant mass of the product
Double_t kin_Wmin(Double_t M)
{ 	return M; }

// maximum invariant mass of the product
Double_t kin_Wmax(Double_t M , Double_t Es , Double_t m)
{ 	return TMath::Sqrt(M*M + 2.0*(Es-m)*(M-m)); }

// minimum energy of scattered electron
Double_t kin_Epmin(Double_t m)
{	return m; }

// maximum energy of scattered electron
Double_t kin_Epmax(Double_t Es , Double_t theta, Double_t M , Double_t m)
{
	Double_t c = TMath::Cos(theta);
	Double_t s = TMath::Sin(theta);
	Double_t s2 = TMath::Sin(theta/2.0);
	Double_t one = TMath::Sqrt(1.0-m*m*s*s/M/M);
	Double_t top = 1.0 + c*one + (1.0+c)*2*m*m/M/Es*s2*s2;
	Double_t bot = 1.0 + c*one + (1.0+c)*2*Es/M*s2*s2;
	return Es*top/bot;
}