//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;
}
//==================================================================================================