//34567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890
// Subroutines "direct" numerical integration of one loop Feynman diagrams.
// D.E. Soper
// See date below.
//
#include "NphotonD.h"
#include <cmath>
#include <iostream>
#include <cassert>
using namespace std;
//
void giveversion()
{
cout << "            Subroutines for NphotonD version 1.0" << endl;
cout << "                      11 December 2008" << endl;
}
//==================================================================================================
// Implementation code for class DESloopmomenta.
//
// Constructor.
DESloopmomenta::DESloopmomenta(int seedIN, int nvertsIN, DESreal4vector pIN[], 
                         int indexAIN, double alphaIN[])
: r(seedIN)     // Construct random number generator instance r.
, indx(nvertsIN)  // Construct index function instance indx.
{
seed = seedIN;
nverts = nvertsIN;
double aI,aII;
DESreal4vector vI,vII,v;
int i, mu;
//
for(i = 1; i <= nverts; i++) p[i] = pIN[i];
indexA = indexAIN;
alpha[0] = 0.0;
beta[0] = 0.0;
for(i = 1; i < maxsize; i++)
{
  alpha[i] = alphaIN[i];
  beta[i] = beta[i-1] + alpha[i-1]; // Thus beta[i] is the probability for method < i.
}
indexAp1 = indx(indexA + 1);
for(mu = 0; mu <= 3; mu++) Q[1].in(mu) = 0.0;  // This is a convenient choice to start.
for (i = 1; i <= nverts - 1; i++)
{
  Q[i+1] = Q[i] + p[i];
}
Pplain = Q[nverts] - Q[1];
Pbar = Q[indexA] - Q[indexAp1];
dotPPbar = dot(Pplain, Pbar);
// Some checks.
if(Pplain(0) <= 0.0) cout << "NONSTANDARD: Pplain(0) <= 0.0" << endl;
if(Pplain(3) <= 0.0) cout << "NONSTANDARD: Pplain(3) <= 0.0" << endl;
if( pow(Pplain(1),2) + pow(Pplain(2),2) > 0.0000001 )
                      cout << "NONSTANDARD: Pplain_T != 0.0" << endl;
if(Pbar(0) <= 0.0) cout << "NONSTANDARD: Pbar(0) <= 0.0" << endl;
if(Pbar(3) >= 0.0) cout << "NONSTANDARD: Pbar(3) >= 0.0" << endl;
if( pow(Pbar(1),2) + pow(Pbar(2),2) > 0.0000001 )
cout << "Pbar_T != 0.0" << endl;
// End checks.
tiny = 1.0e-30;
small = 1.0e-10;
huge = 1.0e30;
aI = dot((Q[nverts] - Q[indexAp1]),Pbar)/dotPPbar;
aII = dot((Q[indexA] - Q[1]),Pplain)/dotPPbar;
vI = aI*Q[1] + (1 - aI)*Q[nverts];
vII =  aII*Q[indexAp1] + (1 - aII)*Q[indexA];
MTsq = - square(vI - vII);
if(indexA == 1 || indexAp1 == nverts) MTsq = 2.0*dotPPbar/double(nverts*nverts);
v = 0.5*(vI + vII);
for (i = 1; i <= nverts; i++) Q[i] = Q[i] - v; // Reset origin of loop momenta.
// Parameters for sampling methods.
// Method 1;
Awide = double(nverts)/4.0;        // Presumably this is a good choice.
if(nverts == 4) Awide = 1.5;       // Here we are helped by renoralization.
// Methods 2 and 9  (near double parton scattering singularity).
// We use a large range for version a and a small range for version b.
msmall4_a = pow(MTsq,2);
mlarge4_a = pow(dotPPbar,2);
if( mlarge4_a < 2.0*msmall4_a ) mlarge4_a = 2.0*msmall4_a;
msmall4_b = 0.5*pow(MTsq,2);
mlarge4_b = 4.0*msmall4_b;
// Methods 3 and 10 (near soft singularities).
// We use a more peaked function for version a, less peaked for version b.
Asoft_a = 0.2;
Bsoft_a = 0.2;
Asoft_b = 0.5;
Bsoft_b = 0.5;
// Methods 4 and 11 (near collinear singularities).
// We use a more peaked function for version a, less peaked for version b.
Acoll_a = 0.2;
Bcoll_a = 0.2;
Acoll_b = 0.5;
Bcoll_b = 0.8;
// Methods 5 and 12 (near the initial state parton collinear singularities).
// We use a more peaked function for version a, less peaked for version b.
AIcoll_a = 0.6;
BIcoll_a = 0.2;
AIcoll_b = 1.0;
BIcoll_b = 0.8;
// Methods 6 and 13 (near intersections of time-like separated cones).
// We use a more peaked function for version a, less peaked for version b.
Atlike_a = 0.3;
Atlike_b = 0.6;
// Methods 7 and 14 (near intersections of space-like separated cones).
// We use a more peaked function for version a, less peaked for version b.
Aslike_a = 0.3;
Aslike_b = 0.6;
// Methods 8 and 15 (near intersections of light-like separated cones).
// We use a more peaked function for version a, less peaked for version b.
Allike_a = 0.8;
Bllike_a = 0.2;
Msqllike_a = 0.001*dotPPbar;
Allike_b = 0.4;
Bllike_b = 0.5;
Msqllike_b = 0.01*dotPPbar;
// Methods 16 and 17 (near intersections of light-like separated cones plus third cone).
// We use a more peaked function for version a, less peaked for version b.
tildeallike3_a = 10.0*MTsq/dotPPbar;
if( tildeallike3_a > 0.1 ) tildeallike3_a = 0.1;
Allike3_a = 0.2;
Bllike3_a = 0.2;
Msqllike3_a = 0.001*dotPPbar;
tildeallike3_b = 0.5;
Allike3_b = 0.6;
Bllike3_b = 0.5;
Msqllike3_b = 0.01*dotPPbar;
//  
}
// There is no default constructor.
//--------------------------------------------------------------------------------------------------
void DESloopmomenta::givenewgraph(DESreal4vector pIN[], int indexAIN)
{
double aI,aII;
DESreal4vector vI,vII,v;
int i, mu;
for(i = 1; i <= nverts; i++) p[i] = pIN[i];
indexA = indexAIN;
indexAp1 = indx(indexA + 1);
for(mu = 0; mu <= 3; mu++) Q[1].in(mu) = 0.0;  // This is a convenient choice to start.
for (i = 1; i <= nverts - 1; i++)
{
  Q[i+1] = Q[i] + p[i];
}
Pplain = Q[nverts] - Q[1];
Pbar = Q[indexA] - Q[indexAp1];
dotPPbar = dot(Pplain, Pbar);
// Some checks.
if(Pplain(0) <= 0.0) cout << "NONSTANDARD: Pplain(0) <= 0.0" << endl;
if(Pplain(3) <= 0.0) cout << "NONSTANDARD: Pplain(3) <= 0.0" << endl;
if( pow(Pplain(1),2) + pow(Pplain(1),2) > 0.0000001 )
                      cout << "NONSTANDARD: Pplain_T != 0.0" << endl;
if(Pbar(0) <= 0.0) cout << "NONSTANDARD: Pbar(0) <= 0.0" << endl;
if(Pbar(3) >= 0.0) cout << "NONSTANDARD: Pbar(3) >= 0.0" << endl;
if( pow(Pbar(1),2) + pow(Pbar(1),2) > 0.0000001 )
cout << "Pbar_T != 0.0" << endl;
// End checks.
aI = dot((Q[nverts] - Q[indexAp1]),Pbar)/dotPPbar;
aII = dot((Q[indexA] - Q[1]),Pplain)/dotPPbar;
vI = aI*Q[1] + (1 - aI)*Q[nverts];
vII =  aII*Q[indexAp1] + (1 - aII)*Q[indexA];
MTsq = - square(vI - vII);
if(indexA == 1 || indexAp1 == nverts) MTsq = 2.0*dotPPbar/double(nverts*nverts);
v = 0.5*(vI + vII);
for (i = 1; i <= nverts; i++) Q[i] = Q[i] - v; // Reset origin of loop momenta.
return;
}
//--------------------------------------------------------------------------------------------------
void DESloopmomenta::givepoint(DESreal4vector ellIN)
{
  ell= ellIN;
  return;
}
//--------------------------------------------------------------------------------------------------
void DESloopmomenta::getpoint(DESreal4vector& ellOUT)
{
  ellOUT = ell;
  return;
}
//--------------------------------------------------------------------------------------------------
int DESloopmomenta::getmethod()
{
  return method;
}
//--------------------------------------------------------------------------------------------------
double DESloopmomenta::getMTsq()
{
  return MTsq;
}
//--------------------------------------------------------------------------------------------------
int DESloopmomenta::getindexA()
{
  return indexA;
}
//--------------------------------------------------------------------------------------------------
void DESloopmomenta::getQs(DESreal4vector Qout[maxsize])
{
  for(int i = 1; i <= nverts; i++) Qout[i] = Q[i];
  return;
}
//--------------------------------------------------------------------------------------------------
void DESloopmomenta::getps(DESreal4vector pout[maxsize])
{
  for(int i = 1; i <= nverts; i++) pout[i] = p[i];
  return;
}
//--------------------------------------------------------------------------------------------------
// Generate new point ell.
// Here ellOUT is a real vector {ell[0], ell[1], ell[2], ell[3]}.
void DESloopmomenta::newell(DESreal4vector& ellOUT)
{
static const double twoPI = 6.283185307179586;
static const double PIover4 = 0.7853981633974483;
if(r.random() < alpha[1])
{
  //
  // Method 1: roughly uniform over a wide range.
  //
  // First construct point at random on unit sphere in 4 dimensions.
  // d^3 \Omega = (1/2) d\theta_A d\theta_B d cos^2(theta_C) where
  // 0 < theta_A < 2 Pi, 0 < theta_A < 2 Pi, 0 < theta_C < Pi/2 and
  // ell[0] = cosC*cosA; ell[1] = cosC*sinA; ell[2] = sinC*cosB; ell[3] = sinC*sinB;
  // The density of points is 1/area = 1/(2 Pi^2).
  //  
  double thetaA = r.random()*twoPI;
  double sinA = sin(thetaA);
  double cosA = cos(thetaA);
  double thetaB = r.random()*twoPI;
  double sinB = sin(thetaB);
  double cosB = cos(thetaB);
  double temp = r.random();
  double cosC = sqrt(temp);
  double sinC = sqrt(1.0 - temp);
  ell.in(0) = cosC*cosA;
  ell.in(1) = cosC*sinA;
  ell.in(2) = sinC*cosB;
  ell.in(3) = sinC*sinB;
  //
  // Now fix the length of ell.
  double ellE4;
  ellE4 = pow(r.random(), -1.0/(Awide - 1.0)) - 1.0;
  ellE4 = ellE4 * pow(dotPPbar,2);
  temp = pow(ellE4,0.25);
  ell = temp*ell;
  ellOUT = ell;
  method = 1;
  return;
} // End method 1.
else if (r.random() < alpha[2]/(1.0 - beta[2] + small))
{
  //
  // Method 2: distributed near double parton scattering singularity.
  //  
  double thetaA = r.random()*twoPI;
  double sinA = sin(thetaA);
  double cosA = cos(thetaA);
  double thetaB = r.random()*twoPI;
  double sinB = sin(thetaB);
  double cosB = cos(thetaB);
  double temp = r.random();
  double cosC = sqrt(temp);
  double sinC = sqrt(1.0 - temp);
  ell.in(0) = cosC*cosA;
  ell.in(1) = cosC*sinA;
  ell.in(2) = sinC*cosB;
  ell.in(3) = sinC*sinB;
  //
  // Now fix the length of ell.
  double ellE4;
  double ratio = pow(mlarge4_a/msmall4_a, r.random());
  ellE4 = (mlarge4_a - msmall4_a*ratio)/(ratio - 1.0);
  temp = pow(ellE4,0.25);
  ell = temp*ell;
  ellOUT = ell;
  method = 2;
  return;
} // End method 2.
else if (r.random() < alpha[3]/(1.0 - beta[3] + small))
{
  //
  // Method 3: distributed near the Q[i] soft singularity.
  //
  int ivert = int(r.random()*double(nverts) + 1.0); // Equal probabilities for ivert = 1, ... , nverts.
  DESreal4vector Pplus  = Q[indx(ivert+1)] - Q[ivert];
  if(Pplus(0) < 0.0)  Pplus  = - Pplus;   // So Pplus is on the positive light cone.
  DESreal4vector Pminus = Q[indx(ivert-1)] - Q[ivert];
  if(Pminus(0) < 0.0) Pminus = - Pminus;  // So Pminus is on the positive light cone.
  // The sign choices for Pplus and Pminus make Pplus + Pminus always timelike and
  // Pplus - Pminus always spacelike, which is convenient for generating ell.
  // Choose x.
  double sx = 1.0;
  if(r.random() < 0.5) sx = -1.0;
  double x = sx*pow(r.random(),1/Asoft_a);
  // Choose y.
  double sy = 1.0;
  if(r.random() < 0.5) sy = -1.0;
  double y = sy*pow(r.random(),1/Asoft_a);
  // Choose phi.
  double phi = r.random()*twoPI;
  // Choose elllperp.
  double xyabs = abs(x*y);
  double temp = pow(1.0 + xyabs, Bsoft_a)/pow(xyabs, Bsoft_a) - 1.0;
  double lambdaperp = pow( temp*r.random() + 1.0, 1/Bsoft_a) - 1.0;
  double Misq = dot(Pplus,Pminus);  // This is positive with our definitions.
  double Psumabs = sqrt(2.0*Misq);  //This is sqrt( square(Pplus + Pminus) ).
  double ellTsq  = Misq*xyabs*lambdaperp;
  double ellTabs = sqrt(ellTsq);
  // Construct ellT in frame where Psum is in 0 direction and Pdiff is in the +/- 3 direction.
  // Subscripts 2 refer to this frame.
  // Subscripts 1 refer to frame with Psum is in 0 direction but Pdiff in another direction.
  // Then no subscripts refer to original frame.  
  DESreal4vector Psum  = Pplus + Pminus;
  DESreal4vector Pdiff = Pplus - Pminus;
  DESreal4vector Evec, Zvec; // This initializes them to zero.
  Evec.in(0) = Psumabs;     // So square(Evec) = square(Psum);
  DESreal4vector Pdiff_1 = boost(Evec, Psum, Pdiff);   // Pdiff in frame 1.
  if( Pdiff_1(3) > 0 )  Zvec.in(3) = Psumabs;
  else  Zvec.in(3) = - Psumabs;
  DESreal4vector ellT_2;  // This initializes it to zero.
  ellT_2.in(1) = ellTabs*sin(phi);
  ellT_2.in(2) = ellTabs*cos(phi);
  DESreal4vector ellT_1 = boost(Pdiff_1, Zvec, ellT_2); // ellT in frame 1.
  DESreal4vector ellT = boost(Psum, Evec, ellT_1);
  ell = Q[ivert] + x*Pplus + y*Pminus + ellT;
//  cout << "Make: x = " << x << " y = " << y << " ellTsq = " << ellTsq << endl;
  ellOUT = ell;
  method = 3;
  return;  
} // End method 3.
else if (r.random() < alpha[4]/(1.0 - beta[4] + small))
{
  //
  // Method 4: distributed near the Q[i] - Q[i+1] collinear singularity.
  //
  bool reverse = false;
  if( r.random() > 0.5 ) reverse = true; // Exchange roles of i and i+1.
  int ivert = int(r.random()*double(nverts) + 1.0); // Equal probabilities for ivert = 1, ... , nverts.
  DESreal4vector Pplus;
  DESreal4vector Pminus;
  if(!reverse)
  {
    Pplus  = Q[indx(ivert+1)] - Q[ivert];
    Pminus = Q[indx(ivert-1)] - Q[ivert];  
  }
  else // if(reverse)
  {
    Pplus = Q[ivert] - Q[indx(ivert+1)];
    Pminus = Q[indx(ivert+1)] - Q[indx(ivert+2)];
  }
  if(Pplus(0)*Pminus(0) < 0.0)  Pminus = - Pminus;   // So dot[Pplus,Pminus] > 0.
  double Misq = dot(Pplus,Pminus);  // This is positive with our definitions.
  double Psumabs = sqrt(2.0*Misq);  // This is sqrt( square(Pplus + Pminus) ).
  // Choose x.
  double x = pow(r.random(), 1.0/Acoll_a);
  // Choose sign of y.
  double sy = 1.0;
  if(r.random() < 0.5) sy = -1.0;
  // Choose phi.
  double phi;
  phi = r.random()*twoPI;
  // Choose random variables rz and rdmod = rd^{1/B}.
  double rz = r.random();
  double rdmod = pow(r.random(), 1.0/Bcoll_a);
  // Choose ellTsq and absy.
  double ellTsq = Misq*rz*rdmod;
  double absy = (1.0 - rz)*rdmod/x;
  double ellTabs = sqrt(ellTsq);
  double y = sy*absy;
  // Construct ellT in frame where Psum is in 0 direction and Pdiff is in the +/- 3 direction.
  // Subscripts 2 refer to this frame.
  // Subscripts 1 refer to frame with Psum is in 0 direction but Pdiff in another direction.
  // Then no subscripts refer to original frame.  
  DESreal4vector Psum  = Pplus + Pminus;
  DESreal4vector Pdiff = Pplus - Pminus;
  DESreal4vector Evec, Zvec; // This initializes them to zero.
  if(Psum(0) > 0) Evec.in(0) = Psumabs;     // So square(Evec) = square(Psum);
  else Evec.in(0) = - Psumabs;              // So Evec(0) has same sign as Psum(0);
  DESreal4vector Pdiff_1 = boost(Evec, Psum, Pdiff);   // Pdiff in frame 1.
  if( Pdiff_1(3) > 0 )  Zvec.in(3) = Psumabs;
  else  Zvec.in(3) = - Psumabs;
  DESreal4vector ellT_2;  // This initializes it to zero.
  ellT_2.in(1) = ellTabs*sin(phi);
  ellT_2.in(2) = ellTabs*cos(phi);
  DESreal4vector ellT_1 = boost(Pdiff_1, Zvec, ellT_2); // ellT in frame 1.
  DESreal4vector ellT = boost(Psum, Evec, ellT_1);
  if(!reverse)
  {
   ell = Q[ivert] + x*Pplus + y*Pminus + ellT; 
  }
  else // if(reverse)
  {
    ell = Q[indx(ivert+1)] + x*Pplus + y*Pminus + ellT;
  }
  ellOUT = ell;
  method = 4;
  return;  
} // End method 4.
else if (r.random() < alpha[5]/(1.0 - beta[5] + small))
{
  //
  // Method 5: distributed near the Q[i] - Q[i+1] collinear singularity
  // for the initial state momenta, i = nverts and i = indexA.
  //
  bool reverse = false;
  if( r.random() > 0.5 ) reverse = true; // Exchange roles of i and i+1.
  int ivert;
  if ( r.random() > 0.5 ) ivert = nverts;
  else ivert = indexA;
  DESreal4vector Pplus;
  DESreal4vector Pminus;
  if(!reverse)
  {
    Pplus  = Q[indx(ivert+1)] - Q[ivert];
    Pminus = Q[indx(ivert-1)] - Q[ivert];  
  }
  else // if(reverse)
  {
    Pplus = Q[ivert] - Q[indx(ivert+1)];
    Pminus = Q[indx(ivert+1)] - Q[indx(ivert+2)];
  }
  if(Pplus(0)*Pminus(0) < 0.0)  Pminus = - Pminus;   // So dot[Pplus,Pminus] > 0.
  double Misq = dot(Pplus,Pminus);  // This is positive with our definitions.
  double Psumabs = sqrt(2.0*Misq);  //This is sqrt( square(Pplus + Pminus) ).
                                    // Choose x.
  double x = pow(r.random(), 1.0/AIcoll_a);
  // Choose sign of y.
  double sy = 1.0;
  if(r.random() < 0.5) sy = -1.0;
  // Choose phi.
  double phi;
  phi = r.random()*twoPI;
  // Choose random variables rz and rdmod = rd^{1/B}.
  double rz = r.random();
  double rdmod = pow(r.random(), 1.0/BIcoll_a);
  // Choose ellTsq and absy.
  double ellTsq = Misq*rz*rdmod;
  double absy = (1.0 - rz)*rdmod/x;
  double ellTabs = sqrt(ellTsq);
  double y = sy*absy;
  // Construct ellT in frame where Psum is in 0 direction and Pdiff is in the +/- 3 direction.
  // Subscripts 2 refer to this frame.
  // Subscripts 1 refer to frame with Psum is in 0 direction but Pdiff in another direction.
  // Then no subscripts refer to original frame.  
  DESreal4vector Psum  = Pplus + Pminus;
  DESreal4vector Pdiff = Pplus - Pminus;
  DESreal4vector Evec, Zvec; // This initializes them to zero.
  if(Psum(0) > 0) Evec.in(0) = Psumabs;     // So square(Evec) = square(Psum);
  else Evec.in(0) = - Psumabs;              // So Evec(0) has same sign as Psum(0);
  DESreal4vector Pdiff_1 = boost(Evec, Psum, Pdiff);   // Pdiff in frame 1.
  if( Pdiff_1(3) > 0 )  Zvec.in(3) = Psumabs;
  else  Zvec.in(3) = - Psumabs;
  DESreal4vector ellT_2;  // This initializes it to zero.
  ellT_2.in(1) = ellTabs*sin(phi);
  ellT_2.in(2) = ellTabs*cos(phi);
  DESreal4vector ellT_1 = boost(Pdiff_1, Zvec, ellT_2); // ellT in frame 1.
  DESreal4vector ellT = boost(Psum, Evec, ellT_1);
  if(!reverse)
  {
    ell = Q[ivert] + x*Pplus + y*Pminus + ellT; 
  }
  else // if(reverse)
  {
    ell = Q[indx(ivert+1)] + x*Pplus + y*Pminus + ellT;
  }
  ellOUT = ell;
  method = 5;
  return;  
} // End method 5.
else if (r.random() < alpha[6]/(1.0 - beta[6] + small))
{
  //
  // Method 6: distributed near forward lightcone from Q[i]
  // and backward lightcone from Q[j] where Q[j] - Q[i] is timelike
  // with positive energy.
  //
  // Choose which case, (i,j), at random.
  int i = 0, j = 0;
  int numberofLcases = 0;
  if (indexA > 2) numberofLcases = (indexA - 1)*(indexA - 2)/2;
  int numberofRcases = 0;
  if (nverts > indexA + 2) numberofRcases = (nverts - indexA - 1)*(nverts - indexA - 2)/2;
  int numberofcases = numberofLcases + numberofRcases;
  int casenumber = 0;
  double randomN = numberofcases*r.random();
  if( randomN <= numberofLcases )
  {
    casenumber = 0;
    for(int itrial = 1; itrial <= indexA - 2; itrial++)
    {
      for(int jtrial = itrial + 2; jtrial <= indexA; jtrial++)
      {
        casenumber += 1;
        if (casenumber > randomN)  // We have found which case.
        {
          i = itrial;
          j = jtrial;
          goto ijfound6;
        }
      }
    }
  }
  else // if( randomN > numberofLcases )
  {
    casenumber = numberofLcases;
    for(int itrial = indexA + 1; itrial <= nverts - 2; itrial++)
    {
      for(int jtrial = itrial + 2; jtrial <= nverts; jtrial++)
      {
        casenumber += 1;
        if (casenumber > randomN)  // We have found which case.
        {
          i = itrial;
          j = jtrial;
          goto ijfound6;
        }
      }
    }
  }
  cout << "(i,j) not found in newell method 6" << endl;
  exit(1);
  ijfound6: ; // null statement: start here when (i,j) is found.
  DESreal4vector K = Q[j] - Q[i];
  double Ksq = square(K);
  if(Ksq <= 0.0)
  {
    cout << "This is horrible, Ksq <= 0 in newell method 6." << endl;
    exit(1);
  }
  double absK = sqrt(Ksq);
  DESreal4vector Kprime(absK,0.0,0.0,0.0);
  DESreal4vector lprime;
  double costheta = 2.0*(r.random() - 0.5);
  double sintheta = sqrt(1.0 - costheta*costheta);
  double phi = twoPI*r.random();
  double lsq = 0.25*Ksq*pow(r.random(),1.0/Atlike_a);
  if( r.random() < 0.5 ) lsq = - lsq;
  double lmKsq = 0.25*Ksq*pow(r.random(),1.0/Atlike_a);
  if( r.random() < 0.5 ) lmKsq = - lmKsq;
  lprime.in(0) = (lsq - lmKsq + Ksq)/2.0/absK;
  double absveclprime = sqrt( pow(lprime(0),2) - lsq );
  lprime.in(1) = absveclprime*sintheta*cos(phi);
  lprime.in(2) = absveclprime*sintheta*sin(phi);
  lprime.in(3) = absveclprime*costheta;
  // We have the result in the special frame.
  // Now we boost it back to the c.m. frame
  // and add Q[i] to get ell.
  ell = boost(K,Kprime,lprime) + Q[i];
  ellOUT = ell;
  method = 6;
  return;  
} // End method 6.
else if (r.random() < alpha[7]/(1.0 - beta[7] + small))
{
  //
  // Method 7: distributed near lightcone from Q[i]
  // and lightcone from Q[j] where Q[j] - Q[i] is spacelike.
  //
  // Choose which case, (i,j), at random.
  int i = 0, j = 0;
  int numberofcases = indexA*(nverts - indexA) - 2;
  int casenumber = 0;
  double randomN = numberofcases*r.random();
  for(int itrial = 1; itrial <= indexA; itrial++)
  {
    for(int jtrial = indexA + 1; jtrial <= nverts; jtrial++)
    {
      if( !( (   (itrial == 1)    && (jtrial == nverts)    )   ||
             ( (itrial == indexA) && (jtrial == indexA + 1)) ) )
      {
        casenumber +=1;
        if (casenumber > randomN)  // We have found which case.
        {
          i = itrial;
          j = jtrial;
          goto ijfound7;
        }
      }
    }
  }
  cout << "(i,j) not found in newell method 7" << endl;
  exit(1);
  ijfound7: ;
  DESreal4vector K = Q[j] - Q[i];
  double Ksq = square(K);
  if(Ksq >= 0.0)
  {
    cout << "This is horrible, Ksq >= 0 in newell method 7." << endl;
    exit(1);
  }
  double absK = sqrt(-Ksq);  // Note the sign.
  DESreal4vector Kprime(0.0,0.0,0.0,absK);
  DESreal4vector lprime;
  double phi = twoPI*r.random();
  double coshomega = tan( PIover4*(1.0 + r.random()) );
  double sinhomega = sqrt(coshomega*coshomega - 1.0);
  double lsq = 0.25*Ksq*pow(r.random(),1.0/Aslike_a);
  if( r.random() < 0.5 ) lsq = - lsq;
  double lmKsq = 0.25*Ksq*pow(r.random(),1.0/Aslike_a);
  if( r.random() < 0.5 ) lmKsq = - lmKsq;
  lprime.in(3) = (lmKsq - lsq - Ksq)/2.0/absK;
  double tildelsq = lsq + pow(lprime(3),2);
  if(tildelsq < 0.0)
  {
    cout << "tildelsq negative in newell method 7." << endl;
    exit(1);
  }
  double abstildel = sqrt(tildelsq);
  double Eprime = abstildel*coshomega;
  if(r.random() < 0.5) Eprime = - Eprime;
  lprime.in(0) = Eprime;
  lprime.in(1) = abstildel*sinhomega*cos(phi);
  lprime.in(2) = abstildel*sinhomega*sin(phi);
  // We have the result in the special frame.
  // Now we boost it back to the c.m. frame
  // and add Q[i] to form ell.
  ell = boost(K,Kprime,lprime) + Q[i];
  ellOUT = ell;
  method = 7;
  return;  
} // End method 7.
else if (r.random() < alpha[8]/(1.0 - beta[8] + small))
{
  //
  // Method 8: distributed near lightcone from Q[i]
  // and lightcone from Q[j] where Q[j] - Q[i] is lightlike.
  //
  // Choose which case, (i,j), at random.
  // Switch i and j if needed to make Q[j](0) - Q[i](0) > 0.
  int i = 0, j = 0;
  int ii = ceil(r.random()*nverts);
  if( ii == nverts )
  {
    i = 1;
    j = nverts;
  }
  else if( ii == indexA )
  {
    i = indexA + 1;
    j = indexA;
  }
  else
  {
    i = ii;
    j = ii + 1;
  }
  // Define vector K.
  DESreal4vector K = Q[j] - Q[i];
  // Let vecK be the 3-vector part of K.
  DESreal4vector vecK = K;
  vecK.in(0) = 0.0;
  // Define vecK rotated so that vec K is in +/- z-direction.
  DESreal4vector vecKprime(0.0,0.0,0.0,K(0));
  if( K(3) < 0.0 ) vecKprime.in(3) = - K(0);
  // We construct l_2, in the "prime" system in which
  // K^+ > 0 and all other components of K vanish.
  double temp = 2.0*r.random() - 1.0;
  double x = 0.5*pow(abs(temp),1.0/(1.0 - Allike_a));
  if( temp < 0.0 ) x = - x;
  if( r.random() < 0.5 ) x = 1.0 - x; // 1/2 < x < 3/2 choice
  temp =  2.0*r.random() - 1.0;
  double lsq = Msqllike_a/( pow(abs(temp),-1.0/Bllike_a) - 1.0 );
  if( temp < 0.0 ) lsq = - lsq;
  temp =  2.0*r.random() - 1.0;
  double lmKsq = Msqllike_a/( pow(abs(temp),-1.0/Bllike_a) - 1.0 );
  if( temp < 0.0 ) lmKsq = - lmKsq;
  double lperpsq = -(1.0 - x)*lsq - x*lmKsq;
  // We need lperpsq > 0. If it is not, we could choose a new point,
  // but it is faster to just reverse lsq and lmKsq.
  if(lperpsq < 0.0)
  {
    lsq = - lsq;
    lmKsq = - lmKsq;
    lperpsq = - lperpsq;
  }
  double lperpabs = sqrt(lperpsq);
  double phi = twoPI*r.random();
  DESreal4vector lprime;
  temp = (lsq - lmKsq)/4.0/K(0);
  lprime.in(0) = x*K(0) + temp;
  lprime.in(1) = lperpabs*cos(phi);
  lprime.in(2) = lperpabs*sin(phi);
  if( K(3) > 0.0 ) lprime.in(3) = x*K(0) - temp;
  else lprime.in(3) = - x*K(0) + temp;
  // Now we rotate it back to the c.m. frame
  // and add Q[i] to form ell.
  ell = boost(vecK,vecKprime,lprime) + Q[i];
  ellOUT = ell;
  method = 8;
  return;
} // End method 8.
else if (r.random() < alpha[9]/(1.0 - beta[9] + small))
{
  //
  // Method 9: distributed near double parton scattering singularity.
  //  
  double thetaA = r.random()*twoPI;
  double sinA = sin(thetaA);
  double cosA = cos(thetaA);
  double thetaB = r.random()*twoPI;
  double sinB = sin(thetaB);
  double cosB = cos(thetaB);
  double temp = r.random();
  double cosC = sqrt(temp);
  double sinC = sqrt(1.0 - temp);
  ell.in(0) = cosC*cosA;
  ell.in(1) = cosC*sinA;
  ell.in(2) = sinC*cosB;
  ell.in(3) = sinC*sinB;
  //
  // Now fix the length of ell.
  double ellE4;
  double ratio = pow(mlarge4_b/msmall4_b, r.random());
  ellE4 = (mlarge4_b - msmall4_b*ratio)/(ratio - 1.0);
  temp = pow(ellE4,0.25);
  ell = temp*ell;
  ellOUT = ell;
  method = 9;
  return;
} // End method 9.
else if (r.random() < alpha[10]/(1.0 - beta[10] + small))
{
  //
  // Method 10: distributed near the Q[i] soft singularity.
  //
  int ivert = int(r.random()*double(nverts) + 1.0); // Equal probabilities for ivert = 1, ... , nverts.
  DESreal4vector Pplus  = Q[indx(ivert+1)] - Q[ivert];
  if(Pplus(0) < 0.0)  Pplus  = - Pplus;   // So Pplus is on the positive light cone.
  DESreal4vector Pminus = Q[indx(ivert-1)] - Q[ivert];
  if(Pminus(0) < 0.0) Pminus = - Pminus;  // So Pminus is on the positive light cone.
  // The sign choices for Pplus and Pminus make Pplus + Pminus always timelike and
  // Pplus - Pminus always spacelike, which is convenient for generating ell.
  // Choose x.
  double sx = 1.0;
  if(r.random() < 0.5) sx = -1.0;
  double x = sx*pow(r.random(),1/Asoft_b);
  // Choose y.
  double sy = 1.0;
  if(r.random() < 0.5) sy = -1.0;
  double y = sy*pow(r.random(),1/Asoft_b);
  // Choose phi.
  double phi = r.random()*twoPI;
  // Choose elllperp.
  double xyabs = abs(x*y);
  double temp = pow(1.0 + xyabs, Bsoft_b)/pow(xyabs, Bsoft_b) - 1.0;
  double lambdaperp = pow( temp*r.random() + 1.0, 1/Bsoft_b) - 1.0;
  double Misq = dot(Pplus,Pminus);  // This is positive with our definitions.
  double Psumabs = sqrt(2.0*Misq);  //This is sqrt( square(Pplus + Pminus) ).
  double ellTsq  = Misq*xyabs*lambdaperp;
  double ellTabs = sqrt(ellTsq);
  // Construct ellT in frame where Psum is in 0 direction and Pdiff is in the +/- 3 direction.
  // Subscripts 2 refer to this frame.
  // Subscripts 1 refer to frame with Psum is in 0 direction but Pdiff in another direction.
  // Then no subscripts refer to original frame.  
  DESreal4vector Psum  = Pplus + Pminus;
  DESreal4vector Pdiff = Pplus - Pminus;
  DESreal4vector Evec, Zvec; // This initializes them to zero.
  Evec.in(0) = Psumabs;     // So square(Evec) = square(Psum);
  DESreal4vector Pdiff_1 = boost(Evec, Psum, Pdiff);   // Pdiff in frame 1.
  if( Pdiff_1(3) > 0 )  Zvec.in(3) = Psumabs;
  else  Zvec.in(3) = - Psumabs;
  DESreal4vector ellT_2;  // This initializes it to zero.
  ellT_2.in(1) = ellTabs*sin(phi);
  ellT_2.in(2) = ellTabs*cos(phi);
  DESreal4vector ellT_1 = boost(Pdiff_1, Zvec, ellT_2); // ellT in frame 1.
  DESreal4vector ellT = boost(Psum, Evec, ellT_1);
  ell = Q[ivert] + x*Pplus + y*Pminus + ellT;
//  cout << "Make: x = " << x << " y = " << y << " ellTsq = " << ellTsq << endl;
  ellOUT = ell;
  method = 10;
  return;  
} // End method 10.
else if (r.random() < alpha[11]/(1.0 - beta[11] + small))
{
  //
  // Method 11: distributed near the Q[i] - Q[i+1] collinear singularity.
  //
  bool reverse = false;
  if( r.random() > 0.5 ) reverse = true; // Exchange roles of i and i+1.
  int ivert = int(r.random()*double(nverts) + 1.0); // Equal probabilities for ivert = 1, ... , nverts.
  DESreal4vector Pplus;
  DESreal4vector Pminus;
  if(!reverse)
  {
    Pplus  = Q[indx(ivert+1)] - Q[ivert];
    Pminus = Q[indx(ivert-1)] - Q[ivert];  
  }
  else // if(reverse)
  {
    Pplus = Q[ivert] - Q[indx(ivert+1)];
    Pminus = Q[indx(ivert+1)] - Q[indx(ivert+2)];
  }
  if(Pplus(0)*Pminus(0) < 0.0)  Pminus = - Pminus;   // So dot[Pplus,Pminus] > 0.
  double Misq = dot(Pplus,Pminus);  // This is positive with our definitions.
  double Psumabs = sqrt(2.0*Misq);  //This is sqrt( square(Pplus + Pminus) ).
  // Choose x.
  double x = pow(r.random(), 1.0/Acoll_b);
  // Choose sign of y.
  double sy = 1.0;
  if(r.random() < 0.5) sy = -1.0;
  // Choose phi.
  double phi;
  phi = r.random()*twoPI;
  // Choose random variables rz and rdmod = rd^{1/B}.
  double rz = r.random();
  double rdmod = pow(r.random(), 1.0/Bcoll_b);
  // Choose ellTsq and absy.
  double ellTsq = Misq*rz*rdmod;
  double absy = (1.0 - rz)*rdmod/x;
  double ellTabs = sqrt(ellTsq);
  double y = sy*absy;
  // Construct ellT in frame where Psum is in 0 direction and Pdiff is in the +/- 3 direction.
  // Subscripts 2 refer to this frame.
  // Subscripts 1 refer to frame with Psum is in 0 direction but Pdiff in another direction.
  // Then no subscripts refer to original frame.  
  DESreal4vector Psum  = Pplus + Pminus;
  DESreal4vector Pdiff = Pplus - Pminus;
  DESreal4vector Evec, Zvec; // This initializes them to zero.
  if(Psum(0) > 0) Evec.in(0) = Psumabs;     // So square(Evec) = square(Psum);
  else Evec.in(0) = - Psumabs;              // So Evec(0) has same sign as Psum(0);
  DESreal4vector Pdiff_1 = boost(Evec, Psum, Pdiff);   // Pdiff in frame 1.
  if( Pdiff_1(3) > 0 )  Zvec.in(3) = Psumabs;
  else  Zvec.in(3) = - Psumabs;
  DESreal4vector ellT_2;  // This initializes it to zero.
  ellT_2.in(1) = ellTabs*sin(phi);
  ellT_2.in(2) = ellTabs*cos(phi);
  DESreal4vector ellT_1 = boost(Pdiff_1, Zvec, ellT_2); // ellT in frame 1.
  DESreal4vector ellT = boost(Psum, Evec, ellT_1);
  if(!reverse)
  {
   ell = Q[ivert] + x*Pplus + y*Pminus + ellT; 
  }
  else // if(reverse)
  {
    ell = Q[indx(ivert+1)] + x*Pplus + y*Pminus + ellT;
  }
  ellOUT = ell;
  method = 11;
  return;  
} // End method 11.
else if (r.random() < alpha[12]/(1.0 - beta[12] + small))
{
  //
  // Method 12: distributed near the Q[i] - Q[i+1] collinear singularity
  // for the initial state momenta, i = nverts and i = indexA.
  //
  bool reverse = false;
  if( r.random() > 0.5 ) reverse = true; // Exchange roles of i and i+1.
  int ivert;
  if ( r.random() > 0.5 ) ivert = nverts;
  else ivert = indexA;
  DESreal4vector Pplus;
  DESreal4vector Pminus;
  if(!reverse)
  {
    Pplus  = Q[indx(ivert+1)] - Q[ivert];
    Pminus = Q[indx(ivert-1)] - Q[ivert];  
  }
  else // if(reverse)
  {
    Pplus = Q[ivert] - Q[indx(ivert+1)];
    Pminus = Q[indx(ivert+1)] - Q[indx(ivert+2)];
  }
  if(Pplus(0)*Pminus(0) < 0.0)  Pminus = - Pminus;   // So dot[Pplus,Pminus] > 0.
  double Misq = dot(Pplus,Pminus);  // This is positive with our definitions.
  double Psumabs = sqrt(2.0*Misq);  //This is sqrt( square(Pplus + Pminus) ).
                                    // Choose x.
  double x = pow(r.random(), 1.0/AIcoll_b);
  // Choose sign of y.
  double sy = 1.0;
  if(r.random() < 0.5) sy = -1.0;
  // Choose phi.
  double phi;
  phi = r.random()*twoPI;
  // Choose random variables rz and rdmod = rd^{1/B}.
  double rz = r.random();
  double rdmod = pow(r.random(), 1.0/BIcoll_b);
  // Choose ellTsq and absy.
  double ellTsq = Misq*rz*rdmod;
  double absy = (1.0 - rz)*rdmod/x;
  double ellTabs = sqrt(ellTsq);
  double y = sy*absy;
  // Construct ellT in frame where Psum is in 0 direction and Pdiff is in the +/- 3 direction.
  // Subscripts 2 refer to this frame.
  // Subscripts 1 refer to frame with Psum is in 0 direction but Pdiff in another direction.
  // Then no subscripts refer to original frame.  
  DESreal4vector Psum  = Pplus + Pminus;
  DESreal4vector Pdiff = Pplus - Pminus;
  DESreal4vector Evec, Zvec; // This initializes them to zero.
  if(Psum(0) > 0) Evec.in(0) = Psumabs;     // So square(Evec) = square(Psum);
  else Evec.in(0) = - Psumabs;              // So Evec(0) has same sign as Psum(0);
  DESreal4vector Pdiff_1 = boost(Evec, Psum, Pdiff);   // Pdiff in frame 1.
  if( Pdiff_1(3) > 0 )  Zvec.in(3) = Psumabs;
  else  Zvec.in(3) = - Psumabs;
  DESreal4vector ellT_2;  // This initializes it to zero.
  ellT_2.in(1) = ellTabs*sin(phi);
  ellT_2.in(2) = ellTabs*cos(phi);
  DESreal4vector ellT_1 = boost(Pdiff_1, Zvec, ellT_2); // ellT in frame 1.
  DESreal4vector ellT = boost(Psum, Evec, ellT_1);
  if(!reverse)
  {
    ell = Q[ivert] + x*Pplus + y*Pminus + ellT; 
  }
  else // if(reverse)
  {
    ell = Q[indx(ivert+1)] + x*Pplus + y*Pminus + ellT;
  }
  ellOUT = ell;
  method = 12;
  return;  
} // End method 12.
else if (r.random() < alpha[13]/(1.0 - beta[13] + small))
{
  //
  // Method 13: distributed near forward lightcone from Q[i]
  // and backward lightcone from Q[j] where Q[j] - Q[i] is timelike
  // with positive energy.
  //
  // Choose which case, (i,j), at random.
  int i = 0, j = 0;
  int numberofLcases = 0;
  if (indexA > 2) numberofLcases = (indexA - 1)*(indexA - 2)/2;
  int numberofRcases = 0;
  if (nverts > indexA + 2) numberofRcases = (nverts - indexA - 1)*(nverts - indexA - 2)/2;
  int numberofcases = numberofLcases + numberofRcases;
  int casenumber = 0;
  double randomN = numberofcases*r.random();
  if( randomN <= numberofLcases )
  {
    casenumber = 0;
    for(int itrial = 1; itrial <= indexA - 2; itrial++)
    {
      for(int jtrial = itrial + 2; jtrial <= indexA; jtrial++)
      {
        casenumber += 1;
        if (casenumber > randomN)  // We have found which case.
        {
          i = itrial;
          j = jtrial;
          goto ijfound13;
        }
      }
    }
  }
  else // if( randomN > numberofLcases )
  {
    casenumber = numberofLcases;
    for(int itrial = indexA + 1; itrial <= nverts - 2; itrial++)
    {
      for(int jtrial = itrial + 2; jtrial <= nverts; jtrial++)
      {
        casenumber += 1;
        if (casenumber > randomN)  // We have found which case.
        {
          i = itrial;
          j = jtrial;
          goto ijfound13;
        }
      }
    }
  }
  cout << "(i,j) not found in newell method 13" << endl;
  exit(1);
  ijfound13: ; // null statement: start here when (i,j) is found.
  DESreal4vector K = Q[j] - Q[i];
  double Ksq = square(K);
  if(Ksq <= 0.0)
  {
    cout << "This is horrible, Ksq <= 0 in newell method 13." << endl;
    exit(1);
  }
  double absK = sqrt(Ksq);
  DESreal4vector Kprime(absK,0.0,0.0,0.0);
  DESreal4vector lprime;
  double costheta = 2.0*(r.random() - 0.5);
  double sintheta = sqrt(1.0 - costheta*costheta);
  double phi = twoPI*r.random();
  double lsq = 0.25*Ksq*pow(r.random(),1.0/Atlike_b);
  if( r.random() < 0.5 ) lsq = - lsq;
  double lmKsq = 0.25*Ksq*pow(r.random(),1.0/Atlike_b);
  if( r.random() < 0.5 ) lmKsq = - lmKsq;
  lprime.in(0) = (lsq - lmKsq + Ksq)/2.0/absK;
  double absveclprime = sqrt( pow(lprime(0),2) - lsq );
  lprime.in(1) = absveclprime*sintheta*cos(phi);
  lprime.in(2) = absveclprime*sintheta*sin(phi);
  lprime.in(3) = absveclprime*costheta;
  // We have the result in the special frame.
  // Now we boost it back to the c.m. frame
  // and add Q[i] to get ell.
  ell = boost(K,Kprime,lprime) + Q[i];
  ellOUT = ell;
  method = 13;
  return;  
} // End method 13.
else if (r.random() < alpha[14]/(1.0 - beta[14] + small))
{
  //
  // Method 14: distributed near lightcone from Q[i]
  // and lightcone from Q[j] where Q[j] - Q[i] is spacelike.
  //
  // Choose which case, (i,j), at random.
  int i = 0, j = 0;
  int numberofcases = indexA*(nverts - indexA) - 2;
  int casenumber = 0;
  double randomN = numberofcases*r.random();
  for(int itrial = 1; itrial <= indexA; itrial++)
  {
    for(int jtrial = indexA + 1; jtrial <= nverts; jtrial++)
    {
      if( !( (   (itrial == 1)    && (jtrial == nverts)    )   ||
             ( (itrial == indexA) && (jtrial == indexA + 1)) ) )
      {
        casenumber +=1;
        if (casenumber > randomN)  // We have found which case.
        {
          i = itrial;
          j = jtrial;
          goto ijfound14;
        }
      }
    }
  }
  cout << "(i,j) not found in newell method 14" << endl;
  exit(1);
  ijfound14: ;
  DESreal4vector K = Q[j] - Q[i];
  double Ksq = square(K);
  if(Ksq >= 0.0)
  {
    cout << "This is horrible, Ksq >= 0 in newell method 14." << endl;
    exit(1);
  }
  double absK = sqrt(-Ksq);  // Note the sign.
  DESreal4vector Kprime(0.0,0.0,0.0,absK);
  DESreal4vector lprime;
  double phi = twoPI*r.random();
  double coshomega = tan( PIover4*(1.0 + r.random()) );
  double sinhomega = sqrt(coshomega*coshomega - 1.0);
  double lsq = 0.25*Ksq*pow(r.random(),1.0/Aslike_b);
  if( r.random() < 0.5 ) lsq = - lsq;
  double lmKsq = 0.25*Ksq*pow(r.random(),1.0/Aslike_b);
  if( r.random() < 0.5 ) lmKsq = - lmKsq;
  lprime.in(3) = (lmKsq - lsq - Ksq)/2.0/absK;
  double tildelsq = lsq + pow(lprime(3),2);
  if(tildelsq < 0.0)
  {
    cout << "tildelsq negative in newell method 14." << endl;
    exit(1);
  }
  double abstildel = sqrt(tildelsq);
  double Eprime = abstildel*coshomega;
  if(r.random() < 0.5) Eprime = - Eprime;
  lprime.in(0) = Eprime;
  lprime.in(1) = abstildel*sinhomega*cos(phi);
  lprime.in(2) = abstildel*sinhomega*sin(phi);
  // We have the result in the special frame.
  // Now we boost it back to the c.m. frame
  // and add Q[i] to form ell.
  ell = boost(K,Kprime,lprime) + Q[i];
  ellOUT = ell;
  method = 14;
  return;  
} // End method 14.
else if (r.random() < alpha[15]/(1.0 - beta[15] + small))
{
  //
  // Method 15: distributed near lightcone from Q[i]
  // and lightcone from Q[j] where Q[j] - Q[i] is lightlike.
  //
  // Choose which case, (i,j), at random.
  // Switch i and j if needed to make Q[j](0) - Q[i](0) > 0.
  int i = 0, j = 0;
  int ii = ceil(r.random()*nverts);
  if( ii == nverts )
  {
    i = 1;
    j = nverts;
  }
  else if( ii == indexA )
  {
    i = indexA + 1;
    j = indexA;
  }
  else
  {
    i = ii;
    j = ii + 1;
  }
  // Define vector K.
  DESreal4vector K = Q[j] - Q[i];
  // Let vecK be the 3-vector part of K.
  DESreal4vector vecK = K;
  vecK.in(0) = 0.0;
  // Define vecK rotated so that vec K is in +/- z-direction.
  DESreal4vector vecKprime(0.0,0.0,0.0,K(0));
  if( K(3) < 0.0 ) vecKprime.in(3) = - K(0);
  // We construct l_2, in the "prime" system in which
  // K^+ > 0 and all other components of K vanish.
  double temp = 2.0*r.random() - 1.0;
  double x = 0.5*pow(abs(temp),1.0/(1.0 - Allike_b));
  if( temp < 0.0 ) x = - x;
  if( r.random() < 0.5 ) x = 1.0 - x; // 1/2 < x < 3/2 choice
  temp =  2.0*r.random() - 1.0;
  double lsq = Msqllike_b/( pow(abs(temp),-1.0/Bllike_b) - 1.0 );
  if( temp < 0.0 ) lsq = - lsq;
  temp =  2.0*r.random() - 1.0;
  double lmKsq = Msqllike_b/( pow(abs(temp),-1.0/Bllike_b) - 1.0 );
  if( temp < 0.0 ) lmKsq = - lmKsq;
  double lperpsq = -(1.0 - x)*lsq - x*lmKsq;
  // We need lperpsq > 0. If it is not, we could choose a new point,
  // but it is faster to just reverse lsq and lmKsq.
  if(lperpsq < 0.0)
  {
    lsq = - lsq;
    lmKsq = - lmKsq;
    lperpsq = - lperpsq;
  }
  double lperpabs = sqrt(lperpsq);
  double phi = twoPI*r.random();
  DESreal4vector lprime;
  temp = (lsq - lmKsq)/4.0/K(0);
  lprime.in(0) = x*K(0) + temp;
  lprime.in(1) = lperpabs*cos(phi);
  lprime.in(2) = lperpabs*sin(phi);
  if( K(3) > 0.0 ) lprime.in(3) = x*K(0) - temp;
  else lprime.in(3) = - x*K(0) + temp;
  // Now we rotate it back to the c.m. frame
  // and add Q[i] to form ell.
  ell = boost(vecK,vecKprime,lprime) + Q[i];
  ellOUT = ell;
  method = 15;
  return;
} // End method 15.
else if (r.random() < alpha[16]/(1.0 - beta[16] + small))
{
  //
  // Method 16: distributed near lightcone from Q[i]
  // and lightcone from Q[j] where Q[j] - Q[i] is lightlike
  // and also near lightcone from Q[k], which intersects
  // the line from Q[i] to Q[j].
  //
  // We are mainly concerned with the double parton scattering
  // singularity. Therefore if 3 <= A <= N - 3 we do the
  // following. We take one of four cases:
  // i = 1,   j = N,   k = A
  // i = N,   j = 1,   k = A+1
  // i = A+1, j = A,   k = N
  // i = A,   j = A+1, k = 1
  // In cases other than 3 <= A <= N - 3, we take one of
  // i = 1,   j = N,   k = 2, 3,..., A
  // i = N,   j = 1,   k = A+1, A+2,..., N-1
  // i = A+1, j = A,   k = A+2, A+3,..., N
  // i = A,   j = A+1, k = 1, 2, ..., A-1
  // That is a total of 2*N - 4 cases.
  // 
  int i,j,k;
  if ( 3 <= indexA && indexA <= nverts - 3)
  {
    double ran = r.random();
    if( ran < 0.25 )
    {
      i = 1;
      j = nverts;
      k = indexA;
    }
    else if( ran < 0.50 )
    {
      i = nverts;
      j = 1;
      k = indexA + 1;
    }
    else if( ran < 0.75 )
    {
      i = indexA + 1;
      j = indexA;
      k = nverts;
    }
    else
    {
      i = indexA;
      j = indexA + 1;
      k = 1;
    }
  }
  else // if not 3 <= A <= N
  {
    // Choose casenumber at random in range from 1 to 2*nverts - 4.
    int casenumber = int( (2*nverts - 4)*r.random() ) + 1;
    if(casenumber <= indexA - 1)
    {
      i = 1;
      j = nverts;
      k = casenumber + 1;
    }
    else if(casenumber <= nverts - 2)
    {
      i = nverts;
      j = 1;
      k = casenumber + 1;
    }
    else if(casenumber <= 2*nverts - indexA - 3)
    {
      i = indexA + 1;
      j = indexA;
      k = casenumber - nverts + indexA + 3;
    }
    else // 2*nverts - indexA - 2 <= casenumber <= 2*nverts - 4
    {
      i = indexA;
      j = indexA + 1;
      k = casenumber - 2*nverts + indexA + 3;
    }
   } // End if ( 3 <= indexA && indexA <= nverts - 3) ... else ...
  // We now have i,j,k.
  ell = ell_ijkA(i,j,k);
  ellOUT = ell;
  method = 16;
  return;
} // End method 16.
else if (r.random() < alpha[17]/(1.0 - beta[17] + small))
{
  //
  // Method 17: same as method 16 but with different parameters,
  // distributed near lightcone from Q[i]
  // and lightcone from Q[j] where Q[j] - Q[i] is lightlike
  // and also near lightcone from Q[k], which intersects
  // the line from Q[i] to Q[j].
  //
  // We are mainly concerned with the double parton scattering
  // singularity. Therefore if 3 <= A <= N - 3 we do the
  // following. We take one of four cases:
  // i = 1,   j = N,   k = A
  // i = N,   j = 1,   k = A+1
  // i = A+1, j = A,   k = N
  // i = A,   j = A+1, k = 1
  // In cases other than 3 <= A <= N - 3, we take one of
  // i = 1,   j = N,   k = 2, 3,..., A
  // i = N,   j = 1,   k = A+1, A+2,..., N-1
  // i = A+1, j = A,   k = A+2, A+3,..., N
  // i = A,   j = A+1, k = 1, 2, ..., A-1
  // That is a total of 2*N - 4 cases.
  // 
  int i,j,k;
  if ( 3 <= indexA && indexA <= nverts - 3)
  {
    double ran = r.random();
    if( ran < 0.25 )
    {
      i = 1;
      j = nverts;
      k = indexA;
    }
    else if( ran < 0.50 )
    {
      i = nverts;
      j = 1;
      k = indexA + 1;
    }
    else if( ran < 0.75 )
    {
      i = indexA + 1;
      j = indexA;
      k = nverts;
    }
    else
    {
      i = indexA;
      j = indexA + 1;
      k = 1;
    }
  }
  else // if not 3 <= A <= N
  {
    // Choose casenumber at random in range from 1 to 2*nverts - 4.
    int casenumber = int( (2*nverts - 4)*r.random() ) + 1;
    if(casenumber <= indexA - 1)
    {
      i = 1;
      j = nverts;
      k = casenumber + 1;
    }
    else if(casenumber <= nverts - 2)
    {
      i = nverts;
      j = 1;
      k = casenumber + 1;
    }
    else if(casenumber <= 2*nverts - indexA - 3)
    {
      i = indexA + 1;
      j = indexA;
      k = casenumber - nverts + indexA + 3;
    }
    else // 2*nverts - indexA - 2 <= casenumber <= 2*nverts - 4
    {
      i = indexA;
      j = indexA + 1;
      k = casenumber - 2*nverts + indexA + 3;
    }
  } // End if ( 3 <= indexA && indexA <= nverts - 3) ... else ...
  // We now have i,j,k.
  ell = ell_ijkB(i,j,k);
  ellOUT = ell;
  method = 17;
  return;
} // End method 17.
else if (r.random() < alpha[18]/(1.0 - beta[18] + small))
{
  //
  // Method 18: distributed near lightcone from Q[i]
  // and lightcone from Q[j] where j = i + 1 
  // so Q[j] - Q[i] is lightlike
  // and also near lightcone from Q[k],
  // k = j + 1 or else k = i - 1.
  // That is a total of 2*N cases.
  // 
  int i,j,k;
    // Choose casenumber at random in range from 1 to 2*nverts.
    int casenumber = int( 2*nverts*r.random() ) + 1;
    if(casenumber <= nverts - 2)
    {
      i = casenumber;
      j = i + 1;
      k = i + 2;
    }
    else if(casenumber == nverts - 1)
    {
      i = nverts - 1;
      j = nverts;
      k = 1;
    }
    else if(casenumber == nverts)
    {
      i = nverts;
      j = 1;
      k = 2;
    }
    else if(casenumber == nverts + 1)
    {
      i = 1;
      j = 2;
      k = nverts;
    }
    else if(casenumber <= 2*nverts - 1)
    {
      i = casenumber - nverts;
      j = i + 1;
      k = i - 1 ;
    }
    else // casenumber = 2*nverts
    {
      i = nverts;
      j = 1;
      k = nverts - 1;
    }
  // We now have i,j,k.
  ell = ell_ijkA(i,j,k);
  ellOUT = ell;
  method = 18;
  return;  
} // End method 18.
else if (r.random() < alpha[19]/(1.0 - beta[19] + small))
{
  //
  // Method 19: distributed near lightcone from Q[i]
  // and lightcone from Q[j] where j = i + 1 
  // so Q[j] - Q[i] is lightlike
  // and also near lightcone from Q[k],
  // k = j + 1 or else k = i - 1.
  // That is a total of 2*N cases. 
  // Same as method 18, but with different parameters. 
  // 
  int i,j,k;
  // Choose casenumber at random in range from 1 to 2*nverts.
  int casenumber = int( 2*nverts*r.random() ) + 1;
  if(casenumber <= nverts - 2)
  {
    i = casenumber;
    j = i + 1;
    k = i + 2;
  }
  else if(casenumber == nverts - 1)
  {
    i = nverts - 1;
    j = nverts;
    k = 1;
  }
  else if(casenumber == nverts)
  {
    i = nverts;
    j = 1;
    k = 2;
  }
  else if(casenumber == nverts + 1)
  {
    i = 1;
    j = 2;
    k = nverts;
  }
  else if(casenumber <= 2*nverts - 1)
  {
    i = casenumber - nverts;
    j = i + 1;
    k = i - 1 ;
  }
  else // casenumber = 2*nverts
  {
    i = nverts;
    j = 1;
    k = nverts - 1;
  }
  // We now have i,j,k.
  ell = ell_ijkB(i,j,k);
  ellOUT = ell;
  method = 19;
  return;    
} // End method 19.
else if (r.random() < alpha[20]/(1.0 - beta[20] + small))
{
  //
  // Method 20: distributed near lightcone from Q[i]
  // and lightcone from Q[j] where Q[j] - Q[i] is lightlike
  // and also near lightcone from Q[k], which intersects
  // the line from Q[i] to Q[j].
  //
  // We take one of (for the generic case 1 < A < nverts - 1)
  // i = 1,   j = N,   k = 3,..., A          (empty for A = 2)
  // i = N,   j = 1,   k = A+1, A+2,..., N-2 (empty for A = N-2)
  // i = A+1, j = A,   k = A+3,..., N        (empty for A = N-2)
  // i = A,   j = A+1, k = 1, 2, ..., A-2    (empty for A = 2)
  // That is a total of 2*N - 8 cases.
  // Special cases (always with a total of 2*N - 8 cases):
  // A = 1:
  // i = 1,   j = N,   k = none
  // i = N,   j = 1,   k = 3,..., N-2
  // i = A+1, j = A,   k = 4,..., N-1
  // i = A,   j = A+1, k = none
  // A = N-1:
  // i = 1,   j = N,   k = 3,..., N-2
  // i = N,   j = 1,   k = none
  // i = A+1, j = A,   k = none
  // i = A,   j = A+1, k = 2, ..., N-3
  // 
  int i,j,k;
  // Choose casenumber at random in range from 1 to 2*nverts - 8.
  int casenumber = int( (2*nverts - 8)*r.random() ) + 1;
  if((casenumber <= indexA - 2)&&(casenumber <= nverts - 4))
  {
    i = 1;
    j = nverts;
    k = casenumber + 2;
  }
  else if(casenumber <= nverts - 4)
  {
    i = nverts;
    j = 1;
    k = casenumber + 2;
  }
  else if((casenumber <= 2*nverts - indexA - 6)&&(casenumber <= 2*nverts - 8))
  {
    i = indexA + 1;
    j = indexA;
    k = casenumber - nverts + indexA + 6;
  }
  else //  casenumber <= 2*nverts - 8
  {
    i = indexA;
    j = indexA + 1;
    k = casenumber - 2*nverts + indexA + 6;
  }
  // We now have i,j,k.
  ell = ell_ijkA(i,j,k);
  ellOUT = ell;
  method = 20;
  return;
} // End method 20.
else //if (r.random() < alpha[21]/(1.0 - beta[21] + small))
{
  //
  // Method 21: distributed near lightcone from Q[i]
  // and lightcone from Q[j] where Q[j] - Q[i] is lightlike
  // and also near lightcone from Q[k], which intersects
  // the line from Q[i] to Q[j]. Same as Method 20 except
  // that it uses different parameters, as contained in 
  // ell_ijkB(i,j,k) instead of ell_ijkA(i,j,k).
  //
  int i,j,k;
  // Choose casenumber at random in range from 1 to 2*nverts - 8.
  int casenumber = int( (2*nverts - 8)*r.random() ) + 1;
  if((casenumber <= indexA - 2)&&(casenumber <= nverts - 4))
  {
    i = 1;
    j = nverts;
    k = casenumber + 2;
  }
  else if(casenumber <= nverts - 4)
  {
    i = nverts;
    j = 1;
    k = casenumber + 2;
  }
  else if((casenumber <= 2*nverts - indexA - 6)&&(casenumber <= 2*nverts - 8))
  {
    i = indexA + 1;
    j = indexA;
    k = casenumber - nverts + indexA + 6;
  }
  else //  casenumber <= 2*nverts - 8
  {
    i = indexA;
    j = indexA + 1;
    k = casenumber - 2*nverts + indexA + 6;
  }
  // We now have i,j,k.
  ell = ell_ijkB(i,j,k);
  ellOUT = ell;
  method = 21;
  return;
} // End method 21.
  
} // End of newell(ellOUT)
//--------------------------------------------------------------------------------------------------
// Return rho(ell) for a point ell.
double DESloopmomenta::rhoell()
{
static const double invtwoPIsq = 0.05066059182116889;  // 1/(2*Pi^2)
static const double invfourPI  = 0.07957747154594767;  // 1/(4*Pi)
static const double invtwoPI   = 0.1591549430918953;   // 1/(2*Pi)
static const double invPIsq    = 0.1013211836423378;   // 1/Pi^2
double MT4 = pow(MTsq,2);
double M04 = pow(dotPPbar,2);
double ellEsq = ell(0)*ell(0) + ell(1)*ell(1) + ell(2)*ell(2) + ell(3)*ell(3);
double ellE4 = pow(ellEsq,2);
double rhoval = 0.0;
double rhotemp = 0.0;
// The various methods follow.
if( alpha[1] > tiny )
{
// Method 1: roughly uniform over a wide range.
rhotemp = alpha[1]*invtwoPIsq*4.0;
rhotemp = rhotemp * (Awide - 1.0) /M04 /pow((1.0 + ellE4/M04),Awide);
// cout << "Method 1, rhotemp = " << rhotemp << endl;
//
rhoval = rhoval + rhotemp;
}
// End of method 1.
//
if( alpha[2] > tiny )
{
// Method 2: distributed near double parton scattering singularity.
rhotemp = alpha[2]*invtwoPIsq*4.0;
rhotemp = rhotemp * (mlarge4_a - msmall4_a) /log(mlarge4_a/msmall4_a) 
          /(mlarge4_a + ellE4) /(msmall4_a + ellE4);
// cout << "Method 2, rhotemp = " << rhotemp << endl;
//
rhoval = rhoval + rhotemp;
}
// End of method 2.
//
if( alpha[3] > tiny )
{
// Method 3: distributed near the soft singularities.
DESreal4vector Pplus, Pminus, ellrel, ellT;
double x, y, dotPplusPminus, Misq, absxy, ellTsq;
for(int ivert = 1; ivert <= nverts; ivert++)
{
  rhotemp = 0.0;  // Initial value. It will survive outside range of integral for this method.
  Pplus  = Q[indx(ivert+1)] - Q[ivert];
  Pminus = Q[indx(ivert-1)] - Q[ivert];
  ellrel = ell - Q[ivert];
  dotPplusPminus = dot(Pplus,Pminus);
  Misq = abs(dotPplusPminus);
  x = dot(ellrel,Pminus)/dotPplusPminus;
  if(abs(x) < 1.0)
  {
    y = dot(ellrel,Pplus) /dotPplusPminus;
    if( abs(y) < 1.0)
    {
      ellT = ellrel - x*Pplus - y*Pminus;
      ellTsq = - square(ellT); // This is positive.
      if( ellTsq < Misq )
      {
        absxy = abs(x*y);
        if(absxy > tiny)  // This is to prevent rhotemp = nan.
        {
          rhotemp = alpha[3]/double(nverts);
          rhotemp *= Asoft_a*Asoft_a*Bsoft_a*invfourPI/Misq/Misq;
          rhotemp /= pow(1.0 + absxy, Bsoft_a) -  pow(absxy, Bsoft_a);
          rhotemp /= pow(absxy, 1.0 - Asoft_a);
          rhotemp *= pow(Misq/(absxy*Misq + ellTsq), 1.0 - Bsoft_a);
        }
        else rhotemp = huge;
      } // Close if( ellTsq < Misq ).
    }   // Close if( abs(y) < 1.0).
  }     // Close if(abs(x) < 1.0).
  //
  rhoval = rhoval + rhotemp;
} // Close for(int ivert = 1; ivert <= nverts; ivert++). 
}
// End of method 3.
//
if( alpha[4] > tiny )
{
// Method 4: distributed near the collinear singularities.
DESreal4vector Pplus, Pminus, ellrel, ellT;
double x, y, dotPplusPminus, Misq, absxy, ellTsq;
for(int ivert = 1; ivert <= nverts; ivert++)
{
  for(int ireverse = 1; ireverse <= 2; ireverse++)  // Whether to go in reverse sense.
    {
    bool reverse = false;
    if(ireverse == 2) reverse = true;
    rhotemp = 0.0;  // Initial value. It will survive outside range of integral for this method.
    if(!reverse)
    {
      Pplus  = Q[indx(ivert+1)] - Q[ivert];
      Pminus = Q[indx(ivert-1)] - Q[ivert];  
    }
    else // if(reverse)
    {
      Pplus = Q[ivert] - Q[indx(ivert+1)];
      Pminus = Q[indx(ivert+1)] - Q[indx(ivert+2)];
    }
    if(Pplus(0)*Pminus(0) < 0.0)  Pminus = - Pminus;   // So dot[Pplus,Pminus] > 0.
    dotPplusPminus = dot(Pplus,Pminus);   // This is positive with our definitions.
    Misq = dotPplusPminus;
    if(!reverse) ellrel = ell - Q[ivert];
    else ellrel = ell - Q[indx(ivert+1)];
    x = dot(ellrel,Pminus)/dotPplusPminus;
    if(x > 0.0 && x < 1.0)
    {
      y = dot(ellrel,Pplus) /dotPplusPminus;
      double absy = abs(y);
      ellT = ellrel - x*Pplus - y*Pminus;
      ellTsq = - square(ellT); // This is positive.
      double denom = x*absy + ellTsq/Misq;
      if(denom < 1.0)
      {
        rhotemp = 0.5*alpha[4]/double(nverts);
        rhotemp *= Acoll_a*Bcoll_a*invtwoPI/Misq/Misq;
        rhotemp *= pow(x,Acoll_a);
        rhotemp /= pow(denom,2.0 - Bcoll_a);
      } // Close if(denom < 1.0).
    }   // Close if(x > 0.0 && x < 1.0).
    //
    rhoval = rhoval + rhotemp;
  }  // Close for(int ireverse = 1; ireverse <= 2; ireverse++).
} // Close for(int ivert = 1; ivert <= nverts; ivert++).
}
// End of method 4.
//
if( alpha[5] > tiny )
{
// Method 5: distributed near the Q[i] - Q[i+1] collinear singularity
// for the initial state momenta, i = nverts and i = indexA.
DESreal4vector Pplus, Pminus, ellrel, ellT;
double x, y, dotPplusPminus, Misq, absxy, ellTsq;
for(int itemp = 1; itemp <= 2; itemp++)
{
  int ivert;
  if(itemp == 1) ivert = nverts;
  else ivert = indexA;
  for(int ireverse = 1; ireverse <= 2; ireverse++)  // Whether to go in reverse sense.
  {
    bool reverse = false;
    if(ireverse == 2) reverse = true;
    rhotemp = 0.0;  // Initial value. It will survive outside range of integral for this method.
    if(!reverse)
    {
      Pplus  = Q[indx(ivert+1)] - Q[ivert];
      Pminus = Q[indx(ivert-1)] - Q[ivert];  
    }
    else // if(reverse)
    {
      Pplus = Q[ivert] - Q[indx(ivert+1)];
      Pminus = Q[indx(ivert+1)] - Q[indx(ivert+2)];
    }
    if(Pplus(0)*Pminus(0) < 0.0)  Pminus = - Pminus;   // So dot[Pplus,Pminus] > 0.
    dotPplusPminus = dot(Pplus,Pminus);   // This is positive with our definitions.
    Misq = dotPplusPminus;
    if(!reverse) ellrel = ell - Q[ivert];
    else ellrel = ell - Q[indx(ivert+1)];
    x = dot(ellrel,Pminus)/dotPplusPminus;
    if(x > 0.0 && x < 1.0)
    {
      y = dot(ellrel,Pplus) /dotPplusPminus;
      double absy = abs(y);
      ellT = ellrel - x*Pplus - y*Pminus;
      ellTsq = - square(ellT); // This is positive.
      double denom = x*absy + ellTsq/Misq;
      if(denom < 1.0)
      {
        rhotemp = 0.25*alpha[5];
        rhotemp *= AIcoll_a*BIcoll_a*invtwoPI/Misq/Misq;
        rhotemp *= pow(x,AIcoll_a);
        rhotemp /= pow(denom,2.0 - BIcoll_a);
      } // Close if(denom < 1.0).
    }   // Close if(x > 0.0 && x < 1.0).
          //
    rhoval = rhoval + rhotemp;
  }  // Close for(int ireverse = 1; ireverse <= 2; ireverse++).
} // Close for(int ivert = 1; ivert <= nverts; ivert++).
}
// End of method 5.
//
if( alpha[6] > tiny )
{
// Method 6: distributed near forward lightcone from Q[i]
// and backward lightcone from Q[j] where Q[j] - Q[i] is timelike
// with positive energy.  
int numberofcases = 0;
if (indexA > 2) numberofcases += (indexA - 1)*(indexA - 2)/2;
if (nverts > indexA + 2) 
   numberofcases += (nverts - indexA - 1)*(nverts - indexA - 2)/2;
for(int i = 1; i <= indexA - 2; i++)
{
  for(int j = i + 2; j <= indexA; j++)
  {
    DESreal4vector K = Q[j] - Q[i];
    double Ksq = square(K);
    if(Ksq <= 0.0)
    {
      cout << "This is horrible, Ksq <= 0 in rhoell method 6." << endl;
      exit(1);
    }
    double absK = sqrt(Ksq);
    DESreal4vector l = ell - Q[i];
    DESreal4vector lmK = ell - Q[j];
    double lsq = square(l);
    double lmKsq = square(lmK);
    if( (abs(lsq) < 0.25*Ksq) && (abs(lmKsq) < 0.25*Ksq) )
    {
      if( abs(lsq*lmKsq) > tiny) // This is to prevent rhotemp = nan.
      {
        double abslvec = sqrt( pow(dot(l,K),2)/Ksq - lsq );
        rhotemp = alpha[6]/numberofcases;
        rhotemp *= invfourPI*absK/abslvec;
        rhotemp *= Atlike_a*Atlike_a/abs(lsq*lmKsq);
        double temp = 16.0*abs(lsq*lmKsq)/Ksq/Ksq;
        rhotemp *= pow(temp,Atlike_a);
      }
      else rhotemp = huge;
      rhoval = rhoval + rhotemp;
    }
  }
}
for(int i = indexA + 1; i <= nverts - 2; i++)
{
  for(int j = i + 2; j <= nverts; j++)
  {
    DESreal4vector K = Q[j] - Q[i];
    double Ksq = square(K);
    if(Ksq < 0.0)
    {
      cout << "This is horrible, Ksq < 0 in rhoell method 6." << endl;
      exit(1);
    }
    double absK = sqrt(Ksq);
    DESreal4vector l = ell - Q[i];
    DESreal4vector lmK = ell - Q[j];
    double lsq = square(l);
    double lmKsq = square(lmK);
    if( (abs(lsq) < 0.25*Ksq) && (abs(lmKsq) < 0.25*Ksq) )
    {
      if( abs(lsq*lmKsq) > tiny) // This is to prevent rhotemp = nan.
      {
        double abslvec = sqrt( pow(dot(l,K),2)/Ksq - lsq );
        rhotemp = alpha[6]/numberofcases;
        rhotemp *= invfourPI*absK/abslvec;
        rhotemp *= Atlike_a*Atlike_a/abs(lsq*lmKsq);
        double temp = 16.0*abs(lsq*lmKsq)/Ksq/Ksq;
        rhotemp *= pow(temp,Atlike_a);
        rhoval = rhoval + rhotemp;
      }
      else rhotemp = huge;
    }
  }
}
}
// End of method 6.
//
if( alpha[7] > tiny )
{
// Method 7: distributed near lightcone from Q[i]
// and lightcone from Q[j] where Q[j] - Q[i] is spacelike.
int numberofcases = indexA*(nverts - indexA) - 2;
for(int i = 1; i <= indexA; i++)
{
  for(int j = indexA + 1; j <= nverts; j++)
  {
    if( !( (   (i == 1)    && (j == nverts)    )   ||
           ( (i == indexA) && (j == indexA + 1)) ) )
    {
      DESreal4vector K = Q[j] - Q[i];
      double Ksq = square(K);
      if(Ksq >= 0.0)
      {
        cout << "This is horrible, Ksq >= 0 in rhoell method 7." << endl;
        exit(1);
      }
      double absK = sqrt(-Ksq);       // Note sign.
      DESreal4vector l = ell - Q[i];
      DESreal4vector lmK = ell - Q[j];
      double lsq = square(l);
      double lmKsq = square(lmK);
      if( (abs(lsq) < 0.25*abs(Ksq)) && (abs(lmKsq) < 0.25*abs(Ksq)) )
      {
        double tildelsq = lsq - pow(dot(l,K),2)/Ksq;
        DESreal4vector Kprime(0.0,0.0,0.0,absK);
        DESreal4vector nprime(1.0,0.0,0.0,0.0);
        DESreal4vector n = boost(K,Kprime,nprime);
        rhotemp = alpha[7]/numberofcases;
        rhotemp *= Aslike_a*Aslike_a/abs(lsq*lmKsq)*invPIsq;
        double temp = 16.0*abs(lsq*lmKsq)/Ksq/Ksq;
        if( temp > tiny) // This protects us against getting nan.
        {
          rhotemp *= pow(temp,Aslike_a);
          rhotemp *= sqrt(tildelsq*abs(Ksq))/(tildelsq + pow(dot(l,n),2));
        }
        else rhotemp = huge;
        rhoval = rhoval + rhotemp;
      }
    }
  }
}
}
// End of method 7.
//
if( alpha[8] > tiny )
{
// Method 8: distributed near lightcone from Q[i]
// and lightcone from Q[j] where Q[j] - Q[i] is lightlike.
int numberofcases = nverts;
for(int ii = 1; ii <= nverts; ii++)
{
  int i,j;
  if( ii == nverts)
  {
    i = 1;
    j = nverts;
  }
  else if( ii == indexA )
  {
    i = indexA + 1;
    j = indexA;
  }
  else
  {
    i = ii;
    j = ii + 1;
  }
  // Define vectors K and Kbar.
  DESreal4vector K = Q[j] - Q[i];
  DESreal4vector Kbar = - K;
  Kbar.in(0) = K(0);
  DESreal4vector l = ell - Q[i];
  DESreal4vector lmK = ell - Q[j];
  double lsq = square(l);
  double lmKsq = square(lmK);
  double x = dot(Kbar,l)/2/pow(K(0),2);
  //
  if(abs(x-0.5) < 1.0)
  {
    if( abs(lsq*lmKsq*x*(1.0-x)) > tiny) // This protects us against getting nan.
    {
      rhotemp = alpha[8]/numberofcases;
      rhotemp *= Bllike_a*Bllike_a*(1.0 - Allike_a)*invtwoPI;
      double xmod;
      if(x < 0.5) xmod = abs(x);
      else xmod = abs(1-x);
      rhotemp *= pow(2*xmod,-Allike_a);
      double temp = (1.0 + Msqllike_a/abs(lsq))*(1.0 + Msqllike_a/abs(lmKsq));
      rhotemp *= pow(temp,-Bllike_a-1);
      rhotemp *= Msqllike_a/pow(lsq,2);
      rhotemp *= Msqllike_a/pow(lmKsq,2);
    }
    else rhotemp = huge;
    rhoval = rhoval + rhotemp;
  }
}
}
// End of method 8.
//
if( alpha[9] > tiny )
{
// Method 9: distributed near double parton scattering singularity.
rhotemp = alpha[9]*invtwoPIsq*4.0;
rhotemp = rhotemp * (mlarge4_b - msmall4_b) /log(mlarge4_b/msmall4_b) 
          /(mlarge4_b + ellE4) /(msmall4_b + ellE4);
// cout << "Method 9, rhotemp = " << rhotemp << endl;
//
rhoval = rhoval + rhotemp;
}
// End of method 9.
if( alpha[10] > tiny )
{
// Method 10: distributed near the soft singularities.
// DESreal4vector Pplus, Pminus, ellrel, ellT;
// double x, y, dotPplusPminus, Misq, absxy, ellTsq;
DESreal4vector Pplus, Pminus, ellrel, ellT;
double x, y, dotPplusPminus, Misq, absxy, ellTsq;
for(int ivert = 1; ivert <= nverts; ivert++)
{
  rhotemp = 0.0;  // Initial value. It will survive outside range of integral for this method.
  Pplus  = Q[indx(ivert+1)] - Q[ivert];
  Pminus = Q[indx(ivert-1)] - Q[ivert];
  ellrel = ell - Q[ivert];
  dotPplusPminus = dot(Pplus,Pminus);
  Misq = abs(dotPplusPminus);
  x = dot(ellrel,Pminus)/dotPplusPminus;
  if(abs(x) < 1.0)
  {
    y = dot(ellrel,Pplus) /dotPplusPminus;
    if( abs(y) < 1.0)
    {
      ellT = ellrel - x*Pplus - y*Pminus;
      ellTsq = - square(ellT); // This is positive.
      if( ellTsq < Misq )
      {
        absxy = abs(x*y);
        if(absxy > tiny)  // This is to prevent rhotemp = nan.
        {
          rhotemp = alpha[10]/double(nverts);
          rhotemp *= Asoft_b*Asoft_b*Bsoft_b*invfourPI/Misq/Misq;
          rhotemp /= pow(1.0 + absxy, Bsoft_b) -  pow(absxy, Bsoft_b);
          rhotemp /= pow(absxy, 1.0 - Asoft_b);
          rhotemp *= pow(Misq/(absxy*Misq + ellTsq), 1.0 - Bsoft_b);
        }
        else rhotemp = huge;
      } // Close if( ellTsq < Misq ).
    }   // Close if( abs(y) < 1.0).
  }     // Close if(abs(x) < 1.0).
  //
  rhoval = rhoval + rhotemp;
} // Close for(int ivert = 1; ivert <= nverts; ivert++).
}
// End of method 10.
//
if( alpha[11] > tiny )
{
// Method 11: distributed near the collinear singularities.
DESreal4vector Pplus, Pminus, ellrel, ellT;
double x, y, dotPplusPminus, Misq, ellTsq;
for(int ivert = 1; ivert <= nverts; ivert++)
{
  for(int ireverse = 1; ireverse <= 2; ireverse++)  // Whether to go in reverse sense.
    {
    bool reverse = false;
    if(ireverse == 2) reverse = true;
    rhotemp = 0.0;  // Initial value. It will survive outside range of integral for this method.
    if(!reverse)
    {
      Pplus  = Q[indx(ivert+1)] - Q[ivert];
      Pminus = Q[indx(ivert-1)] - Q[ivert];  
    }
    else // if(reverse)
    {
      Pplus = Q[ivert] - Q[indx(ivert+1)];
      Pminus = Q[indx(ivert+1)] - Q[indx(ivert+2)];
    }
    if(Pplus(0)*Pminus(0) < 0.0)  Pminus = - Pminus;   // So dot[Pplus,Pminus] > 0.
    dotPplusPminus = dot(Pplus,Pminus);   // This is positive with our definitions.
    Misq = dotPplusPminus;
    if(!reverse) ellrel = ell - Q[ivert];
    else ellrel = ell - Q[indx(ivert+1)];
    x = dot(ellrel,Pminus)/dotPplusPminus;
    if(x > 0.0 && x < 1.0)
    {
      y = dot(ellrel,Pplus) /dotPplusPminus;
      double absy = abs(y);
      ellT = ellrel - x*Pplus - y*Pminus;
      ellTsq = - square(ellT); // This is positive.
      double denom = x*absy + ellTsq/Misq;
      if(denom < 1.0)
      {
        rhotemp = 0.5*alpha[11]/double(nverts);
        rhotemp *= Acoll_b*Bcoll_b*invtwoPI/Misq/Misq;
        rhotemp *= pow(x,Acoll_b);
        rhotemp /= pow(denom,2.0 - Bcoll_b);
      } // Close if(denom < 1.0).
    }   // Close if(x > 0.0 && x < 1.0).
    //
    rhoval = rhoval + rhotemp;
  }  // Close for(int ireverse = 1; ireverse <= 2; ireverse++).
} // Close for(int ivert = 1; ivert <= nverts; ivert++).
}
// End of method 11.
//
if( alpha[12] > tiny )
{
// Method 12: distributed near the Q[i] - Q[i+1] collinear singularity
// for the initial state momenta, i = nverts and i = indexA.
DESreal4vector Pplus, Pminus, ellrel, ellT;
double x, y, dotPplusPminus, Misq, ellTsq;
for(int itemp = 1; itemp <= 2; itemp++)
{
  int ivert;
  if(itemp == 1) ivert = nverts;
  else ivert = indexA;
  for(int ireverse = 1; ireverse <= 2; ireverse++)  // Whether to go in reverse sense.
  {
    bool reverse = false;
    if(ireverse == 2) reverse = true;
    rhotemp = 0.0;  // Initial value. It will survive outside range of integral for this method.
    if(!reverse)
    {
      Pplus  = Q[indx(ivert+1)] - Q[ivert];
      Pminus = Q[indx(ivert-1)] - Q[ivert];  
    }
    else // if(reverse)
    {
      Pplus = Q[ivert] - Q[indx(ivert+1)];
      Pminus = Q[indx(ivert+1)] - Q[indx(ivert+2)];
    }
    if(Pplus(0)*Pminus(0) < 0.0)  Pminus = - Pminus;   // So dot[Pplus,Pminus] > 0.
    dotPplusPminus = dot(Pplus,Pminus);   // This is positive with our definitions.
    Misq = dotPplusPminus;
    if(!reverse) ellrel = ell - Q[ivert];
    else ellrel = ell - Q[indx(ivert+1)];
    x = dot(ellrel,Pminus)/dotPplusPminus;
    if(x > 0.0 && x < 1.0)
    {
      y = dot(ellrel,Pplus) /dotPplusPminus;
      double absy = abs(y);
      ellT = ellrel - x*Pplus - y*Pminus;
      ellTsq = - square(ellT); // This is positive.
      double denom = x*absy + ellTsq/Misq;
      if(denom < 1.0)
      {
        rhotemp = 0.25*alpha[12];
        rhotemp *= AIcoll_b*BIcoll_b*invtwoPI/Misq/Misq;
        rhotemp *= pow(x,AIcoll_b);
        rhotemp /= pow(denom,2.0 - BIcoll_b);
      } // Close if(denom < 1.0).
    }   // Close if(x > 0.0 && x < 1.0).
          //
    rhoval = rhoval + rhotemp;
  }  // Close for(int ireverse = 1; ireverse <= 2; ireverse++).
} // Close for(int ivert = 1; ivert <= nverts; ivert++).
}
// End of method 12.
//
if( alpha[13] > tiny )
{
// Method 13: distributed near forward lightcone from Q[i]
// and backward lightcone from Q[j] where Q[j] - Q[i] is timelike
// with positive energy.
int numberofcases = 0;
if (indexA > 2) numberofcases += (indexA - 1)*(indexA - 2)/2;
if (nverts > indexA + 2) 
   numberofcases += (nverts - indexA - 1)*(nverts - indexA - 2)/2;
for(int i = 1; i <= indexA - 2; i++)
{
  for(int j = i + 2; j <= indexA; j++)
  {
    DESreal4vector K = Q[j] - Q[i];
    double Ksq = square(K);
    if(Ksq <= 0.0)
    {
      cout << "This is horrible, Ksq <= 0 in rhoell method 13." << endl;
      exit(1);
    }
    double absK = sqrt(Ksq);
    DESreal4vector l = ell - Q[i];
    DESreal4vector lmK = ell - Q[j];
    double lsq = square(l);
    double lmKsq = square(lmK);
    if( (abs(lsq) < 0.25*Ksq) && (abs(lmKsq) < 0.25*Ksq) )
    {
      if( abs(lsq*lmKsq) > tiny) // This is to prevent rhotemp = nan.
      {
        double abslvec = sqrt( pow(dot(l,K),2)/Ksq - lsq );
        rhotemp = alpha[13]/numberofcases;
        rhotemp *= invfourPI*absK/abslvec;
        rhotemp *= Atlike_b*Atlike_b/abs(lsq*lmKsq);
        double temp = 16.0*abs(lsq*lmKsq)/Ksq/Ksq;
        rhotemp *= pow(temp,Atlike_b);
      }
      else rhotemp = huge;
      rhoval = rhoval + rhotemp;
    }
  }
}
for(int i = indexA + 1; i <= nverts - 2; i++)
{
  for(int j = i + 2; j <= nverts; j++)
  {
    DESreal4vector K = Q[j] - Q[i];
    double Ksq = square(K);
    if(Ksq < 0.0)
    {
      cout << "This is horrible, Ksq < 0 in rhoell method 13." << endl;
      exit(1);
    }
    double absK = sqrt(Ksq);
    DESreal4vector l = ell - Q[i];
    DESreal4vector lmK = ell - Q[j];
    double lsq = square(l);
    double lmKsq = square(lmK);
    if( (abs(lsq) < 0.25*Ksq) && (abs(lmKsq) < 0.25*Ksq) )
    {
      if( abs(lsq*lmKsq) > tiny) // This is to prevent rhotemp = nan.
      {
        double abslvec = sqrt( pow(dot(l,K),2)/Ksq - lsq );
        rhotemp = alpha[13]/numberofcases;
        rhotemp *= invfourPI*absK/abslvec;
        rhotemp *= Atlike_b*Atlike_b/abs(lsq*lmKsq);
        double temp = 16.0*abs(lsq*lmKsq)/Ksq/Ksq;
        rhotemp *= pow(temp,Atlike_b);
        rhoval = rhoval + rhotemp;
      }
      else rhotemp = huge;
    }
  }
}
}
// End of method 13.
//
if( alpha[14] > tiny )
{
// Method 14: distributed near lightcone from Q[i]
// and lightcone from Q[j] where Q[j] - Q[i] is spacelike.
int numberofcases = indexA*(nverts - indexA) - 2;
for(int i = 1; i <= indexA; i++)
{
  for(int j = indexA + 1; j <= nverts; j++)
  {
    if( !( (   (i == 1)    && (j == nverts)    )   ||
           ( (i == indexA) && (j == indexA + 1)) ) )
    {
      DESreal4vector K = Q[j] - Q[i];
      double Ksq = square(K);
      if(Ksq >= 0.0)
      {
        cout << "This is horrible, Ksq >= 0 in rhoell method 14." << endl;
        exit(1);
      }
      double absK = sqrt(-Ksq);       // Note sign.
      DESreal4vector l = ell - Q[i];
      DESreal4vector lmK = ell - Q[j];
      double lsq = square(l);
      double lmKsq = square(lmK);
      if( (abs(lsq) < 0.25*abs(Ksq)) && (abs(lmKsq) < 0.25*abs(Ksq)) )
      {
        double tildelsq = lsq - pow(dot(l,K),2)/Ksq;
        DESreal4vector Kprime(0.0,0.0,0.0,absK);
        DESreal4vector nprime(1.0,0.0,0.0,0.0);
        DESreal4vector n = boost(K,Kprime,nprime);
        rhotemp = alpha[14]/numberofcases;
        rhotemp *= Aslike_b*Aslike_b/abs(lsq*lmKsq)*invPIsq;
        double temp = 16.0*abs(lsq*lmKsq)/Ksq/Ksq;
        if( temp > tiny) // This protects us against getting nan.
        {
          rhotemp *= pow(temp,Aslike_b);
          rhotemp *= sqrt(tildelsq*abs(Ksq))/(tildelsq + pow(dot(l,n),2));
        }
        else rhotemp = huge;
        rhoval = rhoval + rhotemp;
      }
    }
  }
}
}
// End of method 14.
//
if( alpha[15] > tiny )
{
// Method 15: distributed near lightcone from Q[i]
// and lightcone from Q[j] where Q[j] - Q[i] is lightlike.
int numberofcases = nverts;
for(int ii = 1; ii <= nverts; ii++)
{
  int i,j;
  if( ii == nverts)
  {
    i = 1;
    j = nverts;
  }
  else if( ii == indexA )
  {
    i = indexA + 1;
    j = indexA;
  }
  else
  {
    i = ii;
    j = ii + 1;
  }
  // Define vectors K and Kbar.
  DESreal4vector K = Q[j] - Q[i];
  DESreal4vector Kbar = - K;
  Kbar.in(0) = K(0);
  DESreal4vector l = ell - Q[i];
  DESreal4vector lmK = ell - Q[j];
  double lsq = square(l);
  double lmKsq = square(lmK);
  double x = dot(Kbar,l)/2/pow(K(0),2);
  //
  if(abs(x-0.5) < 1.0)
  {
    if( abs(lsq*lmKsq*x*(1.0-x)) > tiny) // This protects us against getting nan.
    {
      rhotemp = alpha[15]/numberofcases;
      rhotemp *= Bllike_b*Bllike_b*(1.0 - Allike_b)*invtwoPI;
      double xmod;
      if(x < 0.5) xmod = abs(x);
      else xmod = abs(1-x);
      rhotemp *= pow(2*xmod,-Allike_b);
      double temp = (1.0 + Msqllike_b/abs(lsq))*(1.0 + Msqllike_b/abs(lmKsq));
      rhotemp *= pow(temp,-Bllike_b-1);
      rhotemp *= Msqllike_b/pow(lsq,2);
      rhotemp *= Msqllike_b/pow(lmKsq,2);
    }
    else rhotemp = huge;
    rhoval = rhoval + rhotemp;
  }
}
}
// End of method 15.
//
if( alpha[16] > tiny )
{
// Method 16: distributed near lightcone from Q[i]
// and lightcone from Q[j] where Q[j] - Q[i] is lightlike
// and also near lightcone from Q[k], which intersects
// the line from Q[i] to Q[j].  
if ( 3 <= indexA && indexA <= nverts - 3)
{
  int numberofcases = 4;
  double prefactor = alpha[16]/numberofcases;
  rhoval = rhoval + prefactor*rho_ijkA(1,        nverts,   indexA,   ell);
  rhoval = rhoval + prefactor*rho_ijkA(nverts,   1,        indexA+1, ell);
  rhoval = rhoval + prefactor*rho_ijkA(indexA+1, indexA,   nverts,   ell);
  rhoval = rhoval + prefactor*rho_ijkA(indexA,   indexA+1, 1,        ell);
}
else // if not 3 <= A <= N
{
  int numberofcases = 2*nverts - 4;
  double prefactor = alpha[16]/numberofcases;
  for( int k = 2; k <= indexA; k++)
  {
    rhoval = rhoval + prefactor*rho_ijkA(1, nverts, k, ell);
  }
  for( int k = indexA+1; k <= nverts-1; k++)
  {
    rhoval = rhoval + prefactor*rho_ijkA(nverts, 1, k, ell);
  }
  for( int k = indexA+2; k <= nverts; k++)
  {
    rhoval = rhoval + prefactor*rho_ijkA(indexA+1, indexA, k, ell);
  }
  for( int k = 1; k <= indexA-1; k++)
  {
    rhoval = rhoval + prefactor*rho_ijkA(indexA, indexA+1, k, ell);
  }
}
}
// End of method 16.
//
if( alpha[17] > tiny )
{
// Method 17: same as method 16, but with different parameters,
// distributed near lightcone from Q[i]
// and lightcone from Q[j] where Q[j] - Q[i] is lightlike
// and also near lightcone from Q[k], which intersects
// the line from Q[i] to Q[j].
if ( 3 <= indexA && indexA <= nverts - 3)
{
  int numberofcases = 4;
  double prefactor = alpha[17]/numberofcases;
  rhoval = rhoval + prefactor*rho_ijkB(1,        nverts,   indexA,   ell);
  rhoval = rhoval + prefactor*rho_ijkB(nverts,   1,        indexA+1, ell);
  rhoval = rhoval + prefactor*rho_ijkB(indexA+1, indexA,   nverts,   ell);
  rhoval = rhoval + prefactor*rho_ijkB(indexA,   indexA+1, 1,        ell);
}
else // if not 3 <= A <= N
{
  int numberofcases = 2*nverts - 4;
  double prefactor = alpha[17]/numberofcases;
  for( int k = 2; k <= indexA; k++)
  {
    rhoval = rhoval + prefactor*rho_ijkB(1, nverts, k, ell);
  }
  for( int k = indexA+1; k <= nverts-1; k++)
  {
    rhoval = rhoval + prefactor*rho_ijkB(nverts, 1, k, ell);
  }
  for( int k = indexA+2; k <= nverts; k++)
  {
    rhoval = rhoval + prefactor*rho_ijkB(indexA+1, indexA, k, ell);
  }
  for( int k = 1; k <= indexA-1; k++)
  {
    rhoval = rhoval + prefactor*rho_ijkB(indexA, indexA+1, k, ell);
  }
}
}
// End of method 17.  
//
if( alpha[18] > tiny )
{
// Method 18: distributed near lightcone from Q[i]
// and lightcone from Q[j] where j = i + 1 
// so Q[j] - Q[i] is lightlike
// and also near lightcone from Q[k],
// k = j + 1 or else k = i - 1.
// That is a total of 2*N cases.  
int numberofcases = 2*nverts;
double prefactor = alpha[18]/numberofcases;
int i, j, k;
for( i = 1; i <= nverts - 2; i++)
{
  j = i + 1;
  k = i + 2;
  rhoval = rhoval + prefactor*rho_ijkA(i, j, k, ell);
}
for( i = 2; i <= nverts - 1; i++)
{
  j = i + 1;
  k = i - 1;
  rhoval = rhoval + prefactor*rho_ijkA(i, j, k, ell);
}
i = nverts - 1;
j = nverts;
k = 1;
rhoval = rhoval + prefactor*rho_ijkA(i, j, k, ell);
i = nverts;
j = 1;
k = 2;
rhoval = rhoval + prefactor*rho_ijkA(i, j, k, ell);
i = 1;
j = 2;
k = nverts;
rhoval = rhoval + prefactor*rho_ijkA(i, j, k, ell);
i = nverts;
j = 1;
k = nverts - 1;
rhoval = rhoval + prefactor*rho_ijkA(i, j, k, ell);
}
// End of method 18.
//
if( alpha[19] > tiny )
{
// Method 19: distributed near lightcone from Q[i]
// and lightcone from Q[j] where j = i + 1 
// so Q[j] - Q[i] is lightlike
// and also near lightcone from Q[k],
// k = j + 1 or else k = i - 1.
// That is a total of 2*N cases.
// This is the same as method 18, but with different parameters.
int numberofcases = 2*nverts;
double prefactor = alpha[19]/numberofcases;
int i, j, k;
for( i = 1; i <= nverts - 2; i++)
{
  j = i + 1;
  k = i + 2;
  rhoval = rhoval + prefactor*rho_ijkB(i, j, k, ell);
}
for( i = 2; i <= nverts - 1; i++)
{
  j = i + 1;
  k = i - 1;
  rhoval = rhoval + prefactor*rho_ijkB(i, j, k, ell);
}
i = nverts - 1;
j = nverts;
k = 1;
rhoval = rhoval + prefactor*rho_ijkB(i, j, k, ell);
i = nverts;
j = 1;
k = 2;
rhoval = rhoval + prefactor*rho_ijkB(i, j, k, ell);
i = 1;
j = 2;
k = nverts;
rhoval = rhoval + prefactor*rho_ijkB(i, j, k, ell);
i = nverts;
j = 1;
k = nverts - 1;
rhoval = rhoval + prefactor*rho_ijkB(i, j, k, ell);
}
// End of method 19.
//
if( alpha[20] > tiny )
{
  // Method 20: distributed near lightcone from Q[i]
  // and lightcone from Q[j] where Q[j] - Q[i] is lightlike
  // and also near lightcone from Q[k], which intersects
  // the line from Q[i] to Q[j].  
  int numberofcases = 2*nverts - 8;
  double prefactor = alpha[20]/numberofcases;
  if( indexA == 1 )
  {
    for(int k = 3; k <= nverts-2; k++)
    {
      rhoval = rhoval + prefactor*rho_ijkA(nverts, 1, k, ell);
    }
    for(int k = 4; k <= nverts-1; k++)
    {
      rhoval = rhoval + prefactor*rho_ijkA(indexA+1, indexA, k, ell);
    }
  }
  else if(indexA == nverts-1)
  {
    for(int k = 3; k <= nverts-2; k++)
    {
      rhoval = rhoval + prefactor*rho_ijkA(1, nverts, k, ell);
    }
    for(int k = 2; k <= nverts-3; k++)
    {
      rhoval = rhoval + prefactor*rho_ijkA(indexA, indexA+1, k, ell);
    }
  }
  else // indexA is not 1 or nverts-1
  {
    for( int k = 3; k <= indexA; k++)
    {
      rhoval = rhoval + prefactor*rho_ijkA(1, nverts, k, ell);
    }
    for( int k = indexA+1; k <= nverts-2; k++)
    {
      rhoval = rhoval + prefactor*rho_ijkA(nverts, 1, k, ell);
    }
    for( int k = indexA+3; k <= nverts; k++)
    {
      rhoval = rhoval + prefactor*rho_ijkA(indexA+1, indexA, k, ell);
    }
    for( int k = 1; k <= indexA-2; k++)
    {
      rhoval = rhoval + prefactor*rho_ijkA(indexA, indexA+1, k, ell);
    }    
  }
}
// End of method 20.
if( alpha[21] > tiny )
{
  // Method 21: distributed near lightcone from Q[i]
  // and lightcone from Q[j] where Q[j] - Q[i] is lightlike
  // and also near lightcone from Q[k], which intersects
  // the line from Q[i] to Q[j].  Same as method 20, but
  // with different parameters, using rho_ijkB instead
  // of rho_ijkA.
  int numberofcases = 2*nverts - 8;
  double prefactor = alpha[21]/numberofcases;
  if( indexA == 1 )
  {
    for(int k = 3; k <= nverts-2; k++)
    {
      rhoval = rhoval + prefactor*rho_ijkB(nverts, 1, k, ell);
    }
    for(int k = 4; k <= nverts-1; k++)
    {
      rhoval = rhoval + prefactor*rho_ijkB(indexA+1, indexA, k, ell);
    }
  }
  else if(indexA == nverts-1)
  {
    for(int k = 3; k <= nverts-2; k++)
    {
      rhoval = rhoval + prefactor*rho_ijkB(1, nverts, k, ell);
    }
    for(int k = 2; k <= nverts-3; k++)
    {
      rhoval = rhoval + prefactor*rho_ijkB(indexA, indexA+1, k, ell);
    }
  }
  else // indexA is not 1 or nverts-1
  {
    for( int k = 3; k <= indexA; k++)
    {
      rhoval = rhoval + prefactor*rho_ijkB(1, nverts, k, ell);
    }
    for( int k = indexA+1; k <= nverts-2; k++)
    {
      rhoval = rhoval + prefactor*rho_ijkB(nverts, 1, k, ell);
    }
    for( int k = indexA+3; k <= nverts; k++)
    {
      rhoval = rhoval + prefactor*rho_ijkB(indexA+1, indexA, k, ell);
    }
    for( int k = 1; k <= indexA-2; k++)
    {
      rhoval = rhoval + prefactor*rho_ijkB(indexA, indexA+1, k, ell);
    }    
  }
}
// End of method 21.
//
return rhoval;
} // End of function rhoell().
//--------------------------------------------------------------------------------------------------
// Return new 4-vector ell distributed near lightcones i and j, where Q[j] - Q[i] is lightlike,
// and also near the intersection of i and j and a third lightcone k.
// This function uses the ``a'' set of parameters.
DESreal4vector DESloopmomenta::ell_ijkA(int i, int j, int k)
{
  static const double twoPI = 6.283185307179586;
  const double sqrt2 = 1.414213562373095;
  DESreal4vector K = Q[j] - Q[i];
  if(abs(square(K)) > small)
  {
    cout << "Vector K was supposed to be lightlike but K^2 is " << square(K) << endl;
    cout << "i = " << i << " j = " << j << " k = " << k << endl;
    exit(1);     
  }
  DESreal4vector L = Q[k] - Q[i];
  double Lsq = square(L);
  double dotKL = dot(K,L);
  if (abs(dotKL) < tiny)
  {
    cout << "dot(K,L) should have been non-zero." << endl;
    exit(1);
  }
  double a = Lsq/2.0/dotKL;
  // Define vectors in frame 1,
  // related to frame 0 by possible reflections in t and z.
  int sign0 = 1;
  if( K(0) < 0.0 ) sign0 = -1;
  int sign3 = 1;
  if( K(3) < 0.0 ) sign3 = -1;
  DESreal4vector K1(sign0*K(0), K(1), K(2), sign3*K(3));
  DESreal4vector L1(sign0*L(0), L(1), L(2), sign3*L(3));
  // Define vectors in frame 2,
  // related by a rotation such that K2 is in plus direction.
  DESreal4vector K2(K1(0), 0.0, 0.0, K1(0));
  DESreal4vector K1vec(0.0, K1(1), K1(2), K1(3));
  DESreal4vector K2vec(0.0, K2(1), K2(2), K2(3));
  DESreal4vector L2 = boost(K2vec, K1vec, L1);
  // Define vectors in frame 3,
  // related by a boost such that L3(perp) = 0.
  // L3(-) = L2(-); L3(+) = L2(+) - L2(perp)^2/2/L2(-) ; L3(perp) = 0.
  // K3 = K2.
  // The boost is
  // v3(+) = v2(+) - v2(perp)*u(perp) + 0.5*v2(-)*u(perp)^2
  // v3(-) = v2(-)
  // v3(perp) = v2(perp) - v2(-)*u(perp)
  // where u(perp) = L2(perp)/L2(-).
  // We do not actually use these vectors, so we calculate only
  double K3plus = sqrt2*K2(0);
  double L2minus = (L2(0) - L2(3))/sqrt2;
  double u_1 = L2(1)/L2minus;
  double u_2 = L2(2)/L2minus;     
  // We construct l_3, in the 3 system.
  double temp = 2.0*r.random() - 1.0;
  double x = pow( abs(temp), 1.0/Allike3_a );
  if (temp < 0) x = - x;
  x = tildeallike3_a*x + a;
  temp =  2.0*r.random() - 1.0;
  double lsq = Msqllike3_a/( pow(abs(temp),-1.0/Bllike3_a) - 1.0 );
  if( temp < 0.0 ) lsq = - lsq;
  temp =  2.0*r.random() - 1.0;
  double lmKsq = Msqllike3_a/( pow(abs(temp),-1.0/Bllike3_a) - 1.0 );
  if( temp < 0.0 ) lmKsq = - lmKsq;
  double lperpsq = -(1.0 - x)*lsq - x*lmKsq;
  // We need lperpsq > 0. If it is not, we could choose a new point,
  // but it is faster to just reverse lsq and lmKsq.
  if(lperpsq < 0.0)
  {
    lsq = - lsq;
    lmKsq = - lmKsq;
    lperpsq = - lperpsq;
  }
  double lperpabs = sqrt(lperpsq);
  double phi = twoPI*r.random();
  double l3plus = x*K3plus;
  double l3minus = (lsq - lmKsq)/2.0/K3plus;
  double l3_1 = lperpabs*cos(phi);
  double l3_2 = lperpabs*sin(phi);
  // Now we transform back to the c.m. frame
  // and add Q[i] to form ell.
  double l2minus = l3minus;
  double l2plus = l3plus + (l3_1*u_1 + l3_2*u_2) + 0.5*l3minus*(u_1*u_1 + u_2*u_2);
  double l2_1 = l3_1 + l3minus*u_1;
  double l2_2 = l3_2 + l3minus*u_2;
  DESreal4vector l2((l2plus + l2minus)/sqrt2, l2_1, l2_2, (l2plus - l2minus)/sqrt2);
  DESreal4vector l1 = boost(K1vec,K2vec,l2);
  DESreal4vector l(sign0*l1(0), l1(1), l1(2), sign3*l1(3));
  DESreal4vector ell = l + Q[i];
  return ell;
}  
// End of function ell_ijkA(i, j, k)
//--------------------------------------------------------------------------------------------------
// Return new 4-vector ell distributed near lightcones i and j, where Q[j] - Q[i] is lightlike,
// and also near the intersection of i and j and a third lightcone k.
// This function uses the ``b'' set of parameters.
DESreal4vector DESloopmomenta::ell_ijkB(int i, int j, int k)
{
  static const double twoPI = 6.283185307179586;
  const double sqrt2 = 1.414213562373095;
  DESreal4vector K = Q[j] - Q[i];
  if(abs(square(K)) > small)
  {
    cout << "Vector K was supposed to be lightlike but K^2 is " << square(K) << endl;
    cout << "i = " << i << " j = " << j << " k = " << k << endl;
    exit(1);     
  }
  DESreal4vector L = Q[k] - Q[i];
  double Lsq = square(L);
  double dotKL = dot(K,L);
  if (abs(dotKL) < tiny)
  {
    cout << "dot(K,L) should have been non-zero." << endl;
    exit(1);
  }
  double a = Lsq/2.0/dotKL;
  // Define vectors in frame 1,
  // related to frame 0 by possible reflections in t and z.
  int sign0 = 1;
  if( K(0) < 0.0 ) sign0 = -1;
  int sign3 = 1;
  if( K(3) < 0.0 ) sign3 = -1;
  DESreal4vector K1(sign0*K(0), K(1), K(2), sign3*K(3));
  DESreal4vector L1(sign0*L(0), L(1), L(2), sign3*L(3));
  // Define vectors in frame 2,
  // related by a rotation such that K2 is in plus direction.
  DESreal4vector K2(K1(0), 0.0, 0.0, K1(0));
  DESreal4vector K1vec(0.0, K1(1), K1(2), K1(3));
  DESreal4vector K2vec(0.0, K2(1), K2(2), K2(3));
  DESreal4vector L2 = boost(K2vec, K1vec, L1);
  // Define vectors in frame 3,
  // related by a boost such that L3(perp) = 0.
  // L3(-) = L2(-); L3(+) = L2(+) - L2(perp)^2/2/L2(-) ; L3(perp) = 0.
  // K3 = K2.
  // The boost is
  // v3(+) = v2(+) - v2(perp)*u(perp) + 0.5*v2(-)*u(perp)^2
  // v3(-) = v2(-)
  // v3(perp) = v2(perp) - v2(-)*u(perp)
  // where u(perp) = L2(perp)/L2(-).
  // We do not actually use these vectors, so we calculate only
  double K3plus = sqrt2*K2(0);
  double L2minus = (L2(0) - L2(3))/sqrt2;
  double u_1 = L2(1)/L2minus;
  double u_2 = L2(2)/L2minus;     
  // We construct l_3, in the 3 system.
  double temp = 2.0*r.random() - 1.0;
  double x = pow( abs(temp), 1.0/Allike3_b );
  if (temp < 0) x = - x;
  x = tildeallike3_b*x + a;
  temp =  2.0*r.random() - 1.0;
  double lsq = Msqllike3_b/( pow(abs(temp),-1.0/Bllike3_b) - 1.0 );
  if( temp < 0.0 ) lsq = - lsq;
  temp =  2.0*r.random() - 1.0;
  double lmKsq = Msqllike3_b/( pow(abs(temp),-1.0/Bllike3_b) - 1.0 );
  if( temp < 0.0 ) lmKsq = - lmKsq;
  double lperpsq = -(1.0 - x)*lsq - x*lmKsq;
  // We need lperpsq > 0. If it is not, we could choose a new point,
  // but it is faster to just reverse lsq and lmKsq.
  if(lperpsq < 0.0)
  {
    lsq = - lsq;
    lmKsq = - lmKsq;
    lperpsq = - lperpsq;
  }
  double lperpabs = sqrt(lperpsq);
  double phi = twoPI*r.random();
  double l3plus = x*K3plus;
  double l3minus = (lsq - lmKsq)/2.0/K3plus;
  double l3_1 = lperpabs*cos(phi);
  double l3_2 = lperpabs*sin(phi);
  // Now we transform back to the c.m. frame
  // and add Q[i] to form ell.
  double l2minus = l3minus;
  double l2plus = l3plus + (l3_1*u_1 + l3_2*u_2) + 0.5*l3minus*(u_1*u_1 + u_2*u_2);
  double l2_1 = l3_1 + l3minus*u_1;
  double l2_2 = l3_2 + l3minus*u_2;
  DESreal4vector l2((l2plus + l2minus)/sqrt2, l2_1, l2_2, (l2plus - l2minus)/sqrt2);
  DESreal4vector l1 = boost(K1vec,K2vec,l2);
  DESreal4vector l(sign0*l1(0), l1(1), l1(2), sign3*l1(3));
  DESreal4vector ell = l + Q[i];
  return ell;
}  
// End of function ell_ijkB(i, j, k)
//--------------------------------------------------------------------------------------------------
// Return density for ell chosen near lightcones i and j, where Q[j] - Q[i] is lightlike,
// and also near the intersection of i and j and a third lightcone k.
// This function uses the ``a'' set of parameters.
double DESloopmomenta::rho_ijkA(int i, int j, int k, DESreal4vector ell)
{
  static const double twoPI = 6.283185307179586;
  static const double zero = 0.0;
  double rho;
  DESreal4vector K = Q[j] - Q[i];
  DESreal4vector L = Q[k] - Q[i];
  double Lsq = square(L);
  double dotKL = dot(K,L);
  if (abs(dotKL) < tiny)
  {
    cout << "dot(K,L) should have been non-zero." << endl;
    exit(1);
  }
  double a = Lsq/2.0/dotKL;
  double lsq = square(ell - Q[i]);
  if( abs(lsq) < tiny ) return huge;
  double lmKsq = square(ell - Q[j]);
  if( abs(lmKsq) < tiny ) return huge;
  double x = ( dot(ell - Q[i],Q[k] - Q[i]) - a * dot(ell - Q[i],Q[j] - Q[i]) )/dotKL;
  if ( x > a + tildeallike3_a + small) return zero;
  if ( x < a - tildeallike3_a - small) return zero;
  if ( abs(x-a) < tiny ) return huge;
  rho = Allike3_a*pow(Bllike3_a,2)/twoPI/tildeallike3_a;
  rho *= pow(abs(tildeallike3_a/(x - a)),1.0 - Allike3_a);
  rho *= pow(1.0 + Msqllike3_a/abs(lsq), - Bllike3_a - 1.0);
  rho *= Msqllike3_a/pow(lsq,2);
  rho *= pow(1.0 + Msqllike3_a/abs(lmKsq), - Bllike3_a - 1.0);
  rho *= Msqllike3_a/pow(lmKsq,2);
  return rho;
}
// End of function rho_ijkA(i, j, k, ell)
//--------------------------------------------------------------------------------------------------
// Return density for ell chosen near lightcones i and j, where Q[j] - Q[i] is lightlike,
// and also near the intersection of i and j and a third lightcone k.
// This function uses the ``b'' set of parameters.
double DESloopmomenta::rho_ijkB(int i, int j, int k, DESreal4vector ell)
{
  static const double twoPI = 6.283185307179586;
  static const double zero = 0.0;
  double rho;
  DESreal4vector K = Q[j] - Q[i];
  DESreal4vector L = Q[k] - Q[i];
  double Lsq = square(L);
  double dotKL = dot(K,L);
  if (abs(dotKL) < tiny)
  {
    cout << "dot(K,L) should have been non-zero." << endl;
    exit(1);
  }
  double a = Lsq/2.0/dotKL;
  double lsq = square(ell - Q[i]);
  if( abs(lsq) < tiny ) return huge;
  double lmKsq = square(ell - Q[j]);
  if( abs(lmKsq) < tiny ) return huge;
  double x = ( dot(ell - Q[i],Q[k] - Q[i]) - a * dot(ell - Q[i],Q[j] - Q[i]) )/dotKL;
  if ( x > a + tildeallike3_b + small) return zero;
  if ( x < a - tildeallike3_b - small) return zero;
  if ( abs(x-a) < tiny ) return huge;
  rho = Allike3_b*pow(Bllike3_b,2)/twoPI/tildeallike3_b;
  rho *= pow(abs(tildeallike3_b/(x - a)),1.0 - Allike3_b);
  rho *= pow(1.0 + Msqllike3_b/abs(lsq), - Bllike3_b - 1.0);
  rho *= Msqllike3_b/pow(lsq,2);
  rho *= pow(1.0 + Msqllike3_b/abs(lmKsq), - Bllike3_b - 1.0);
  rho *= Msqllike3_b/pow(lmKsq,2);
  return rho;
}
// End of function rho_ijkB(i, j, k, ell)
// End of implementation code for DESloopmomenta.
//==================================================================================================
// Implementation code for class DESdeform.
// Constructor.
DESdeform::DESdeform(const int& nvertsIN, const int& indexAIN, const DESreal4vector QIN[maxsize])
  :indx(nvertsIN)  // Construct instance of index function.
{
  for(int mu = 0; mu<=3; mu++) zero4vector.in(mu) = 0.0;
  small = 1.0e-10;  // A small number. Use to avoid 1/0.
  nverts = nvertsIN;
  indexA = indexAIN;
  indexAp1 = indexA + 1;             // We know indexA < nverts;
  for(int i = 0; i < maxsize; i++)
  {
    for(int mu = 0; mu<=3; mu++)
    {
      Q[i].in(mu) = 0.0;
    }
  }
  for(int i = 1; i <= nverts; i++) Q[i] = QIN[i];
  P = Q[nverts] - Q[1];
  Pbar = Q[indexA] - Q[indexAp1];
  PdotPbar = dot(P,Pbar);
  if (PdotPbar < 0.0) cout << "Warning: PdotPbar < 0 ." << endl;
  M1 = 0.05*sqrt(PdotPbar);  // Parameter.
  M2 = sqrt(PdotPbar);     // Parameter.
  M3 = sqrt(PdotPbar);     // Parameter.
  gamma1 = 0.7;            // Parameter.
  gamma2 = 1.0;            // Parameter.
  Pc.in(0) = P(0);         // The covariant components of P.
  Pbarc.in(0) = Pbar(0);   // The covariant components of Pbar.
  for(int mu = 1; mu <= 3; mu++)
  {
    Pc.in(mu) = - P(mu);
    Pbarc.in(mu) = - Pbar(mu);
  }
  // Some checks.
  if (abs(P(3) + Pbar(3)) > 1.0e-10)
  {
    cout << "Warning: we were supposed to be working in the c.m. frame." << endl;
  }
  if( abs(P(1)) + abs(Pbar(1)) + abs(P(2)) + abs(Pbar(2)) > 1.0e-10)
  {
    cout << "Warning: we were supposed to allign the beams with the z-axis." << endl;
  }
  if( P(0) < 0.0 || Pbar(0) < 0.0)
  {
    cout << "Warning: we were supposed have positive beam energies." << endl;
  }
  for(int i = 1; i <= nverts; i++)
  {
    if( abs(square(Q[i]- Q[indx(i+1)])) >  1.0e-10 )
    {
      cout << "Warning: we should have massless photons." << endl;
    }
  }
  // End checks.
}
// There is no default constructor.
//--------------------------------------------------------------------------------------------------
// Compute complex lc for a point ell, also compute jacobian;
void DESdeform::deform(DESreal4vector& ellIN, DEScmplx4vector& lcOUT, complex<double>& jacobian)
{
  static const complex<double> ii(0.0,1.0);
  double cj;
  DESreal4vector kappaj;
  ell = ellIN;
  double x    = dot(ell - Q[indexAp1], Pbar) /PdotPbar;
  double xbar = dot(ell - Q[1], P) /PdotPbar;
  kappa = zero4vector;
  DESreal4vector gradlogcj;
  double gradkappa[4][4] = {{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},
                            {0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}};
  //
  double gval = g(ell);
  DESreal4vector gradloggval = gradlogg(ell);
  double cjsum = 0.0;
  DESreal4vector gradcjsum(0.0,0.0,0.0,0.0);
  for (int j=1; j <= nverts; j++)
  {
    // Calculate c_j and gradlogc_j.
    c_j(j, gval, gradloggval, cj, gradlogcj);
    cjsum += cj;
    gradcjsum = gradcjsum + cj*gradlogcj;
    // Use c_j and gradlogc_j to calculate additions 
    // to kappa and its gradient.
    kappaj = - cj * (ell - Q[j]);
    kappa = kappa + kappaj;
    for(int mu = 0; mu <= 3; mu++)
    {
      gradkappa[mu][mu] -= cj;
      for(int nu = 0; nu <= 3; nu++)
      {
        gradkappa[mu][nu] -= (ell(mu) - Q[j](mu))*cj*gradlogcj(nu);
      }
    }
  }
  //
  for (int indexpm = -1; indexpm <= 1; indexpm += 2)
  {
    // Calculate c_pm and gradlogc_pm, called here cj and gradlogcj.
    c_pm(indexpm, x, xbar, cj, gradlogcj);
    // Use c_pm and gradlogc_pm to calculate additions 
    // to kappa and its gradient.
    kappaj = (indexpm*cj) * (P + Pbar);
    kappa = kappa + kappaj;
    for(int mu = 0; mu <= 3; mu++)
    {
      for(int nu = 0; nu <= 3; nu++)
      {
        gradkappa[mu][nu] += (P(mu) + Pbar(mu))*(indexpm*cj)*gradlogcj(nu);
      }
    }     
  }
// ============ End of calculation of kappa0 and grad kappa0 ============
// So far kappa and gradkappa are really kappa_0 and its gradient.
// Now we need the lambda_i.
  double lambda_i;
  // Set our maximum possible value for lambda. One could change the value.
  double lambda = 1.0;
  int ifound = - 777;
  if( 0.25/cjsum < lambda)
  {
    lambda = 0.25/cjsum;
    ifound = 0;           // lambda_0 is the smallest so far.
  }
  double kappasq, lisq, dotlikappa;
  for(int i = 1; i <= nverts; i++)
  {
    kappasq = square(kappa);
    dotlikappa = dot(kappa,ell-Q[i]);
    lisq = square(ell-Q[i]);
    if( kappasq*lisq > 2.0*pow(dotlikappa,2) ) // Region 1.
    {
      lambda_i = 0.5*sqrt(lisq/kappasq);
    }
    else if( kappasq*lisq > 0.0 ) // Region 2, 0.0 < kappasq*lisq < 2.0*pow(dotlikappa,2).
    {
      lambda_i = sqrt(4.0*pow(dotlikappa,2) - kappasq*lisq);
      lambda_i *= 0.5/abs(kappasq);
    }
    else // Region 3, kappasq*lisq < 0.0.
    {
      lambda_i = sqrt(4.0*pow(dotlikappa,2) - 2.0*kappasq*lisq);
      lambda_i *= 0.5/abs(kappasq);    
    }
    // We now have lambda_i.
    if(lambda_i < lambda)
    {
      lambda = lambda_i;
      ifound = i;
    }
  }
  // We have lambda.
  // Now calculate its gradient.
  DESreal4vector gradlambda;
  if(ifound < 0) gradlambda = zero4vector;
  else if(ifound == 0)
  {
    gradlambda = (- 0.25/pow(cjsum,2))*gradcjsum;
  }
  else if(ifound <= nverts)
  {
    kappasq = square(kappa);
    dotlikappa = dot(kappa,ell-Q[ifound]);
    lisq = square(ell-Q[ifound]);
    // We need the gradient of kappasq.
    DESreal4vector gradkappasq;
    for(int nu = 0; nu <= 3; nu++)
    {
      double temp = 2.0*kappa(0)*gradkappa[0][nu];
      for(int mu = 1; mu <= 3; mu++)
      {
        temp -= 2.0*kappa(mu)*gradkappa[mu][nu];
      }
      gradkappasq.in(nu) = temp;
    }
    // We need the gradient of lisq.
    DESreal4vector gradlisq;
    gradlisq.in(0) = 2.0*(ell(0) - Q[ifound](0));
    for(int nu = 1; nu <= 3; nu++)
    {
      gradlisq.in(nu) = - 2.0*(ell(nu) - Q[ifound](nu));
    }
    // We need the gradient of dotlikappa, for which we need
    // also kappac, the covariant components of kappa.
    DESreal4vector graddotlikappa;
    DESreal4vector kappacov(kappa(0), -kappa(1), -kappa(2), -kappa(3));
    for(int nu = 0; nu <= 3; nu++)
    {
      double temp = kappacov(nu);
      temp += (ell(0) - Q[ifound](0))*gradkappa[0][nu];
      for(int mu = 1; mu <= 3; mu++)
      {
        temp -= (ell(mu) - Q[ifound](mu))*gradkappa[mu][nu];
      }
      graddotlikappa.in(nu) = temp;
    }
    // Now we have the ingredients, so we calculate gradlambda.
    DESreal4vector gradlambdasq, gradnumerator;
    double numerator;
    if( kappasq*lisq > 2.0*pow(dotlikappa,2) ) // Region 1.
    {
      //lambda_i = 0.5*sqrt(lisq/kappasq);
      numerator = kappasq*lisq;
      gradnumerator = kappasq*gradlisq + gradkappasq*lisq;
    }
    else if( kappasq*lisq > 0.0 ) // Region 2, 0.0 < kappasq*lisq < 2.0*pow(dotlikappa,2).
    {
      //lambda_i = sqrt(4.0*pow(dotlikappa,2) - kappasq*lisq);
      //lambda_i *= 0.5/abs(kappasq);
      numerator = 4.0*pow(dotlikappa,2) - kappasq*lisq;
      gradnumerator = 8.0*dotlikappa*graddotlikappa 
                     - kappasq*gradlisq - gradkappasq*lisq;
    }
    else // Region 3, kappasq*lisq < 0.0.
    {
      //lambda_i = sqrt(4.0*pow(dotlikappa,2) - 2.0*kappasq*lisq);
      //lambda_i *= 0.5/abs(kappasq);
      numerator = 4.0*pow(dotlikappa,2) - 2.0*kappasq*lisq;
      gradnumerator = 8.0*dotlikappa*graddotlikappa
                     - 2.0*kappasq*gradlisq - 2.0*gradkappasq*lisq;
    }
    gradlambdasq = gradnumerator - (2.0*numerator/kappasq)*gradkappasq;
    gradlambdasq = (1.0/4.0/pow(kappasq,2))*gradlambdasq;
    gradlambda = (1.0/2.0/lambda)*gradlambdasq;
  }
  // We now have lambda and gradlambda. 
  // Use this to get the new kappa and gradkappa.
  for(int nu = 0; nu <= 3; nu++)
  {
    for(int mu = 0; mu <= 3; mu++)
    {
      gradkappa[mu][nu] = lambda*gradkappa[mu][nu] + kappa(mu)*gradlambda(nu);
    }
  }
  kappa = lambda*kappa;
//
// ============ End of calculation of kappa and grad kappa ============
//
  DEScmplx4vector kappac, ellc;
  for(int mu = 0; mu <= 3; mu++)  // We need this for the type conversion.
  {
    kappac.in(mu) = kappa(mu);
    ellc.in(mu) = ell(mu);
  }
  lc = ellc + ii * kappac;
  lcOUT = lc;
  complex<double> gradlc[maxsize][maxsize]; 
       // We need this dimension for Determinant(gradlc,dimension).
  for(int mu = 0; mu <= 3; mu++)
  {
    for(int nu = 0; nu <= 3; nu++)
    {
      if(mu == nu) gradlc[mu][nu] = 1.0 + ii * gradkappa[mu][nu];
      else         gradlc[mu][nu] = ii * gradkappa[mu][nu];
    }
  }
  int dimension = 4;
  jacobian = Determinant(gradlc,dimension);
  return;
} // End of DESdeform::deform(ellIN, lcOUT, jacobian)
//--------------------------------------------------------------------------------------------------
void DESdeform::c_j(int& j, double& gval, DESreal4vector& gradloggval, 
                    double& cj, DESreal4vector& gradlogcj)
{
cj = 0.0;                    // This will remain unless changed.
gradlogcj = zero4vector;     // This will remain unless changed.
//
// There are three parts: 
//   A = 1,
//   A = N-1,
//   1 < A < N-1.
//
//------------------------------
//         A = 1
//------------------------------
if( indexA == 1 )
{
  //  ---------  3 <= j <= N-1  ---------
  if(j >= 3 && j <= nverts-1)
  {
    if( (! ellinCminus(j-1)) && (! ellinCplus(j+1)) &&
        (! ellinCminus(1)) && (! ellinCplus(1)) )
    {
      cj  = hminus(ell - Q[j-1]);
      cj *= hplus(ell - Q[j+1]);
      cj *= hminus(ell - Q[1]);
      cj *= hplus(ell - Q[1]);
      cj *= gval;
      //
      gradlogcj = gradloghminus(ell - Q[j-1]);
      gradlogcj = gradlogcj + gradloghplus(ell - Q[j+1]);
      gradlogcj = gradlogcj + gradloghminus(ell - Q[1]);
      gradlogcj = gradlogcj + gradloghplus(ell - Q[1]);
      gradlogcj = gradlogcj + gradloggval;
    } // Close if( (! ellinCminus(j-1)) ... )
  }
  //  ---------      j = 1      ---------
  else if(j == 1)
  {
    if( (! ellinCminus(nverts-1)) && (! ellinCplus(3)) )
    {
      cj = hminus(ell - Q[nverts-1]);
      cj *= hplus(ell - Q[3]);
      cj *= gval;
      //
      gradlogcj = gradloghminus(ell - Q[nverts-1]);
      gradlogcj = gradlogcj + gradloghplus(ell - Q[3]);
      gradlogcj = gradlogcj + gradloggval;
    } // Close if( (! ellinCminus(nverts-1)) ... )
  }
  //  ---------      j = 2      ---------
  else if(j == 2)
  {
    if( (! ellinCplus(3)) && (! ellinCplus(1)) )
    {
      cj  = hplus(ell - Q[3]);
      cj *= hplus(ell - Q[1]);
      cj *= gval;
      //
      gradlogcj = gradloghplus(ell - Q[3]);
      gradlogcj = gradlogcj + gradloghplus(ell- Q[1]);
      gradlogcj = gradlogcj + gradloggval;
    } // Close if( (! ellinCplus(3)) ... )
  }
  //  ---------      j = N      ---------
  else if( j == nverts )
  {
    if( (! ellinCminus(nverts-1)) && (! ellinCminus(1)) )
    {
      cj  = hminus(ell - Q[nverts-1]);
      cj *= hminus(ell - Q[1]);
      cj *= gval;
      //
      gradlogcj = gradloghminus(ell - Q[nverts-1]);
      gradlogcj = gradlogcj + gradloghminus(ell- Q[1]);
      gradlogcj = gradlogcj + gradloggval;
    } // Close if( (! ellinCminus(nverts-1)) ... ) 
  } 
  else
  {
    cout << "Bad value of j in c_j for A = 1: " << j << endl;
    exit(1);
  }
}
//------------------------------
//         A = N-1
//------------------------------
else if (indexA == nverts-1)
{
  //  ---------  3 <= j <= N-1  ---------
  if(j >= 2 && j <= nverts-2)
  {
    if( (! ellinCminus(j-1)) && (! ellinCplus(j+1)) &&
        (! ellinCminus(nverts)) && (! ellinCplus(nverts)) )
    {
      cj  = hminus(ell - Q[j-1]);
      cj *= hplus(ell - Q[j+1]);
      cj *= hminus(ell - Q[nverts]);
      cj *= hplus(ell - Q[nverts]);
      cj *= gval;
      //
      gradlogcj = gradloghminus(ell - Q[j-1]);
      gradlogcj = gradlogcj + gradloghplus(ell - Q[j+1]);
      gradlogcj = gradlogcj + gradloghminus(ell - Q[nverts]);
      gradlogcj = gradlogcj + gradloghplus(ell - Q[nverts]);
      gradlogcj = gradlogcj + gradloggval;
    } // Close if( (! ellinCminus(j-1)) ... )
  }
  //  ---------      j = 1      ---------
  else if(j == 1)
  {
    if( (! ellinCplus(2)) && (! ellinCplus(nverts)) )
    {
      cj = hplus(ell - Q[2]);
      cj *= hplus(ell - Q[nverts]);
      cj *= gval;
      //
      gradlogcj = gradloghplus(ell - Q[2]);
      gradlogcj = gradlogcj + gradloghplus(ell - Q[nverts]);
      gradlogcj = gradlogcj + gradloggval;
    } // Close if( (! ellinCplus(2)) ... )
  }
  //  ---------      j = N-1      ---------
  else if(j == nverts-1)
  {
    if( (! ellinCminus(nverts-2)) && (! ellinCminus(nverts)) )
    {
      cj  = hminus(ell - Q[nverts-2]);
      cj *= hminus(ell - Q[nverts]);
      cj *= gval;
      //
      gradlogcj = gradloghminus(ell - Q[nverts-2]);
      gradlogcj = gradlogcj + gradloghminus(ell- Q[nverts]);
      gradlogcj = gradlogcj + gradloggval;
    } // Close if( (! ellinCplus(nverts-2) ... )
  }
  //  ---------      j = N      ---------
  else if( j == nverts )
  {
    if( (! ellinCplus(2)) && (! ellinCminus(nverts-2)) )
    {
      cj  = hplus(ell - Q[2]);
      cj *= hminus(ell - Q[nverts-2]);
      cj *= gval;
      //
      gradlogcj = gradloghplus(ell - Q[2]);
      gradlogcj = gradlogcj + gradloghminus(ell- Q[nverts-2]);
      gradlogcj = gradlogcj + gradloggval;
    } // Close if( (! ellinCminus(2)) ... ) 
  }
  else
  {
    cout << "Bad value of j in c_j for A = N-1: " << j << endl;
    exit(1);
  }
}
//------------------------------
//          1 < A < N-1
//------------------------------
else
{
  //  --------- 2 <= j <= A-1 ---------
  if(j >=2 && j <= indexA-1)
  {
    if( (! ellinCminus(j-1)) && (! ellinCplus(j+1)) && 
        (! ellinCminus(nverts)) && (! ellinCplus(indexAp1)) )
    {
      cj  = hminus(ell - Q[j-1]);
      cj *= hplus(ell - Q[j+1]);
      cj *= hminus(ell - Q[nverts]);
      cj *= hplus(ell - Q[indexAp1]);
      cj *= gval;
      //
      gradlogcj = gradloghminus(ell - Q[j-1]);
      gradlogcj = gradlogcj + gradloghplus(ell - Q[j+1]);
      gradlogcj = gradlogcj + gradloghminus(ell - Q[nverts]);
      gradlogcj = gradlogcj + gradloghplus(ell - Q[indexAp1]);
      gradlogcj = gradlogcj + gradloggval;
    } // Close if( (! ellinCminus(j-1)) ... )
  }
  //  --------- A+2 <= j <= N-1 ---------
  else if(j >= indexA+2 && j <= nverts-1)
  {
    if( (! ellinCminus(j-1)) && (! ellinCplus(j+1)) &&
        (! ellinCminus(indexA)) && (! ellinCplus(1)) )
    {
      cj  = hminus(ell - Q[j-1]);
      cj *= hplus(ell - Q[j+1]);
      cj *= hminus(ell - Q[indexA]);
      cj *= hplus(ell - Q[1]);
      cj *= gval;
      //
      gradlogcj = gradloghminus(ell - Q[j-1]);
      gradlogcj = gradlogcj + gradloghplus(ell - Q[j+1]);
      gradlogcj = gradlogcj + gradloghminus(ell - Q[indexA]);
      gradlogcj = gradlogcj + gradloghplus(ell - Q[1]);
      gradlogcj = gradlogcj + gradloggval;
    } // Close if( (! ellinCminus(j-1)) ... )
  }
  //  ---------      j = 1      ---------
  else if(j == 1)
  {
    if( (! ellinCplus(2)) && (! ellinCminus(nverts-1)) && (! ellinCplus(indexAp1)) )
    {
      cj  = hplus(ell - Q[2]);
      cj *= hminus(ell - Q[nverts-1]);
      cj *= hplus(ell - Q[indexAp1]);
      cj *= gval;
      //
      gradlogcj = gradloghplus(ell - Q[2]);
      gradlogcj = gradlogcj + gradloghminus(ell - Q[nverts-1]);
      gradlogcj = gradlogcj + gradloghplus(ell - Q[indexAp1]);
      gradlogcj = gradlogcj + gradloggval;
    } // Close if( (! ellinCplus(2)) ... )
  }
  //  ---------      j = A      ---------
  else if(j == indexA)
  {
    if( (! ellinCminus(indexA-1)) && (! ellinCplus(indexA+2)) && (! ellinCminus(nverts)) )
    {
      cj  = hminus(ell - Q[indexA-1]);
      cj *= hplus(ell - Q[indexA+2]);
      cj *= hminus(ell - Q[nverts]);
      cj *= gval;
      //
      gradlogcj = gradloghminus(ell - Q[indexA-1]);
      gradlogcj = gradlogcj + gradloghplus(ell - Q[indexA+2]);
      gradlogcj = gradlogcj + gradloghminus(ell- Q[nverts]);
      gradlogcj = gradlogcj + gradloggval;
    } // Close if( (! ellinCminus(indexA-1)) ... )
  }
  //  ---------     j = A+1     ---------
  else if(j == indexAp1)
  {
    if( (! ellinCplus(indexA+2)) && (! ellinCminus(indexA-1)) && (! ellinCplus(1)) )
    {
      cj  = hplus(ell - Q[indexA+2]);
      cj *= hminus(ell - Q[indexA-1]);
      cj *= hplus(ell - Q[1]);
      cj *= gval;
      //
      gradlogcj = gradloghplus(ell - Q[indexA+2]);
      gradlogcj = gradlogcj + gradloghminus(ell - Q[indexA-1]);
      gradlogcj = gradlogcj + gradloghplus(ell- Q[1]);
      gradlogcj = gradlogcj + gradloggval;
    } // Close if( (! ellinCplus(indexA+2) ... )
  }
  //  ---------      j = N      ---------
  else if( j == nverts )
  {
    if( (! ellinCminus(nverts-1)) && (! ellinCplus(2)) && (! ellinCminus(indexA)) )
    {
      cj  = hminus(ell - Q[nverts-1]);
      cj *= hplus(ell - Q[2]);
      cj *= hminus(ell - Q[indexA]);
      cj *= gval;
      //
      gradlogcj = gradloghminus(ell - Q[nverts-1]);
      gradlogcj = gradlogcj + gradloghplus(ell - Q[2]);
      gradlogcj = gradlogcj + gradloghminus(ell- Q[indexA]);
      gradlogcj = gradlogcj + gradloggval;
    } // Close if( (! ellinCminus(nverts-1)) ... ) 
  }
  else
  {
    cout << "Bad value of j in c_j for generic A: " << j << endl;
    exit(1);
  }
}
return;
} // End of c_j(j, gval, gradloggval, cj, gradlogcj).
//--------------------------------------------------------------------------------------------------
void DESdeform::c_pm(int& indexpm, double& x, double& xbar, double& cj, DESreal4vector& gradlogcj)
{
  cj = 0.0;                    // This will remain unless changed.
  gradlogcj = zero4vector;     // This will remain unless changed.
                               //
// There are three parts: 
//   A = 1,
//   A = N-1,
//   1 < A < N-1.
//
//------------------------------
//         A = 1
//------------------------------
if( indexA == 1 )
{
  if(indexpm == 1)
  {
    if( (x + xbar > 0) && (! ellinCminus(nverts)) )
    {
      cj = hminus(ell - Q[nverts]);
      cj *= (x + xbar);
      cj *= gminus(ell);
      //
      gradlogcj = gradloghminus(ell - Q[nverts]);
      gradlogcj = gradlogcj + (1.0/(x + xbar)/PdotPbar)*(Pc + Pbarc);
      gradlogcj = gradlogcj + gradloggminus(ell);
    } // Close if( (x + xbar > 0) ... )    
  }
  else if( indexpm == -1)
  {
    if( (x + xbar < 0) && (! ellinCplus(2)) )
    {
      cj = hplus(ell - Q[2]);
      cj *= -(x + xbar);
      cj *= gplus(ell);
      //
      gradlogcj = gradloghplus(ell - Q[2]);
      gradlogcj = gradlogcj + (1.0/(x + xbar)/PdotPbar)*(Pc + Pbarc);
      gradlogcj = gradlogcj + gradloggplus(ell);
    } // Close if( (x + xbar < 0) ... )    
  }
  else
  {
    cout << "This cannot happen, in c_pm, indexpm is not +1 or -1. " << indexpm << endl;
    exit(1);
  }
}
//------------------------------
//         A = N-1
//------------------------------
else if(indexA == nverts - 1)
{
  if(indexpm == 1)
  {
    if( (x + xbar > 0) && (! ellinCminus(nverts-1)) )
    {
      cj = hminus(ell - Q[nverts-1]);
      cj *= (x + xbar);
      cj *= gminus(ell);
      //
      gradlogcj = gradloghminus(ell - Q[nverts-1]);
      gradlogcj = gradlogcj + (1.0/(x + xbar)/PdotPbar)*(Pc + Pbarc);
      gradlogcj = gradlogcj + gradloggminus(ell);
    } // Close if( (x + xbar > 0) ... )  return;    
  }
  else if( indexpm == -1)
  {
    if( (x + xbar < 0) && (! ellinCplus(1)) )
    {
      cj = hplus(ell - Q[1]);
      cj *= -(x + xbar);
      cj *= gplus(ell);
      //
      gradlogcj = gradloghplus(ell - Q[1]);
      gradlogcj = gradlogcj + (1.0/(x + xbar)/PdotPbar)*(Pc + Pbarc);
      gradlogcj = gradlogcj + gradloggplus(ell);
    } // Close if( (x + xbar < 0) ... )    
  }
  else
  {
    cout << "This cannot happen, in c_pm, indexpm is not +1 or -1. " << indexpm << endl;
    exit(1);
  }
}
//------------------------------
//        1 < A < N-1
//------------------------------
else
{
  if(indexpm == 1)
  {
    if( (x + xbar > 0) && (! ellinCminus(indexA)) && (! ellinCminus(nverts)) )
    {
      cj  = hminus(ell - Q[indexA]);
      cj *= hminus(ell - Q[nverts]);
      cj *= (x + xbar);
      cj *= gminus(ell);
      //
      gradlogcj = gradloghminus(ell - Q[indexA]);
      gradlogcj = gradlogcj + gradloghminus(ell - Q[nverts]);
      gradlogcj = gradlogcj + (1.0/(x + xbar)/PdotPbar)*(Pc + Pbarc);
      gradlogcj = gradlogcj + gradloggminus(ell);
    } // Close if( (x + xbar > 0) ... )    
  }
  else if( indexpm == -1)
  {
    if( (x + xbar < 0) && (! ellinCplus(1)) && (! ellinCplus(indexAp1)) )
    {
      cj  = hplus(ell - Q[1]);
      cj *= hplus(ell - Q[indexAp1]);
      cj *= -(x + xbar);
      cj *= gplus(ell);
      //
      gradlogcj = gradloghplus(ell - Q[1]);
      gradlogcj = gradlogcj + gradloghplus(ell - Q[indexAp1]);
      gradlogcj = gradlogcj + (1.0/(x + xbar)/PdotPbar)*(Pc + Pbarc);
      gradlogcj = gradlogcj + gradloggplus(ell);   
    } // Close if( (x + xbar < 0) ... )      
  }
  else
  {
    cout << "This cannot happen, in c_pm, indexpm is not +1 or -1. " << indexpm << endl;
    exit(1);
  }  
}
return;
}// End of c_pm(index, x, xbar, cj, gradlogcj)
//--------------------------------------------------------------------------------------------------
// Check if ell is inside the positive light cone starting from Q[j];
bool DESdeform::ellinCplus(int j)
{
  return (square(ell - Q[j]) > 0.0 && ell(0) > Q[j](0));
}
//--------------------------------------------------------------------------------------------------
// Check if ell is inside the negative light cone starting from Q[j];
bool DESdeform::ellinCminus(int j)
{
  return (square(ell - Q[j]) > 0.0 && ell(0) < Q[j](0));
}
//--------------------------------------------------------------------------------------------------
// Function that smoothly vanishes if k is inside the backward lightcone;
double DESdeform::hminus(const DESreal4vector k)
{
  double h;
  double E = k(0);
  double absk = sqrt( k(1)*k(1) + k(2)*k(2) + k(3)*k(3));
  if( absk + E < small ) return 0.0;
  else
  {
    double temp = pow(absk + E,2);
    h = temp/(temp + pow(M1,2));
    return h;
  }
}
//--------------------------------------------------------------------------------------------------
// Function that smoothly vanishes if k is inside the forward lightcone;
double DESdeform::hplus(const DESreal4vector k)
{
  double h;
  double E = k(0);
  double absk = sqrt( k(1)*k(1) + k(2)*k(2) + k(3)*k(3));
  if( absk - E < small) return 0.0;
  else
  {
    double temp = pow(absk - E,2);
    h = temp/(temp + pow(M1,2));
    return h;
  }
}
//--------------------------------------------------------------------------------------------------
// Gradient of log(hminus(k));
 DESreal4vector DESdeform::gradloghminus(const DESreal4vector k)
{
  DESreal4vector grad;
  double E = k(0);
  double absk = sqrt( k(1)*k(1) + k(2)*k(2) + k(3)*k(3));
  if( absk + E < small ) return zero4vector;
  else
  {
    double w = absk + E;
    DESreal4vector gradw(1.0, k(1)/absk, k(2)/absk, k(3)/absk);
    grad = (2.0*pow(M1,2)/w/(pow(w,2) + pow(M1,2)))*gradw;
    return grad;
  }
}
//--------------------------------------------------------------------------------------------------
// Gradient of log(hplus(k));
DESreal4vector DESdeform::gradloghplus(const DESreal4vector k)
{
  DESreal4vector grad;
  double E = k(0);
  double absk = sqrt( k(1)*k(1) + k(2)*k(2) + k(3)*k(3));
  if( absk - E < small) return zero4vector;
  else
  {
    double w = absk - E;
    DESreal4vector gradw(-1.0, k(1)/absk, k(2)/absk, k(3)/absk);
    grad = (2.0*pow(M1,2)/w/(pow(w,2) + pow(M1,2)))*gradw;
    return grad;
  }
}
//--------------------------------------------------------------------------------------------------
// Function that smoothly vanishes for large k.
double DESdeform::g(const DESreal4vector k)
{
  double kEsq = k(0)*k(0) + k(1)*k(1) + k(2)*k(2) + k(3)*k(3);
  double M2sq = M2*M2;
  double g = gamma1 * M2sq/(kEsq + M2sq);
  return g;
}
//--------------------------------------------------------------------------------------------------
// Function emphasizing backward light cones.
double DESdeform::gplus(const DESreal4vector k)
{
  double omega = sqrt( k(1)*k(1) + k(2)*k(2) + k(3)*k(3) + M3*M3 );
  double E = k(0);
  double gplus = gamma2/(1.0 + pow(1.0 + E/omega, 2));
  return gplus;
}
//--------------------------------------------------------------------------------------------------
// Function emphasizing forward light cones.
double DESdeform::gminus(const DESreal4vector k)
{
  double omega = sqrt( k(1)*k(1) + k(2)*k(2) + k(3)*k(3) + M3*M3 );
  double E = k(0);
  double gminus = gamma2/(1.0 + pow(1.0 - E/omega, 2));
  return gminus;
}
//--------------------------------------------------------------------------------------------------
// Gradient of log(g(k));
// Note that grad propto the contravariant components of k since we differentiate 
// the euclidian square of k.
DESreal4vector DESdeform::gradlogg(const DESreal4vector k)
{
  DESreal4vector grad;
  double kEsq = k(0)*k(0) + k(1)*k(1) + k(2)*k(2) + k(3)*k(3);
  grad = (-2.0/( kEsq + M2*M2 ))*k;
  return grad;
}
//--------------------------------------------------------------------------------------------------
// Gradient of log(gplus(k));
DESreal4vector DESdeform::gradloggplus(const DESreal4vector k)
{
  DESreal4vector grad;
  double omega = sqrt( k(1)*k(1) + k(2)*k(2) + k(3)*k(3) + M3*M3 );
  double E = k(0);
  double ratio = E/omega;
  double factor = 2.0*(1 + ratio)/(1.0 + pow((1.0 + ratio), 2))/omega;
  grad.in(0) = - factor;
  for(int mu = 1; mu<=3; mu++) grad.in(mu) = factor*ratio/omega * k(mu);
  return grad;
}
//--------------------------------------------------------------------------------------------------
// Gradient of log(gminus(k));
DESreal4vector DESdeform::gradloggminus(const DESreal4vector k)
{
  DESreal4vector grad;
  double omega = sqrt( k(1)*k(1) + k(2)*k(2) + k(3)*k(3) + M3*M3 );
  double E = k(0);
  double ratio = E/omega;
  double factor = 2.0*(1 - ratio)/(1.0 + pow((1.0 - ratio), 2))/omega;
  grad.in(0) =  factor;
  for(int mu = 1; mu<=3; mu++) grad.in(mu) = - factor*ratio/omega * k(mu);
  return grad;
}
// End of implementation code for class DESdeform.

