/* =============================================================
   The 1-loop and 2-loop continuum gg -> gamma gamma amplitudes.
   with light quark masses and top quark contribution neglected
   ============================================================= */

#include "ATOOLS/Math/MathTools.H"

using namespace ATOOLS;

#include "gggamgam1l.C"  // gggamgam 1-loop
                         // amplitudes, with and without phases.
                         // as another contributing process, at NNLO
#include "gggamgam2l.C"  // gggamgam 2-loop amplitudes (finite parts)
                         // with phases removed.  Also reads in 1-loop
                         // amplitudes, with and without phases.
                         // as another contributing process, at NNLO

/* =============================================================
   The 1-loop continuum gg -> gamma gamma amplitudes.
   The light quark (u,d,c,s,b) masses included
   ============================================================= */

// NEW helicity LABELLING:
// hel = 1 means ++ -> ++ (or equivalently -- -> --) or "mmpp"
// hel = 2 means ++ -> -- (or equivalently -- -> ++) or "pppp" 
// Other helicity configurations are irrelevant for Higgs interference

// OLD helicity LABELLING:   pppp = 1, mppp = 2, mmpp = 3, mpmp = 4, mppm = 5
// Unless specified defaulting to new helicity labelling

/* The one-loop light-by-light amplitudes */

// Basic scalar integrals, with a uniform mass m in the loop:

Complex B0(double m, double s){
 Complex temp; 
 double x;
 if (Abs(s/4./m/m) < 0.01)  // use large mass expansion here:
     { return s/6./m/m * ( 1. + 1./10. * s/m/m * (
			   1. + 1./7. * s/m/m * ( 
			   1. + 1./6. * s/m/m ))) ; }
 if (Abs(s/4./m/m) > 300.0) // use small mass expansion here:
   { if (s < 0) { return 2. - log(-s/m/m)
		    + 2.*m*m/s * ( 1. + log(-s/m/m) )
                    + m*m*m*m/s/s * ( -1. + 2. * log(-s/m/m) ) ; }
     else { return 2. - log(s/m/m)
		    + 2.*m*m/s * ( 1. + log(s/m/m) )
	            + m*m*m*m/s/s * ( -1. + 2. * log(s/m/m) ) 
	    + I * PI_DEF * ( 1. - 2.*m*m/s - 2.*m*m*m*m/s/s ) ; }
   }
 if (s < 0) 
    { x = sqrt(1.-4.*m*m/s); 
      temp = 2. + x * log((x-1.)/(x+1.)) ; }
 if (s >= 0 && s < 4.*m*m) 
      { x = sqrt(-1.+4.*m*m/s) ;
        temp = real( 2. + I * x * log((I*x-1.)/(I*x+1.)) ) ; }
 if (s >= 4.*m*m)
    { x = sqrt(1.-4.*m*m/s); 
      temp = 2. + x * ( log((1.-x)/(x+1.)) + I * PI_DEF ) ; }
 return temp; 
}

Complex C0(double m, double s){
 Complex temp, temp1; 
 double x;
 if (Abs(s/4./m/m) < 0.01)    // use large mass expansion here:
     { return -0.5/m/m * ( 1. + 1./12. * s/m/m * (
			   1. + 2./15. * s/m/m * ( 
			   1. + 9./56. * s/m/m ))) ; }
 if (Abs(s/4./m/m) > 300.0) // use small mass expansion here:
   { if (s < 0) { temp1 = log(-s/m/m) - 2.*m*m/s - 3.*m*m*m*m/s/s 
			       - 20./3.*m*m*m*m*m*m/s/s/s ; 
                  return 0.5/s * temp1 * temp1 ; }
     else { temp1 = log(s/m/m) - 2.*m*m/s - 3.*m*m*m*m/s/s 
			       - 20./3.*m*m*m*m*m*m/s/s/s - I * PI_DEF ;
            return 0.5/s * temp1 * temp1 ; }
   }
 if (s < 0) 
    { x = sqrt(1.-4.*m*m/s); 
      temp1 = log((x+1.)/(x-1.));
      temp = 0.5/s * temp1 * temp1 ; }
 if (s > 0 && s < 4.*m*m) 
      { x = sqrt(-1.+4.*m*m/s) ;
        temp1 = log((I*x-1.)/(I*x+1.)) ;
        temp = 0.5/s * temp1 * temp1 ; }
 if (s >= 4.*m*m)
    { x = sqrt(1-4.*m*m/s); 
      temp1 = log((1.-x)/(x+1.)) + I * PI_DEF ;
      temp = 0.5/s * temp1 * temp1 ; }
 return temp; 
}

