//////////////////////////////////////////////////////////////////////////
//									//
// Pairs.C								//
//									//
// Calculates the e+e- pair production rate including polarization	//
// effects for a given kinematics.  The kinematics are set by the       //
// initial photon energy kin, the mass of the e+e- pair M, the recoil	//
// momentum squared qR^2, the azimuthal angle of the recoil momentum	//
// phiR, the azimuthal angle of the outgoing electron about the pair	//
// momentum axis phi-, and energy of the outgoing electron E-.  The	//
// returned value is the differential cross section measured in		//
// microbarns/GeV^4/r, differential in (d^3 qR  dphi- dE-).  Another 	//
// useful expression for the differential measure is	 		//
//    (d^3 qR dphi- dE-) = (M / 2 kin) (dM dqR^2 dphiR dphi- dE-)	//
// The cross section contains the form factor squared that is supposed  //
// to represent the carbon atom.					//
//									//
//////////////////////////////////////////////////////////////////////////

// #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 posXhat(1,0,0);
const TThreeVectorReal negXhat(-1,0,0);
const TThreeVectorReal posYhat(0,1,0);
const TThreeVectorReal negYhat(0,-1,0);
const TThreeVectorReal posZhat(0,0,1);
const TThreeVectorReal negZhat(0,0,-1);

const Double_t PI_=2*atan2(1,0);

TRandom2 random(0);

Double_t Pairs(Double_t *var, Double_t *par)
{
   Double_t kin=par[0];
   Double_t Eele=par[1]=var[0];
   Double_t phie=par[2];
   Double_t Mpair=par[3];
   Double_t qR2=par[4];
   Double_t phiR=par[5];

   Double_t qRmin=sqr(Mpair)/(2*kin);
   Double_t Ppair=kin-qRmin;
   Double_t theta2=(qR2-sqr(qRmin))/(kin*Ppair);
   if (theta2 < 0) {
      cout << "<";
      return 0;
   }
   Double_t theta=sqrt(theta2);
   TThreeVectorReal qRecoil(cos(phiR)*Ppair*theta,
                            sin(phiR)*Ppair*theta,
                            qRmin+theta2*Ppair/2 );

   TPhoton gIn;
   TLepton eOut(mElectron), pOut(mElectron);

   // Solve for the rest of the kinematics
   Double_t Epos = kin-Eele;
   Double_t qmin=sqr(mElectron)*kin/(2*Eele*Epos);
   if (qRecoil[3] < qmin) {
      cout << ".";
      return 0;
   }

   TThreeVectorReal p;
   gIn.SetMom(p.SetPolar(kin,0,0));
   Double_t pStar2=sqr(Mpair/2)-sqr(mElectron);
   if (pStar2 < 0) {
      cout << "*";
      return 0;
   }
   Double_t pStar=sqrt(pStar2);
   Double_t gamma=kin/Mpair;
   Double_t beta=sqrt(1-1/sqr(gamma));
   Double_t costhetaStar=((Eele/gamma)-Mpair/2)/(beta*pStar);
   if (fabs(costhetaStar) > 1) {
      cout << "!";
      return 0;
   }
   Double_t sinthetaStar=sqrt(1-sqr(costhetaStar));
   Double_t p1Long=gamma*(beta*Mpair/2+pStar*costhetaStar);
   Double_t p2Long=gamma*(beta*Mpair/2-pStar*costhetaStar);
   Double_t pTrans=pStar*sinthetaStar;
   p[1] = pTrans*cos(phie)-qRecoil[1]*p1Long/Ppair;
   p[2] = pTrans*sin(phie)-qRecoil[2]*p1Long/Ppair;
   p[3] = p1Long;
   eOut.SetMom(p);
   p = gIn.Mom()-qRecoil-p;
   pOut.SetMom(p);

   // Set the initial,final polarizations
   gIn.SetPol(posXhat);
   eOut.AllPol();
   pOut.AllPol();

   // Multiply the basic cross section by the carbon atomic form factor
   // Here I chose a simple parameterization for the form factor
   const Double_t Z=6;
   const Double_t beta=111*pow(Z,-1/3.)/mElectron;	// ff cutoff in /GeV
   const Double_t Fff=1/(1+sqr(beta)*qRecoil.LengthSqr());
   Double_t result = TCrossSection::PairProduction(gIn,eOut,pOut);
   result *= sqr(Z*(1-Fff));
   return result;

   // The unpolarized Bethe-Heitler cross section is given here for comparison
   Double_t delta=136*mElectron/pow(Z,0.33333) * kin/(Eele*Epos);
   Double_t aCoul=sqr(alphaQED*Z);
   Double_t fCoul=aCoul*(1/(1+aCoul)+0.20206-0.0369*aCoul
                       +0.0083*pow(aCoul,2)-0.002*pow(aCoul,3));
   Double_t xsi=log(1440/pow(Z,0.66667))/(log(183/pow(Z,0.33333)-fCoul));
   Double_t FofZ=(8./3.)*log(Z) + (kin < 0.05)? 0 : 8*fCoul;
   Double_t Phi1=20.867-3.242*delta+0.625*sqr(delta);
   Double_t Phi2=20.209-1.930*delta-0.086*sqr(delta);
   Double_t Phi0=21.12-4.184*log(delta+0.952);
   if (delta > 1) {
      Phi1=Phi2=Phi0;
   }
   result = hbarcSqr/sqr(mElectron)*pow(alphaQED,3)/kin
            * Z*(Z+xsi)
            * (
                 (sqr(Eele)+sqr(Epos))/sqr(kin)*(Phi1-FofZ/2)
                 + (2./3.)*(Eele*Epos)/sqr(kin)*(Phi2-FofZ/2)
              );
}