//--------------------------------------------------------------------------------------------------
//--------------------------------------------------------------------------------------------------
// Compute complex lc for a point ell, also compute jacobian;
// This is a version for testing: we get just the deformation for index jIN, namely c_j;
// We can also get \tilde c_+ using jIN = -1 and we can get \tilde c_- using jIN = -2.
// For the moment, this is defined for indexA not equal to 1 or nverts -1.
// Also, this does not include lambda -- that is lambda is set to 1.
//
void DESdeform::deformtest(DESreal4vector& ellIN, DEScmplx4vector& lcOUT, complex<double>& jacobian,
                           int jIN)
{
  static const complex<double> ii(0.0,1.0);
  double cj;
  DESreal4vector kappaj;
  ell = ellIN;
  double x    = dot(ell - Q[indexAp1], Pbar) /PdotPbar;
  double xbar = dot(ell - Q[1], P) /PdotPbar;
  kappa = zero4vector;
  DESreal4vector gradlogcj;
  double gradkappa[4][4] = {{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},
                            {0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}};
  //
  double gval = g(ell);
  DESreal4vector gradloggval = gradlogg(ell);
  // =========== General case: 1 < A < N - 1 ===========
  if( (1 < indexA) && (indexA < nverts - 1) )
  {
    // Vertices on left:
    if( (! ellinCminus(nverts)) && (! ellinCplus(indexAp1)) )
    {
      for(int j = 2; j<= indexA - 1; j++)
      {
        if( (! ellinCminus(j-1)) && (! ellinCplus(j+1)) )
        {
          cj  = hminus(ell - Q[j-1])*hplus(ell - Q[j+1]);
          cj *= hminus(ell - Q[nverts]);
          cj *= hplus(ell - Q[indexAp1]);
          cj *= gval;
          //
          if(j != jIN) cj = 0.0; // We demand that j = jIN to count this.
          if(j == jIN)
          {
            cout << "cj = " <<  cj << endl;
            cout << "h- = " << hminus(ell - Q[j-1]);
            cout << " h+ = " << hplus(ell - Q[j+1]);
            cout << " h- = " << hminus(ell - Q[nverts]);
            cout << " h+ = " << hplus(ell - Q[indexAp1]);
            cout << " g = " << gval << endl;
          }
          //
          kappaj = - cj * (ell - Q[j]);
          kappa = kappa + kappaj;
          //
          gradlogcj = gradloghminus(ell - Q[j-1]);
          gradlogcj = gradlogcj + gradloghplus(ell - Q[j+1]);
          gradlogcj = gradlogcj + gradloghminus(ell - Q[nverts]);
          gradlogcj = gradlogcj + gradloghplus(ell - Q[indexAp1]);
          gradlogcj = gradlogcj + gradloggval;
          for(int mu = 0; mu <= 3; mu++)
          {
            gradkappa[mu][mu] -= cj;
            for(int nu = 0; nu <= 3; nu++)
            {
              gradkappa[mu][nu] -= (ell(mu) - Q[j](mu))*cj*gradlogcj(nu);
            }
          }
        } // Close if( (! ellinCminus(j-1)) && (! ellinCplus(j+1)) )
      } // Close for(int j = 2; j<= indexA - 1; j++)
    } // Close if( (! ellinCminus(nverts)) && (! ellinCplus(indexAp1)) )
    // Vertices on right:
    if( (! ellinCminus(indexA)) && (! ellinCplus(1)) )
    {
      for(int j = indexA + 2; j<= nverts - 1; j++)
      {
        if( (! ellinCminus(j-1)) && (! ellinCplus(j+1)) )
        {
          cj  = hminus(ell - Q[j-1])*hplus(ell - Q[j+1]);
          cj *= hminus(ell - Q[indexA]);
          cj *= hplus(ell - Q[1]);
          cj *= gval;
          //
          if(j != jIN) cj = 0.0; // We demand that j = jIN to count this.
          if(j == jIN)
          {
            cout << "cj = " <<  cj << endl;
            cout << "h- = " << hminus(ell - Q[j-1]);
            cout << " h+ = " << hplus(ell - Q[j+1]);
            cout << " h- = " << hminus(ell - Q[indexA]);
            cout << " h+ = " << hplus(ell - Q[1]);
            cout << " g = " << gval << endl;
          }          
          //
          kappaj = - cj * (ell - Q[j]);
          kappa = kappa + kappaj;
          //
          gradlogcj = gradloghminus(ell - Q[j-1]);
          gradlogcj = gradlogcj + gradloghplus(ell - Q[j+1]);
          gradlogcj = gradlogcj + gradloghminus(ell - Q[indexA]);
          gradlogcj = gradlogcj + gradloghplus(ell - Q[1]);
          gradlogcj = gradlogcj + gradloggval;
          for(int mu = 0; mu <= 3; mu++)
          {
            gradkappa[mu][mu] -= cj;
            for(int nu = 0; nu <= 3; nu++)
            {
              gradkappa[mu][nu] -= (ell(mu) - Q[j](mu))*cj*gradlogcj(nu);
            }
          }
        } // Close if( (! ellinCminus(j-1)) && (! ellinCplus(j+1)) )
      } // Close for(int j = indexA + 2; j<= nverts - 1; j++)
    } // Close if( (! ellinCminus(indexA)) && (! ellinCplus(1)) )
    // Vertex at j = 1:
    if( (! ellinCplus(2)) && (! ellinCminus(nverts-1)) && (! ellinCplus(indexAp1)) )
    {
      cj  = hplus(ell - Q[2]);
      cj *= hminus(ell - Q[nverts-1]);
      cj *= hplus(ell - Q[indexAp1]);
      cj *= gval;
      //
      if(1 != jIN) cj = 0.0; // We demand that j = jIN to count this.
      if(1 == jIN)
      {
        cout << "cj = " <<  cj << endl;
        cout << " h+ = " << hplus(ell - Q[2]);
        cout << " h- = " << hminus(ell - Q[nverts-1]);
        cout << " h+ = " << hplus(ell - Q[indexAp1]);
        cout << " g = " << gval << endl;
      }
      //
      kappaj = - cj * (ell - Q[1]);
      kappa = kappa + kappaj;
      //
      gradlogcj = gradloghplus(ell - Q[2]);
      gradlogcj = gradlogcj + gradloghminus(ell - Q[nverts-1]);
      gradlogcj = gradlogcj + gradloghplus(ell - Q[indexAp1]);
      gradlogcj = gradlogcj + gradloggval;
      for(int mu = 0; mu <= 3; mu++)
      {
        gradkappa[mu][mu] -= cj;
        for(int nu = 0; nu <= 3; nu++)
        {
          gradkappa[mu][nu] -= (ell(mu) - Q[1](mu))*cj*gradlogcj(nu);
        }
      }      
    } // Close if( (! ellinCplus(2))  ... )
    // Vertex at j = A:
    if( (! ellinCminus(indexA-1)) && (! ellinCplus(indexA+2)) && (! ellinCminus(nverts)) )
    {
      cj  = hminus(ell - Q[indexA-1]);
      cj *= hplus(ell - Q[indexA+2]);
      cj *= hminus(ell - Q[nverts]);
      cj *= gval;
      //
      if(indexA != jIN) cj = 0.0; // We demand that j = jIN to count this.
      if(indexA == jIN)
      {
        cout << "cj = " <<  cj << endl;
        cout << " h- = " << hminus(ell - Q[indexA-1]);
        cout << " h+ = " << hplus(ell - Q[indexA+2]);
        cout << " h- = " << hminus(ell - Q[nverts]);
        cout << " g = " << gval << endl;
      }      
      //      
      kappaj = - cj * (ell - Q[indexA]);
      kappa = kappa + kappaj;
      //
      gradlogcj = gradloghminus(ell - Q[indexA-1]);
      gradlogcj = gradlogcj + gradloghplus(ell - Q[indexA+2]);
      gradlogcj = gradlogcj + gradloghminus(ell- Q[nverts]);
      gradlogcj = gradlogcj + gradloggval;
      for(int mu = 0; mu <= 3; mu++)
      {
        gradkappa[mu][mu] -= cj;
        for(int nu = 0; nu <= 3; nu++)
        {
          gradkappa[mu][nu] -= (ell(mu) - Q[indexA](mu))*cj*gradlogcj(nu);
        }
      }      
    } // Close if( (! ellinCminus(indexA-1))  ... )
    // Vertex at j = A + 1:
    if( (! ellinCplus(indexA+2)) && (! ellinCminus(indexA-1)) && (! ellinCplus(1)) )
    {
      cj  = hplus(ell - Q[indexA+2]);
      cj *= hminus(ell - Q[indexA-1]);
      cj *= hplus(ell - Q[1]);
      cj *= gval;
      //
      if(indexA+1 != jIN) cj = 0.0; // We demand that j = jIN to count this.
      if(indexA+1 == jIN)
      {
        cout << "cj = " <<  cj << endl;
        cout << " h+ = " << hplus(ell - Q[indexA+2]);
        cout << " h- = " << hminus(ell - Q[indexA-1]);
        cout << " h+ = " << hplus(ell - Q[1]);
        cout << " g = " << gval << endl;
      }            
      //      
      kappaj = - cj * (ell - Q[indexAp1]);
      kappa = kappa + kappaj;
      //
      gradlogcj = gradloghplus(ell - Q[indexA+2]);
      gradlogcj = gradlogcj + gradloghminus(ell - Q[indexA-1]);
      gradlogcj = gradlogcj + gradloghplus(ell- Q[1]);
      gradlogcj = gradlogcj + gradloggval;
      for(int mu = 0; mu <= 3; mu++)
      {
        gradkappa[mu][mu] -= cj;
        for(int nu = 0; nu <= 3; nu++)
        {
          gradkappa[mu][nu] -= (ell(mu) - Q[indexAp1](mu))*cj*gradlogcj(nu);
        }
      }      
    } // Close if( (! ellinCplus(indexA+2)  ... )
    // Vertex at j = N:
    if( (! ellinCminus(nverts-1)) && (! ellinCplus(2)) && (! ellinCminus(indexA)) )
    {
      cj  = hminus(ell - Q[nverts-1]);
      cj *= hplus(ell - Q[2]);
      cj *= hminus(ell - Q[indexA]);
      cj *= gval;
      //
      if(nverts != jIN) cj = 0.0; // We demand that j = jIN to count this.
      if(nverts == jIN)
      {
        cout << "cj = " <<  cj << endl;
        cout << " h- = " << hminus(ell - Q[nverts-1]);
        cout << " h+ = " << hplus(ell - Q[2]);
        cout << " h- = " << hminus(ell - Q[indexA]);
        cout << " g = " << gval << endl;
      }                  
      //      
      kappaj = - cj * (ell - Q[nverts]);
      kappa = kappa + kappaj;
      //
      gradlogcj = gradloghminus(ell - Q[nverts-1]);
      gradlogcj = gradlogcj + gradloghplus(ell - Q[2]);
      gradlogcj = gradlogcj + gradloghminus(ell- Q[indexA]);
      gradlogcj = gradlogcj + gradloggval;
      for(int mu = 0; mu <= 3; mu++)
      {
        gradkappa[mu][mu] -= cj;
        for(int nu = 0; nu <= 3; nu++)
        {
          gradkappa[mu][nu] -= (ell(mu) - Q[nverts](mu))*cj*gradlogcj(nu);
        }
      }      
    } // Close if( (! ellinCminus(nverts-1))  ... )
    // Now \tilde c_+ contribution:
    if( (x + xbar > 0) && (! ellinCminus(indexA)) && (! ellinCminus(nverts)) )
    {
      cj = hminus(ell - Q[indexA]) * hminus(ell - Q[nverts]);
      cj *= (x + xbar);
      cj *= gminus(ell);
      //
      if(-1 != jIN) cj = 0.0; // We demand that j = jIN to count this.
      //      
      kappaj = cj * (P + Pbar);
      kappa = kappa + kappaj;
      //
      gradlogcj = gradloghminus(ell - Q[indexA]);
      gradlogcj = gradlogcj + gradloghminus(ell - Q[nverts]);
      gradlogcj = gradlogcj + (1.0/(x + xbar)/PdotPbar)*(Pc + Pbarc);
      gradlogcj = gradlogcj + gradloggminus(ell);
      for(int mu = 0; mu <= 3; mu++)
      {
        for(int nu = 0; nu <= 3; nu++)
        {
          gradkappa[mu][nu] += (P(mu) + Pbar(mu))*cj*gradlogcj(nu);
        }
      }      
    } // Close if( (x + xbar > 0) ... )
    // Now \tilde c_- contribution:
    if( (x + xbar < 0) && (! ellinCplus(1)) && (! ellinCplus(indexAp1)) )
    {
      cj = hplus(ell - Q[1]) * hplus(ell - Q[indexAp1]);
      cj *= -(x + xbar);
      cj *= gplus(ell);
      //
      if(-2 != jIN) cj = 0.0; // We demand that j = jIN to count this.
      //      
      kappaj = - cj * (P + Pbar);
      kappa = kappa + kappaj;
      //
      gradlogcj = gradloghplus(ell - Q[1]);
      gradlogcj = gradlogcj + gradloghplus(ell - Q[indexAp1]);
      gradlogcj = gradlogcj + (1.0/(x + xbar)/PdotPbar)*(Pc + Pbarc);
      gradlogcj = gradlogcj + gradloggplus(ell);
      for(int mu = 0; mu <= 3; mu++)
      {
        for(int nu = 0; nu <= 3; nu++)
        {
          gradkappa[mu][nu] -= (P(mu) + Pbar(mu))*cj*gradlogcj(nu);
        }
      }      
    } // Close if( (x + xbar < 0) ... )
  } // Close if( (1 < indexA) && (indexA < nverts - 1) )
  // ====================== A = 1 ======================
  else if( indexA == 1 )
  {
    // Vertices on right:
    if( ! ellinCplus(1) )
    {
      for(int j = 3; j<= nverts - 1; j++)
      {
        if( (! ellinCminus(j-1)) && (! ellinCplus(j+1)) )
        {
          cj  = hminus(ell - Q[j-1])*hplus(ell - Q[j+1]);
          cj *= hplus(ell - Q[1]);
          cj *= gval;
          kappaj = - cj * (ell - Q[j]);
          kappa = kappa + kappaj;
          //
          gradlogcj = gradloghminus(ell - Q[j-1]);
          gradlogcj = gradlogcj + gradloghplus(ell - Q[j+1]);
          gradlogcj = gradlogcj + gradloghplus(ell - Q[1]);
          gradlogcj = gradlogcj + gradloggval;
          for(int mu = 0; mu <= 3; mu++)
          {
            gradkappa[mu][mu] -= cj;
            for(int nu = 0; nu <= 3; nu++)
            {
              gradkappa[mu][nu] -= (ell(mu) - Q[j](mu))*cj*gradlogcj(nu);
            }
          }
        } // Close if( (! ellinCminus(j-1)) && (! ellinCplus(j+1)) )
      } // Close for(int j = 3; j<= nverts - 1; j++)
    } // Close if( ! ellinCplus(1) )
    // Vertex at j = 1:
    if( (! ellinCminus(nverts-1)) && (! ellinCplus(3)) )
    {
      cj  = hminus(ell - Q[nverts-1]);
      cj *= hplus(ell - Q[3]);
      cj *= gval;
      kappaj = - cj * (ell - Q[1]);
      kappa = kappa + kappaj;
      //
      gradlogcj = gradloghminus(ell - Q[nverts-1]);
      gradlogcj = gradlogcj + gradloghplus(ell- Q[3]);
      gradlogcj = gradlogcj + gradloggval;
      for(int mu = 0; mu <= 3; mu++)
      {
        gradkappa[mu][mu] -= cj;
        for(int nu = 0; nu <= 3; nu++)
        {
          gradkappa[mu][nu] -= (ell(mu) - Q[1](mu))*cj*gradlogcj(nu);
        }
      }      
    } // Close if( (! ellinCminus(nverts-1))  ... )
    // Vertex at j = 2:
    if( (! ellinCplus(3)) && (! ellinCplus(1)) )
    {
      cj  = hplus(ell - Q[3]);
      cj *= hplus(ell - Q[1]);
      cj *= gval;
      kappaj = - cj * (ell - Q[2]);
      kappa = kappa + kappaj;
      //
      gradlogcj = gradloghplus(ell - Q[3]);
      gradlogcj = gradlogcj + gradloghplus(ell- Q[1]);
      gradlogcj = gradlogcj + gradloggval;
      for(int mu = 0; mu <= 3; mu++)
      {
        gradkappa[mu][mu] -= cj;
        for(int nu = 0; nu <= 3; nu++)
        {
          gradkappa[mu][nu] -= (ell(mu) - Q[2](mu))*cj*gradlogcj(nu);
        }
      }      
    } // Close if( (! ellinCplus(indexA+2)  ... )
    // Vertex at j = N:
    if( (! ellinCminus(nverts-1)) && (! ellinCminus(1)) )
    {
      cj  = hminus(ell - Q[nverts-1]);
      cj *= hminus(ell - Q[1]);
      cj *= gval;
      kappaj = - cj * (ell - Q[nverts]);
      kappa = kappa + kappaj;
      //
      gradlogcj = gradloghminus(ell - Q[nverts-1]);
      gradlogcj = gradlogcj + gradloghminus(ell- Q[1]);
      gradlogcj = gradlogcj + gradloggval;
      for(int mu = 0; mu <= 3; mu++)
      {
        gradkappa[mu][mu] -= cj;
        for(int nu = 0; nu <= 3; nu++)
        {
          gradkappa[mu][nu] -= (ell(mu) - Q[nverts](mu))*cj*gradlogcj(nu);
        }
      }      
    } // Close if( (! ellinCminus(nverts-1))  ... )
    // Now \tilde c_+ contribution:
    if( (x + xbar > 0) && (! ellinCminus(nverts)) )
    {
      cj = hminus(ell - Q[nverts]);
      cj *= (x + xbar);
      cj *= gminus(ell);
      kappaj = cj * (P + Pbar);
      kappa = kappa + kappaj;
      //
      gradlogcj = gradloghminus(ell - Q[nverts]);
      gradlogcj = gradlogcj + (1.0/(x + xbar)/PdotPbar)*(Pc + Pbarc);
      gradlogcj = gradlogcj + gradloggminus(ell);
      for(int mu = 0; mu <= 3; mu++)
      {
        for(int nu = 0; nu <= 3; nu++)
        {
          gradkappa[mu][nu] += (P(mu) + Pbar(mu))*cj*gradlogcj(nu);
        }
      }
    } // Close if( (x + xbar > 0) ... )
    // Now \tilde c_- contribution:
    if( (x + xbar < 0) && (! ellinCplus(2)) )
    {
      cj = hplus(ell - Q[2]);
      cj *= -(x + xbar);
      cj *= gplus(ell);
      kappaj = - cj * (P + Pbar);
      kappa = kappa + kappaj;
      //
      gradlogcj = gradloghplus(ell - Q[2]);
      gradlogcj = gradlogcj + (1.0/(x + xbar)/PdotPbar)*(Pc + Pbarc);
      gradlogcj = gradlogcj + gradloggplus(ell);
      for(int mu = 0; mu <= 3; mu++)
      {
        for(int nu = 0; nu <= 3; nu++)
        {
          gradkappa[mu][nu] -= (P(mu) + Pbar(mu))*cj*gradlogcj(nu);
        }
      }      
    } // Close if( (x + xbar < 0) ... )
  } // Close else if( indexA == 1 )
  // ==================== A = N - 1 ====================
  else if( indexA == nverts - 1 )
  {
    // Vertices on left:
    if( (! ellinCminus(nverts)) && (! ellinCplus(nverts)) )
    {
      for(int j = 2; j<= nverts - 2; j++)
      {
        if( (! ellinCminus(j-1)) && (! ellinCplus(j+1)) )
        {
          cj  = hminus(ell - Q[j-1])*hplus(ell - Q[j+1]);
          cj *= hminus(ell - Q[nverts]);
          cj *= hplus(ell - Q[nverts]);
          cj *= gval;
          kappaj = - cj * (ell - Q[j]);
          kappa = kappa + kappaj;
          //
          gradlogcj = gradloghminus(ell - Q[j-1]);
          gradlogcj = gradlogcj + gradloghplus(ell - Q[j+1]);
          gradlogcj = gradlogcj + gradloghminus(ell - Q[nverts]);
          gradlogcj = gradlogcj + gradloghplus(ell - Q[nverts]);
          gradlogcj = gradlogcj + gradloggval;
          for(int mu = 0; mu <= 3; mu++)
          {
            gradkappa[mu][mu] -= cj;
            for(int nu = 0; nu <= 3; nu++)
            {
              gradkappa[mu][nu] -= (ell(mu) - Q[j](mu))*cj*gradlogcj(nu);
            }
          }
        } // Close if( (! ellinCminus(j-1)) && (! ellinCplus(j+1)) )
      } // Close for(int j = 2; j<= indexA - 1; j++)
    } // Close if( (! ellinCminus(nverts)) && (! ellinCplus(nverts)) )
    // Vertex at j = 1:
    if( (! ellinCplus(2)) && (! ellinCplus(nverts)) )
    {
      cj  = hplus(ell - Q[2]);
      cj *= hplus(ell - Q[nverts]);
      cj *= gval;
      kappaj = - cj * (ell - Q[1]);
      kappa = kappa + kappaj;
      //
      gradlogcj = gradloghplus(ell - Q[2]);
      gradlogcj = gradlogcj + gradloghplus(ell - Q[nverts]);
      gradlogcj = gradlogcj + gradloggval;
      for(int mu = 0; mu <= 3; mu++)
      {
        gradkappa[mu][mu] -= cj;
        for(int nu = 0; nu <= 3; nu++)
        {
          gradkappa[mu][nu] -= (ell(mu) - Q[1](mu))*cj*gradlogcj(nu);
        }
      }      
    } // Close if( (! ellinCplus(2))  ... )
    // Vertex at j = N - 1:
    if( (! ellinCminus(nverts-2)) && (! ellinCminus(nverts)) )
    {
      cj  = hminus(ell - Q[nverts-2]);
      cj *= hminus(ell - Q[nverts]);
      cj *= gval;
      kappaj = - cj * (ell - Q[nverts-1]);
      kappa = kappa + kappaj;
      //
      gradlogcj = gradloghminus(ell - Q[nverts-2]);
      gradlogcj = gradlogcj + gradloghminus(ell- Q[nverts]);
      gradlogcj = gradlogcj + gradloggval;
      for(int mu = 0; mu <= 3; mu++)
      {
        gradkappa[mu][mu] -= cj;
        for(int nu = 0; nu <= 3; nu++)
        {
          gradkappa[mu][nu] -= (ell(mu) - Q[nverts-1](mu))*cj*gradlogcj(nu);
        }
      }      
    } // Close if( (! ellinCminus(nverts-2))  ... )
    // Vertex at j = N:
    if( (! ellinCplus(2)) && (! ellinCminus(nverts-2)) )
    {
      cj = hplus(ell - Q[2]);
      cj *= hminus(ell - Q[nverts-2]);
      cj *= gval;
      kappaj = - cj * (ell - Q[nverts]);
      kappa = kappa + kappaj;
      //
      gradlogcj = gradloghplus(ell - Q[2]);
      gradlogcj = gradlogcj + gradloghminus(ell- Q[nverts-2]);
      gradlogcj = gradlogcj + gradloggval;
      for(int mu = 0; mu <= 3; mu++)
      {
        gradkappa[mu][mu] -= cj;
        for(int nu = 0; nu <= 3; nu++)
        {
          gradkappa[mu][nu] -= (ell(mu) - Q[nverts](mu))*cj*gradlogcj(nu);
        }
      }      
    } // Close if( (! ellinCplus(2))  ... )
    // Now \tilde c_+ contribution:
    if( (x + xbar > 0) && (! ellinCminus(nverts-1)) )
    {
      cj = hminus(ell - Q[nverts-1]);
      cj *= (x + xbar);
      cj *= gminus(ell);
      kappaj = cj * (P + Pbar);
      kappa = kappa + kappaj;
      //
      gradlogcj = gradloghminus(ell - Q[nverts-1]);
      gradlogcj = gradlogcj + (1.0/(x + xbar)/PdotPbar)*(Pc + Pbarc);
      gradlogcj = gradlogcj + gradloggminus(ell);
      for(int mu = 0; mu <= 3; mu++)
      {
        for(int nu = 0; nu <= 3; nu++)
        {
          gradkappa[mu][nu] += (P(mu) + Pbar(mu))*cj*gradlogcj(nu);
        }
      }      
    } // Close if( (x + xbar > 0) ... )
    // Now \tilde c_- contribution:
    if( (x + xbar < 0) && (! ellinCplus(1)) )
    {
      cj = hplus(ell - Q[1]);
      cj *= -(x + xbar);
      cj *= gplus(ell);
      kappaj = - cj * (P + Pbar);
      kappa = kappa + kappaj;
      //
      gradlogcj = gradloghplus(ell - Q[1]);
      gradlogcj = gradlogcj + (1.0/(x + xbar)/PdotPbar)*(Pc + Pbarc);
      gradlogcj = gradlogcj + gradloggplus(ell);
      for(int mu = 0; mu <= 3; mu++)
      {
        for(int nu = 0; nu <= 3; nu++)
        {
          gradkappa[mu][nu] -= (P(mu) + Pbar(mu))*cj*gradlogcj(nu);
        }
      }      
    } // Close if( (x + xbar < 0) ... )
  } // Close else if( indexA == nverts - 1 )
  else 
  {
    cout << "Bad value for index A." << endl;
    exit(1);
  }
  // ============ End of calculation of kappa and grad kappa ============
  DEScmplx4vector kappac, ellc;
  for(int mu = 0; mu <= 3; mu++)  // We need this for the type conversion.
  {
    kappac.in(mu) = kappa(mu);
    ellc.in(mu) = ell(mu);
  }
  lc = ellc + ii * kappac;
  lcOUT = lc;
  complex<double> gradlc[maxsize][maxsize]; 
       // We need this dimension for Determinant(gradlc,dimension).
  for(int mu = 0; mu <= 3; mu++)
  {
    for(int nu = 0; nu <= 3; nu++)
    {
      if(mu == nu) gradlc[mu][nu] = 1.0 + ii * gradkappa[mu][nu];
      else         gradlc[mu][nu] = ii * gradkappa[mu][nu];
    }
  }
  int dimension = 4;
  jacobian = Determinant(gradlc,dimension);
  return;
}
//--------------------------------------------------------------------------------------------------
//==================================================================================================
// Implementation code for class DESreal4vector.