Complex D0(double m, double s, double t){
 Complex bu, temp; 
 double u, v, a, bv, temp_re, temp_im;
 // Complex u, v, a, bv, temp_re;
 if (t > s) { return D0(m,t,s); }
 u = s/4./m/m;   
 v = t/4./m/m;
 if ((Abs(u) < 0.01) && (Abs(v) < 0.01))   // use large mass expansion here:
   { return 1./m/m/m/m * ( 1./6. + 1./15.*(u+v) 
			   + 2./105. * (2.*(u*u+v*v)+u*v)
                           + 8./945. * (u+v)*(3*(u*u+v*v)-2*u*v) ); }
 if ((Abs(u) > 10000.0) && (Abs(v) > 10000.0)) // small mass expansion here:
 { if (u < 0.) { return 
         ( - (log(s/t)*log(s/t) + PI_DEF*PI_DEF) 
           + 2.*s * C0(m,s) + 2.*t * C0(m,t) )/ (s*t-4.*m*m*(s+t)) ; }
   else return ( - ( log(-s/t)*log(-s/t) - 2. * I * PI_DEF * log(-s/t) ) 
		 + 2.*s * C0(m,s) + 2.*t * C0(m,t) )/ (s*t-4.*m*m*(s+t)) ; 
 }
 a = sqrt(1.-(u+v)/u/v);
 if ((u<0.) || (u>1.)) { bu = sqrt((u-1.)/u); }
 else { bu = I * sqrt((1.-u)/u); }
 bv = sqrt((v-1.)/v);
 temp_re = 2./s/t/a * ( 
              - ReLi2((a+1.)/(a+bu)) - ReLi2((a+1.)/(a-bu)) 
              + ReLi2((a-1.)/(a+bu)) + ReLi2((a-1.)/(a-bu)) 
              - ReLi2((a+1.)/(a+bv)) - ReLi2((a+1.)/(a-bv)) 
              + ReLi2((a-1.)/(a+bv)) + ReLi2((a-1.)/(a-bv)) ); 
 if (s > 4.*m*m) { 
   temp_im = real( -2./s/t * PI_DEF/a * log(v*(a+bu)*(a+bu)) );
   temp = temp_re + I * temp_im ; 
   return temp; }
 else return temp_re;
}

/*=====================================================================
   One-loop light-by-light helicity amplitudes, 
   from Gounaris, Porfyriadis, Renard, Eur. Phys. J C9, 673 (1999).   
   We use our all-outgoing helicity convention, though, and note that
   t and u should be exchanged vs. GPR conventions.
   Here x = cos(theta).  t = -s/2*(1-x), u = -s/2*(1+x). 
======================================================================*/

// W loop: (old hel. label used)

