//////////////////////////////////////////////////////////////////////////
//									//
// Compton.C								//
//									//
// Calculates the Compton differential cross section in the rest frame	//
// of the initial electron, assumed to be free.  Units are microbarns   //
// per unit solid angle of the scattered photon.			//
//									//
//////////////////////////////////////////////////////////////////////////

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

#include "constants.h"

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

Double_t Compton(Double_t *var, Double_t *par)
{
   Double_t theta=var[0];
   Double_t kin=par[0];
   Double_t phi=par[1];
   Int_t siggIn=par[3];
   Int_t sigeIn=par[4];

   TLepton eIn(mElectron), eOut(mElectron);
   TPhoton gIn, gOut;

   // Solve for the rest of the kinematics
   Double_t kout = kin/(1+(kin/mElectron)*(1-cos(theta)));
   TThreeVectorReal p;
   gIn.SetMom(p.SetPolar(kin,0,0));
   eIn.SetMom(zeroVector);
   gOut.SetMom(p.SetPolar(kout,theta,phi));
   eOut.SetMom(gIn.Mom()+eIn.Mom()-gOut.Mom());

   // Set the initial,final polarizations
   switch (siggIn) {
    case -1:
      gIn.SetPol(p.SetPolar(1,PI_,0));
      break;
    case  0:
      gIn.SetPol(zeroVector);
      break;
    case +1:
      gIn.SetPol(p.SetPolar(1,0,0));
      break;
    default:
      cout << "Compton.C :"
           << "bad photon helicity specified, default to zero!"
           << endl;
      gIn.SetPol(zeroVector);
   }
   switch (sigeIn) {
    case -1:
      eIn.SetPol(p.SetPolar(1,PI_,0));
      break;
    case  0:
      eIn.SetPol(zeroVector);
      break;
    case +1:
      eIn.SetPol(p.SetPolar(1,0,0));
      break;
    default:
      cout << "Compton.C :"
           << "bad electron helicity specified, default to zero!"
           << endl;
      eIn.SetPol(zeroVector);
   }
   gOut.AllPol();
   eOut.AllPol();

   Double_t result=TCrossSection::Compton(gIn,eIn,gOut,eOut);
   return result;
}

Double_t ComptonAsym(Double_t *var, Double_t *par)
{
   Double_t sigpp,sigpm,sigmp,sigmm;
   Double_t sigp0,sigm0,sig0p,sig0m;
   par[3] = +1;
   par[4] = +1;
   sigpp = Compton(var,par);
   par[3] = +1;
   par[4] = -1;
   sigpm = Compton(var,par);
   par[3] = -1;
   par[4] = +1;
   sigmp = Compton(var,par);
   par[3] = -1;
   par[4] = -1;
   sigmm = Compton(var,par);
   return (sigpm - sigmm) / (sigpm + sigmm);
}

Int_t demoCompton()
{
   c1 = new TCanvas("c1","Compton Cross Section",200,10,700,500);
   TF1 *asym = new TF1("asym",ComptonAsym,0,3.1416,5);
   asym->SetParameter(0,2e-5);
   asym->SetParameter(1,0);
   asym->SetParameter(2,1.2);
   asym->Draw();
   c1->Update();
}

Int_t demoCompton(Double_t Ephot)
{
   TF1 *asym2 = new TF1("asym2",ComptonAsym,0,3.1416,5);
   asym2->SetParameter(0,Ephot);
   asym2->SetParameter(1,0);
   asym2->SetParameter(2,1.2);
   asym2->DrawCopy("SAME");
   c1->Update();
}

Double_t ComptonBackScatter(Double_t *var, Double_t *par)
{
   Double_t kf=var[0];
   Double_t ki=par[0];
   Double_t phi=par[1];
   Double_t E0=par[2];

   TLepton eIn(mElectron), eOut(mElectron);
   TPhoton gIn, gOut;

   // Solve for the rest of the kinematics - in rest frame of reaction
   Double_t P0 = sqrt(E0*E0-mElectron*mElectron);
   Double_t S = mElectron*mElectron+2*ki*(E0+P0);
   Double_t rootS = sqrt(S);
   Double_t kstar = (S-mElectron*mElectron)/(2*rootS);
   TThreeVectorReal k(0,0,-kstar);
   gIn.SetMom(k);
   TThreeVectorReal p(0,0,kstar);
   eIn.SetMom(p);
   Double_t gamma = (E0+ki)/rootS;
   Double_t eta = (P0-ki)/rootS;
   Double_t costheta = (1-kf/(gamma*kstar))*gamma/eta;
   if (fabs(costheta) > 1.0) { return 0; }
   Double_t theta = acos(costheta);
   gOut.SetMom(k.SetPolar(kstar,theta+PI_,phi));
   eOut.SetMom(p.SetPolar(kstar,theta,phi));

   // Set the initial,final polarizations
   gIn.SetPol(zeroVector);
   eIn.SetPol(zeroVector);
   gOut.AllPol();
   eOut.AllPol();

   Double_t result=TCrossSection::Compton(gIn,eIn,gOut,eOut);

   result*=(2*PI_)/(eta*kstar);	// convert cross section from d(sigma)/d(Omega*)
                               	// to d(sigma)/d(kf)

   Double_t P=600;	  // laser power in Watts (peak times duty factor)
   Double_t G=250;	  // laser cavity gain factor
   Double_t tau=2e-12;	  // laser pulse length (s) times crossing factor
   Double_t rC=10e-6;	  // neck radius of cavity beam (m)
   Double_t rB=10e-6;	  // electron beam radius (m)
   Double_t I=1.0e-6;	  // electron beam current (A)
   Double_t L=0.0450;	  // effective length of cavity (m)
   Int_t    N=2; 	  // number of passes through beam
   Int_t    pulsed=1;	  // indicate whether laser is pulsed

   result*=1e-34;		// convert from microbarns to m^2
   result*=N*P*G/(ki*1.6e-10);	// ki is converted from GeV to Joules
   result/=2*PI_*(rB*rB+rC*rC);	// area of overlap region
   result*=I/1.6e-19;		// rate of electrons in beam
   if (pulsed)
     result*=tau;		// duration of pulse crossing (lab frame)
   else
     result*=L/3e8;		// time spent in crossing region (lab frame)

   return result;
}

Int_t demoComptonBackScatter()
{
   c2 = new TCanvas("c2","Compton Backscatter Cross Section",200,10,700,500);
   f2 = new TF1("f2",ComptonBackScatter,0,12,3);
   f2->SetParameter(0,4.8e-9);
   f2->SetParameter(1,0);
   f2->SetParameter(2,12);
   f2->Draw();
   c2->Update();
}