// Default constructor
DESreal4vector::DESreal4vector()
{
for(int mu = 0; mu <= 3; mu++)
  {
  vector[mu] = 0.0;
  }
}
//--------------------------------------------------------------------------------------------------
// Constructor to initialize to given value.
DESreal4vector::DESreal4vector(double v0, double v1, double v2, double v3)
{
  vector[0] = v0;
  vector[1] = v1;
  vector[2] = v2;
  vector[3] = v3;
}
//--------------------------------------------------------------------------------------------------
// v.in(mu) for input, v.in(mu) = r
double&  DESreal4vector::in(const int& mu)
{
return vector[mu];
}
//--------------------------------------------------------------------------------------------------
// v(mu) for output,  r = v(mu)
const double&  DESreal4vector::operator()(const int& mu) const
{
return vector[mu];
}
//--------------------------------------------------------------------------------------------------
// Sum of two vectors: v3 = v1 + v2
DESreal4vector DESreal4vector::operator+(const DESreal4vector& v2) const
{
DESreal4vector sum;
for(int mu = 0; mu <= 3; mu++)
  {
  sum.vector[mu] =  vector[mu] + v2.vector[mu];
  }
return sum;
}
//--------------------------------------------------------------------------------------------------
// Difference of two vectors: v3 = v1 - v2
DESreal4vector DESreal4vector::operator-(const DESreal4vector& v2) const
{
DESreal4vector difference;
for(int mu = 0; mu <= 3; mu++)
  {
  difference.vector[mu] =  vector[mu] - v2.vector[mu];
  }
return difference;
}
//--------------------------------------------------------------------------------------------------
// Vector times scalar: v2 = v1*c
DESreal4vector  DESreal4vector::operator*(const double& a) const
{
DESreal4vector product;
for(int mu = 0; mu <= 3; mu++)
  {
  product.vector[mu] = a * vector[mu];
  }
return product;
}
//--------------------------------------------------------------------------------------------------
// -1 times vector: v2 = - v1
DESreal4vector  DESreal4vector::operator-() const
{
DESreal4vector negative;
for(int mu = 0; mu <= 3; mu++)
  {
  negative.vector[mu] = - vector[mu];
  }
return negative;
}

