// Main program for numerical integration of one loop Feynman diagrams.
// This version is for 2 photons -> (nverts - 2) photons with a fermion loop.
// D. E. Soper
// See date below.

using namespace std;
#include <iostream>
#include <cmath>
#include <complex>
#include <string>
#include "Nphoton.h"

//----------

//----------
int main()
{
char* date = "25 September 2006";
int nverts;
enum VertexTypes {VECTOR, LEFTHANDED, RIGHTHANDED};
VertexTypes vertextype;
enum PolarizationTypes {NULLPLANE, COULOMB};
PolarizationTypes polarizationtype;
double xcutoff = 1.0e-6;
double lambdadeform = 1.3;
double m0factor = 0.5;
const double zero = 0.0;
const double pivalue = 3.14159265358979;
const double em4over3 = 0.2635971381157268;  // exp(-4/3)
const double em3over2 = 0.2231301601484298;  // exp(-3/2)
const complex<double>ii(0.0,1.0);
const complex<double>sqrtii(0.7071067811865475,0.7071067811865475);
const complex<double>sqrtminusii(0.7071067811865475,-0.7071067811865475);
bool altway = false;
//
int i,j,k,n,mu,graphnumber,ip1;
int numberofgraphs,onlygraphnumber,graphnumbermin,graphnumbermax;
int seed;
int nreno,npoint;
char response[100];
DESxarray x;
DESzarray z;
bool xOK;
double m0sq;
DEScmplx4vector ell;
DEScmplx4vector lvec,lvec1,lvec2;
double ellsqE;
double rts,rtssq;
DESreal4vector p[maxsize];  // outgoing momenta, p[1],...,p[nverts] where we throw away p[0]
DEScmplx4vector Quv[maxsize];  // Equals zero. For UV subtraction if nverts = 4.
DESreal4vector Q[maxsize];
DES_linfo linfo;
double commonfactor;
double rhox,rhoell,prefactor;
complex<double> f,det_dzdx;
complex<double> sum,value,result;
double sumsqR,sumsqI,errorR,errorI;
double sumRplus, sumsqRplus, resultRplus, errorRplus;
double sumRminus,sumsqRminus,resultRminus,errorRminus;
double sumIplus, sumsqIplus, resultIplus, errorIplus;
double sumIminus,sumsqIminus,resultIminus,errorIminus;
double etasum;
double temp,fell;
complex<double> tempC;
complex<double> lambdasqval,numeratorval;
complex<double> numeratorUV1,numeratorUV2A,numeratorUV2B;
complex<double> lambdasqUV1,lambdasqUV2,fUV1,fUV2;
double musq;
complex<double> valplain,valUV;
complex<double> numeratorUV;
complex<double> lamsqplain,lamsqUV;
complex<double> zdetplain;
double etasumplain;
complex<double> Lamsq;
int pi[maxsize];
int PermList[maxperms][maxsize];
int helicity[maxsize];
double npoints;
int const maxmethods = 11;
double alpha[2*maxsize],beta[2*maxsize]; // Probabilities for given distributions \rho_i(x).
// Here alpha[i] are the input probabilities, beta[i] are after redistribution.
int method;
int Ntries[2*maxsize],Nbad[2*maxsize];
double fluct[2*maxsize];
double fluctnorm;
double fluctG[10000];
int triesG[10000];
int totaltriesG;
int factnvertsm1,factnvertsm2;
complex<double> testsumell,testsumuv,testsumsoft[maxsize],testsumcoll[maxsize];
double testsumsqell,testsumsquv,testsumsqsoft[maxsize],testsumsqcoll[maxsize];
complex<double> testvalue;
complex<double> testresultell,testresultuv,testresultsoft[maxsize],testresultcoll[maxsize];
double testerrorell,testerroruv,testerrorsoft[maxsize],testerrorcoll[maxsize];
double theta;     // Angle by which final state is to be rotated.
DESreal4vector v; // Just a generic real 4 vector.
double mass = 0.0;     // The fermion mass. User input can change this.
bool zeromass = true;  // True if the mass is zero.  User input can change this.
bool photon3hto0 = false; // True if polarization of photon 3 should be proportional to p[3].
bool doublemusq = false; // True doubles renormalization scale musq.
//
cout << "Please give number of photons (4, 5, 6, 7, or 8)."<< endl;
cin >> nverts;
factnvertsm1 = factorial(nverts - 1);
factnvertsm2 = factorial(nverts - 2);
cout << "There are "<< factnvertsm1 << " graphs."<< endl;
cout << "Please enter 0 if you want all of them. "
     << "Otherwise enter the number of the graph you want."<< endl;
cin >> onlygraphnumber;
if (onlygraphnumber == 0)
  {
  numberofgraphs = factnvertsm1;
  for(i = 1; i <= factnvertsm1; i++) triesG[i] = 1;  // This gives equal weights to all graphs.
  totaltriesG = 0;
  for(i = 1; i <= factnvertsm1; i++) totaltriesG = totaltriesG + triesG[i];
  }
else
  {
  numberofgraphs = 1;
  for(i = 1; i <= factnvertsm1; i++) triesG[i] = 0;
  triesG[onlygraphnumber] = 1;
  totaltriesG = 1;
  }
cout << "What kind of vertices would you like "<< 
         "(vector or left or right)?" << endl;
cin >> response;
if(strcmp(response,"vector") == 0)
  {
  vertextype = VECTOR;
  cout << "Using vector vertices." << endl;
  }
else if (strcmp(response,"left") == 0)
  {
  vertextype = LEFTHANDED;
  cout << "Using left-handed vertices." << endl;
  }
else if (strcmp(response,"right") == 0)
  {
  vertextype = RIGHTHANDED;
  cout << "Using right-handed vertices." << endl;
  }
else
  {
  cout << "I did not recognize the vertex type. " << response << endl;
  exit(1);
  }
//
cout << "Please polarization type, null or coulomb."<< endl;
cin >> response;
if(strcmp(response,"null") == 0)
  {
  polarizationtype = NULLPLANE;
  cout << "Using null-plane gauge polarizations." << endl;
  }
else if (strcmp(response,"coulomb") == 0)
  {
  polarizationtype = COULOMB;
  cout << "Using Coulomb gauge polarizations." << endl;
  }
else
  {
  cout << "I did not recognize the polarization type. " << response << endl;
  exit(1);
  }
//
if (vertextype == VECTOR)
  {
  cout << "With vector vertices, you can have a fermion mass." << endl;
  cout << "If you want no mass, please enter 0. Otherwise enter the mass you want." << endl;
  cin >> mass;
  if(mass == 0) zeromass = true;
  else zeromass = false;
  }
cout << "For a default calculation, please enter 0. Otherwise enter"<< endl;
cout << " 1 to cut deformation in half"<< endl;
cout << " 2 to double m0"<< endl;
cout << " 3 to double give photon 3 a polarization proportional to its momentum"<< endl;
cout << " 4 to use the alternative formulation"<< endl;
cout << " 5 to double the renormalization scale (for N = 4)"<< endl;
cin >> i;
if(i == 1)
  {
  cout << "Cutting deformation in half." << endl;
  lambdadeform = lambdadeform/2.0;
  }
else if(i == 2)
  {
  cout << "Doubling m0." << endl;
  m0factor = m0factor*2.0;
  }
else if (i==3)
  {
  cout << "Giving photon 3 a polarization proportional to its momentum." << endl;
  photon3hto0 = true;
  }
else if (i==4)
  {
  cout << "Using alternative formulation." << endl;
  altway = true;
  }
else if (i==5)
  {
  cout << "Doubling the renormalization scale (for N = 4)." << endl;
  doublemusq = true;
  }
else cout << "Using default values." << endl;
//
cout << "Please give angle theta for rotation of final state about y-axis."<< endl;
cin >> theta;
cout << "Please give seed."<< endl;
cin >> seed;
cout << "Please give number of points." << endl;
cin >> nreno;
//
cout << endl << endl;
cout << "                       Program Nphoton" << endl;
cout << "                  Z. Nagy and D. E. Soper" << endl;
cout << "                        Version 1.0" << endl;
cout << "                      " << date << endl;
cout << endl;
cout << "This program calculates the scattering amplitude"<<endl;
cout << "for 2 photons to (N-2) photons with an electron loop." << endl;
cout << "Reports the integral equal to the scattering amplitude M" << endl;
cout << "times Sqrt(s)^((N-4)/2) divided by alpha_em^(N/2)." << endl;
cout << "The number of vertices is " << nverts << "." << endl;
if(vertextype == VECTOR) cout << "Using vector vertices." << endl;
else if (vertextype == LEFTHANDED) cout << "Using left-handed vertices." << endl;
else if (vertextype == RIGHTHANDED) cout << "Using right-handed vertices." << endl;
else cout << "I did not recognize the vertex type. " << vertextype << endl;
if(polarizationtype == NULLPLANE) cout << "Using null-plane gauge polarizations." << endl;
else if (polarizationtype == COULOMB) cout << "Using Coulomb gauge polarizations." << endl;
else cout << "I did not recognize the polarizatin type. " << vertextype << endl;
cout << "The seed is " << seed << "." << endl;
if(onlygraphnumber == 0)
  {
  cout << "Using all graphs." << endl;
  }
else
  {
  cout << "Using only graph " << onlygraphnumber << "." <<endl;
  }
if(zeromass) cout << "Fermion mass is zero." << endl;
else cout << "Fermion mass = " << mass << "." << endl;
if(altway) cout << "Using alternative formulation."<<endl;
cout << "Deformation parameter is " << lambdadeform << ". (Default is 1.3.)" << endl;
cout << "Parameter m0 is " << m0factor << " * Sqrt(s). (Default is 0.5 * Sqrt(s).)" << endl;
if(nverts == 4) 
  {
  if (doublemusq) cout << "Renormalization scale is mu^2 = 2 * m0^2." << endl;
  else cout << "Renormalization scale is mu^2 = m0^2." << endl;
  }
//
// ------- Construct a set of q vectors around loop, then matrix S. -------
//
// We simply use an arbitrary choice of the external outgoing momenta.
// The user could specify any desired (lightlike) momenta.
if (nverts == 8)
  {
  p[3].in(1) =  33.5;
  p[3].in(2) =   0.9;
  p[3].in(3) =  25.0;
  p[4].in(1) =   1.5;
  p[4].in(2) =  24.3;
  p[4].in(3) =   0.3;
  p[5].in(1) = -19.1;
  p[5].in(2) = -35.1;
  p[5].in(3) = - 3.3;
  p[6].in(1) =  28.2;
  p[6].in(2) = - 6.6;
  p[6].in(3) =   8.2;
  p[7].in(1) = -12.2;
  p[7].in(2) = - 8.6;
  p[7].in(3) =   8.2;
  }
else if (nverts == 7)
  {
  p[3].in(1) =  33.5;
  p[3].in(2) =  15.9;
  p[3].in(3) =  25.0;
  p[4].in(1) = -12.5;
  p[4].in(2) =  15.3;
  p[4].in(3) =   0.3;
  p[5].in(1) = -10.0;
  p[5].in(2) = -18.0;
  p[5].in(3) = -3.3;
  p[6].in(1) =  18.2;
  p[6].in(2) = - 6.6;
  p[6].in(3) =   8.2;
  }
else if (nverts == 6)
  {
  p[3].in(1) =  33.5;
  p[3].in(2) =  15.9;
  p[3].in(3) =  25.0;
  p[4].in(1) = -12.5;
  p[4].in(2) =  15.3;
  p[4].in(3) =   0.3;
  p[5].in(1) = -10.0;
  p[5].in(2) = -18.0;
  p[5].in(3) = -3.3;
  }
else if (nverts == 5)
  {
  p[3].in(1) =  33.5;
  p[3].in(2) =  15.9;
  p[3].in(3) =  25.0;
  p[4].in(1) = -12.5;
  p[4].in(2) =  15.3;
  p[4].in(3) =   0.3;
  }
else if (nverts == 4)
  {
  p[3].in(1) =  50.0;
  p[3].in(2) =   0.0;
  p[3].in(3) =   0.0;
  }
else
  {
  cout << "Unfortunately, I am not programmed for nverts = " << nverts << "." << endl;
  exit(1);
  }
// Alternative choice that is close to dps singularity. 
// Comment this out if you want the nicer point.
/*if (nverts == 6)
  {
  cout << "Using external momenta close to a dps singularity." << endl;
  p[3].in(1) =  33.5;
  p[3].in(2) =  15.9;
  p[3].in(3) =  25.0;
  p[4].in(1) =  12.5;
  p[4].in(2) =  15.3;
  p[4].in(3) =   0.3;
  p[5].in(1) = -10.0;
  p[5].in(2) = -18.0;
  p[5].in(3) = -3.3;
  }
*/
// Now set \vec p_i for i = nverts by momentum conservation.
p[nverts].in(1) = 0.0;
p[nverts].in(2) = 0.0;
p[nverts].in(3) = 0.0;
for (i = 3; i < nverts; i++)
  {
  p[nverts].in(1) = p[nverts].in(1) - p[i](1);
  p[nverts].in(2) = p[nverts].in(2) - p[i](2);
  p[nverts].in(3) = p[nverts].in(3) - p[i](3);
  }
// Now calculate rts and put final state particles on shell.
rts = 0;
for(i = 3; i <= nverts; i++)
  {
  p[i].in(0) = sqrt(pow(p[i](1),2) + pow(p[i](2),2) + pow(p[i](3),2));
  rts = rts + p[i].in(0);
  }
rtssq = pow(rts,2);
// Finally, set (outgoing) momenta for the incoming particles.
p[1].in(1) =  0.0;
p[1].in(2) =  0.0;
p[1].in(3) = -0.5*rts;
p[1].in(0) = -0.5*rts;
p[2].in(1) =  0.0;
p[2].in(2) =  0.0;
p[2].in(3) =  0.5*rts;
p[2].in(0) = -0.5*rts;
//
// -------------
// Rotate.
for(i = 3; i <= nverts; i++)
  {
  v = p[i];
  p[i].in(1) =    cos(theta)*v(1) + sin(theta)*v(3);
  p[i].in(3) =  - sin(theta)*v(1) + cos(theta)*v(3);
  }
// -------------
// Write out momenta.
cout << "Outgoing momenta, with s = "<< rtssq  << " and theta = "<< theta << "." << endl;
for(i = 1; i <= nverts; i++)
  {
  cout << i<<": "<<p[i](0)<<"  "<<p[i](1)<<"  "<<p[i](2)<<"  "<<p[i](3)<<endl;
  }
  cout << endl;
//----------
if(nverts == 4)
  {
  cout << "theta = "<< theta <<
          "; t/s = "<< square(p[1] + p[3])/square(p[1] + p[2]) <<
          "; u/s = "<< square(p[1] + p[4])/square(p[1] + p[2]) << endl<<endl;
  }
//----------
// Check for double parton scattering singularities.
if(nverts > 5)
  {
  cout << "Check for double parton scattering singularities." << endl;
  for(i = 3; i <= nverts; i++)
    {
    for(j = i+1; j <= nverts; j++)
      {
      temp = pow(p[i](1) + p[j](1),2) + pow(p[i](2) + p[j](2),2) ;
      cout << "(P_T["<< i <<"] + P_T["<< j <<"])^2/s = "<<  temp/rtssq << endl;
      }
    }
  if(nverts > 6)
    {
    for(i = 3; i <= nverts; i++)
      {
      for(j = i+1; j <= nverts; j++)
        {
        for(k = j+1; k <= nverts; k++)
          {
          temp = pow(p[i](1) + p[j](1) + p[k](1),2) + pow(p[i](2) + p[j](2) + p[k](2),2) ;
          cout << "(P_T["<< i <<"] + P_T["<< j <<"] + P_T["<< k <<"])^2/s = "<<  temp/rtssq << endl;
          }
        }
      }
    }
  cout << endl;
  }
//----------
// We specify the helicities we want. This can be changed.
if (nverts == 8)
  {
  helicity[1] =  1;
  helicity[2] = -1;
  helicity[3] = -1;
  helicity[4] =  1;
  helicity[5] =  1;
  helicity[6] = -1;
  helicity[7] =  1;
  helicity[8] = -1;
  }
else if (nverts == 7)
  {
  helicity[1] =  1;
  helicity[2] = -1;
  helicity[3] =  1;
  helicity[4] = -1;
  helicity[5] = -1;
  helicity[6] =  1;
  helicity[7] = -1;
  }
else if (nverts == 6)
  {
  helicity[1] =  1;
  helicity[2] = -1;
  helicity[3] = -1;
  helicity[4] =  1;
  helicity[5] =  1;
  helicity[6] = -1;
  }
else if (nverts == 5)
  {
  helicity[1] =  1;
  helicity[2] = -1;
  helicity[3] = -1;
  helicity[4] =  1;
  helicity[5] =  1;
  }
else if (nverts == 4)
  {
  helicity[1] =  1;
  helicity[2] =  1;
  helicity[3] =  1;
  helicity[4] =  1;
  }
if( photon3hto0 ) helicity[3] = 0;
cout << "External photon helicities: ";
for (i = 1; i <= nverts ; i++) cout << helicity[i] << " ";
cout << endl << "The physical helicities of the incoming photons are the negative of these."<<endl;
//----------
cout << "We will integrate with x[i] > " << xcutoff <<"."<<endl;
//----------
// Probabilities to choose points x with various distributions.
alpha[0]  = 0.00;  // Not used.
alpha[1]  = 0.02;  // Probability for uniform distribution.
alpha[2]  = 0.00;  // ADJUSTED BELOW. Probability for one of "soft-n" distributions.
alpha[3]  = 0.05;  // Probability for one of the modified collinear distributions.
alpha[4]  = 0.02;  // Probability for one of the modified soft distributions.
alpha[5]  = 0.02;  // Probability for one of the "x0" collinear distritutions.
alpha[6]  = 0.04;  // Probability for one of the "x^0 x^n x^np1" collinear distributions.
alpha[7]  = 0.02;  // Probability for one of the collinear-soft distributions.
alpha[8]  = 0.30;  // Probability for the double parton scattering distribution.
alpha[9]  = 0.10;  // Probability for one of the collinear distribtions.
alpha[10] = 0.02;  // Probability for the small x^0, uniform x^i distribution.
alpha[11] = 0.02;  // Probability for the large x^0, uniform x^i distribution.
if(nverts < 6) alpha[8] = 0.1; // Don't need this.
if(nverts < 5)
  { 
  alpha[1]  = 0.30; // More points with uniform distribution for 4 vertices.
  alpha[5]  = 0.00; // Collinear x0 method doesn't work for 4 vertices.
  alpha[6]  = 0.05; // Less of this.
  alpha[8]  = 0.0;  // Don't need this.
  alpha[9]  = 0.0;  // Don't need this.
  alpha[11] = 0.20; // More of this.
  }
if(nverts >= 7)
  {
  alpha[0]  = 0.00;  // Not used.
  alpha[1]  = 0.02;  // Probability for uniform distribution.
  alpha[2]  = 0.00;  // ADJUSTED BELOW. Probability for one of "soft-n" distributions.
  alpha[3]  = 0.10;  // Probability for one of the modified collinear distributions.
  alpha[4]  = 0.00;  // Probability for one of the modified soft distributions.
  alpha[5]  = 0.00;  // Probability for one of the "x0" collinear distritutions.
  alpha[6]  = 0.02;  // Probability for one of the "x^0 x^n x^np1" collinear distributions.
  alpha[7]  = 0.00;  // Probability for one of the collinear-soft distributions.
  alpha[8]  = 0.20;  // Probability for the double parton scattering distribution.
  alpha[9]  = 0.40;  // Probability for one of the collinear distribtions.
  alpha[10] = 0.02;  // Probability for the small x^0, uniform x^i distribution.
  alpha[11] = 0.02;  // Probability for the large x^0, uniform x^i distribution.
  }
alpha[2] = 1.0 - alpha[1];
for(n = 3; n <= maxmethods; n++)
  {
  alpha[2] = alpha[2] - alpha[n];
  }
if(alpha[2] < 0.0)
  {
  cout << "The probabilities entered result in negative alpha[2]." << endl;
  exit(1);
  }
//----------
FermionLoop fermion(nverts);             // Create class for numerator calculation
DES_Smatrix Splain;                      // Create class for matrix S.
DES_Smatrix Smassless;                   // Create class for matrix S fermion mass = 0.
DESxdeform xdeform(nverts, rtssq, lambdadeform);  // Create class for contour deformation.
DESellpoint ellpoint(seed, nverts);               // Create class for numerator loop momentum.
DESxpoint xpoint(seed, nverts, xcutoff, alpha);   // Create class for x point.
//-------
// Set initial S variables and one linfo variable.
Splain.setnverts(nverts);
Smassless.setnverts(nverts);
linfo.nverts = nverts;
m0sq = pow(m0factor*rts,2);
musq = m0sq;  // This is the renormalization scale for nverts = 4.
if(doublemusq) musq = 2.0*musq;
Splain.setm0sq(m0sq);
Smassless.setm0sq(m0sq);
Splain.setrtssq(rtssq);
Smassless.setrtssq(rtssq);
Splain.settype(PLAIN);
Smassless.settype(PLAIN);
Splain.setmass(mass);
Smassless.setmass(zero);
for(i = 1; i <= nverts; i++)
  {
  for(mu = 0; mu <= 3; mu++)
    {
    Q[i].in(mu) = 0.0;         // Modified for PLAIN contribution for each graph. (Q is real.)
    Quv[i].in(mu) = 0.0;       // For UV counterterm in case nverts = 4. (Quv is complex.)
    }
  }
getperms(nverts - 1,PermList); // Get permutations of vertices, always with pi(nverts) = nverts.
//
// Set up integral.
sum = 0.0;
sumsqR = 0.0;
sumsqI = 0.0;
sumRplus = 0.0;
sumRminus = 0.0;
sumIplus = 0.0;
sumIminus = 0.0;
sumsqRplus = 0.0;
sumsqRminus = 0.0;
sumsqIplus = 0.0;
sumsqIminus = 0.0;
testsumell = 0.0;
testsumsqell = 0.0;
testsumuv = 0.0;
testsumsquv = 0.0;
for(i = 1; i <= nverts; i++)
  {
  testsumsoft[i] = 0.0;
  testsumsqsoft[i] = 0.0;
  testsumcoll[i] = 0.0;
  testsumsqcoll[i] = 0.0;
  }
int ngood = 0;
int nbad = 0;
double worst = 300000;
if(nverts == 6) worst = 1.0e10;
if(nverts == 7) worst = 1.0e12;
for(i = 1; i <= maxmethods; i++)
  {
  fluct[i] = 0.0;
  Ntries[i] = 0;
  Nbad[i] = 0;
  }
//
// Loop over choice of graphs.
if(numberofgraphs == factnvertsm1)
  {
  graphnumbermin = 1;
  graphnumbermax = factnvertsm1;
  }
else
  {
  graphnumbermin = onlygraphnumber;
  graphnumbermax = onlygraphnumber;
  }
// Common normalization factor, - (m_0^2/rts^2)  N! (4\pi)^{N/2} /(2\pi)^4.
commonfactor = - m0sq/rtssq*factorial(nverts)*pow(4.0*pivalue,0.5*nverts)/pow(2.0*pivalue,4);
//
for (graphnumber = graphnumbermin; graphnumber <= graphnumbermax; graphnumber ++)
{
// Pick permutation of external legs.
for(i = 1; i <= nverts; i++) pi[i] = PermList[graphnumber][i];
//
// Note that we have previously defined for S:
// m0sq, rtssq, nverts;
// We need to define Q and linfo.Q.
for( mu = 0; mu <= 3; mu++) Q[1].in(mu) = 0.0;  // This is a convenient choice.
for (i = 1; i <= nverts - 1; i++)
  {
  Q[i+1] = Q[i] + p[pi[i]];
  }
for(i = 1; i <= nverts; i++)
  {
  for( mu = 0; mu <= 3; mu++) linfo.Q[i].in(mu) = Q[i](mu);
  }
// -----------------
// We can now make the matrix S and inform other classes:
//
Splain.make(Q);           // Use Q to create new matrix S.
Smassless.make(Q);        // This one with mass = 0.
xpoint.newS(Smassless);   // Initialize xpoint with new matrix S (with mass = 0).
xdeform.newS(Splain);     // Also send new information to xdeform.
//
// Construct the external photon polarization vectors.
for (i = 1; i <= nverts; i++)
  {
  if(polarizationtype == COULOMB) fermion.CalEpsilonCoulomb(i, helicity[pi[i]], p[pi[i]]);
  else if(polarizationtype == NULLPLANE) fermion.CalEpsilonNull(i, helicity[pi[i]], p[pi[i]]);
  else
    {
    cout << "Oops, wrong polarization type." << endl;
    exit(1);
    }
  }
// -------------------------------------------------------------
// The integration routine, integrates f \propto 1/[Lambda^2]^{nverts - 2}
// times 1/(1 + ell(0)^2 + ell(1)^2 + ell(2)^2 + ell(3)^2)^{nverts + 1}
//
for (npoint = 1; npoint <= nreno*triesG[graphnumber]; npoint++)
{
  xpoint.neu(x,xOK,method);            // New point x.
  Ntries[method] = Ntries[method] + 1;
  if (xOK)
    {
    rhox = xpoint.rho();                // rho_x(x).
    ellpoint.neu(ell);                  // New point ell.
    rhoell = ellpoint.rho();            // rho_l(ell).
    ellsqE = 0.0;                       // Euclidean square of ell.
    for( mu = 0; mu <= 3; mu++) ellsqE = ellsqE + pow(real(ell(mu)),2);
    // Factor common to main term and UV subtraction. Uses
    // commonfactor = 
    //  - m0sq/rtssq*factorial(nverts)*pow(4.0*pivalue,0.5*nverts)/pow(2.0*pivalue,4);
    prefactor = commonfactor*double(totaltriesG)/double(triesG[graphnumber]);
    if(altway) // This is an alternative formulation, used as a check.
      {
      prefactor = -prefactor*pow(m0sq/rtssq,2)/rhox/rhoell;  // Extra - m_0^4 with no fell.
      }
    else // This is the default formulation.
      {
      fell = pow(1.0 + ellsqE, - nverts - 1);  // fell = 1/[1+l*P*l]^{N+1}
      prefactor = prefactor*fell/rhox/rhoell;
      }
    xdeform.deform(x,z,det_dzdx,etasum);       // Calculates z and det_dzdx.
    lambdasqval = Lambdasq(nverts,Splain,z);
    if(altway) // This is an alternative formulation, used as a check.
      {
      f = pow(1.0 - z(0),nverts - 1)*pow(z(0),2);
      f = f*pow(rtssq/(lambdasqval + ii*z(0)*(1.0 - z(0))*m0sq*ellsqE),nverts + 1);
      }
    else // This is the default formulation.
      {
      f =  pow(1.0 - z(0),nverts - 3)*pow(rtssq/lambdasqval, nverts - 1);
      }
    // Collect information for finding transformed loop momentum. 
    // linfo.nverts and linfo.Q remain fixed.
    linfo.Lambdasq = lambdasqval;
    linfo.z = z;
    linfo.etasum = etasum;
    linfo.type = PLAIN;
    linfo.ell = ell;
    if(altway) // This is an alternative formulation, used as a check.
      {
      // Transformed loop momentum for numerator, lvec1 with ell and lvec2 with - ell.
      tempC = sqrt(z(0)*(1.0 - z(0))*m0sq);
      lvec1.in(0) =   tempC*sqrtii*ell(0);
      lvec2.in(0) = - tempC*sqrtii*ell(0);
      for(mu = 1; mu <=3; mu++)
        {
        lvec1.in(mu) =   tempC*sqrtminusii*ell(mu);
        lvec2.in(mu) = - tempC*sqrtminusii*ell(mu);
        }
      for(i = 1; i <= nverts; i++)
        {
        lvec1 = lvec1 + z(i)*linfo.Q[i];
        lvec2 = lvec2 + z(i)*linfo.Q[i];
        }
      lvec1 = (1.0/(1.0 - z(0)))*lvec1;
      lvec2 = (1.0/(1.0 - z(0)))*lvec2;
      }
    else // This is the default formulation.
      {
      lvec1 = lmod(linfo); // Transformed loop momentum for numerator function.
      linfo.ell = - ell;
      lvec2 = lmod(linfo); // Transformed loop momentum for - ell.
      }
    if(vertextype == VECTOR)
      {
      if (zeromass)
        {
        numeratorval = (  fermion.NumeratorV(lvec1, linfo.Q) 
                        + fermion.NumeratorV(lvec2, linfo.Q))
                   * 0.5/pow(rts,nverts);  // The numerator function.
        }
      else
        {
        numeratorval = (  fermion.NumeratorM(lvec1, linfo.Q, mass) 
                        + fermion.NumeratorM(lvec2, linfo.Q, mass))
                   * 0.5/pow(rts,nverts);  // The numerator function.
        }      
      }
    else if (vertextype == LEFTHANDED)
      {
      numeratorval = (fermion.NumeratorL(lvec1, linfo.Q) + fermion.NumeratorL(lvec2, linfo.Q))
                   * 0.5/pow(rts,nverts);  // The numerator function.  
      }
    else if (vertextype == RIGHTHANDED)
      {
      numeratorval = (fermion.NumeratorR(lvec1, linfo.Q) + fermion.NumeratorR(lvec2, linfo.Q))
                   * 0.5/pow(rts,nverts);  // The numerator function.  
      }
    else
      {
      cout << "Problem with the vertex type." << endl;
      exit(1);
      }
    value = prefactor*det_dzdx*f*numeratorval;
    valplain = value;         // For diagnostic purposes.
    lamsqplain = lambdasqval; // For diagnostic purposes.
    zdetplain = det_dzdx;     // For diagnostic purposes.
    etasumplain = etasum;     // For diagnostic purposes.
    // Test integrals.
    testvalue = nverts*(nverts - 1)/pivalue/pivalue*fell/rhoell;
    testsumell = testsumell + testvalue;
    testsumsqell = testsumsqell + norm(testvalue);
    testvalue = double(factnvertsm1)/pow(1.0 - z(0),nverts-1);
    testvalue = testvalue*det_dzdx/rhox;
    testsumuv = testsumuv + testvalue;
    testsumsquv = testsumsquv + norm(testvalue);
    for(i = 1; i <= nverts; i++)
      {
      ip1 = i + 1;
      if(ip1 > nverts) ip1  = ip1 - nverts;
      testvalue = double(factnvertsm1)/pow(1.0 - z(i),nverts-1);
      testvalue = testvalue*det_dzdx/rhox;
      testsumsoft[i] = testsumsoft[i] + testvalue;
      testsumsqsoft[i] = testsumsqsoft[i] + norm(testvalue);
      testvalue = double(factnvertsm2)/(z(i) + z(ip1))/pow(1.0 - z(i) - z(ip1),nverts - 2);
      testvalue = testvalue*det_dzdx/rhox;
      testsumcoll[i] = testsumcoll[i] + testvalue;
      testsumsqcoll[i] = testsumsqcoll[i] + norm(testvalue);
      }
    // Now UV subtraction, if needed.
    if (nverts == 4)
      {
      // musq is the renormalization scale.
      lambdasqUV1 = - pow(1.0 - z(0),2)*musq*em4over3 + ii*z(0)*(1.0 - z(0))*m0sq;
      lambdasqUV2 = - pow(1.0 - z(0),2)*musq*em3over2 + ii*z(0)*(1.0 - z(0))*m0sq;
      if(altway) // This is an alternative formulation, used as a check.
        {
        f = pow(1.0 - z(0),nverts - 1)*pow(z(0),2);
        fUV1 = f*pow(rtssq/(lambdasqUV1 + ii*z(0)*(1.0 - z(0))*m0sq*ellsqE),nverts + 1);
        fUV2 = f*pow(rtssq/(lambdasqUV2 + ii*z(0)*(1.0 - z(0))*m0sq*ellsqE),nverts + 1);
        }
      else // This is the default formulation.
        {
        fUV1 =  pow(1.0 - z(0),nverts - 3)*pow(rtssq/lambdasqUV1, nverts - 1);
        fUV2 =  pow(1.0 - z(0),nverts - 3)*pow(rtssq/lambdasqUV2, nverts - 1);
        }
      if(altway) // This is an alternative formulation, used as a check.
        {
        tempC = sqrt(z(0)*(1.0 - z(0))*m0sq);
        lvec.in(0) =   tempC*sqrtii*ell(0);
        for(mu = 1; mu <=3; mu++) lvec.in(mu) = tempC*sqrtminusii*ell(mu);
        lvec = (1.0/(1.0 - z(0)))*lvec;
        numeratorUV1 = fermion.NumeratorV(lvec, Quv)/pow(rts,nverts);
        numeratorUV2A = fermion.NumeratorUV(lvec)/pow(rts,nverts);
        numeratorUV2B = numeratorUV2A;
        }
      else // This is the default formulation.
        {
        linfo.ell = ell;
        linfo.type = UV;    // Temporarily switch type.
        linfo.Lambdasq = lambdasqUV1;
        lvec = lmod(linfo);
        numeratorUV1 = fermion.NumeratorV(lvec, Quv)/pow(rts,nverts); 
        numeratorUV2A = fermion.NumeratorUV(lvec)/pow(rts,nverts);
        linfo.Lambdasq = lambdasqUV2;
        lvec = lmod(linfo);
        numeratorUV2B = fermion.NumeratorUV(lvec)/pow(rts,nverts); 
        linfo.type = PLAIN; // Restore type.    
        }
      // The numerator function.
      // Note that for the UV subtraction with nverts = 4, the numerator is an
      // even function of ell, so we do not need to average over ell and - ell.
      if(vertextype != VECTOR)
        {
        cout << "We can do nverts = 4 only for vector vertices." << endl;
        exit(1);
        }
      valUV = prefactor*det_dzdx*(fUV1*(numeratorUV1 - numeratorUV2A)
                                 + fUV2*numeratorUV2B);
      value = value - valUV;
      lamsqUV = lambdasqval;      // For diagnostic purposes.
      numeratorUV = numeratorval; // For diagnostic purposes.
      }
    // Done with contributions to value.
    sum = sum + value;
    sumsqR = sumsqR + pow(real(value),2);
    sumsqI = sumsqI + pow(imag(value),2);
    ngood = ngood + 1;
    fluct[method] = fluct[method] + norm(value);
    fluctG[graphnumber] = fluctG[graphnumber] + norm(value);
    temp = real(value);
    if (temp > 0.0)
      {
      sumRplus = sumRplus + temp;
      sumsqRplus = sumsqRplus + pow(temp,2);
      }
    else
      {
      sumRminus = sumRminus - temp;
      sumsqRminus = sumsqRminus +  pow(temp,2);
      }
    temp = imag(value);
    if (temp > 0.0)
      {
      sumIplus = sumIplus + temp;
      sumsqIplus = sumsqIplus + pow(temp,2);
      }
    else
      {
      sumIminus = sumIminus - temp;
      sumsqIminus = sumsqIminus +  pow(temp,2);
      }
    // Diagnostic report for selected points.
    if (abs(value) > worst || (npoint < 2 && graphnumber < 2) )
      {
      if (abs(value) > worst)
        {
        worst = abs(value);
        cout << endl << "worst abs(value) so far is "<< worst << 
                " for graph " << graphnumber << " point " << npoint << endl;
        }
      else 
        {
        cout << endl << "for point " << npoint << " for graph " << graphnumber << endl;
        }
      for (i = 0; i <= nverts; i++)
        {
        cout << " x("<<i<<") = " << x(i) << endl;
        }
      cout << endl;
      for (i = 0; i <= nverts; i++)
        {
        cout << " z("<<i<<") = " << z(i) << endl;
        }
      cout << endl;
      cout << " rhox = " << rhox << "  fell = " << fell  << " rhoell = " << rhoell << endl;
      cout << " prefactor = " << prefactor << endl << endl;
      cout << "valplain = "<< valplain << endl;
      cout << "Lambdasqplain = "<< lamsqplain << endl;
      if(nverts == 4)
        {
        cout << "fermion.NumeratorV(lvec, Quv)/pow(rts,nverts) is " 
           << fermion.NumeratorV(lvec, Quv)/pow(rts,nverts) << "." << endl;
        cout << "fermion.AltNumeratorV(lvec, Quv)/pow(rts,nverts) is " 
           << fermion.NumeratorV(lvec, Quv)/pow(rts,nverts) << "." << endl;
        }
      if(vertextype == VECTOR)
        {
        cout << "N(l)/rts^nverts  = " << fermion.NumeratorV(lvec1, linfo.Q)/pow(rts,nverts) << endl;
        cout << "N(-l)/rts^nverts = " << fermion.NumeratorV(lvec2, linfo.Q)/pow(rts,nverts) << endl;
        if(nverts < 6)
          {
          cout << "Compare alternative N(l)/rts^nverts: " 
               << fermion.AltNumeratorV(lvec1, linfo.Q)/pow(rts,nverts) << endl;
          }
        }
      else if (vertextype == LEFTHANDED)
        {
        cout << "N(l)/rts^nverts  = " << fermion.NumeratorL(lvec1, linfo.Q)/pow(rts,nverts) << endl;
        cout << "N(-l)/rts^nverts = " << fermion.NumeratorL(lvec2, linfo.Q)/pow(rts,nverts) << endl;
        }
      else if (vertextype == RIGHTHANDED)
        {
        cout << "N(l)/rts^nverts  = " << fermion.NumeratorR(lvec1, linfo.Q)/pow(rts,nverts) << endl;
        cout << "N(-l)/rts^nverts = " << fermion.NumeratorR(lvec2, linfo.Q)/pow(rts,nverts) << endl;
        }
      else
        {
        cout << "Problem with the vertex type." << endl;
        exit(1);
        }
//
      cout << "  Numerator values for vertices of different types" << endl;
      cout << "  N_V(l)/rts^nverts  = " << fermion.NumeratorV(lvec1, linfo.Q)/pow(rts,nverts) << endl;
      cout << "  N_L(l)/rts^nverts  = " << fermion.NumeratorL(lvec1, linfo.Q)/pow(rts,nverts) << endl;
      cout << "  N_R(l)/rts^nverts  = " << fermion.NumeratorR(lvec1, linfo.Q)/pow(rts,nverts) << endl;
      cout << "  N_V(-l)/rts^nverts = " << fermion.NumeratorV(lvec2, linfo.Q)/pow(rts,nverts) << endl;
      cout << "  N_L(-l)/rts^nverts = " << fermion.NumeratorL(lvec2, linfo.Q)/pow(rts,nverts) << endl;
      cout << "  N_R(-l)/rts^nverts = " << fermion.NumeratorR(lvec2, linfo.Q)/pow(rts,nverts) << endl;
//
      cout << "10*(1-z^0)^nverts*numerator/rts^{nverts-2}/Lambdasq = " 
              << 10.0*pow(1.0 - z(0),nverts)*numeratorval*rts*rts/lamsqplain << endl;
      cout << "zdetplain = "<< zdetplain << endl;
      cout << "Sum of the eta^i values, sumplain = "<< etasumplain << endl;
      cout << "Cancellation: ratio of sum of absolute values of Lambdasq(x) pieces" << endl;
      cout << "to abs(Lambdasq(z)) " << Lambdasqabs(nverts, Splain, x)/abs(lamsqplain) << endl ;
      if (nverts == 4)
        {
        cout << "UV subtraction:" << endl;
        cout << "valUV = "<< valUV << endl;
        cout << "LambdasqUV = "<< lamsqUV << endl;
        cout << "Nuv(l)/rts^nverts  = " << numeratorUV << endl;
        }
      cout << endl;
  //=========================
  //=========================
      cout << "--------Reports Completed---------" << endl;
      } // Close if (abs(value) > worst || (npoint < 2 && graphnumber < 20) ).
        // Done with diagnostic report for selected points.
    }  
  else
    {
    nbad = nbad + 1;
    Nbad[method] = Nbad[method] + 1;
    } // Close if (xOK) ... else ...
} // Close for (npoint = 1; npoint <= nreno; npoint++)
} // Close for (graphnumber = graphnumbermin; graphnumber <= graphnumbermax; graphnumber ++)
//
cout << endl << "fraction of bad points = "<< double(nbad)/double(ngood + nbad) << endl << endl;
npoints = nreno*totaltriesG;
result = sum/npoints;
errorR = sqrt( (sumsqR/npoints - pow(real(result),2))/(npoints - 1) );
errorI = sqrt( (sumsqI/npoints - pow(imag(result),2))/(npoints - 1) );
resultRplus = sumRplus/npoints;
resultRminus = sumRminus/npoints;
resultIplus = sumIplus/npoints;
resultIminus = sumIminus/npoints;
errorRplus =  sqrt( (sumsqRplus/npoints  - pow(resultRplus,2) )  /(npoints - 1));
errorRminus = sqrt( (sumsqRminus/npoints - pow(resultRminus,2))  /(npoints - 1));
errorIplus =  sqrt( (sumsqIplus/npoints  - pow(resultIplus,2) )  /(npoints - 1));
errorIminus = sqrt( (sumsqIminus/npoints - pow(resultIminus,2))  /(npoints - 1));
// Also the test integrals.
testresultell = testsumell/double(ngood);
testerrorell = sqrt( (testsumsqell/double(ngood) - norm(testresultell)) /double(ngood-1));
testresultuv = testsumuv/double(npoints);
testerroruv = sqrt( (testsumsquv/double(npoints) - norm(testresultuv)) /double(npoints - 1));
for(i = 1; i <= nverts; i++)
  {
  testresultsoft[i] = testsumsoft[i]/double(npoints);
  testerrorsoft[i] = 
    sqrt( (testsumsqsoft[i]/double(npoints) - norm(testresultsoft[i])) /double(npoints - 1));
  testresultcoll[i] = testsumcoll[i]/double(npoints);
  testerrorcoll[i] =
    sqrt( (testsumsqcoll[i]/double(npoints) - norm(testresultcoll[i])) /double(npoints - 1));
  }
// The fluctuations.
for(i = 1; i <= maxmethods; i++) fluct[i] = fluct[i]/npoints;
// Normalized fluctuations by point selection method.
for(i = 1; i <= maxmethods; i++) beta[i] = double(Ntries[i])/npoints;
fluctnorm = 0.0;
j = 0;
for(i = 1; i <= maxmethods; i++) 
  {
  if (beta[i] > 0.0) 
    {
    fluctnorm = fluctnorm + fluct[i]/beta[i];
    j = j + 1;
    }
  }
fluctnorm = fluctnorm/j;
// Fluctuations by graph number.
temp = 0.0;
for(i = 1; i <= factnvertsm1; i++) 
  {
  fluctG[i] = fluctG[i]/nreno/double(triesG[i]);
  temp = temp + fluctG[i];
  }
temp = temp/double(factnvertsm1);
for(i = 1; i <= factnvertsm1; i++) 
  {
  fluctG[i] = fluctG[i]/temp;
  }
// Output.
if(onlygraphnumber == 0)
  {
  cout << "For " << factnvertsm1 << " graphs, after "<< nreno  
     << " * "<< totaltriesG << " points, we have a result." << endl;
  }
else
  {
  cout << "For the single graph " << onlygraphnumber << ", after "<< nreno  
     << " * "<< totaltriesG << " points, we have a result." << endl;
  }
cout << "integral = " << result << " +/- ("<< errorR << "," << errorI << ")" <<endl;
cout << "|integral| = " << abs(result) << " +/- " << 
        (abs(real(result))*errorR + abs(imag(result))*errorI)/abs(result) << endl<< endl;
if(nverts == 4)
  {
  cout << "Here  is what Binoth, Glover, Marquard and van der Bij have for the amplitude." << endl;
  cout << "External photon helicities: ";
  for (i = 1; i <= nverts ; i++) cout << helicity[i] << " ";
  cout << endl << "The physical helicities of the incoming photons are the negative of these."<<endl;
  cout << "Amplitude = " << fourphoton(helicity, p) << endl << endl;
  }
cout << "positive real contributions:      " << resultRplus  << " +/- " << errorRplus  << endl;
cout << "negative real contributions:      " << resultRminus << " +/- " << errorRminus << endl;
cout << "positive imaginary contributions: " << resultIplus  << " +/- " << errorIplus  << endl;
cout << "negative imaginary contributions: " << resultIminus << " +/- " << errorIminus << endl <<endl;
cout << "Sampling methods." << endl;
cout << "methods in: "<< endl;
for(i = 1; i <= maxmethods; i++)
  {
  if (alpha[i] > 0.0) cout << i << ": "<<alpha[i] << " " ;
  }
cout << endl << "methods actual: "<< endl;
for(i = 1; i <= maxmethods; i++)
  {
  if (beta[i] > 0.0) cout << i << ": "<< beta[i] << " " ;
  }
cout << endl << "fraction of bad points: "<< endl;
for(i = 1; i <= maxmethods; i++)
  {
  if (beta[i] > 0.0) cout << i << ": "<< double(Nbad[i])/double(Ntries[i]) << " " ;
  }
cout <<endl<< "Relative fluctuations per point: "<< endl;
for(i = 1; i <= maxmethods; i++)
  {
  if (beta[i] > 0.0) cout << i << ": "<< fluct[i]/beta[i]/fluctnorm << " " ;
  }
cout << endl;
// Report test integrals.
cout << "Here are the test integrals."<< endl;
cout << "Here npoints = " << npoints << " while ngood = " << ngood << endl;
cout << "Integral over ell: " << testresultell << " +/- " << testerrorell << endl;
cout << "Integral over z, UV region: " << testresultuv << " +/- " << testerroruv << endl;
for(i = 1; i <= nverts; i++)
  {
  cout << "Integral over z, soft region "<< i <<": " 
       << testresultsoft[i] << " +/- " << testerrorsoft[i] << endl;
  }
for(i = 1; i <= nverts; i++)
  {
  cout << "Integral over z, coll region " << i << ": " 
       << testresultcoll[i] << " +/- " << testerrorcoll[i] << endl;
  }
// Report of fluctuations by graph number.
/*
cout << endl << "Fluctuations by graph number:" << endl;
j = 0;
for(i = 1; i <=factnvertsm1; i++)
  {
  j = j + 1;
  cout << i << "  " << fluctG[i] << "   ";
  if(j >= 5)
    {
    j = 0;
    cout << endl;
    }
  }
cout << endl;
*/
//
}  //end of main()
// ---------------------------------------------------