// Specification file for classes etc. 
// used for "direct" numerical integration of one loop Feynman diagrams.
// D.E. Soper
// 11 December 2008

using namespace std;
#include <complex>
const int maxsize = 30;
const int maxperms = 5041; // This is 7! + 1.

void giveversion(); // Prints version and date of subroutines.

class DESrandom
{
public:
  DESrandom();                    // Default constructor. (Sets seed = 0.)
  DESrandom(int seed);            // Constructor. Sets seed.   
  bool reSetSeed(/*in*/int seed); // Precondition : 0<=seed<259200
                                  // Post Condition :seed is reset, return 0 if seed is invalid.
  const double& random();         // return random value.

private:
  int lbit;
  double fac;
  int j;
  int ir[250];
  double rr[250];
  void newran();  
};

class DESreal4vector
{
public:
  DESreal4vector();				                                   // Default constructor
  DESreal4vector(double v0, double v1, double v2, double v3);  // Initializes to (v0,v1,v2,v3).
  double& in(const int& mu);						                     // v.in(mu) for input, v.in(mu) = r
  const double& operator()(const int& mu) const;	           // v(mu) for output,  r = v(mu)
  DESreal4vector operator+ (const DESreal4vector& v2) const; // Sum of two vectors: v3 = v1 + v2
  DESreal4vector operator- (const DESreal4vector& v2) const; // Difference of two vectors: v3 = v1 - v2
  DESreal4vector operator* (const double& a) const;		       // Vector times scalar: v2 = v1*c
  DESreal4vector operator- () const;		                     // -1 times vector: v2 = - v1
private:
  double vector[4];				                                   // The vector, v[0], v[1], v[2], v[3]
};

// Scalar times vector: v2 = c*v1
DESreal4vector operator*(const double& a, const DESreal4vector& v);
// Dot product of two 4-vectors: dot(v1,v2)
double dot(const DESreal4vector& v1, const DESreal4vector& v2);
// Square of one 4-vector: square(v) = dot(v,v) 
double square(const DESreal4vector& v); 
// Boost a vector v with boost that transforms vold to vnew.
DESreal4vector boost(const DESreal4vector& vnew, const DESreal4vector& vold, const DESreal4vector& v);

class DESindex
{
public:
  DESindex();                         // Default constructor with nverts = 1 000 000.
  DESindex(int nverts);               // Constructor, needs nverts.
  int operator()(int i) const;        // Returns i mod nverts.
  int getnverts() const;              // Returns  nverts.
private:
  int nverts;                         // Number of vertices in the loop;
};

class DEScmplx4vector
{
public:
  DEScmplx4vector();				                                   // Default (and only) constructor
  complex<double>& in(const int& mu);						               // v.in(mu) for input, v.in(mu) = r
  const complex<double>& operator()(const int& mu) const;	     // v(mu) for output,  r = v(mu)
  DEScmplx4vector operator+ (const DEScmplx4vector& v2) const; // Sum of two vectors: v3 = v1 + v2
  DEScmplx4vector operator- (const DEScmplx4vector& v2) const; // Difference of two vectors: v3 = v1 - v2
  DEScmplx4vector operator* (const complex<double>& a) const;  // Vector times scalar: v2 = v1*c
  DEScmplx4vector operator- () const;		                       // -1 times vector: v2 = - v1
private:
  complex<double> vector[4];				                           // The vector, v[0], v[1], v[2], v[3]
};

// Scalar times vector: v2 = c*v1
DEScmplx4vector operator*(const complex<double>& a, const DEScmplx4vector& v); 
// Dot product of two 4-vectors: dot(v1,v2) 
complex<double> dot(const DEScmplx4vector& v1, const DEScmplx4vector& v2); 
// Square of one 4-vector: square(v) = dot(v,v) 
complex<double> square(const DEScmplx4vector& v); 

// Function to generate the trace of v1slash*v2slash*...*vnslash.
complex<double> DiracTrace(DEScmplx4vector v[2*maxsize], int n);