Complex lbyl_W_loop(int hel, double x, double m_W, double s){
 double t, u, m_W_2, m_W_4;
 Complex temp;
 t = -s/2.*(1.-x);    
 u = -s/2.*(1.+x);
 m_W_2 = m_W*m_W;
 m_W_4 = m_W_2*m_W_2;
 if (hel==1) { temp = // pppp
       - 12. + 24. * m_W_4 * ( D0(m_W,s,t) + D0(m_W,s,u) + D0(m_W,t,u) ) ; }
 else if (hel==2) { temp = 
       - 12. + 24. * m_W_4 * ( D0(m_W,s,t) + D0(m_W,s,u) + D0(m_W,t,u) )
           + 12. * m_W_2 * s * t * u 
                * ( D0(m_W,s,t)/u/u + D0(m_W,s,u)/t/t + D0(m_W,t,u)/s/s )
           - 24. * m_W_2 * ( 1./s + 1./t + 1./u )
                      * ( t * C0(m_W,t) + u * C0(m_W,u) + s * C0(m_W,s) ) ; } 
 else if (hel==3) { temp = // mmpp
         12. - 12. * (1. + 2.*u/s) * B0(m_W,u) 
             - 12. * (1. + 2.*t/s) * B0(m_W,t)
         + 24. * m_W_2 * t * u/s * D0(m_W,u,t)
         + 16. * (1. - 1.5 * m_W_2/s - 0.75 * t * u/s/s )
              * ( 2. * t * C0(m_W,t) + 2. * u * C0(m_W,u) 
                - t * u * D0(m_W,t,u) )
         + 8. * (s-m_W_2) * (s-3.*m_W_2)
                      * (D0(m_W,s,t) + D0(m_W,s,u) + D0(m_W,t,u)) ; }
 else if (hel==4)      // cross s <-> u from mmpp
   { temp = 
          12. - 12. * (1. + 2.*s/u) * B0(m_W,s) 
              - 12. * (1. + 2.*t/u) * B0(m_W,t)
         + 24. * m_W_2 * t * s/u * D0(m_W,s,t)
         + 16. * (1. - 1.5 * m_W_2/u - 0.75 * t * s/u/u )
              * ( 2. * t * C0(m_W,t) + 2. * s * C0(m_W,s) 
                - t * s * D0(m_W,t,s) )
         + 8. * (u-m_W_2) * (u-3.*m_W_2)
                      * (D0(m_W,u,t) + D0(m_W,u,s) + D0(m_W,t,s)) ; }
 else if (hel==5)      // exchange t <-> u, i.e. x <-> -x, in mpmp=4.
   { return lbyl_W_loop(4,-x,m_W,s) ; }
 return temp;
}

// fermion loop: (old hel. label used)

Complex lbyl_f_loop(int hel, double x, double m_f, double s){
 double t, u, m_f_2, m_f_4;
 Complex temp;
 t = -s/2.*(1.-x);    
 u = -s/2.*(1.+x);
 m_f_2 = m_f*m_f;
 m_f_4 = m_f_2*m_f_2;
 if (hel==1) { return - 2./3. * lbyl_W_loop(1,x,m_f,s) ; } // pppp
 else if (hel==2) { return - 2./3. * lbyl_W_loop(2,x,m_f,s) ; }
 else if (hel==3) { temp = // mmpp
           - 8. + 8. * (1. + 2.*u/s) * B0(m_f,u)  
                + 8. * (1. + 2.*t/s) * B0(m_f,t)
           - 8. * ( (t*t+u*u)/s/s - 4.*m_f_2/s ) 
               * ( t * C0(m_f,t) + u * C0(m_f,u) )
           + 8. * m_f_2 * (s-2.*m_f_2) * (D0(m_f,s,t) + D0(m_f,s,u))
           - 4. * ( 4. * m_f_4 - (2.*s*m_f_2 + t*u) * (t*t+u*u)/s/s 
                  + 4. * m_f_2 * t*u/s ) * D0(m_f,t,u) ; }
 else if (hel==4)   // cross s <-> u from mmpp
   { temp = - 8. + 8. * (1. + 2.*s/u) * B0(m_f,s) 
                + 8. * (1. + 2.*t/u) * B0(m_f,t)
           - 8. * ( (t*t+s*s)/u/u - 4.*m_f_2/u ) 
               * ( t * C0(m_f,t) + s * C0(m_f,s) )
           + 8. * m_f_2 * (u-2.*m_f_2) * (D0(m_f,u,t) + D0(m_f,u,s))
           - 4. * ( 4. * m_f_4 - (2.*u*m_f_2 + t*s) * (t*t+s*s)/u/u 
                    + 4. * m_f_2 * t*s/u ) * D0(m_f,t,s) ; }
 else if (hel==5)      // exchange t <-> u, i.e. x <-> -x, in mpmp=4.
   { return lbyl_f_loop(4,-x,m_f,s) ; }
 return temp;
}