// Scalar times vector: v2 = c*v1
DESreal4vector operator*(const double& a, const DESreal4vector& v)
{
DESreal4vector product;
for(int mu = 0; mu <= 3; mu++)
  {
  product.in(mu) = a * v(mu);  // We use public operator here.
  }
return product;
}
//--------------------------------------------------------------------------------------------------
// Dot product of two 4-vectors: dot(v1,v2)
double dot(const DESreal4vector& v1, const DESreal4vector& v2)
{
double product;
product = v1(0)*v2(0);
for(int mu = 1; mu <= 3; mu++)
  {
  product = product - v1(mu)*v2(mu); // We use public operator here.
  }
return product;
}
//--------------------------------------------------------------------------------------------------
// Square of one 4-vector: square(v) = dot(v,v) 
double square(const DESreal4vector& v)
{
double product;
product = v(0)*v(0);
for(int mu = 1; mu <= 3; mu++)
  {
  product = product - v(mu)*v(mu); // We use public operator here.
  }
return product;
}
// End of implementation code for DESreal4vector.
//--------------------------------------------------------------------------------------------------
// Boost a vector v with boost that transforms vold to vnew.
DESreal4vector boost(const DESreal4vector& vnew, const DESreal4vector& vold, const DESreal4vector& v)
{
  DESreal4vector vout,vsum;
  vout = v;
  vsum = vold + vnew;
  double voldsq = square(vold);
  double vsumsq = square(vsum);
  // Some checks that one could omit.
  if (abs(voldsq) < 10e-10*abs(vsumsq))
  {
    cout << "Oops, vold should not be lightlike in function boost." << endl;
    exit(1);
  }
  if (abs(vsumsq) < 10e-10*abs(voldsq))
  {
    cout << "Oops, vold + vnew should not be lightlike in function boost." << endl;
    exit(1);
  }
  if (abs(voldsq - square(vnew)) > 10e-8*abs(vsumsq))
  {
    cout << "Oops, vold squared and vnew squared should be equal in function boost." << endl;
    exit(1);
  }
  // End of checks.
  vout = vout - (2.0/vsumsq)*vsum*dot(vsum,v);
  vout = vout + (2.0/voldsq)*vnew*dot(vold,v);
  return vout;
}
// End of function boost(vnew, vold, v).
//==================================================================================================
// Implementation code for class DEScmplx4vector.