Int_t demoPairs(Double_t E0=9.,
                Double_t Eele=4.5,
                Double_t phie=0,
                Double_t Mpair=2e-3,
                Double_t qR2=1e-6,
                Double_t phiR=0.)
{
   c1 = new TCanvas("c1","Pair Production Rate",200,10,700,500);
   f1 = new TF1("f1",Pairs,0,Eele,6);
   Double_t params[6];
   params[0] = E0;
   params[1] = Eele;
   params[2] = phie;
   params[3] = Mpair;
   params[4] = qR2;
   params[5] = phiR;
   f1->SetParameter(0,params[0]);
   f1->SetParameter(1,params[1]);
   f1->SetParameter(2,params[2]);
   f1->SetParameter(3,params[3]);
   f1->SetParameter(4,params[4]);
   f1->SetParameter(5,params[5]);
   //f1->DrawCopy("same");
   f1->Draw();
   c1->Update();
}

Int_t genPairs(Int_t N, Double_t kin=9., TFile *hfile=0, TTree *tree=0, Int_t prescale=1000)
{
   struct event_t {
      Double_t E0;
      Double_t Eele;
      Double_t phie;
      Double_t Mpair;
      Double_t qR2;
      Double_t phiR;
      Double_t diffXS;
      Double_t weight;
      Double_t weightedXS;
   } event;
   TString leaflist("E0/D:Eele/D:phie/D:Mpair/D:qR2/D:phiR/D:diffXS/D:weight/D:weightedXS");
   event.E0 = kin;

   if (hfile != 0) {
      if (tree == 0) {
         TString title;
         title.Form("e+e- pair production data, Egamma=%f",event.E0);
         tree = new TTree("epairXS",title);
      }
      tree->Branch("event",&event,leaflist,65536);
   }

   Double_t sum=0;
   for (int n=1; n<=N; n++) { 
      event.weight = 1;

      // generate Eele uniform on [0,E0]
      event.Eele = random.Uniform(event.E0);
      event.weight *= event.E0;
   
      // generate phie uniform on [0,2pi]
      event.phie = random.Uniform(2*PI_);
      event.weight *= 2*PI_;

      // generate phiR uniform on [0,2pi]
      event.phiR = random.Uniform(2*PI_);
      event.weight *= 2*PI_;
   
      // generate Mpair with weight (M0/M)^3
      Double_t M0=2*mElectron;
      event.Mpair = M0/sqrt(random.Uniform(1.));
      event.weight *= pow(event.Mpair,3)/(2*M0*M0);

      // generate qR2 with weight 1/(q02 + qR2)^2
      Double_t q02=sqr(5e-5); // 50 keV/c cutoff parameter
      Double_t u=random.Uniform(1);
      event.qR2 = q02*(1-u)/(u+1e-50);
      event.weight *= sqr(event.qR2+q02)/q02;

      // overall measure Jacobian factor
      event.weight *= event.Mpair/(2*event.E0);

      Double_t *par=&event;
      Double_t *var=&event.Eele;
      event.diffXS = Pairs(var,par);
      event.weightedXS = event.diffXS*event.weight;

      if (tree != 0) {
         tree->Fill();
      }

      sum += event.weightedXS;
      if (n/prescale*prescale == n) {
         cout << sum/n << endl;
      }
   }

   if (tree != 0) {
      tree->FlushBaskets();
   }
   if (hfile != 0) {
      hfile->Write();
   }
}