// 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 "NphotonD.h"

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

//----------
int main()
{
char* date = "11 December 2008";
const double pivalue = 3.14159265358979;
//
int i,j,k,mu;
double temp;
//
// ------------------------- Set number of vertices. -------------------------------
// int nverts = 6;
int nverts = 8;
//* cout << "Please give number of photons (4, 5, 6, 7, or 8)."<< endl;
//* cin >> nverts;
// ---------------------------- Set graph number. ----------------------------------
int factnvertsm1 = factorial(nverts - 1);
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;
int onlygraphnumber;
cin >> onlygraphnumber;
int numberofgraphs, graphnumbermin, graphnumbermax;
if (onlygraphnumber == 0)
  {
  numberofgraphs = factnvertsm1;
  graphnumbermin = 1;
  graphnumbermax = factnvertsm1;
  }
else
  {
  numberofgraphs = 1;
  graphnumbermin = onlygraphnumber;
  graphnumbermax = onlygraphnumber;
  }
// -------------------------------- Set theta. -------------------------------------
  double theta = 0.0;
  cout << "Please give rotation angle. " 
       << "(2.32 is close to double parton scattering singuarity.)" << endl;
  cin >> theta;
// -------------------------------- Set seed. --------------------------------------
int seed= 1;
cout << "Please give seed."<< endl;
cin >> seed;
// ---------------------- Set number of integration points. ------------------------
cout << "Please give number of points." << endl;
int nreno;
cin >> nreno;
// --------------------------------- Set mass. -------------------------------------
bool zeromass = true;
double mass = 0.0;
// ----------------------------- Set vertex type. ----------------------------------
enum VertexTypes {VECTOR, LEFTHANDED, RIGHTHANDED};
VertexTypes vertextype;
vertextype = VECTOR;
cout << "Using vector vertices." << endl;
// --------------------------- Set polarization type. ------------------------------
enum PolarizationTypes {NULLPLANE, COULOMB};
PolarizationTypes polarizationtype;
polarizationtype = COULOMB;
cout << "Using Coulomb gauge polarizations." << endl;
// -------------------------- Set method probabilities. ----------------------------
// Probabilities for methods for choosing new loop momentum points.
double alpha[maxsize];
for(i = 0; i < maxsize; i++) alpha[i] = 0.0;
alpha[1]  = 1.0; // Method 1: roughly uniform over a wide range. Number changed below.
alpha[2]  = 0.10; // Method 2: distributed near double parton scattering singularity.
alpha[3]  = 0.05; // Method 3: distributed near the soft singularities. 
alpha[4]  = 0.03; // Method 4: distributed near the collinear singularities.
alpha[5]  = 0.01; // Method 5: distributed near the initial state collinear singularities.
alpha[6]  = 0.005; // Method 6: distributed near timelike separated cones.
alpha[7]  = 0.02; // Method 7: distributed near spacelike separated cones.
alpha[8]  = 0.03; // Method 8: distributed near lightlike separated cones.
alpha[9]  = 0.10; // Method 9: distributed near double parton scattering singularity.
alpha[10] = 0.06; // Method 10: distributed near the soft singularities. 
alpha[11] = 0.005; // Method 11: distributed near the collinear singularities.
alpha[12] = 0.005; // Method 12: distributed near the initial state collinear singularities.
alpha[13] = 0.02; // Method 13: distributed near timelike separated cones.
alpha[14] = 0.005; // Method 14: distributed near spacelike separated cones.
alpha[15] = 0.06; // Method 15: distributed near lightlike separated cones.
alpha[16] = 0.13; // Method 16: distributed near lightlike separated cones plus third cone.
alpha[17] = 0.23; // Method 17: distributed near lightlike separated cones plus third cone.
alpha[18] = 0.05; // Method 18: distributed near lightlike separted cones plus next cone.
alpha[19] = 0.07; // Method 19: distributed near lightlike separted cones plus next cone.
alpha[20] = 0.0;  // Method 20: like 16, but gets more intersections when nverts = 8.
alpha[21] = 0.0;  // Method 21: like 17, but gets more intersections when nverts = 8.
//
if( nverts == 8)
{
  alpha[1]  = 1.0; // Method 1: roughly uniform over a wide range. Number changed below.
  alpha[2]  = 0.03; // Method 2: distributed near double parton scattering singularity.
  alpha[3]  = 0.02; // Method 3: distributed near the soft singularities. 
  alpha[4]  = 0.03; // Method 4: distributed near the collinear singularities.
  alpha[5]  = 0.10; // Method 5: distributed near the initial state collinear singularities.
  alpha[6]  = 0.01; // Method 6: distributed near timelike separated cones.
  alpha[7]  = 0.005; // Method 7: distributed near spacelike separated cones.
  alpha[8]  = 0.05; // Method 8: distributed near lightlike separated cones.
  alpha[9]  = 0.03; // Method 9: distributed near double parton scattering singularity.
  alpha[10] = 0.04; // Method 10: distributed near the soft singularities. 
  alpha[11] = 0.005; // Method 11: distributed near the collinear singularities.
  alpha[12] = 0.005; // Method 12: distributed near the initial state collinear singularities.
  alpha[13] = 0.03; // Method 13: distributed near timelike separated cones.
  alpha[14] = 0.005; // Method 14: distributed near spacelike separated cones.
  alpha[15] = 0.03; // Method 15: distributed near lightlike separated cones.
  alpha[16] = 0.05; // Method 16: distributed near lightlike separated cones plus third cone.
  alpha[17] = 0.06; // Method 17: distributed near lightlike separated cones plus third cone.
  alpha[18] = 0.08; // Method 18: distributed near lightlike separted cones plus next cone.
  alpha[19] = 0.07; // Method 19: distributed near lightlike separted cones plus next cone.
  alpha[20] = 0.16;  // Method 20: like 16, but gets more intersections when nverts = 8.
  alpha[21] = 0.18;  // Method 21: like 17, but gets more intersections when nverts = 8.
}
int maxindex = 21;
if( maxindex > maxsize - 1)
{
  cout << "We need to enlarge space for alpha[i]." << endl;
  exit(1);
}
for(i = 2; i <= maxindex; i++) alpha[1] = alpha[1] - alpha[i];
// ---------------------------- Set external momenta. ------------------------------
DESreal4vector p[maxsize];  // outgoing momenta p[1],...,p[nverts] ,throw away p[0].
double rts,rtssq;
// We simply use an arbitrary choice of the external outgoing momenta.
if (nverts == 8)
{
  p[3].in(1) =  33.5;
  p[3].in(2) =   5.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) = -20.0;
  p[5].in(3) =  -3.3;
  p[6].in(1) =  18.2;
  p[6].in(2) =  11.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
{
  cout << "We need nverts > 4." << endl;
  exit(1);
}
// 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);
//
double virtualitycut = 1.0e-8*rtssq;  // This cut could be changed.
// Rotate the final state momenta.
DESreal4vector v;
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);
}
// 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;
//
// ---------------------- Finished setting external momenta. -----------------------
// ---------------------------- Set photon helicities. -----------------------------
int helicity[maxsize];
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;
  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 
{
  cout << "We need nverts > 4." << endl;
  exit(1);
}
// ---------------------- Finished setting photon helicities. ----------------------
//
cout << endl << endl;
cout << "                     Program NphotonD" << endl;
cout << "              N Photon Scattering, Direct Method" << endl;
cout << "                W. Gong, Z. Nagy, D. E. Soper" << endl;
cout << "                        Version 1.0" << endl;
cout << "                      " << date << endl;
giveversion(); // Prints version and date of subroutines.
cout << endl;
cout << "The number of vertices is " << nverts << "." << endl;
if(onlygraphnumber == 0)
  {
  cout << "Using all graphs." << endl;
  }