// Default constructor
DEScmplx4vector::DEScmplx4vector()
{
for(int mu = 0; mu <= 3; mu++)
  {
  vector[mu] = 0.0;
  }
}
//--------------------------------------------------------------------------------------------------
// v.in(mu) for input, v.in(mu) = r
complex<double>&  DEScmplx4vector::in(const int& mu)
{
return vector[mu];
}
//--------------------------------------------------------------------------------------------------
// v(mu) for output,  r = v(mu)
const complex<double>&  DEScmplx4vector::operator()(const int& mu) const
{
return vector[mu];
}
//--------------------------------------------------------------------------------------------------
// Sum of two vectors: v3 = v1 + v2
DEScmplx4vector DEScmplx4vector::operator+(const DEScmplx4vector& v2) const
{
DEScmplx4vector sum;
for(int mu = 0; mu <= 3; mu++)
  {
  sum.vector[mu] =  vector[mu] + v2.vector[mu];
  }
return sum;
}
//--------------------------------------------------------------------------------------------------
// Difference of two vectors: v3 = v1 - v2
DEScmplx4vector DEScmplx4vector::operator-(const DEScmplx4vector& v2) const
{
DEScmplx4vector difference;
for(int mu = 0; mu <= 3; mu++)
  {
  difference.vector[mu] =  vector[mu] - v2.vector[mu];
  }
return difference;
}
//--------------------------------------------------------------------------------------------------
// Vector times scalar: v2 = v1*c
DEScmplx4vector  DEScmplx4vector::operator*(const complex<double>& a) const
{
DEScmplx4vector product;
for(int mu = 0; mu <= 3; mu++)
  {
  product.vector[mu] = a * vector[mu];
  }
return product;
}
//--------------------------------------------------------------------------------------------------
// -1 times vector: v2 = - v1
DEScmplx4vector  DEScmplx4vector::operator-() const
{
DEScmplx4vector negative;
for(int mu = 0; mu <= 3; mu++)
  {
  negative.vector[mu] = - vector[mu];
  }
return negative;
}
//--------------------------------------------------------------------------------------------------
// Scalar times vector: v2 = c*v1
DEScmplx4vector operator*(const complex<double>& a, const DEScmplx4vector& v)
{
DEScmplx4vector product;
for(int mu = 0; mu <= 3; mu++)
  {
  product.in(mu) = a * v(mu);  // We use public operator here.
  }
return product;
}
//--------------------------------------------------------------------------------------------------
// Dot product of two 4-vectors: dot(v1,v2)
complex<double> dot(const DEScmplx4vector& v1, const DEScmplx4vector& v2)
{
complex<double> product;
product = v1(0)*v2(0);
for(int mu = 1; mu <= 3; mu++)
  {
  product = product - v1(mu)*v2(mu); // We use public operator here.
  }
return product;
}
//--------------------------------------------------------------------------------------------------
// Square of one 4-vector: square(v) = dot(v,v) 
complex<double> square(const DEScmplx4vector& v)
{
complex<double> product;
product = v(0)*v(0);
for(int mu = 1; mu <= 3; mu++)
  {
  product = product - v(mu)*v(mu); // We use public operator here.
  }
return product;
}