// scalar loop: (old hel. label used)

Complex lbyl_s_loop(int hel, double x, double m_sc, double s){
 double t, u, m_s_2, m_s_4;
 Complex temp;
 t = -s/2.*(1.-x);    
 u = -s/2.*(1.+x);
 m_s_2 = m_sc*m_sc;
 m_s_4 = m_s_2*m_s_2;
 if (hel==1) { return 1./3. * lbyl_W_loop(1,x,m_sc,s) ; } // pppp
 else if (hel==2) { return 1./3. * lbyl_W_loop(2,x,m_sc,s) ; }
 else if (hel==3) { temp = // mmpp
             4. - 4. * (1. + 2.*u/s) * B0(m_sc,u) 
                - 4. * (1. + 2.*t/s) * B0(m_sc,t)
           + 8. * m_s_2 * t*u/s * D0(m_sc,t,u)
           - 8. * m_s_2/s * ( 1. + u*t/2./m_s_2/s )
               * ( 2. * t * C0(m_sc,t) + 2. * u * C0(m_sc,u) 
                 - t * u * D0(m_sc,t,u) )
          + 8. * m_s_4 * (D0(m_sc,s,t) + D0(m_sc,s,u) + D0(m_sc,t,u)) ; }
 else if (hel==4)    // cross s <-> u from mmpp
   { temp = 4. - 4. * (1. + 2.*s/u) * B0(m_sc,s) 
               - 4. * (1. + 2.*t/u) * B0(m_sc,t)
           + 8. * m_s_2 * t*s/u * D0(m_sc,t,s)
           - 8. * m_s_2/u * ( 1. + s*t/2./m_s_2/u )
               * ( 2. * t * C0(m_sc,t) + 2. * s * C0(m_sc,s) 
                 - t * s * D0(m_sc,t,s) )
         + 8. * m_s_4 * (D0(m_sc,u,t) + D0(m_sc,u,s) + D0(m_sc,t,s)) ; }
 else if (hel==5)      // exchange t <-> u, i.e. x <-> -x, in mpmp=4.
   { return lbyl_s_loop(4,-x,m_sc,s) ; }
 return temp;
}

/*==================================================
 Add up all Standard Model one-loop contributions to 
 gg -> gamma gamma amplitude.  s is in GeV^2
 (Formula is not reliable at hadronic thresholds, of course.) */

// old hel. label used

Complex A_cont_1l(int hel, double x, double s, double muR){
 // for the heaviest quark we'll use running masses (because we have them);
 // for u,d,s it doesn't matter, so we'll just leave them as m_q(m_q).
 double mt = M_t(muR);  double mb = M_b(muR); double mc = M_c(muR);
 return 1./2. * alpha0 * alpha_s(muR) 
     * ( 4./9. * ( lbyl_f_loop(hel,x,m_u,s) + lbyl_f_loop(hel,x,mc,s)
                 + lbyl_f_loop(hel,x,mt,s) )
       + 1./9. * ( lbyl_f_loop(hel,x,m_d,s) + lbyl_f_loop(hel,x,m_s,s)
	        + lbyl_f_loop(hel,x,mb,s) ) ) ;
}

/* ================================================================
   The continuum g g -> gamma gamma amplitude with only
   two interference relevent helicity configurations */

// new hel. label used

Complex A_c_1l(double M_H, double cth, int hel, double muR){
  int oldhel;
  if (hel==2) oldhel=1; else if (hel==1) oldhel=3;
  return A_cont_1l(oldhel, cth, M_H*M_H, muR);
}