else
  {
  cout << "Using only graph " << onlygraphnumber << "." <<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 polarization type. " << vertextype << endl;
cout << "The seed is " << seed << "." << endl;  
cout << "We will cut |(l - Q_i)^2| < " << virtualitycut/rtssq << " times s." << endl;
// -------------
// 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;
//----------
// Check for double parton scattering singularities.
double MTsq;    // Scale for closeness to double parton scattering singularity.
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;
//----------
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 << "Probabilities for methods for choosing new loop momenta." << endl;
//for(i = 1; i <= maxindex; i++) cout << i << ": " << alpha[i] << endl;
//cout << endl;
//----------
// Get permutations of vertices, always with pi(N) = 1.
// Also pi(indexA) = 2.
int indexAlist[maxperms];
int PermList[maxperms][maxsize];
getperms(nverts, PermList, indexAlist);
// -------------------------------------------------------------
// Instance of class loopmomenta set using graph number 1.
int graphnumber = 1;
int pi[maxsize];
for(i = 1; i <= nverts; i++) pi[i] = PermList[graphnumber][i];
int indexA;
indexA = indexAlist[graphnumber];  // pi[indexA] = 2.
DESreal4vector pnew[maxsize];
for(i = 1; i <= nverts; i++) pnew[i] = p[pi[i]];
// Create class member for numerator loop momentum.
DESloopmomenta loopmomenta(seed, nverts, pnew, indexA, alpha);
// -------------------------------------------------------------
// The integration routine.
cout << endl << "We will integrate with " << nreno << " points per graph. " << endl << endl;
// Set up integral.
complex<double> sum, value;
double sumsqR, sumsqI;
//
sum = 0.0;
sumsqR = 0.0;
sumsqI = 0.0;
//
complex<double> sumtest, valuetest, ftest, lcsqE;
double sumsqRtest, sumsqItest;
double Mtestsq = 1000.0; // This could be changed.
sumtest = 0.0;
sumsqRtest = 0.0;
sumsqItest = 0.0;
//
double worst = 1.0e14;
// ----- Set up how many times to evaluate each graph ------------
int pointsG[10000];
int totalpoints = 0;
double weightG[10000];
double totalweightG = 0.0;
if(factorial(nverts - 1) >= 10000)
{
  cout << "We need to enlarge array pointsG[...]." << endl;
  exit(1);
}
for(graphnumber = graphnumbermin; graphnumber <= graphnumbermax; graphnumber++)
{
  // Pick permutation of external legs.
  for(i = 1; i <= nverts; i++) pi[i] = PermList[graphnumber][i];
  indexA = indexAlist[graphnumber];  // pi[indexA] = 2.
  // Permuted external momenta.
  for(i = 1; i <= nverts; i++) pnew[i] = p[pi[i]];
  // Provide new data for this graph to our instance of
  // class DESloopmomenta.
  loopmomenta.givenewgraph(pnew, indexA);
  double PtsqScale = loopmomenta.getMTsq();
  if(indexA == 1 | indexA == nverts-1) PtsqScale = rtssq;
  double smallestdot = abs(dot(pnew[nverts],pnew[1]));
  double temp;
  for (i = 1; i < nverts; i++)
  { 
    temp = abs(dot(pnew[i],pnew[i+1]));
    if(temp < smallestdot) smallestdot = temp;
  }
  if(smallestdot < PtsqScale) PtsqScale = smallestdot;
  weightG[graphnumber] = pow(PtsqScale,-0.3); // This weight factor could be changed.
  totalweightG += weightG[graphnumber];
}
for(graphnumber = graphnumbermin; graphnumber <= graphnumbermax; graphnumber++)
{
  double temp = nreno*weightG[graphnumber]/totalweightG*(graphnumbermax - graphnumbermin + 1); 
  pointsG[graphnumber] = int(temp+1);
  // cout << graphnumber << " pointsG = " << pointsG[graphnumber] << endl;
  totalpoints += pointsG[graphnumber];
}
cout << "total points = " << totalpoints << endl << endl;
  
// ----- Variables for recording fluctions according to point selection method. ----
int Ntries[maxsize], Nbad[maxsize];
double fluct[maxsize], beta[maxsize];
int method = 0;
for(i = 1; i < maxsize; i++)
  {
  fluct[i] = 0.0;
  Ntries[i] = 0;
  Nbad[i] = 0;
  }
//
// Looop over graphs.
for(graphnumber = graphnumbermin; graphnumber <= graphnumbermax; graphnumber++)
{
// Prepare calculation for this graph number.
  complex<double> sumGraph = 0.0;
  double sumsqRGraph = 0.0, sumsqIGraph = 0.0;
// Pick permutation of external legs.
for(i = 1; i <= nverts; i++) pi[i] = PermList[graphnumber][i];
indexA = indexAlist[graphnumber];  // pi[indexA] = 2.
// Permuted external momenta.
for(i = 1; i <= nverts; i++) pnew[i] = p[pi[i]];
// Class member for numerator calculation
FermionLoop fermion(nverts); 
// 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);
  }
}
// Provide new data for this graph to our instance of
// class DESloopmomenta.
loopmomenta.givenewgraph(pnew, indexA); 
//
DESreal4vector Q[maxsize];
loopmomenta.getQs(Q);
DEScmplx4vector Qc[maxsize]; // Complex version.
for(i = 1; i <= nverts; i++)
{
  for(mu = 0; mu <= 3; mu++) Qc[i].in(mu) = Q[i](mu); // Type convesion.
}
DESreal4vector ell;     // This is the loop momentum before deformation.
DESdeform elldeform(nverts, indexA, Q);  // Class member for deforming ell.
DEScmplx4vector lc;
DEScmplx4vector Qi,ki;
complex<double> deformjacobian;
double rhoell;
complex<double> numerator, denominator;
// Preparation for this graph number completed.
double prefactor = 16.0 * pow(rts,nverts-4) * pow(4.0*pivalue, 0.5*nverts - 4.0);
double weightG = double(pointsG[graphnumber])/double(totalpoints);
// Loop over points for a given graph.
for (int npoint = 1; npoint <= pointsG[graphnumber]; npoint++)
  {
  loopmomenta.newell(ell);                   // New point ell.
  method = loopmomenta.getmethod();          // The method used to get this point.
  Ntries[method]++ ;
  bool goodpoint = true;
  if( ell(0)*ell(0) + ell(1)*ell(1) + ell(2)*ell(2) + ell(3)*ell(3) > 10.0e10 )
  {
    goodpoint = false;
  }
  for(i = 1; i <= nverts; i++)
  {
    if(abs(square(ell-Q[i])) < virtualitycut) goodpoint = false;
  }
  if( goodpoint )
  { // Execute main code if point is good.
  rhoell = loopmomenta.rhoell();             // rho_l(ell).
  elldeform.deform(ell, lc, deformjacobian); // Deformed ell = lc, with jacobian.
  // The numerator.
  if(vertextype == VECTOR)
  {
    if (zeromass)
    {
      numerator = fermion.NumeratorV(lc, Qc);
    }
    else
    {
      numerator = fermion.NumeratorM(lc, Qc, mass); 
    }      
  }
  else if (vertextype == LEFTHANDED)
  {
    numerator = fermion.NumeratorL(lc, Qc); 
  }
  else if (vertextype == RIGHTHANDED)
  {
    numerator = fermion.NumeratorR(lc, Qc);
  }
  else
  {
    cout << "Problem with the vertex type." << endl;
    exit(1);
  }
  // The denominator.
  denominator = 1.0;
  for (i = 1; i <= nverts; i++)
  {
    ki = lc - Qc[i]; // The momuntum in leg i.
    denominator *= dot(ki,ki);
  }
  value = prefactor /rhoell /weightG * deformjacobian * numerator / denominator;
  //
  for (i = 1; i <= nverts; i++)
  {
    ki = lc - Qc[i]; // The momuntum in leg i.
    complex<double> kisq;
    kisq = dot(ki,ki);
    // Check for bad deformation.
    if( (imag(kisq) < - 10*abs(real(kisq))) && (abs(kisq) < 0.00001*rtssq))
    {
      cout << "i = " << i << " imag(kisq) = "<< imag(kisq) 
           << " for real(kisq) = " << real(kisq) << endl;
      cout << ell(0) << " " << lc(0) << endl;
      cout << ell(1) << " " << lc(1) << endl;
      cout << ell(2) << " " << lc(2) << endl;
      cout << ell(3) << " " << lc(3) << endl << endl;
      cout << "i = 1:   (ell - Q[i])^2 = " << square(ell - Q[1]) 
           << "  (lc - Q[i])^2 = " << square(lc - Qc[1]) << endl;
      cout << "i = N:   (ell - Q[i])^2 = " << square(ell - Q[nverts]) 
           << "  (lc - Q[i])^2 = " << square(lc - Qc[nverts]) << endl;
      cout << "i = A:   (ell - Q[i])^2 = " << square(ell - Q[indexA]) 
           << "  (lc - Q[i])^2 = " << square(lc - Qc[indexA]) << endl;
      cout << "i = A+1: (ell - Q[i])^2 = " << square(ell - Q[indexA+1]) 
           << "  (lc - Q[i])^2 = " << square(lc - Qc[indexA+1]) << endl;
      cout << "|value| is " << abs(value) << endl;
      cout << "Graph number " << graphnumber<< " index A = "<< indexA
           <<", event number " << npoint << endl << endl;
    }
  }    
  //
  if (! (abs(value) < worst))
  {
    worst = abs(value);
    cout << "Worst |value| so far." << endl;
    cout << "For graph number " << graphnumber<< ", event number "
         << npoint << " |value| = "<< abs(value) << endl;
    cout << "ell     deformed lc " << endl;
    cout << ell(0) << " " << lc(0) << endl;
    cout << ell(1) << " " << lc(1) << endl;
    cout << ell(2) << " " << lc(2) << endl;
    cout << ell(3) << " " << lc(3) << endl;
    cout << "deform jacobian is " << deformjacobian <<  endl;
    cout << "rhoell*rts^4 is " << rhoell*pow(rts,4) <<  endl;
    cout << "numerator/rts^nverts is " << numerator/pow(rts,nverts) <<  endl;
    cout << "denominator/rts^(2*nverts) is " << denominator/pow(rts,2*nverts) <<  endl;
    for (i = 1; i <= nverts; i++)
    {
      ki = lc - Qc[i]; // The momuntum in leg i.
      complex<double> kiEsq;
      kiEsq = ki(0)*ki(0) + ki(1)*ki(1) +ki(2)*ki(2) + ki(3)*ki(3) ; 
      cout << i << ":  (lc - Qc[i])^2 = " << dot(ki,ki)
           << "   (lc - Qc[i])_E^2 = " << kiEsq << endl;
    }
    cout << endl;
    //
    for (i = 1; i <= nverts; i++)
    {
      cout << i << ":  (ell - Q[i])^2 = " << dot(ell - Q[i],ell - Q[i]) << endl;
    }
    cout << endl;
    //
    for (j = 1; j <= 0; j++) // This code not used.
    //for (j = 1; j <= nverts; j++)
    {
      DEScmplx4vector lcX;
      complex<double> deformjacobianX;
      // Deformed ell = lcX, with jacobianX for contribution c_j.
      cout << "Contribution from c_j for j = "<< j << endl;
      elldeform.deformtest(ell, lcX, deformjacobianX,j);
      for (i = 1; i <= nverts; i++)
      {
        ki = lcX - Qc[i]; // The momuntum in leg i.
        complex<double> kiEsq;
        kiEsq = ki(0)*ki(0) + ki(1)*ki(1) +ki(2)*ki(2) + ki(3)*ki(3) ; 
        cout << i << ":  (lc - Qc[i])^2 = " << dot(ki,ki)
          << "   (lc - Qc[i])_E^2 = " << kiEsq << endl;
      }
      cout << endl;        
    }
  }
  // Also, we calculate a test integral.
  lcsqE = 0.0;
  for( mu = 0; mu <= 3; mu++) lcsqE = lcsqE + pow(lc(mu),2);
  ftest = exp(- lcsqE/Mtestsq)/pow(pivalue*Mtestsq,2);  // This integrates to 1.
  valuetest = ftest/rhoell /weightG * deformjacobian
               /double(graphnumbermax - graphnumbermin + 1);
  } // End of execution of main code for goodpoint = true.
  else
  {
    value = 0.0;
    valuetest = 0.0;
    Nbad[method]++ ;
  }
  //  
  // Now we have value.
  sum = sum + value;
  sumsqR = sumsqR + pow(real(value),2);
  sumsqI = sumsqI + pow(imag(value),2);
  fluct[method] = fluct[method] + norm(value);
  sumGraph = sumGraph + value*weightG;
  sumsqRGraph = sumsqRGraph + pow(real(value*weightG),2);
  sumsqIGraph = sumsqIGraph + pow(imag(value*weightG),2);
  //
  sumtest += valuetest;
  sumsqRtest += pow(real(valuetest),2);
  sumsqItest += pow(imag(valuetest),2);
 //
} // Close for (npoint = 1; npoint <= nreno; npoint++)
  complex<double> resultGraph = sumGraph/double(pointsG[graphnumber]);
  double errorRGraph = sqrt( (sumsqRGraph/pointsG[graphnumber] 
                           - pow(real(resultGraph),2)) /(pointsG[graphnumber] - 1) );
  double errorIGraph = sqrt ( (sumsqIGraph/pointsG[graphnumber] 
                           - pow(imag(resultGraph),2)) /(pointsG[graphnumber] - 1) );
  double smallestdot = abs(dot(pnew[nverts],pnew[1]));
  double temp;
  for (i = 1; i < nverts; i++)
  { 
    temp = abs(dot(pnew[i],pnew[i+1]));
    if(temp < smallestdot) smallestdot = temp;
  }
