// Specification file for classes etc. used for numerical integration of one loop Feynman diagrams.
// D.E. Soper
// 25 September 2006

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

class DESrandom
{
public:
  DESrandom(/*in*/int seed);      // Constructor. Sets seed.   
  DESrandom();                    // Default constructor. (puts seed as 0)
  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 (and only) constructor
  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); 

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

enum SubtrTypes {PLAIN,UV};

struct DES_Sinfo
{
int nverts;                 // Number of vertices in loop.
double m0sq;                // Square of mass parameter m0.
double rtssq;               // Square of c.m. energy, rts.
DESreal4vector Q[maxsize];  // Momenta around the loop, Q[1],...,Q[nverts]. Throw away Q[0].
SubtrTypes type;            // Type of subtraction term S is used for.
int n;                      // Index of subtraction term S is used for.
};

struct DES_dpsinfo
{
int n;                     // Number of choices for A & B, = 0 or 1.
int A;                     // Index A: lines A and A+1 may be collinear.
int Ap1;                   // A + 1, cyclically defined.
int B;                     // Index B: lines B and B+1 may be collinear.
int Bp1;                   // B + 1, cyclically defined.
double fA;                 // Momentum fraction for lines A and A+1.
double fB;                 // Momentum fraction for lines B and B+1.
double oneminusfA;         // 1 - fA.
double oneminusfB;         // 1 - fB.
double df;                 // Estimated difference ~ |fA - f_B|.
};

struct DES_collinfo
{
int N;                     // Index n; lines n and n+1 may be collinear.
int Np1;                   // N + 1, cyclically defined.
int A;                     // Index A: \xi_A may be not so tiny.
int B;                     // Index B: \xi_B may be not so tiny.
double f;                  // Momentum fraction for lines n and n+1.
double df;                 // Estimated difference ~ |fA - f_B|.
};

class DES_Smatrix
{
public:
  DES_Smatrix();                           // Default constructor.
  void make(DESreal4vector Q[maxsize]);    // Makes matrix S of type needed.
  double& in(const int& i, const int& j);	 // s.in(i,j) for input, s.in(i,j) = ... .
  const double& operator()(const int& i, const int& j) const;
                                           // s(i,j) for output, ... = s(i,j).
  void setnverts(int nverts);              // Set number of vertices in loop.
  const int& getnverts();                  // Report number of vertices in loop.
  void setrtssq(double rtssq);             // Set square of c.m. energy.
  const double& getrtssq();                // Report square of c.m. energy.
  void setm0sq(double m0sq);               // Set squared mass parameter for UV convergence.
  const double& getm0sq();                 // Report squared mass parameter for UV convergence.
  void setmass(double mass);               // Set fermion mass.
  const double& getmass();                 // Report fermion mass.
  void settype(SubtrTypes type);           // Set type of term S is used for (PLAIN or UV).
  const SubtrTypes& gettype();             // Report type of term S is used for (PLAIN or UV).
private:
  double Smatrix[maxsize][maxsize];        // The matrix S.
  double rtssq;                            // Square of c.m. energy.
  double m0sq;                             // Squared mass parameter for UV convergence.
  double mass;                             // Fermion mass.
  double twomsq;                           // Two times square of fermion mass.
  SubtrTypes type;                         // Type of subtraction term S is used for.
  int nverts;                              // Number of vertices in loop.
};

class DESxarray
{
public:
  DESxarray();                                  // Default constructor
  double& in(const int& i);                     // x.in(i) for input, x.in(i) = ... .
  const double& operator()(const int& i) const; // x(i) for output, ... = x(i).
private:
  double x[maxsize];                            // The array x.
};

class DESzarray
{
public:
  DESzarray();                                           // Default constructor
  complex<double>& in(const int& i);                     // z.in(i) for input, z.in(i) = ... .
  const complex<double>& operator()(const int& i) const; // z(i) for output, ... = z(i).
private:
  complex<double> z[maxsize];                            // The array z.
};

