//////////////////////////////////////////////////////////////////////////
//									//
// 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 final electron phi-, and energy of	//
// the final electron E-.  The returned value is the differential cross	//
// section 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-) = 1/2 (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);

Double_t Pairs(Double_t *var, Double_t *par)
{
   Double_t kin=par[0];
   Double_t Eele=par[1];
   Double_t phie=par[2];
   Double_t Mpair=par[3];
   Double_t qR2=par[4]=var[0];
   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 << "<" << endl;
      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 << "." << endl;
      return 0;
   }

   TThreeVectorReal p;
   gIn.SetMom(p.SetPolar(kin,0,0));
   Double_t Aquad=Eele*(kin-qRecoil[3])*sqr(mElectron/kin);
   Double_t Bquad=2*Eele*qRecoil.Rho()*cos(phie-phiR)*(mElectron/kin);
   Double_t Cquad=qRecoil.LengthSqr()-2*Epos*qRecoil[3]
                                     +kin*sqr(mElectron)/Eele;
   Double_t Discrim=sqr(Bquad)-4*Aquad*Cquad;
   if (Discrim < 0) {
      cout << "*" << endl;
      return 0;
   }
   Double_t thetae_1=(-Bquad+sqrt(Discrim))/(2*Aquad);
   Double_t thetae_2=(-Bquad-sqrt(Discrim))/(2*Aquad);
   Double_t thetae=thetae_1;
   if (thetae < 0) {
      thetae = (-Bquad-sqrt(Discrim))/(2*Aquad);
      if (thetae < 0) {
         cout << "!" << endl;
         return 0;
      }
   }
   thetae *= (mElectron/kin);
   eOut.SetMom(p.SetPolar(Eele,thetae,phie));
   p = gIn.Mom()-qRecoil-p;
   pOut.SetMom(p);

   // Set the initial,final polarizations
   gIn.SetPol(zeroVector);
   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));

   // Instead of dressing the cross section up with the form factor, here I
   // actually integrate the form factor squared times the differential cross
   // section over recoil momentum d^3 q.  The range of the integral is over
   // the range (qmin,infinity) for the longitudinal component, and over the
   // infinite plane for the transverse components.  I found that the pair
   // production cross section is very flat in the transverse components of
   // qRecoil if I rescale it with the empirical factor
   //      sqr(qRecoil.LengthSqr()) / (sqr(qRecoil[1]) + sqr(qRecoil[2]))
   // This allows me to do the integral over the transverse components of
   // qRecoil analytically, as follows.
   //
   //           /  2              2    2   -4
   //     FFI = | d q  | 1-FF(q) |  ( q   q  )
   //           /                      T
   //             
   // where q is the momentum transfer 3-vector (time-like component is zero),
   // qL is its longitudinal component, and Q0 = m*m*E0/(E1*E1+E2*E2) with m
   // the mass of the electron, E0 the initial gamma energy, and E1,E2 the
   // final electron and positron energies.
   //
   // This integral is readily calculated analytically using the above dipole
   // form factor, leading to the following result.  

   //result *= sqr(qRecoil.LengthSqr()/qRecoil.Rho());
   //result *= sqr(qRecoil[3]);
   Double_t r=pOut.Mom()[0]/eOut.Mom()[0];
   Double_t weight=(qR2*r*2e-6)+sqr(sqr(Mpair)-qR2*r);
   weight *= qR2-sqr(qRmin);
   if (par[6] == 0) {
      return result*weight;
   }
   else {
      return 1/(weight+1e-30);
   }
   result *= sqr(qRecoil.LengthSqr());
   result /= 1/qRecoil[3]+Q0/sqr(qRecoil[3]);
   result *= FFI;
   result *= 2*PI_;

   // 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.,
                Double_t wgt=0)
{
   c1 = new TCanvas("c1","Pair Production Rate",200,10,700,500);
   f1 = new TF1("f1",Pairs,0,qR2,7);
   Double_t params[7];
   params[0] = E0;
   params[1] = Eele;
   params[2] = phie;
   params[3] = Mpair;
   params[4] = qR2;
   params[5] = phiR;
   params[6] = wgt;
   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->SetParameter(6,params[6]);
   f1->Draw();
   c1->Update();
//   cout << "Total cross section is " << f1->Integral(0,E0,params,0.1)
//        << " microbarns" << endl;
}