// Printout for each graph, if desired.
//  cout << graphnumber 
//       << " A = " << indexA << " MTsq = " << loopmomenta.getMTsq()
//       << " Min p[i]*p[i+1] = " << smallestdot
//       << "  integral = " << resultGraph 
//       << " +/- ("<< errorRGraph << "," << errorIGraph << ")" << endl;
}// Close for(int graphnumber = graphnumbermin; graphnumber <= graphnumbermax; graphnumber++)
// Calculate results.
complex<double> result;
double errorR, errorI;
//
double npoints = double(totalpoints);
result = sum/npoints;
errorR = sqrt( (sumsqR/npoints - pow(real(result),2))/(npoints - 1) );
errorI = sqrt( (sumsqI/npoints - pow(imag(result),2))/(npoints - 1) );
//
complex<double> resulttest = sumtest/npoints;
double errorRtest = sqrt( (sumsqRtest/npoints - pow(real(resulttest),2))/(npoints - 1) );
double errorItest = sqrt( (sumsqItest/npoints - pow(imag(resulttest),2))/(npoints - 1) );
// Fluctuation report.
double fluctnorm = 0.0;
for (i = 1; i <= maxindex; i++)
{
  fluct[i] = fluct[i]/npoints;
  beta[i] = double(Ntries[i])/npoints;
  if (beta[i] > 0.0) fluctnorm += fluct[i];
}
cout << " " << endl;
cout << "Sampling methods." << endl;
cout << "methods in: "<< endl;
for(i = 1; i <= maxindex; i++)
  {
  if (alpha[i] > 0.0) cout << i << ": "<<alpha[i] << " " ;
  }
cout << endl << "methods actual: "<< endl;
for(i = 1; i <= maxindex; i++)
  {
  if (beta[i] > 0.0) cout << i << ": "<< beta[i] << " " ;
  }
cout << endl << "fraction of bad points: "<< endl;
for(i = 1; i <= maxindex; 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 <= maxindex; i++)
  {
  if (beta[i] > 0.0) cout << i << ": "<< fluct[i]/beta[i]/fluctnorm << " " ;
  }
cout << endl;
cout << " " << endl;
// End fluctuation report.
//
cout << "test integral = " << resulttest << " +/- ("<< errorRtest << "," << errorItest << ")" 
     << endl << endl;
//
cout << "integral = " << result << " +/- ("<< errorR << "," << errorI << ")" << endl;
cout << "|integral| = " << abs(result) << " +/- "
     << (abs(real(result)*errorR) + abs(imag(result)*errorI))/abs(result)
     << endl << endl;
//
return 0;
}  //end of main()
// ---------------------------------------------------