// End of implementation code for DEScmplx4vector.
//==================================================================================================
// Implementation code for class DESindex.
// Default constructor.
DESindex::DESindex()
{
  nverts = 1000000;  // This default value is deliberately absurd.
}
//--------------------------------------------------------------------------------------------------
// Constructor.
DESindex::DESindex(int nvertsIN)
{
  nverts = nvertsIN;
}
//--------------------------------------------------------------------------------------------------
int DESindex::getnverts() const
{
  return nverts;
}
//--------------------------------------------------------------------------------------------------
// Operator to return i mod nverts.
int DESindex::operator()(int i) const
{
  int theindex;
  if(i > 0)
  {
    if(i <= nverts)
    {
      theindex = i;   // 0 < i <= nverts, the most common case
    }
    else if(i <= 2*nverts)
    {
      theindex = i - nverts;  // nverts < i <= 2*nverts
    }
    else
    {
      theindex = (i-1)%nverts + 1; // 2*nverts < i , use the general formula
    }
  }
  else if(i > -nverts)
  {
    theindex = i + nverts;  // -nverts < i <= 0
  }
  else
  {
    theindex = (i-1)%nverts + 1; // i <= - nverts, use the general formula
  }
  return theindex;
}
//==================================================================================================
//==================================================================================================
//==================================================================================================
// Implementation code for class DESrandom.
// input : seed (0<=seed<259200)                                            
// output: random numbers (between 2^(-32) and 1-2^(-32) ) 
// Code by D. Sung and DES.

