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