class DESloopmomenta
{
public:
  DESloopmomenta(int seedIN, int nvertsIN, DESreal4vector pIN[], 
              int indexAIN, double alphaIN[]);      // Constructor. (There is no default constructor.)
  void givenewgraph(DESreal4vector pIN[], int indexAIN) ; // Give data for new graph.
  void newell(DESreal4vector& ellOUT);           // New point ell.
  double rhoell();                               // Density rho(ell).
  double getMTsq();                           // Scale for double parton scattering singularity.
  int getindexA();                            // Initial state vertices are between propagators  
                                              //    nverts and 1 and between indexA and indexA + 1.
  void getQs(DESreal4vector Q[]);             // The vectors Q[i].
  void getps(DESreal4vector p[]);             // The vectors p[i].
  void givepoint(DESreal4vector ellIN);       // Supply a new point ell.
  void getpoint(DESreal4vector& ellOUT);      // Get the value of the current point ell.
  int getmethod();                            // Method used to get this point.
  
private:
  int seed;                      // Seed for random numbers.
  int nverts;                    // Number of vertices in loop.
  double alpha[maxsize];         // Probabilities for different newell methods.
  double beta[maxsize];          // Cumulative probabilities for different newell methods.
  DESreal4vector p[maxsize];     // The external momenta.
  DESreal4vector Q[maxsize];     // The fixed momenta in the denominators for each propagator.
  int indexA;                    // Initial state vertices are between (nverts,1) and (indexA,indexA + 1).
  int indexAp1;                  // Equals indexA + 1 cyclically defined.
  double MTsq;                   // Scale for closeness to double parton scattering singularity.
  DESreal4vector Pplain,Pbar;    // vectors P and Pbar.
  double tiny;                   // A very small number.
  double small;                  // A pretty small number.
  double huge;                   // A very big number.
  double dotPPbar;               // P * Pbar.
  double Awide;                  // Parameter for distribution of ell, method 1.
  double msmall4_a;              // Parameter for distribution of ell, method 2.
  double mlarge4_a;              // Parameter for distribution of ell, method 2.
  double msmall4_b;              // Parameter for distribution of ell, method 9.
  double mlarge4_b;              // Parameter for distribution of ell, method 9.
  double Asoft_a;                // Parameter for distribution of ell, method 3.
  double Bsoft_a;                // Parameter for distribution of ell, method 3.
  double Asoft_b;                // Parameter for distribution of ell, method 10.
  double Bsoft_b;                // Parameter for distribution of ell, method 10.
  double Acoll_a;                // Parameter for distribution of ell, method 4.
  double Bcoll_a;                // Parameter for distribution of ell, method 4.
  double Acoll_b;                // Parameter for distribution of ell, method 11.
  double Bcoll_b;                // Parameter for distribution of ell, method 11.
  double AIcoll_a;               // Parameter for distribution of ell, method 5.
  double BIcoll_a;               // Parameter for distribution of ell, method 5.
  double AIcoll_b;               // Parameter for distribution of ell, method 12.
  double BIcoll_b;               // Parameter for distribution of ell, method 12.
  double Atlike_a;               // Parameter for distribution of ell, method 6.
  double Atlike_b;               // Parameter for distribution of ell, method 13.
  double Aslike_a;               // Parameter for distribution of ell, method 7.
  double Aslike_b;               // Parameter for distribution of ell, method 14.
  double Allike_a;               // Parameter for distribution of ell, method 8.
  double Bllike_a;               // Parameter for distribution of ell, method 8.
  double Msqllike_a;             // Parameter for distribution of ell, method 8.
  double Allike_b;               // Parameter for distribution of ell, method 15.
  double Bllike_b;               // Parameter for distribution of ell, method 15.
  double Msqllike_b;             // Parameter for distribution of ell, method 15.
  double tildeallike3_a;         // Parameter for distribution of ell, method 16.
  double Allike3_a;              // Parameter for distribution of ell, method 16.
  double Bllike3_a;              // Parameter for distribution of ell, method 16.
  double Msqllike3_a;            // Parameter for distribution of ell, method 16.
  double tildeallike3_b;         // Parameter for distribution of ell, method 17.
  double Allike3_b;              // Parameter for distribution of ell, method 17.
  double Bllike3_b;              // Parameter for distribution of ell, method 17.
  double Msqllike3_b;            // Parameter for distribution of ell, method 17.
  DESrandom r;                   // Random number generator, r.random() returns a random number.
  DESindex indx;                 // Index function, indx(i) returns i mod nverts.
  DESreal4vector ell;            // The point ell.
  int method;                    // The method used to select the point.
  DESreal4vector ell_ijkA(int i, int j, int k); // Function for choosing ell, method 16;
  double rho_ijkA(int i, int j, int k, DESreal4vector ell); // Density for method 16.
  DESreal4vector ell_ijkB(int i, int j, int k); // Function for choosing ell, method 17;
  double rho_ijkB(int i, int j, int k, DESreal4vector ell); // Density for method 17.
};