class DESxpoint
{
public:
  DESxpoint(int seedIN, int nvertsIN, DES_Smatrix sIN, double xcutoffIN, double alphaIN[]); 
  DESxpoint(int seedIN, int nvertsIN, double xcutoffIN, double alphaIN[]); 
                                               // Constructors. (There is no default constructor.)
  void newS(DES_Smatrix sIN);                  // Change matrix S.
  void neu(DESxarray& xOUT, bool& xOK, int& method);
          // New point x; xOK is true if point is OK; method is the index of the method used.
  double rho();                                // Density rho(x).
  double altrho(const DESxarray& xIN);  // Reset x and get rho(x) from rho().

private:
  int seed;                      // Seed for random numbers.
  int nverts;                    // Number of vertices in loop.
  int nvertsp1;                  // nverts + 1.
  double xcutoff;                // We demand x[i] > xcutoff.
  double alpha[2*maxsize];       // Probabilities for given distributions.
  double Ncollx0;                // Number of cases for the "x0" collinear distritutions.
  double capA[maxsize];          // Limit on g_n(x). Must be <= 1.
  double cparam;                 // Parameter to control power: rho_n ~ 1/g_n(x)^{N-1-cparam}
  double c10;                    // Parameter to control power in method 10.
  double c11;                    // Parameter to control power in method 11.
  double s[maxsize][maxsize];    // Matrix: \Lambda = (-1/2) S_ij x^i x^j.
  double s0;                     // Parameter with dimensions momentum^2, cancels from physical results.
  double m0sq;                   // Squared mass parameter used to control UV convergence.
  double xsmax;                  // Maximum xs in "x^0 x^n x^np1" collinear distributions.
  double x[maxsize];             // The point x.
  bool OK;                       // Point is OK. All x_i > 0.
  DESrandom r;                   // Random number generator, specified above.
  double a[maxsize][maxsize];    // Coefficients a_{ni} in function g_n(x).
  double b[maxsize];             // Coefficients b_{n} in function g_n(x).
  double g(const int& n);        // Scaled approximation to \Lambda in soft region n, g_n(x).
  double lambda[maxsize];        // Equals s0 / min {S_{n-1,n+1}, S_ni}.
  double partialprob[maxsize];   // This is sum_{j=1}^n w^(s)_j/Sum_i w^(s)_i in notes.
  double wsoftsum;               // This is sum_n w^(s)_n in notes.
  double xAminusmax[maxsize];    // Limit on xA for cs-minus method of choosing points.
  double xBminusmax[maxsize];    // Limit on xB for cs-minus method of choosing points.
  double xnminusmax[maxsize];    // Limit on xn for cs-minus method of choosing points.
  double xAplusmax[maxsize];	   // Limit on xA for cs-plus method of choosing points.
  double xBplusmax[maxsize];	   // Limit on xB for cs-plus method of choosing points.
  double xnp1plusmax[maxsize];   // Limit on xnp1 for cs-plus method of choosing points.
  DES_dpsinfo dps;               // Information for sampling for double parton scattering singularity.
  DES_collinfo coll[maxsize];    // Information for sampling for collinear singularity.
  int collN;                     // Number of choices of n,A,B for collinear sampling.
};

class DESxdeform
{
public:
  DESxdeform(int nvertsIN, double rtssq, double lambdadeformIN);    // Constructor (S = 0)
  void newS(DES_Smatrix sIN);                                       // Change matrix S.
  void deform(/*in*/ DESxarray& x, 
              /*out*/ DESzarray& z, complex<double>& jacobian, double& etasum);
            // Returns deformed point z = x + i y, deformation jacobian dz/dx, & sum of eta^i.

private:
  int nverts;                      // Number of vertices in loop.
  double rtssq;                    // Square of c.m. energy rts.
  double lambdadeform;             // Parameter to control amount of deformation.
  double s[maxsize][maxsize];      // Matrix: \Lambda = (-1/2) S_ij x^i x^j.
};

class DESellpoint
{
public:
  DESellpoint(int seedIN, int nvertsIN);  // Constructor. (There is no default constructor.)   
  void neu(DEScmplx4vector& ellOUT);      // New point ell, returned as complex.
  double rho();                           // Density rho(ell).
  
private:
  int seed;                      // Seed for random numbers.
  int nverts;                    // Number of vertices in loop.
  double Aradial;                // Parameter controlling distribution of |ell|.
  DESrandom r;                   // Random number generator, specified above.
  DESreal4vector ell;            // The point ell.
  double ellEsq;                 // Euclidian ell^2 = ell(0)^2 + ell(1)^2 + ell(2)^2 + ell(3)^2.
};

void DESsort(double x[], const int& nmax);  // Sorts an array x of length nmax and returns sorted array
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
complex<double> Lambdasq(const int nverts, DES_Smatrix s, DESzarray z);
                                            // Function Lambda^2(z).
double Lambdasqabs(const int nverts, DES_Smatrix s, DESxarray x);
                                            // Function Lambda^2_abs(x), using |s(i,j)|.

struct DES_linfo
{
int nverts;                  // Number of vertices in loop.
complex<double> Lambdasq;    // The denominator factor Lambda^2(z).
DEScmplx4vector Q[maxsize];  // Momenta around the loop, q[1],...,q[nverts]. Throw away q[0].
DESzarray z;                 // Complex Feynman parameters z[0],z[1],...,z[nverts].
double etasum;               // Sum of eta[0], eta[1],..., eta[nverts].
DEScmplx4vector ell;         // Loop momentum before shifts
SubtrTypes type;             // Type of subtraction term S is used for.
};

DEScmplx4vector lmod(const DES_linfo& linfo);
                      // Transformed loop momentum for numerator function.

complex<double> numerator(const DEScmplx4vector& lvec, const DEScmplx4vector Q[]);
                     // Numerator function N(lvec,Q).

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 n, int PermList[maxperms][maxsize]);
// This calls
void perm(int n, int pi[], int& count, bool up, int PermList[maxperms][maxsize]);
//==================================================================================================
// Function to produce the four-photon scattering amplitude.
double fourphoton(int helicity[maxsize],DESreal4vector p[maxsize]);
//==================================================================================================