//////////////////////////////////////////////////////////////////////////
//									//
// Qpolar.C								//
//									//
// Calculates the quantum amplitude for the scattering of a polarized 	//
// virtual photon from a polarized quark.  The final state consists of 	//
// only the scattered quark.  The photon energy is adjusted to keep the	//
// final quark on-shell. The mass of the quark is adjustable.		//
//									//
//////////////////////////////////////////////////////////////////////////

// #include "TPhoton.h"
// #include "TLepton.h"
// #include "TCrossSection.h"

#include "constants.h"
#include <complex>

inline Double_t sqr(Double_t x) { return x*x; }
inline Complex_t sqr(Complex_t x) { return x*x; }

const TThreeVectorReal zeroVector(0,0,0);
const TThreeVectorReal posZhat(0,0,1);
const TThreeVectorReal negZhat(0,0,-1);

Complex_t Qpolar(Double_t *var, Double_t *par)
{
   Double_t thetai=var[0];  // angle of incident quark wrt the photon axis (r)
   Double_t pin=par[0];     // momentum of the incident quark [GeV/c]
   Double_t Q2=par[1];      // Q^2 of the virtual photon [GeV^2]
   Double_t sigeIn=par[2];  // polarization state of initial quark

   const double mQuark(0);

   TLepton qIn(mQuark), qFi(mQuark);
   TPhoton gIn;

   // Solve for the rest of the kinematics
   Double_t Ein=sqrt(pin*pin+mQuark*mQuark);
   TFourVectorReal pIn(Ein,pin*sin(thetai),0,-pin*cos(thetai));
   Double_t A=pIn[0]*pIn[0]-pIn[3]*pIn[3];
   Double_t B=-Q2*pIn[3];
   Double_t C=-Q2*(pIn[0]*pIn[0]+Q2/4);
   Double_t D=B*B-4*A*C;
   Double_t qin1=(D > 0)? (-B-sqrt(D))/(2*A) : 0;
   Double_t qin2=(D > 0)? (-B+sqrt(D))/(2*A) : 0;
   TFourVectorReal phIn(0,0,0,0);
   if ((qin1 > 1e-6) && (qin1*qin1 > Q2)) {
     phIn = TFourVectorReal(sqrt(qin1*qin1-Q2),0,0,qin1);
   }
   else if ((qin2 > 1e-6) && (qin2*qin2 > Q2)) {
     phIn = TFourVectorReal(sqrt(qin2*qin2-Q2),0,0,qin2);
   }
   TFourVectorReal pFi(pIn+phIn);
   pFi[0] = sqrt(mQuark*mQuark+pFi[1]*pFi[1]+pFi[2]*pFi[2]+pFi[3]*pFi[3]);
   qIn.SetMom(pIn);
   qFi.SetMom(pFi);
   gIn.SetMom(pFi-pIn);

   // Set the initial,final polarizations
   switch (sigeIn) {
    case -1:
      qIn.SetPol(negZhat);
      qFi.SetPol(negZhat);
      break;
    case  0:
      qIn.SetPol(zeroVector);
      qFi.SetPol(zeroVector);
      break;
    case +1:
      qIn.SetPol(posZhat);
      qFi.SetPol(posZhat);
      break;
    default:
      cout << "Qpolar.C :"
           << "bad quark helicity specified, default to zero!"
           << endl;
      qIn.SetPol(zeroVector);
      qFi.SetPol(zeroVector);
   }

   // Compute the current amplitude
// qIn.Print();
// qFi.Print();
// gIn.Print();
   TDiracMatrix qmat1,qmat2;
   qmat1.SetUUbar(qIn.Mom(),qIn.SDM());
   qmat2.SetUUbar(qFi.Mom(),qFi.SDM());
   TDiracMatrix Aslash,AslashStar;
   AslashStar.Slash(gIn.EpsStar(1));
   Aslash.Slash(gIn.Eps(1));
   TDiracMatrix cur;
   cur = AslashStar*qmat2*Aslash*qmat1;
// qmat1.Print();
// qmat2.Print();
// Aslash.Print();
// AslashStar.Print();
// cur.Print();

// double ans;
// cout << "continue? ";
// cin >> ans;

   return cur.Trace();
}

Double_t reQpolar(Double_t *var, Double_t *par)
{
   Complex_t amp = Qpolar(var,par);
   return amp.real();
}

Double_t imQpolar(Double_t *var, Double_t *par)
{
   Complex_t amp = Qpolar(var,par);
   return amp.imag();
}

Double_t absQpolar(Double_t *var, Double_t *par)
{
   Complex_t amp = Qpolar(var,par);
   return abs(amp);
}

Int_t demoQpolar()
{
   c1 = new TCanvas("c1","quark scattering amplitude",200,10,700,500);
   f1 = new TF1("f1",absQpolar,0,2,3);
   f1->SetParameter(0,1.0);
   f1->SetParameter(1,4.0);
   f1->SetParameter(2,+1);
   f1->Draw();
   c1->Update();
   c2 = new TCanvas("c2","quark scattering amplitude",200,10,700,500);
   f2 = new TF1("f2",reQpolar,0,2,3);
   f2->SetParameter(0,1.0);
   f2->SetParameter(1,4.0);
   f2->SetParameter(2,+1);
   f2->SetLineColor(2);
   f2->Draw();
   c2->Update();
   c3 = new TCanvas("c3","quark scattering amplitude",200,10,700,500);
   f3 = new TF1("f3",imQpolar,0,2,3);
   f3->SetParameter(0,1.0);
   f3->SetParameter(1,4.0);
   f3->SetParameter(2,+1);
   f3->SetLineColor(3);
   f3->Draw();
   c3->Update();
   c4 = new TCanvas("c4","quark scattering amplitude",200,10,700,500);
   f4 = new TF1("f4",absQpolar,0,2,3);
   f4->SetParameter(0,1.0);
   f4->SetParameter(1,4.0);
   f4->SetParameter(2,-1);
   f4->Draw();
   c4->Update();
   c5 = new TCanvas("c5","quark scattering amplitude",200,10,700,500);
   f5 = new TF1("f5",reQpolar,0,2,3);
   f5->SetParameter(0,1.0);
   f5->SetParameter(1,4.0);
   f5->SetParameter(2,-1);
   f5->SetLineColor(2);
   f5->Draw();
   c5->Update();
   c6 = new TCanvas("c6","quark scattering amplitude",200,10,700,500);
   f6 = new TF1("f6",imQpolar,0,2,3);
   f6->SetParameter(0,1.0);
   f6->SetParameter(1,4.0);
   f6->SetParameter(2,-1);
   f6->SetLineColor(3);
   f6->Draw();
   c6->Update();
}