class DESdeform
{
public:
  DESdeform(const int& nverts, const int& indexA, const DESreal4vector Q[maxsize]);    
  // Constructor.
  void deform(DESreal4vector& ellIN, DEScmplx4vector& lcOUT, complex<double>& jacobian);
  // Returns deformed point lc = ell + i kappa, deformation jacobian dl/dell.
  void deformtest(DESreal4vector& ellIN, DEScmplx4vector& lcOUT, complex<double>& jacobian,
                  int jIN);
  // Test version. As long as indexA is not 1 or nverts -1, this returns deformed point
  // only the contribution from j = jIN to lc = ell + i kappa, deformation jacobian dl/dell.
  // For the \tilde c_+ contribution, use jIN = -1 and for the \tilde c_- contribution,
  // use jIN = -2.
  
private:
  DESreal4vector zero4vector;    // The zero 4-vector;
  double small;                  // A pretty small number.
  int nverts;                    // Number of vertices in loop.
  DESindex indx;                 // indx(i) returns i mod nverts.
  int indexA;                    // Initial state vertices are between propagators nverts and 1       
  int indexAp1;                  // and between propagators indexA and indexAp1.
  DESreal4vector Q[maxsize];     // Momentum vectors for propagators.
  DESreal4vector P, Pbar;        // The incoming momenta
  DESreal4vector Pc, Pbarc;      // The incoming momenta, covariant components.
  double PdotPbar;               // Their dot product.
  double M1;                     // Parameter, default = 0.05 Sqrt(PdotPbar).
  double M2;                     // Parameter, default = Sqrt(PdotPbar).
  double M3;                     // Parameter, default = Sqrt(PdotPbar).
  double gamma1;                 // Parameter, default = 0.7.
  double gamma2;                 // Parameter, default = 1.
  DESreal4vector ell;            // The real part of lc.
  DESreal4vector kappa;          // The imaginary part of lc.
  DEScmplx4vector lc;            // ell + i kappa.
  bool ellinCplus(int j);   // ell is inside the positive light cone starting from Q[j].
  bool ellinCminus(int j);   // ell is inside the negative light cone starting from Q[j].
  double hplus(const DESreal4vector k);  // Function that smoothly vanishes if k in C+.
  double hminus(const DESreal4vector k); // Function that smoothly vanishes if k in C-.
  DESreal4vector gradloghplus(const DESreal4vector k); //  \partial_\mu log(hplus(k)).
  DESreal4vector gradloghminus(const DESreal4vector k); // \partial_\mu log(hminus(k)).
  double g(const DESreal4vector k);  // Function that smoothly vanishes for large k.
  double gplus(const DESreal4vector k);  // Function emphasizing backward light cones.
  double gminus(const DESreal4vector k);  // Function emphasizing forward light cones.
  DESreal4vector gradlogg(const DESreal4vector k); //  \partial_\mu log(g(k)).
  DESreal4vector gradloggplus(const DESreal4vector k); //  \partial_\mu log(gplus(k)).
  DESreal4vector gradloggminus(const DESreal4vector k); //  \partial_\mu log(gminus(k)).
  void c_j(int& j, double& gval, DESreal4vector& gradloggval, 
           double& cj, DESreal4vector& gradlogcj); // Compute c_j and its gradient.
  void c_pm(int& indexpm, double& x, double& xbar, 
            double& cj, DESreal4vector& gradlogcj); // Compute c_pm and its gradient.
};