DESrandom::DESrandom(/*in*/int seed) //constructor(initializes seed)   
{
  reSetSeed(seed);
}      
//--------------------------------------------------------------------------------------------------
DESrandom::DESrandom()  //default constructor. (puts seed as 0)
{
  reSetSeed(0);
}      
//--------------------------------------------------------------------------------------------------
bool DESrandom::reSetSeed(/*in*/int seed) 
              //preCondition 0<=seed<259200
              //postCond. : seed is reset, return 0 if seed is invalid.
{
  lbit=31;
  fac=pow(2.0,-lbit);

  int i, ifac1, isk, idi, i1s;
  int im=259200;
  int ia=7141;
  int ic=54773;

  assert(seed >=0);
  assert(seed < im);

  for(int i=0; i<10 ; i++)
    seed=(seed*ia + ic)%im;
  ifac1=((int(pow(2.0,lbit-1))-1)*2+1)/im;
  for(int i=1;i<250;i++)
    {
    seed=(seed*ia+ic)%im;
    ir[i]=ifac1*seed;
    }

  ir[0]=1;
  i=0; 
  isk=250/lbit;
  idi=1;
  i1s=1;

  for(int lb=1;lb<lbit;lb++)
  {
    i=i+isk;
    idi=idi*2;
    i1s=i1s+idi;
    ir[i]=ir[i]|idi;        //bitwise operator. | : or
    ir[i]=ir[i]&i1s;        // & : and
  }

  newran();
  return true;
}
//--------------------------------------------------------------------------------------------------
const double& DESrandom::random()   //return random value.
{
  if(j>=249)
    newran();
  j=j+1;
  return rr[j];
}
//--------------------------------------------------------------------------------------------------   
void DESrandom::newran()
{
  for(int n=0;n<103;n++)
  {
  ir[n]=ir[n+147]^ir[n];
  rr[n]=fac*(double(ir[n])+0.5);
  }
  for(int n=103;n<250;n++)
  {
  ir[n]=ir[n-103]^ir[n];       //^ : bitwise op. exclusive or.
  rr[n]=fac*(double(ir[n])+0.5);
  }

j=-1;
}
// End of implementation code for DESrandom.
//==================================================================================================
// Dummy numerator function N(ell,Q).
complex<double> numerator(const DEScmplx4vector& lvec, const DEScmplx4vector Q[])
{
complex<double> thenumerator;
thenumerator = 1.0;
return thenumerator;
}
//==================================================================================================
// Implementation code for class FermionLoop. Authors: W. Gong & D. E. Soper
//
FermionLoop::FermionLoop(int nvertsIN) // Constructor. Initializes gamma.
{
int i,j;
//
nverts = nvertsIN;
// Initialize polarization vectors to zero. They will be changed later.
for(i = 0; i < maxsize; i++)
  {
  for(j = 0; j <= 3; j++)
    {
    epsln[i].in(j) = 0;
    }
  }
}
//--------------------------------------------------------------------------------------------------
void FermionLoop::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".
{
epsln[nv] = epsilon;
}
//--------------------------------------------------------------------------------------------------
void FermionLoop::CalEpsilonNull(const int nv, const int hIN, const DESreal4vector& pIN)
               // 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.
               // This version has null plane gauge.
{
DESreal4vector p;
DEScmplx4vector epsilon;
static const complex<double>ii(0.0,1.0);
complex<double> e1,e2,eps0;
int h;

p = pIN;
if(p(0) < 0) p = - p; // This gives the physical momentum for the incoming particles.
h = - hIN;  // Gives epsilon^* for outgoing particles and - h convention for incoming particles.
//
if(h == 1)
  {
  e1 =  1.0/sqrt(2.0);   // plus helicity
  e2 =   ii/sqrt(2.0);   // plus helicity
  }
else if(h == -1)
  {
  e1 = 1.0/sqrt(2.0);   // minus helicity
  e2 = - ii/sqrt(2.0);  // minus helicity
  }
else if (h == 0) // return \epsilon^\mu \propto p^\mu.
  {
  epsilon.in(0) = 1.0;
  epsilon.in(1) = p(1)/p(0);
  epsilon.in(2) = p(2)/p(0);
  epsilon.in(3) = p(3)/p(0);
  epsln[nv] = epsilon;
  return;
  }
else
 {
 cout << "Wrong helicity input!" << endl;
 exit(1);
 }
if(p(3) > 0.0) // Use epsilon^+ = 0 vectors.
  {
  eps0 = (p(1)*e1 + p(2)*e2)/(p(0) + p(3));
  epsilon.in(0) = eps0;
  epsilon.in(1) = e1;
  epsilon.in(2) = e2;
  epsilon.in(3) = - eps0;
  }
else                // Use epsilon^+ = 0 vectors.
  {
  e2 = - e2;        // eps_T(h) -> - eps_T(-h).
  eps0 = (p(1)*e1 + p(2)*e2)/(p(0) - p(3));
  epsilon.in(0) = eps0;
  epsilon.in(1) = e1;
  epsilon.in(2) = e2;
  epsilon.in(3) = eps0;
  }
epsln[nv] = epsilon;
}
//--------------------------------------------------------------------------------------------------
void FermionLoop::CalEpsilonCoulomb(const int nv, const int hIN, const DESreal4vector& pIN)
               // 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.
               // This version has Coulomb gauge.
{
DESreal4vector p;
DEScmplx4vector epsilon;
double E;
double r;
double temp;
static const complex<double>ii(0.0,1.0);
complex<double> e1,e2;
int h;

p = pIN;
if(p(0) < 0) p = - p; // This gives the physical momentum for the incoming particles.
h = - hIN;  // This gives epsilon^* for outgoing particles and - h convention for incoming particles.

if(h == 1)
  {
  e1 = 1.0/sqrt(2.0);   // plus helicity
  e2 = ii/sqrt(2.0);   // plus helicity
  }
else if(h == -1)
  {
  e1 = 1.0/sqrt(2.0);   // minus helicity
  e2 = - ii/sqrt(2.0);  // minus helicity
  }
else if (h == 0) // return \epsilon^\mu \propto p^\mu.
  {
  epsilon.in(0) = 1.0;
  epsilon.in(1) = p(1)/p(0);
  epsilon.in(2) = p(2)/p(0);
  epsilon.in(3) = p(3)/p(0);
  epsln[nv] = epsilon;
  return;
  }
else
 {
 cout << "Wrong helicity input!" << endl;
 exit(1);
 }

epsilon.in(0) = 0.0;
epsilon.in(1) = e1;
epsilon.in(2) = e2;
epsilon.in(3) = 0.0;
E = sqrt( p(1)*p(1) + p(2)*p(2) + p(3)*p(3) );
r = sqrt( p(1)*p(1) + p(2)*p(2) );
if(r/E > 1.0e-10)
  {
  epsilon.in(1) = p(2)*(e1*p(2)-e2*p(1))/(r*r) + p(1)*p(3)*(e1*p(1)+e2*p(2))/(E*r*r);
  epsilon.in(2) = -p(1)*(e1*p(2)-e2*p(1))/(r*r) + p(2)*p(3)*(e1*p(1)+e2*p(2))/(E*r*r);
  temp = 1.0 - pow((p(3)/E), 2);
  epsilon.in(3) = -sqrt(temp)*(e1*p(1)+e2*p(2))/r;
  }
else if(p(3) < 0)        // Momentum along the negative z-axis.
  {
  epsilon.in(1) = - e1;  // eps_T(h) -> - eps_T(-h).
  epsilon.in(2) =   e2;
  }
epsln[nv] = epsilon;
}
//--------------------------------------------------------------------------------------------------
complex<double> FermionLoop::NumeratorV(const DEScmplx4vector& l, const DEScmplx4vector Q[])
{
static complex<double> ii(0.0,1.0);
complex<double> numerator;
int nv,i,j,mu;
complex<double> eslash[5][5];
complex<double> kslash[5][5];
complex<double> product[5][5];
complex<double> product1[5][5];
DEScmplx4vector k,eps;
//
// Compute the trace of the product of vertex factors and propagators:
// Trace{ epsilonslash[nverts]*kslash[nverts]*...*epsilonslash[1]*kslash[1] }
// 
// Start with the unit matrix.
for(i = 1; i <= 4; i++)
  {
  for(j = 1; j <= 4; j++)
    {
    product[i][j] = 0.0;
    }
  }
for(i = 1; i <= 4; i++) product[i][i] = 1.0;
// Now loop over the propagators and vertices
for(nv = 1; nv <= nverts; nv++)
  {
  // Construct kslash.
  for(mu = 0; mu <= 3; mu++) k.in(mu) = l(mu) - Q[nv](mu);
  kslash[1][3] = k(0) - k(3);
  kslash[1][4] = -k(1) + ii*k(2);
  kslash[2][3] = -k(1) - ii*k(2);
  kslash[2][4] = k(0) + k(3);
  //
  kslash[3][1] = k(0) + k(3);
  kslash[3][2] = k(1) - ii*k(2);
  kslash[4][1] = k(1) + ii*k(2);
  kslash[4][2] = k(0) - k(3);
  // Left multiply by kslash.
  product1[1][3] = kslash[1][3]*product[3][3] + kslash[1][4]*product[4][3];
  product1[1][4] = kslash[1][3]*product[3][4] + kslash[1][4]*product[4][4];
  product1[2][3] = kslash[2][3]*product[3][3] + kslash[2][4]*product[4][3];
  product1[2][4] = kslash[2][3]*product[3][4] + kslash[2][4]*product[4][4];
  //
  product1[3][1] = kslash[3][1]*product[1][1] + kslash[3][2]*product[2][1];
  product1[3][2] = kslash[3][1]*product[1][2] + kslash[3][2]*product[2][2];
  product1[4][1] = kslash[4][1]*product[1][1] + kslash[4][2]*product[2][1];
  product1[4][2] = kslash[4][1]*product[1][2] + kslash[4][2]*product[2][2];
  // Construct epsilonslash.
  for(mu = 0; mu <= 3; mu++) eps.in(mu) = epsln[nv](mu);
  eslash[1][3] = eps(0) - eps(3);
  eslash[1][4] = -eps(1) + ii*eps(2);
  eslash[2][3] = -eps(1) - ii*eps(2);
  eslash[2][4] = eps(0) + eps(3);
  //
  eslash[3][1] = eps(0) + eps(3);
  eslash[3][2] = eps(1) - ii*eps(2);
  eslash[4][1] = eps(1) + ii*eps(2);
  eslash[4][2] = eps(0) - eps(3);
  // Left multiply by epsilonslash.
  product[1][1] = eslash[1][3]*product1[3][1] + eslash[1][4]*product1[4][1];
  product[1][2] = eslash[1][3]*product1[3][2] + eslash[1][4]*product1[4][2];
  product[2][1] = eslash[2][3]*product1[3][1] + eslash[2][4]*product1[4][1];
  product[2][2] = eslash[2][3]*product1[3][2] + eslash[2][4]*product1[4][2];
  //
  product[3][3] = eslash[3][1]*product1[1][3] + eslash[3][2]*product1[2][3];
  product[3][4] = eslash[3][1]*product1[1][4] + eslash[3][2]*product1[2][4];
  product[4][3] = eslash[4][1]*product1[1][3] + eslash[4][2]*product1[2][3];
  product[4][4] = eslash[4][1]*product1[1][4] + eslash[4][2]*product1[2][4];
  } // Close for(nv = 1; nv <= nverts; nv++).
// Now take the trace.
numerator = 0.0;
for(i = 1; i <= 4; i++) numerator = numerator + product[i][i];
// Return the value of numerator.
return numerator;
}
//--------------------------------------------------------------------------------------------------
complex<double> FermionLoop::NumeratorM(const DEScmplx4vector& l, const DEScmplx4vector Q[], const double m)
{
static const complex<double> ii(0.0,1.0);
complex<double> numerator;
int nv,i,j,mu;
complex<double> eslash[5][5];
complex<double> kslash[5][5];
complex<double> product[5][5];
complex<double> product1[5][5];
DEScmplx4vector k,eps;
//
// Compute the trace of the product of vertex factors and propagators, with mass:
// Trace{ epsilonslash[nverts]*(kslash[nverts] + m)*...*epsilonslash[1]*(kslash[1] + m) }
// 
// Start with the unit matrix.
for(i = 1; i <= 4; i++)
  {
  for(j = 1; j <= 4; j++)
    {
    product[i][j] = 0.0;
    }
  }
for(i = 1; i <= 4; i++) product[i][i] = 1.0;
// Now loop over the propagators and vertices
for(nv = 1; nv <= nverts; nv++)
  {
  // Construct kslash.
  for(mu = 0; mu <= 3; mu++) k.in(mu) = l(mu) - Q[nv](mu);
  kslash[1][3] = k(0) - k(3);
  kslash[1][4] = -k(1) + ii*k(2);
  kslash[2][3] = -k(1) - ii*k(2);
  kslash[2][4] = k(0) + k(3);
  //
  kslash[3][1] = k(0) + k(3);
  kslash[3][2] = k(1) - ii*k(2);
  kslash[4][1] = k(1) + ii*k(2);
  kslash[4][2] = k(0) - k(3);
  // Left multiply by kslash + m.
  product1[1][1] = m*product[1][1] + kslash[1][3]*product[3][1] + kslash[1][4]*product[4][1];
  product1[1][2] = m*product[1][2] + kslash[1][3]*product[3][2] + kslash[1][4]*product[4][2];
  product1[1][3] = m*product[1][3] + kslash[1][3]*product[3][3] + kslash[1][4]*product[4][3];
  product1[1][4] = m*product[1][4] + kslash[1][3]*product[3][4] + kslash[1][4]*product[4][4];
  //
  product1[2][1] = m*product[2][1] + kslash[2][3]*product[3][1] + kslash[2][4]*product[4][1];
  product1[2][2] = m*product[2][2] + kslash[2][3]*product[3][2] + kslash[2][4]*product[4][2];
  product1[2][3] = m*product[2][3] + kslash[2][3]*product[3][3] + kslash[2][4]*product[4][3];
  product1[2][4] = m*product[2][4] + kslash[2][3]*product[3][4] + kslash[2][4]*product[4][4];
  //
  product1[3][1] = kslash[3][1]*product[1][1] + kslash[3][2]*product[2][1] + m*product[3][1];
  product1[3][2] = kslash[3][1]*product[1][2] + kslash[3][2]*product[2][2] + m*product[3][2];
  product1[3][3] = kslash[3][1]*product[1][3] + kslash[3][2]*product[2][3] + m*product[3][3];
  product1[3][4] = kslash[3][1]*product[1][4] + kslash[3][2]*product[2][4] + m*product[3][4];
//
  product1[4][1] = kslash[4][1]*product[1][1] + kslash[4][2]*product[2][1] + m*product[4][1];
  product1[4][2] = kslash[4][1]*product[1][2] + kslash[4][2]*product[2][2] + m*product[4][2];
  product1[4][3] = kslash[4][1]*product[1][3] + kslash[4][2]*product[2][3] + m*product[4][3];
  product1[4][4] = kslash[4][1]*product[1][4] + kslash[4][2]*product[2][4] + m*product[4][4];
  //
  // Construct epsilonslash.
  for(mu = 0; mu <= 3; mu++) eps.in(mu) = epsln[nv](mu);
  eslash[1][3] = eps(0) - eps(3);
  eslash[1][4] = -eps(1) + ii*eps(2);
  eslash[2][3] = -eps(1) - ii*eps(2);
  eslash[2][4] = eps(0) + eps(3);
  //
  eslash[3][1] = eps(0) + eps(3);
  eslash[3][2] = eps(1) - ii*eps(2);
  eslash[4][1] = eps(1) + ii*eps(2);
  eslash[4][2] = eps(0) - eps(3);
  // Left multiply by epsilonslash.
  product[1][1] = eslash[1][3]*product1[3][1] + eslash[1][4]*product1[4][1];
  product[1][2] = eslash[1][3]*product1[3][2] + eslash[1][4]*product1[4][2];
  product[1][3] = eslash[1][3]*product1[3][3] + eslash[1][4]*product1[4][3];
  product[1][4] = eslash[1][3]*product1[3][4] + eslash[1][4]*product1[4][4];
  //
  product[2][1] = eslash[2][3]*product1[3][1] + eslash[2][4]*product1[4][1];
  product[2][2] = eslash[2][3]*product1[3][2] + eslash[2][4]*product1[4][2];
  product[2][3] = eslash[2][3]*product1[3][3] + eslash[2][4]*product1[4][3];
  product[2][4] = eslash[2][3]*product1[3][4] + eslash[2][4]*product1[4][4];
  //
  product[3][1] = eslash[3][1]*product1[1][1] + eslash[3][2]*product1[2][1];
  product[3][2] = eslash[3][1]*product1[1][2] + eslash[3][2]*product1[2][2];
  product[3][3] = eslash[3][1]*product1[1][3] + eslash[3][2]*product1[2][3];
  product[3][4] = eslash[3][1]*product1[1][4] + eslash[3][2]*product1[2][4];
  //
  product[4][1] = eslash[4][1]*product1[1][1] + eslash[4][2]*product1[2][1];
  product[4][2] = eslash[4][1]*product1[1][2] + eslash[4][2]*product1[2][2];
  product[4][3] = eslash[4][1]*product1[1][3] + eslash[4][2]*product1[2][3];
  product[4][4] = eslash[4][1]*product1[1][4] + eslash[4][2]*product1[2][4];
  } // Close for(nv = 1; nv <= nverts; nv++).
// Now take the trace.
numerator = 0.0;
for(i = 1; i <= 4; i++) numerator = numerator + product[i][i];
// Return the value of numerator.
return numerator;
}
//--------------------------------------------------------------------------------------------------
complex<double> FermionLoop::NumeratorUV(const DEScmplx4vector& l)
{
// This one if for the "extra" part of the UV subtraction.
// This is only for nverts = 4 with vector vertices.
int i;
complex<double> numerator;
numerator = 32.0;
for(i = 1; i <= 4; i++) numerator = numerator * dot(epsln[i],l);
// Return the value of numerator.
return numerator;
}
//--------------------------------------------------------------------------------------------------
complex<double> FermionLoop::AltNumeratorV(const DEScmplx4vector& l, const DEScmplx4vector Q[])
{
DEScmplx4vector v[2*maxsize];
int i;
complex<double> numerator;
//
// Compute the trace of the product of vertex factors and propagators:
// Trace{ epsilonslash[nverts]*kslash[nverts]*...*epsilonslash[1]*kslash[1] }
// In this case, we expand the trace in dot products of 4-vectors.
//
for(i = 1; i<= nverts; i++)
  {
  v[2*(nverts - i + 1)]     = l - Q[i];
  v[2*(nverts - i + 1) - 1] = epsln[i];
  }
numerator = DiracTrace(v,2*nverts); // Expands the trace in dot products of 4 vectors.
return numerator;
}
//--------------------------------------------------------------------------------------------------
complex<double> FermionLoop::NumeratorL(const DEScmplx4vector& l, const DEScmplx4vector Q[])
{
static const complex<double> ii(0.0,1.0);
complex<double> numerator;
int nv,i,j,mu;
complex<double> eslash[5][5];
complex<double> kslash[5][5];
complex<double> product[5][5];
complex<double> product1[5][5];
DEScmplx4vector k,eps;
//
// Compute the trace of the product of vertex factors and propagators:
// Trace{ epsilonslash[nverts]*kslash[nverts]*...*epsilonslash[1]*kslash[1] }
// 
// Start with the matrix diag(0,0,1,1).
for(i = 1; i <= 4; i++)
  {
  for(j = 1; j <= 4; j++)
    {
    product[i][j] = 0.0;
    }
  }
for(i = 3; i <= 4; i++) product[i][i] = 1.0;
// Now loop over the propagators and vertices
for(nv = 1; nv <= nverts; nv++)
  {
  // Construct kslash.
  for(mu = 0; mu <= 3; mu++) k.in(mu) = l(mu) - Q[nv](mu);
  kslash[1][3] = k(0) - k(3);
  kslash[1][4] = -k(1) + ii*k(2);
  kslash[2][3] = -k(1) - ii*k(2);
  kslash[2][4] = k(0) + k(3);
  // We don't need the following part of kslash.
  //kslash[3][1] = k(0) + k(3);
  //kslash[3][2] = k(1) - ii*k(2);
  //kslash[4][1] = k(1) + ii*k(2);
  //kslash[4][2] = k(0) - k(3);
  // Left multiply by kslash.
  product1[1][3] = kslash[1][3]*product[3][3] + kslash[1][4]*product[4][3];
  product1[1][4] = kslash[1][3]*product[3][4] + kslash[1][4]*product[4][4];
  product1[2][3] = kslash[2][3]*product[3][3] + kslash[2][4]*product[4][3];
  product1[2][4] = kslash[2][3]*product[3][4] + kslash[2][4]*product[4][4];
  //
  //product1[3][1] = kslash[3][1]*product[1][1] + kslash[3][2]*product[2][1];
  //product1[3][2] = kslash[3][1]*product[1][2] + kslash[3][2]*product[2][2];
  //product1[4][1] = kslash[4][1]*product[1][1] + kslash[4][2]*product[2][1];
  //product1[4][2] = kslash[4][1]*product[1][2] + kslash[4][2]*product[2][2];
  // Construct epsilonslash.
  for(mu = 0; mu <= 3; mu++) eps.in(mu) = epsln[nv](mu);
  // The following part of epsilon slash is removed by the left-handed projection.
  //eslash[1][3] = eps(0) - eps(3);
  //eslash[1][4] = -eps(1) + ii*eps(2);
  //eslash[2][3] = -eps(1) - ii*eps(2);
  //eslash[2][4] = eps(0) + eps(3);
  //
  eslash[3][1] = eps(0) + eps(3);
  eslash[3][2] = eps(1) - ii*eps(2);
  eslash[4][1] = eps(1) + ii*eps(2);
  eslash[4][2] = eps(0) - eps(3);
  // Left multiply by epsilonslash.
  //product[1][1] = eslash[1][3]*product1[3][1] + eslash[1][4]*product1[4][1];
  //product[1][2] = eslash[1][3]*product1[3][2] + eslash[1][4]*product1[4][2];
  //product[2][1] = eslash[2][3]*product1[3][1] + eslash[2][4]*product1[4][1];
  //product[2][2] = eslash[2][3]*product1[3][2] + eslash[2][4]*product1[4][2];
  //
  product[3][3] = eslash[3][1]*product1[1][3] + eslash[3][2]*product1[2][3];
  product[3][4] = eslash[3][1]*product1[1][4] + eslash[3][2]*product1[2][4];
  product[4][3] = eslash[4][1]*product1[1][3] + eslash[4][2]*product1[2][3];
  product[4][4] = eslash[4][1]*product1[1][4] + eslash[4][2]*product1[2][4];
  } // Close for(nv = 1; nv <= nverts; nv++).
// Now take the trace of the {3,4} diagonal elements.
numerator = 0.0;
for(i = 3; i <= 4; i++) numerator = numerator + product[i][i];
// Return the value of numerator.
return numerator;
}
//--------------------------------------------------------------------------------------------------
complex<double> FermionLoop::NumeratorR(const DEScmplx4vector& l, const DEScmplx4vector Q[])
{
static const complex<double> ii(0.0,1.0);
complex<double> numerator;
int nv,i,j,mu;
complex<double> eslash[5][5];
complex<double> kslash[5][5];
complex<double> product[5][5];
complex<double> product1[5][5];
DEScmplx4vector k,eps;
//
// Compute the trace of the product of vertex factors and propagators:
// Trace{ epsilonslash[nverts]*kslash[nverts]*...*epsilonslash[1]*kslash[1] }
// 
// Start with the matrix diag(1,1,0,0).
for(i = 1; i <= 4; i++)
  {
  for(j = 1; j <= 4; j++)
    {
    product[i][j] = 0.0;
    }
  }
for(i = 1; i <= 2; i++) product[i][i] = 1.0;
// Now loop over the propagators and vertices
for(nv = 1; nv <= nverts; nv++)
  {
  // Construct kslash.
  for(mu = 0; mu <= 3; mu++) k.in(mu) = l(mu) - Q[nv](mu);
  //  We don't need the following part of kslash.
  //kslash[1][3] = k(0) - k(3);
  //kslash[1][4] = -k(1) + ii*k(2);
  //kslash[2][3] = -k(1) - ii*k(2);
  //kslash[2][4] = k(0) + k(3);
  //
  kslash[3][1] = k(0) + k(3);
  kslash[3][2] = k(1) - ii*k(2);
  kslash[4][1] = k(1) + ii*k(2);
  kslash[4][2] = k(0) - k(3);
  // Left multiply by kslash.
  //product1[1][3] = kslash[1][3]*product[3][3] + kslash[1][4]*product[4][3];
  //product1[1][4] = kslash[1][3]*product[3][4] + kslash[1][4]*product[4][4];
  //product1[2][3] = kslash[2][3]*product[3][3] + kslash[2][4]*product[4][3];
  //product1[2][4] = kslash[2][3]*product[3][4] + kslash[2][4]*product[4][4];
  //
  product1[3][1] = kslash[3][1]*product[1][1] + kslash[3][2]*product[2][1];
  product1[3][2] = kslash[3][1]*product[1][2] + kslash[3][2]*product[2][2];
  product1[4][1] = kslash[4][1]*product[1][1] + kslash[4][2]*product[2][1];
  product1[4][2] = kslash[4][1]*product[1][2] + kslash[4][2]*product[2][2];
  // Construct epsilonslash.
  for(mu = 0; mu <= 3; mu++) eps.in(mu) = epsln[nv](mu);
  //
  eslash[1][3] = eps(0) - eps(3);
  eslash[1][4] = -eps(1) + ii*eps(2);
  eslash[2][3] = -eps(1) - ii*eps(2);
  eslash[2][4] = eps(0) + eps(3);
  // The following part of epsilon slash is removed by the right-handed projection.
  //eslash[3][1] = eps(0) + eps(3);
  //eslash[3][2] = eps(1) - ii*eps(2);
  //eslash[4][1] = eps(1) + ii*eps(2);
  //eslash[4][2] = eps(0) - eps(3);
  // Left multiply by epsilonslash.
  product[1][1] = eslash[1][3]*product1[3][1] + eslash[1][4]*product1[4][1];
  product[1][2] = eslash[1][3]*product1[3][2] + eslash[1][4]*product1[4][2];
  product[2][1] = eslash[2][3]*product1[3][1] + eslash[2][4]*product1[4][1];
  product[2][2] = eslash[2][3]*product1[3][2] + eslash[2][4]*product1[4][2];
  //
  //product[3][3] = eslash[3][1]*product1[1][3] + eslash[3][2]*product1[2][3];
  //product[3][4] = eslash[3][1]*product1[1][4] + eslash[3][2]*product1[2][4];
  //product[4][3] = eslash[4][1]*product1[1][3] + eslash[4][2]*product1[2][3];
  //product[4][4] = eslash[4][1]*product1[1][4] + eslash[4][2]*product1[2][4];
  } // Close for(nv = 1; nv <= nverts; nv++).
// Now take the trace of the {1,2} diagonal elements.
numerator = 0.0;
for(i = 1; i <= 2; i++) numerator = numerator + product[i][i];
// Return the value of numerator.
return numerator;
}
// End of implementation code for class FermionLoop.
//==================================================================================================
// Function to generate the trace of v1slash*v2slash*...*vnslash.
complex<double> DiracTrace(DEScmplx4vector v[2*maxsize], int n)
{
static const complex<double> ii(0.0,1.0);
complex<double> thetrace;
int i,j;
DEScmplx4vector u[10];
complex<double> sign;
//
if(n == 1) thetrace = 0.0;
else if(n == 2) thetrace = 4.0*dot(v[1],v[2]);
else
{
for(i = 1; i <= n - 2; i++) u[i] = v[i+2];
sign = 1.0;
thetrace = 0.0;
for(i = 2; i <= n; i++)
  {
  thetrace = thetrace + sign*dot(v[1],v[i])*DiracTrace(u, n-2);
  for(j = 1; j <= n - 3; j++) u[j] = u[j + 1];
  u[n-2] = v[i];
  sign = - sign;
  }
}
return thetrace;
}
//==================================================================================================
// Function getperms(nverts, permlist, indexA).
// Generates a list of permutations of the integers 1,2,...,nverts.
// The ith permutation is pi[j] = PermList[i][j].
// For j > nverts, getperms gives pi[j] = j.
// We generate only the (nverts - 1)! permutations with  pi[1] = 1.
// We define indexA by pi[indexA] = 2.
//
void getperms(int nverts, int PermList[maxperms][maxsize], int indexAlist[maxperms])
{
int count;
bool up;
int i;
int nperm;
int pi[maxsize];
int indexA;
int n, factnvertsm1;
if(nverts > 8)
  {
  cout << "Actually, getperms is not set for nverts > 8." << endl;
  exit(1);
  }
for(i = 0; i < maxsize; i++) pi[i] = i; // Initialization.
for(nperm = 0; nperm < maxperms; nperm++)
{
  for(i = 0; i < maxsize; i++) PermList[nperm][i] = 0; // Initialization.
}
pi[0] = 0;  // We never use this one.
n = nverts - 1;  // We actually permute nverts - 1 objects.
up = false;
count = 0;
perm(n,pi,count,up,PermList); // This function does the work.
// We now have the PermList. Put it in standard form.
DESindex indx(nverts); // Then indx(i) will return i mod nverts.
factnvertsm1 = factorial(nverts-1);
if(count != factnvertsm1)
{
  cout << "Oops, we didn't get the right number of permuations in getperms." << endl;
  exit(1);
}
for(nperm = 1; nperm <= factnvertsm1; nperm++)
{
  for(i = 1; i <= nverts; i++) pi[i] = PermList[nperm][i];
  int i1 = 0;
  int i2 = 0;
  for(i = 1; i <= nverts; i++)
  {
    if(pi[i]==1) i1 = i;
    if(pi[i]==2) i2 = i;
  }
  if(i1 == 0 || i2 == 0) 
  {
    cout << "Oops, i1 and i2 should have been changed in getperms." << endl;
    exit(1);
  }
  indexA = indx(i2-i1);  // pi[indexA] = 2.
  for(i = 1; i <= nverts; i++)
  {
    PermList[nperm][i] = pi[indx(i+i1)];
  }
  indexAlist[nperm] = indexA;
}
return;
}
//--------------
// This calls the recursively defined function perm(...).
void perm(int n, int pi[], int& count, bool up, int PermList[maxperms][maxsize])
{
int i,j,k,temp;
// If called from "higher n" version of perm and here n > 2,
// call perm() to construct permutations of n-1 objects.
if(!up && n > 2) perm(n-1,pi,count,up,PermList);
// Now that we have constructed permuations of (n - 1) objects, we want to
// rotate the first n elements and then construct permutations of new list of
// (n - 1) objects. We will do this for all (n-1) rotated version of our
// n objects (in addition to the unrotated version).
for (i = 1; i < n; i++)
  {
  if (n == 2)
    {
    // Increment count and record this permutation.
    count++;                                
    PermList[count][0] = count;
    for(k = 1; k < maxsize; k++) PermList[count][k] = pi[k]; 
    // Rotate first n elements, here for n = 2.
    temp = pi[1];                            
    pi[1] = pi[2];
    pi[2] = temp;
    count++;
    // Increment count and record this permutation.
    PermList[count][0] = count;
    for(k = 1; k < maxsize; k++) PermList[count][k] = pi[k]; 
    }
  if(n > 2)
    {
    // Rotate first n elements.
    temp = pi[1];
    for (j = 1; j < n; j++) pi[j] = pi[j+1];
    pi[n] = temp;
    // Call perm() to construct permutations of n-1 objects.
    up = false;
    perm(n-1,pi,count,up,PermList);
    }
  }
// We have rotated first n elements (n - 1) times.
// Rotate first n elements once more.
temp = pi[1];
for (j = 1; j < n; j++) pi[j] = pi[j+1];
pi[n] = temp;
up = true;                             // Now we are back at original state.
return;                                // We return "up" to perm with one greater n.
}
// End of function perm(n, pi[], count, up, PermList)
//==================================================================================================
// Function n!.
int factorial( const int& n)
{
  int result = 1;
  for(int i = 2; i <= n; i++ )
    result = result * i;
  return result;
}
// End of function factorial(n).
//==================================================================================================
// Function to give cyclic index in the range 1,...,n.
int cyclic(const int& i, const int& n)
{
  int iout;
  if(i > 0)
  {
    if(i <= n)
    {
      iout = i;   // 0 < i <= n, the most common case
    }
    else if(i <= 2*n)
    {
      iout = i - n;  // n < i <= 2*n
    }
    else
    {
      iout = (i-1)%n + 1; // 2*n < i , use the general formula
    }
  }
  else if(i > -n)
  {
    iout = i + n;  // -n < i <= 0
  }
  else
  {
    iout = (i-1)%n + 1; // i <= -n, use the general formula
  }
  return iout;
}
// End of function cyclic(i,  n).
//==================================================================================================
// Function Determinant(M,d).
// This function calculates the determinant of any input matrix using LU decomposition.
// Code by W. Gong and D.E. Soper
complex<double> Determinant(const complex<double> bb[maxsize][maxsize], const int& dimension)
{
  // Define matrix related variables.
  static const complex<double> one(1.0, 0.0);
  complex<double> determinant;
  complex<double> aa[dimension][dimension];
  int indx[dimension], d;
  // Initialize the determinant.
  determinant = one;
  // Inintialize the matrix to be decomposed with the transferred matrix b.
  for(int i=0; i < dimension; i++)
  {
    for(int j=0; j < dimension; j++)
    {
      aa[i][j] = bb[i][j];
    }
  }
  // Define parameters used in decomposition.
  int i, imax, j, k, flag=1;
  complex <double> dumc,sum;
  double aamax, dumr;
  double vv[maxsize];
  // Initialize the parity parameter.
  d = 1;
  // Get the implicit scaling information.
  for(i = 0; i < dimension; i++)
  {
    aamax=0;
    for(j = 0; j < dimension; j++)
    {
      if(abs(aa[i][j]) > aamax) aamax=abs(aa[i][j]);
    }
    // Set a flag to check if the determinant is zero.
    if(aamax == 0) flag=0;
    // Save the scaling.
    vv[i]=1.0/aamax;
  }
  if(flag == 1)
  {
    for(j = 0; j < dimension; j++)
    {
      for(i = 0; i < j; i++)
      {
        sum = aa[i][j];
        for(k = 0; k < i; k++) sum = sum - aa[i][k]*aa[k][j];
        aa[i][j] = sum;
      }
      //Initialize for the search for largest pivot element.
      aamax=0;
      for(i = j; i < dimension; i++)
      {
        sum = aa[i][j];
        for(k = 0; k < j; k++) sum = sum - aa[i][k]*aa[k][j];
        aa[i][j]=sum;
        // Figure of merit for the pivot.
        dumr = vv[i] * abs(sum);
        // Is it better than the best so far?
        if(dumr >= aamax)
        {
          imax=i;
          aamax=dumr;
        }
      }  // End for(i = j; i < dimension; i++)
         // See if we need to interchange rows.
      if(j != imax)
      {
        for(k = 0; k < dimension; k++)
        {
          dumc = aa[imax][k];
          aa[imax][k] = aa[j][k];
          aa[j][k] = dumc;
        }
        // Change the parity of d.
        d = -d;
        // Interchange the scale factor.
        vv[imax] = vv[j];
      }  // End if(j != imax)
      indx[j]=imax;
      if(j != dimension - 1)
      {
        dumc = 1.0/aa[j][j];
        for(i = j+1; i < dimension; i++) aa[i][j] = aa[i][j]*dumc;
      }
    } // End for(j = 0; j < dimension; j++)
  } // END if(flag == 1)
    // Calculate the determinant using the decomposed matrix.
  if(flag == 0) determinant = 0.0;
  else
  {    // Multiply the diagonal elements.
    for(int diagonal = 0; diagonal < dimension; diagonal++)
    {
      determinant = determinant * aa[diagonal][diagonal];
    }
    determinant=(double)d*determinant;
  } // End if(flag == 0) ... else
  return determinant;
}
// End of function Determinant(M,d).
//==================================================================================================
// Function fourphoton(helicities, pvectors).
// Produces the four-photon scattering amplitude.
// Result is from
//   T.~Binoth, E.~W.~N.~Glover, P.~Marquard and J.~J.~van der Bij,
//   %``Two-loop corrections to light-by-light scattering in supersymmetric  QED,''
//   JHEP {\bf 0205}, 060 (2002)
//   [arXiv:hep-ph/0202266].
//
double fourphoton(int helicity[maxsize],DESreal4vector p[maxsize])
{
double s,t,u;
double X,Y;
double calM;
const double pi = 3.14159265358979;
//-------
s = 2.0*dot(p[1],p[2]);
t = 2.0*dot(p[1],p[3]);
u = 2.0*dot(p[1],p[4]);
X = log(-t/s);
Y = log(-u/s);
if(helicity[1] == 1 && helicity[2] == 1 && helicity[3] == 1 && helicity[4] == 1)
  {
  calM = -8.0;
  }
else if(helicity[1] == 1 && helicity[2] == 1 && helicity[3] == 1 && helicity[4] == -1)
  {
  calM = 8.0;
  }
else if(helicity[1] == 1 && helicity[2] == 1 && helicity[3] == -1 && helicity[4] == -1)
  {
  calM =      - 4.0 - 4.0*((X-Y)*(X-Y) + pi*pi - 2*(X-Y))*t*t/s/s;
  calM = calM - 4.0 - 4.0*((Y-X)*(Y-X) + pi*pi - 2*(Y-X))*u*u/s/s;
  }
else
  {
  cout << "Actually, I am not prepared for that helicity combination." << endl;
  calM = 0.0;
  }
return calM;
}
//==================================================================================================