int factorial(const int& n);             // Function n!.
int cyclic(const int& i, const int& n);     // Function to give cyclic index in the range 1,...,n.
complex<double> Determinant(const complex<double> b[maxsize][maxsize], const int& dimension);
                                            // Determinant of matrix b of size dimension*dimension


class FermionLoop
{
public:
  FermionLoop(int nverts); // Constructor.
               // Postcondition:
               //    Class object is constructed & polarization vectors are
               //    initialized to be zero.
  void FeedEpsilon(const int nv, const DEScmplx4vector& epsilon);
               // Precondition:
               //    nv is the index for the target vertex.
               //    epsilon is a complex 4-vector.
               // Postcondition:
               //    The polarization vector at the given vertex will be set
               //    equal to the given "epsilon".
  void CalEpsilonNull(const int nv, const int h, const DESreal4vector& p);
  void CalEpsilonCoulomb(const int nv, const int h, const DESreal4vector& p);
               // There are two versions, for null-plane gauge or Coulomb gauge.
               // Precondition:
               //    nv is the index for the target vertex.
               //    h could only be either "+1" or "-1", which separately stands
               //    for plus helicity and minus helicity.
               //    l is a 4-vector, q is an array of 4-vectors, where the
               //    external momemtum p[i] = q[i-1] - q[i], and loop momemtum
               //    are l-q[i].
               // Postcondition:
               //    Calculate the polarization vector at the given vertex
               //    with the chosen helicity.
  complex<double> NumeratorV(const DEScmplx4vector& l, const DEScmplx4vector Q[]);
               // Precondition:
               //    Input the polarization vectors and loop momemtum.
               // Postcondition:
               //    Calculate the numerator of the loop for vector couplings.
  complex<double> NumeratorM(const DEScmplx4vector& l, const DEScmplx4vector Q[], const double m);
               // Vector current with nonzero mass.
  complex<double> NumeratorUV(const DEScmplx4vector& l);
               // This one if for the "extra" part of the UV subtraction.
  complex<double> NumeratorL(const DEScmplx4vector& l, const DEScmplx4vector Q[]);
               // Precondition:
               //    Input the polarization vectors and loop momemtum.
               // Postcondition:
               //    Calculate the numerator of the loop for left-handed couplings.
  complex<double> NumeratorR(const DEScmplx4vector& l, const DEScmplx4vector Q[]);
               // Precondition:
               //    Input the polarization vectors and loop momemtum.
               // Postcondition:
               //    Calculate the numerator of the loop for right-handed couplings.
  complex<double> AltNumeratorV(const DEScmplx4vector& l, const DEScmplx4vector Q[]);
               // This version expands the trace in dot products of 4 vectors.
private:
  int nverts;
  DEScmplx4vector epsln[maxsize];
};
//==================================================================================================
// Function to generate a list of permutations of the integers 1,2,...,n.
// The ith permutation is pi[j] = PermList[i][j].
void getperms(int nverts, int PermList[maxperms][maxsize], int indexAlist[maxperms]);
// This calls
void perm(int n, int pi[], int& count, bool up, int PermList[maxperms][maxsize]);
//==================================================================================================
// Function n!.
int factorial(const int& n); 
//==================================================================================================
// Function to give cyclic index in the range 1,...,n. E.g. cyclic(7,5) = 2.
int cyclic(const int& i, const int& n);     // Function to give cyclic index in the range 1,...,n.
//==================================================================================================
complex<double> Determinant(const complex<double> b[maxsize][maxsize], const int& dimension);
// Determinant of matrix b of size dimension*dimension
//==================================================================================================
// Function to produce the four-photon scattering amplitude.
double fourphoton(int helicity[maxsize],DESreal4vector p[maxsize]);
//==================================================================================================