Page MenuHomeHEPForge
Files F10664324

No OneTemporary
Actions

View Options

Index: branches/bbz/src/4F0_LO.C
===================================================================
--- branches/bbz/src/4F0_LO.C (revision 37)
+++ branches/bbz/src/4F0_LO.C (revision 38)
@@ -1,120 +1,119 @@
///
// Calculate the massless limit at LO and NLO
///
#include "Massless_lim.H"
#include "Massless_Kernels.H"
#include <iostream>
#include <vector>
#include <cmath>
using namespace LHAPDF;
using namespace std;
double Massless_Limit::dsigmagg_lo(double y, double z)
{
double ds(0.), m2b(p_lumi->m2bext()), L(log(m2h/m2b)/(2.*M_PI));
double asr(2.*M_PI*p_lumi->ass(m2r));
ds = asr*asr*(Agg_2_2(y,z)*L*L + Agg_2_1(y,z)*L + Agg_2_0(y)); // Eq (29) of FONLL-A paper
return ds; // A20 is not convoluted with anything, z variable of convolution
}
-double Massless_Limit::dsigmaqqb_lo(double y, double z) // This term should be the same as the one in the 5F but at the moment it doe not work
+double Massless_Limit::dsigmaqqb_lo(double y, double z)
{
double ds(0.);
double asr(2.*M_PI*p_lumi->ass(m2r));
double mylfh=0., mylfr=0.;
- ds = asr*asr*deltaqqba_(&y,&mylfh,&mylfr)/pow(4.*M_PI,2)/(1. - y); // Tutte le cose costanti in z sono divise per il Jacobiano
+ ds = asr*asr*deltaqqba_(&y,&mylfh,&mylfr)/(1. - y);
return ds;
}
-double Massless_Limit::Agg_2_0(double y) // This term should be the same as the one in the 5F but at the moment it doe not work
+double Massless_Limit::Agg_2_0(double y)
{
// constant in z!
double mylfh=0., mylfr=0.;
- return deltagga_(&y,&mylfh,&mylfr)/pow(4.*M_PI,2)/(1.-y); // deltagga e' una funzione NNLO di Harlander
+ return deltagga_(&y,&mylfh,&mylfr)/(1.-y);
}
double Massless_Limit::Agg_2_1(double y, double z)
{
double pqg = p_lumi->Pqg(z);
double yz = y/z;
- double mylfh=0., mylfr=0.;
- return 4.*pqg*sbg1_(&yz,&mylfh)/M_PI/4.;
+ double mylfh=0.;
+ return 4.*pqg*sbg1_(&yz,&mylfh);
}
double Massless_Limit::Agg_2_2(double y, double z)
{
double pqg_1 = p_lumi->Pqg(z);
double pqg_2 = p_lumi->Pqg(y/z);
return 2.*y*pqg_1*pqg_2/z; // 2 is b=b~
}
double Massless_Limit::L_qq(const double x, const double y)
{
double lqq(0.);
lqq += p_lumi->Lumi(4,6,x,y,1.,1.,10);//d
lqq += p_lumi->Lumi(6,4,x,y,1.,1.,10);
lqq += p_lumi->Lumi(3,7,x,y,1.,1.,10);//u
lqq += p_lumi->Lumi(7,3,x,y,1.,1.,10);
lqq += p_lumi->Lumi(2,8,x,y,1.,1.,10);//s
lqq += p_lumi->Lumi(8,2,x,y,1.,1.,10);
lqq += p_lumi->Lumi(1,9,x,y,1.,1.,10);//c
lqq += p_lumi->Lumi(9,1,x,y,1.,1.,10);
return lqq;
}
double Massless_Limit::L_gq(const double x, const double y)
{
double lgq(0.);
lgq += p_lumi->Lumi(5,6,x,y,1.,1.,10);
lgq += p_lumi->Lumi(5,4,x,y,1.,1.,10);
lgq += p_lumi->Lumi(5,7,x,y,1.,1.,10);
lgq += p_lumi->Lumi(5,3,x,y,1.,1.,10);
lgq += p_lumi->Lumi(5,8,x,y,1.,1.,10);
lgq += p_lumi->Lumi(5,2,x,y,1.,1.,10);
lgq += p_lumi->Lumi(5,9,x,y,1.,1.,10);
lgq += p_lumi->Lumi(5,1,x,y,1.,1.,10);
return lgq;
}
double Massless_Limit::dsigma0LO(double z[], size_t dim, void* p)
{
double dsg(0.), x,y;
double tau, jac(1.), xi;
double t, tmin;
double lgg,lqq;
double sgg(0.), sqq(0.);
tau = tauh + (1. - tauh)*z[0];
jac *= (1. - tauh);
double ximax = - 0.5*log(tau);
xi = -ximax + 2.*ximax*z[1];
jac *= 2.*ximax;
//standard variables x, y ;
x = sqrt(tau)*exp(xi);
y = sqrt(tau)*exp(-xi);
tmin = y;
t = tmin + (1. - tmin)*z[2];
jac *= (1. - tmin);
// gg
sgg = dsigmagg_lo(y,t);
// q q~
sqq = dsigmaqqb_lo(y,t);
lgg = p_lumi->Lumi(5,5,x,tauh/tau,1.,1.,10);
lqq = L_qq(x,tauh/tau);
- dsg = jac*s0*( lqq * sqq
- + lgg * sgg )/tau;
+ dsg = jac*s0*( lqq * sqq + lgg * sgg )/tau;
return dsg;
}
Index: branches/bbz/src/bbznnlo/bits.f
===================================================================
--- branches/bbz/src/bbznnlo/bits.f (revision 37)
+++ branches/bbz/src/bbznnlo/bits.f (revision 38)
@@ -1,1944 +1,1957 @@
***********************************************
* delta term at nlo
***
double precision function delta1(lfr)
implicit none
double precision pi,cf,z2
double precision lfr,lfac
c lfr here is actually -lfh
lfac = -lfr
cf = 4d0 / 3d0
pi = 3.1415926535897932385d0
z2 = 1.6449340668482264365d0
delta1 = 8.0d0 * z2 - 16.0d0 + 6.0d0 * lfac
- delta1 = cf * delta1
+ delta1 = cf * delta1/4.d0/pi
return
end
***
* subtraction term associated with + prescription at nlo
***
double precision function dtsub1(lfh,tauh)
implicit none
- double precision cf,lfac
+ double precision cf,lfac,pi
double precision ltmin,lfh,tauh
c lfr here is actually -lfh
lfac = -lfh
-
+ pi = 3.1415926535897932385d0
cf = 4.d0/3d0
+
ltmin = dlog(1.0d0 - tauh)
dtsub1 = 8.0d0 * ltmin**2 + 8.0d0*ltmin * lfac
- dtsub1 = cf * dtsub1
+ dtsub1 = cf * dtsub1/4.d0/pi
return
end
***
* a term in front of delta function at nlo
***
double precision function dterms1(xt,lfh)
implicit none
- double precision ca,cf,nf
+ double precision ca,cf,nf,pi
double precision xt,lfh,lfac
double precision dd0,dd1
c lfh here is actually -lfh
lfac = -lfh
-
+ pi = 3.1415926535897932385d0
+
ca = 3.d0
cf = 4d0/3d0
nf = 5d0
dd0 = 1.d0/(1.0d0-xt)
dd1 = dd0 * dlog(1.0d0-xt)
dterms1 = 16.0d0 * dd1
dterms1 = dterms1 + 8.0d0 * lfac * dd0
- dterms1 = cf * dterms1
+ dterms1 = cf * dterms1/4.d0/pi
return
end
***
* regular terms at nlo in qq
***
double precision function sbbq1(xt,lfh)
implicit none
double precision ca,cf,nf
double precision xt,lfh,lfac
- double precision dd0,dd1
+ double precision dd0,dd1,pi
c lfr here is actually -lfh
lfac = -lfh
-
+ pi = 3.1415926535897932385d0
+
ca = 3.d0
cf = 4d0/3d0
nf = 5d0
sbbq1 = -8.0d0*( 1.0d0 + xt )*dlog(1.0d0-xt)
. -4.0d0*( 1.0d0 + xt**2 )*dlog(xt)/(1.0d0-xt)
sbbq1 = sbbq1 - 4.0d0*( 1.0d0 + xt )*lfac
- sbbq1 = xt* cf * sbbq1
+ sbbq1 = xt* cf * sbbq1/4.d0/pi
return
end
***
* term in front of qg function at nlo
***
double precision function sbg1(xt,lfh)
implicit none
double precision ca,tf
double precision xt,lfh,lfac
- double precision z,zmin,log1mz
+ double precision z,zmin,log1mz,pi
c lfh here is actually -lfh
lfac = -lfh
-
+ pi = 3.1415926535897932385d0
+
tf = 0.500000000000000000000d0
z = xt
zmin = 1.0d0-z
log1mz = dlog(zmin)
sbg1 = 2.0d0*( z**2 + zmin**2 )*( 2.0d0*log1mz - dlog(z) ) +
. 1.d0 + 6.0d0*z - 7.0d0*z**2
sbg1 = sbg1 + 2.d0*( z**2 + zmin**2 )*lfac
- sbg1 = xt * tf * sbg1
+ sbg1 = xt * tf * sbg1/4.d0/pi
return
end
c.......................................................................
c............ NNLO .................................................
c.......................................................................
***
* Regular term in bb~
**
double precision function deltabbqa(xx,lfh,lfr)
implicit none
double precision xx,lfh,lfr,z,lfac,lfren
double precision CQQIV,CQQ,CQQIA,CQBBB,CQBACD
double precision CQBBCD,CQBAB,CQBCAR,CQBNFR,CQBCFR
- double precision cax,cv,cd,sw
+ double precision cax,cv,cd,sw,pi
external CQQIV,CQQ,CQQIA,CQBBB,CQBACD
external CQBBCD,CQBAB,CQBCAR,CQBNFR,CQBCFR
sw = 0.2228972225239183d0
cax = 1.d0
cv = -1.d0 + 4.d0/3.d0*sw
cd = 1.d0 + cv**2
-
+ pi = 3.1415926535897932385d0
+
c lfh here is actually -lfh,while lfr = -lfh+lfr
lfren = -lfh+lfr
lfac = -lfh
deltabbqa = xx*(cd*(CQBBB(xx,lfac,lfren) + 2.d0*CQQ(xx,lfac,lfren)
? + 2.d0*CQBBCD(xx,lfac,lfren) )
? - 2.d0*(cv**2)*CQQIV(xx,lfac,lfren)
? + 2.d0*(cax**2)*CQQIA(xx,lfac,lfren)
? + CQBAB(xx,lfac,lfren)
? + cd*(CQBCFR(xx,lfac,lfren) + CQBCAR(xx,lfac,lfren)
? + CQBNFR(xx,lfac,lfren) + 2.d0*CQBACD(xx,lfac,lfren)))
- deltabbqa = deltabbqa/cd
+ deltabbqa = deltabbqa/cd/((4.d0*pi)**2)
return
end
c.......................................................................
double precision function deltabbqf(xx,lfh,lfr)
implicit none
- double precision xx,lfh,lfr
- deltabbqf = 0d0
+ double precision xx,lfh,lfr,pi
+ pi = 3.1415926535897932385d0
+
+ deltabbqf = 0d0/((4.d0*pi)**2)
return
end
***
* Regular term in bg
**
double precision function deltabga(xx,lfh,lfr)
implicit none
double precision xx,lfh,lfr,lfac,lfren
double precision xm1,xp1,dlnx,dlxm1,dlxp1
double precision dli2a,dli3a,zcom,s1mz,dimz,trimz
double precision trimco,tripco
double precision z2,z3,z4
- double precision ca, cf, nf
+ double precision ca, cf, nf,pi
double precision partCA, partCF,scaleCA,scaleCF,scaleNF
double precision dilog,trilog,wgplg
external dilog,trilog,wgplg
c lfh here is actually -lfh,while lfr = -lfh+lfr
lfren = -lfh+lfr
lfac = -lfh
-
+ pi = 3.1415926535897932385d0
ca = 3.d0
cf = 4.d0/3.d0
nf = 5.d0
z2 = 1.6449340668482264365d0
z3 = 1.2020569031595942854d0
z4 = 1.0823232337111381915d0
xm1 = 1.d0 - xx
dlnx = dlog(xx)
dlxm1 = dlog(xm1)
xp1 = 1.0d0 + xx
dlxp1 = dlog(xp1)
dli2a = dilog(xm1)
dli3a = trilog(xm1)
dimz = wgplg(1,1,-xx)
trimz = wgplg(2,1,-xx)
zcom = xm1/xp1
tripco= wgplg(2,1,zcom)
trimco= wgplg(2,1,-zcom)
s1mz = wgplg(1,2,xm1)
partCA =
. dlnx**3 * ( - 3.d0 - 20.d0/3.d0*xx )
. + dlnx**2*dlxm1 * ( 6.d0 - 4.d0*xx**2 + 28.d0*xx )
. + dlnx**2*dlxp1 * ( - 6 - 12*xx**2 - 12*xx )
. + dlnx**2 * ( 5.d0/2.d0 + 173.d0/3.d0*xx**2 + 2.0d0*xx )
. + dlnx*dlxm1**2 * ( - 2.d0 + 12.d0*xx**2 - 44.d0*xx )
. + dlnx*dlxm1*dlxp1 * ( 8.d0 + 16.d0*xx**2 + 16.d0*xx )
. + dlnx*dlxm1 * ( - 10.d0 - 130.d0*xx**2 + 16.d0*xx )
. + dlnx*dlxp1 * ( - 4.d0 - 16.d0*xx**2 - 20.d0*xx )
. + dlnx*dli2a * ( 8.d0*xx**2 - 28.d0*xx )
. + dlnx*dimz * ( - 8.d0 - 16.d0*xx**2 - 16.d0*xx )
. + dlnx*z2 * ( - 16.d0*xx**2 + 40.d0*xx )
. + dlnx * ( - 118.d0/3.d0 + 457.d0/9.d0*xx**2 + 16.d0/3.d0*xx )
. + dlxm1**3 * ( - 13.d0/3.d0 - 26.d0/3.d0*xx**2 + 26.d0/3.d0*xx )
. + dlxm1**2 * ( - 4.d0 + 154.d0/3.d0*xx**2 - 42.d0*xx
. - 16.d0/3.d0/xx )
. + dlxm1*dli2a * ( - 28.d0 - 20.d0*xx**2 - 40.d0*xx )
. + dlxm1*dimz * ( 8.d0 + 16.d0*xx**2 + 16.d0*xx )
. + dlxm1*z2 * ( 16.d0 + 32.d0*xx**2 - 16.d0*xx )
partCA = partCA
. + dlxm1 * ( 70.d0/3.d0 - 74.d0/9.d0*xx**2
. - 25.d0/3.d0*xx - 88.d0/9.d0/xx )
. + dli2a * ( - 22.d0 - 88.d0/3.d0*xx**2 - 64.d0*xx
. - 32.d0/3.d0/xx )
. + dimz * ( - 4.d0 - 16.d0*xx**2 - 20.d0*xx )
. + z2 * ( - 10.d0 - 214.d0/3.d0*xx**2 + 56.d0*xx
. + 16.d0/3.d0/xx )
. + dli3a * ( 30.d0 + 24.d0*xx**2 + 68.d0*xx )
. + trimz * ( 4 + 8*xx**2 + 8*xx )
. + trimco * ( 8.d0 + 16.d0*xx**2 + 16.d0*xx )
. + tripco * ( - 8.d0 - 16.d0*xx**2 - 16.d0*xx )
. + s1mz * ( - 36.d0 - 16.d0*xx**2 - 64.d0*xx )
. + z3 * ( 2.d0 + 4.d0*xx**2 + 8.d0*xx )
. - 539.d0/18.d0 - 1837.d0/54.d0*xx**2 + 613.d0/9.d0*xx
. - 58.d0/27.d0/xx
partCA = - ca * partCA
* -- CF
partCF =
. dlnx**3 * ( 17.d0/6.d0 + 26.d0/3.d0*xx**2 - 17.d0/3.d0*xx )
. + dlnx**2*dlxm1 * ( - 12.d0 - 40.d0*xx**2 + 24.d0*xx )
. + dlnx**2 * ( 11.d0/4.d0 + xx**2 - 15.d0*xx )
. + dlnx*dlxm1**2 * ( 21.d0 + 66.d0*xx**2 - 42.d0*xx )
. + dlnx*dlxm1 * ( - 14.d0 - 96.d0*xx**2 + 96.d0*xx )
. + dlnx*dlxp1 * ( - 8.d0 + 24.d0*xx**2 + 16.d0*xx )
. + dlnx*dli2a * ( - 2.d0 + 4.d0*xx )
. + dlnx*dimz * ( 8.d0 + 16.d0*xx**2 - 16.d0*xx )
. + dlnx*z2 * ( - 12.d0 - 48.d0*xx**2 + 24.d0*xx )
. + dlnx * ( 31.d0/2.d0 + 87.d0*xx**2 - 201.d0/2.d0*xx )
. + dlxm1**3 * ( - 35.d0/3.d0 - 70.d0/3.d0*xx**2 + 70.d0/3.d0*xx )
. + dlxm1**2 * ( 23.d0 + 63.d0*xx**2 - 80.d0*xx )
. + dlxm1*dli2a * ( 6.d0 + 52.d0*xx**2 - 12.d0*xx )
. + dlxm1*z2 * ( 8.d0 + 16.d0*xx**2 - 16.d0*xx )
partCF = partCF
. + dlxm1 * ( - 26.d0 - 88.d0*xx**2 + 135.d0*xx )
. + dli2a * ( 9.d0 - 40.d0*xx**2 + 24.d0*xx )
. + dimz * ( - 8.d0 + 24.d0*xx**2 + 16.d0*xx )
. + z2 * ( - 10.d0 + 24.d0*xx**2 - 4.d0*xx )
. + dli3a * ( 2.d0 - 36.d0*xx**2 - 4.d0*xx )
. + trimz * ( - 16.d0 - 32.d0*xx**2 + 32.d0*xx )
. + s1mz * ( 22.d0 + 68.d0*xx**2 - 44.d0*xx )
. + z3 * ( - 50.d0 - 100.d0*xx**2 + 100.d0*xx )
. + 157.d0/4.d0 + 305.d0/4.d0*xx**2 - 221.d0/2.d0*xx
partCF = - cf * partCF
scaleCA =
. lfac**2*( dlnx * ( 2.d0 + 8.d0*xx )
. + dlxm1 * ( 2.d0 + 4.d0*xx**2 - 4.d0*xx )
. + ( 14.d0 - 9.d0*xx**2 + 2.d0*xx + 4.d0/xx ) / 3.d0 ) +
. lfac*lfren*( -11.d0/3.d0 - 22.d0/3.d0*xx**2 + 22.d0/3.d0*xx )
. + lfac *( dlnx**2 * ( - 4.d0 - 12.d0*xx )
. + dlnx*dlxm1 * ( 4.d0 - 8.d0*xx**2 + 40.d0*xx )
. + dlnx*dlxp1 * ( - 4.d0 - 8.d0*xx**2 - 8.d0*xx )
. + dlnx * ( 7.d0 + 146.d0*xx**2 + 10.d0*xx ) / 3.d0
. + dlxm1**2 * ( 6.d0 + 12.d0*xx**2 - 12.d0*xx )
. + dlxm1 * ( 40.d0/3.d0 - 98.d0/3.d0*xx**2 +
. 64.d0/3.d0*xx + 16.d0/3.d0/xx )
. + dli2a * ( 12.d0 + 8.d0*xx**2 + 24.d0*xx )
. + dimz * ( - 4.d0 - 8.d0*xx**2 - 8.d0*xx )
. + z2 * ( - 8.d0 - 16.d0*xx**2 + 8.d0*xx )
. - 47.d0/6.d0 - 85.d0/18.d0*xx**2 + 29.d0/3.d0*xx
. + 44.d0/9.d0/xx ) +
. lfren *( dlnx * ( 11.d0 + 22.d0*xx**2 - 22.d0*xx ) / 3.d0
. + dlxm1 * ( -22.d0 - 44.d0*xx**2 + 44.d0*xx ) / 3.d0
. - 11.d0/6.d0 + 77.d0/6.d0*xx**2 - 11.d0*xx )
scaleCA = CA * scaleCA
scaleCF =
. lfac**2*( dlnx * ( - 3.d0 - 12.d0*xx**2 + 6.d0*xx )
. + dlxm1 * ( 6.d0 + 12.d0*xx**2 - 12.d0*xx )
. - 3.d0/2.d0 + 6.d0*xx ) +
. lfac *( dlnx**2 * ( 4.d0 + 16.d0*xx**2 - 8.d0*xx )
. + dlnx*dlxm1 * ( - 20.d0 - 64.d0*xx**2 + 40.d0*xx )
. + dlnx * ( 5.d0 + 46.d0*xx**2 - 40.d0*xx )
. + dlxm1**2 * ( 18.d0 + 36.d0*xx**2 - 36.d0*xx )
. + dlxm1 * ( - 16.d0 - 46.d0*xx**2 + 68.d0*xx )
. + dli2a * ( - 24.d0*xx**2 )
. + z2 * ( - 4.d0 - 8.d0*xx**2 + 8.d0*xx )
. + 12.d0 + 11.d0*xx**2 - 34.d0*xx )
scaleCF = CF * scaleCF
scaleNF =
. lfac**2*( - 2.d0/3.d0 - 4.d0/3.d0*xx**2 + 4.d0/3.d0*xx ) +
. lfac*lfren*( 2.d0/3.d0 + 4.d0/3.d0*xx**2 - 4.d0/3.d0*xx ) +
. lfac *( dlnx * ( 2.d0 + 4.d0*xx**2 - 4.d0*xx ) / 3.d0
. + dlxm1 * ( - 4.d0 - 8.d0*xx**2 + 8.d0*xx ) / 3.d0
. - 1.d0/3.d0 + 7.d0/3.d0*xx**2 - 2.d0*xx ) +
. lfren *( dlnx * ( - 2.d0 - 4.d0*xx**2 + 4.d0*xx ) / 3.d0
. + dlxm1 * ( 4.d0 + 8.d0*xx**2 - 8.d0*xx ) / 3.d0
. + 1.d0/3.d0 - 7.d0/3.d0*xx**2 + 2.d0*xx )
scaleNF = NF * scaleNF
deltabga = xx*(partCA + partCF + scaleCA + scaleCF + scaleNF)
+ . /((4.d0*pi)**2)
return
end
c.......................................................................
double precision function deltabgf(xx,lfh,lfr)
implicit none
double precision xx,lfh,lfr
deltabgf = 0d0
return
end
***
* Regular term in gg - CGG
**
double precision function deltagga(xx,lfh,lfr)
implicit none
double precision xx,lfh,lfr,z,lfac,lfren
double precision xm1,xp1,dlnx,dlxm1,dlxp1
double precision dli2a,dli3a,zcom,s1mz,dimz,trimz,smz
double precision trimco,tripco
double precision z2,z3,z4
double precision delCA,delCF1,delCF2,scale
double precision part1, part2
- double precision ca,cf,nf
+ double precision ca,cf,nf,pi
double precision dilog,trilog,wgplg
external dilog,trilog,wgplg
c lfh here is actually -lfh,while lfr = -lfh+lfr
lfren = -lfh+lfr
lfac = -lfh
-
+ pi = 3.1415926535897932385d0
+
ca=3.d0
cf=4.d0/ca
nf=5d0
z2 = 1.6449340668482264365d0
z3 = 1.2020569031595942854d0
z4 = 1.0823232337111381915d0
xm1 = 1.d0 - xx
dlnx = dlog(xx)
dlxm1 = dlog(xm1)
xp1 = 1.0d0 + xx
dlxp1 = dlog(xp1)
dli2a = dilog(xm1)
dli3a = trilog(xm1)
dimz = wgplg(1,1,-xx)
trimz = wgplg(2,1,-xx)
zcom = xm1/xp1
tripco= wgplg(2,1,zcom)
trimco= wgplg(2,1,-zcom)
s1mz = wgplg(1,2,xm1)
smz = wgplg(1,2,-xx)
delCF1 = dli3a * ( 16.0d0 + 64.0d0*xx + 64.0d0*xx**2 ) +
. trimz * ( - 8.0d0 - 16.0d0*xx + 8.0d0*xx**2 ) +
. s1mz * ( - 8.0d0 - 80.0d0*xx - 56.0d0*xx**2 ) +
. smz * ( - 16.0d0 - 32.0d0*xx - 16.0d0*xx**2 ) +
. z3 * ( - 4.0d0 - 8.0d0*xx + 8.0d0*xx**2 ) +
. dlxm1*dli2a * ( - 16.0d0 - 64.0d0*xx - 64.0d0*xx**2 ) +
. dlnx*dli2a * ( - 4.0d0 - 16.0d0*xx - 16.0d0*xx**2 ) +
. dlxp1*dimz * ( - 16.0d0 - 32.0d0*xx - 16.0d0*xx**2 ) +
. dlnx*dimz * ( 16.0d0 + 32.0d0*xx + 8.0d0*xx**2 ) +
. dlnx*z2 * ( 12.0d0 + 40.0d0*xx + 40.0d0*xx**2 ) +
. dlxp1*z2 * ( - 8.0d0 - 16.0d0*xx - 8.0d0*xx**2 ) +
. dlxm1**2*dlnx * ( - 8.0d0 - 32.0d0*xx - 32.0d0*xx**2 ) +
. dlxm1*dlnx**2 * ( 4.0d0 + 16.0d0*xx + 16.0d0*xx**2 ) +
. dlnx**3 * ( - 2.0d0 - 16.0d0/3.0d0*xx
. - 16.0d0/3.0d0*xx**2 ) +
. dlnx**2*dlxp1 * ( 12.0d0 + 24.0d0*xx + 12.0d0*xx**2 ) +
. dlnx*dlxp1**2 * ( - 8.0d0 - 16.0d0*xx - 8.0d0*xx**2 )
delCF2 = dli2a * ( - 20.0d0 - 16.0d0*xx + 56.0d0*xx**2 ) +
. dimz * ( 8.0d0 + 8.0d0*xx ) +
. z2 * ( 20.0d0 + 36.0d0*xx - 48.0d0*xx**2 ) +
. dlxm1**2 * ( - 16.0d0 - 32.0d0*xx + 48.0d0*xx**2 ) +
. dlxm1*dlnx * ( 4.0d0 + 32.0d0*xx - 16.0d0*xx**2 ) +
. dlnx**2 * ( - 6.0d0 - 14.0d0*xx - 8.0d0*xx**2 ) +
. dlnx*dlxp1 * ( 8.0d0 + 8.0d0*xx ) +
. dlxm1 * ( 14.0d0 + 120.0d0*xx - 134.0d0*xx**2 ) +
. dlnx * ( - 23.0d0 - 64.0d0*xx + 105.0d0*xx**2 )
. - 32.0d0 - 66.0d0*xx + 98.0d0*xx**2
delCA = s1mz * ( -8.0d0*xx**2 + 16.0d0*xx - 8.0d0 ) +
. smz * ( 16.0d0*xx**2 + 32.0d0*xx + 16.0d0 ) +
. trimz * ( 24.0d0*xx**2 + 48.0d0*xx + 24.0d0 ) +
. z3 * ( 16.0d0*xx**2 + 32.0d0*xx + 16.0d0 ) +
. dimz * dlnx * ( -24.0d0*xx**2 - 48.0d0*xx - 24.0d0 ) +
. dimz * dlxp1 * ( 16.0d0*xx**2 + 32.0d0*xx + 16.0d0 ) +
. dimz * ( 16.0d0*xx**2 + 32.0d0*xx +16.0d0 )/3.0d0 +
. z2 * dlxp1 * ( 8.0d0*xx**2 + 16.0d0*xx + 8.0d0 ) +
. z2 * ( 8.0d0*xx**2 + 16.0d0*xx +8.0d0 )/3.0d0 +
. dlnx**2 * dlxp1 * ( -12.0d0*xx**2 - 24.0d0*xx - 12.0d0 ) +
. dlnx**2 * ( 50.0d0*xx**2 + 4.0d0*xx - 4.0d0 )/3.0d0 +
. dlnx * dlxp1**2 * ( 8.0d0*xx**2 + 16.0d0*xx + 8.0d0 ) +
. dlnx * dlxp1 * ( 16.0d0*xx**2 + 32.0d0*xx +
. 16 .0d0)/3.0d0 +
. dlnx * ( -50.0d0*xx**2 - 76.0d0/3.0d0*xx - 4.0d0 ) +
. 191.0d0/3.0d0*xx**2 - 48.0d0*xx - 47.0d0/3.0d0
delCA = ca**2 / ( ca**2 - 1.0d0 ) * delCA
scale = ( dlnx * ( - 2.0d0 - 8.0d0*xx - 8.0d0*xx**2 )
. - 4.0d0 - 8.0d0*xx + 12.0d0*xx**2 ) * lfac**2 +
. ( dli2a * ( - 8.0d0 - 32.0d0*xx - 32.0d0*xx**2 ) +
. dlxm1*dlnx * ( - 8.0d0 - 32.0d0*xx - 32.0d0*xx**2 ) +
. dlnx**2 * ( 2.0d0 + 8.0d0*xx + 8.0d0*xx**2 ) +
. dlxm1 * ( - 16.0d0 - 32.0d0*xx + 48.0d0*xx**2 ) +
. dlnx * ( 2.0d0 + 16.0d0*xx - 8.0d0*xx**2 ) +
. 7.0d0 + 60.0d0*xx - 67.0d0*xx**2) * lfac
- deltagga = xx*(delCF1 + delCF2 + delCA + scale)
+ deltagga = xx*(delCF1 + delCF2 + delCA + scale)/(4.d0*pi)**2
return
end
c.......................................................................
double precision function deltaggf(xx,lfh,lfr)
implicit none
double precision xx,lfh,lfr
deltaggf = 0.d0
return
end
***
* Regular term in qq - CQQ
**
double precision function CQQIV(xx,lfh,lfr)
implicit none
double precision partIV
double precision xx,lfh,lfr,z
double precision xm1,xp1,dlnx,dlxm1,dlxp1
double precision dli2a,dli3a,zcom,s1mz,dimz,trimz,smz
double precision trimco,tripco
double precision part1,part2,scale
double precision z2,z3,z4
double precision ca,cf,nf
double precision dilog,trilog,wgplg
double precision cax,cv,cd,sw
external dilog,trilog,wgplg
sw = 0.2228972225239183d0
cax = 1.d0
cv = -1.d0 + 4.d0/3.d0*sw
cd = 1.d0 + cv**2
ca=3.d0
cf=4.d0/ca
nf=5d0
z2 = 1.6449340668482264365d0
z3 = 1.2020569031595942854d0
z4 = 1.0823232337111381915d0
xm1 = 1.d0 - xx
dlnx = dlog(xx)
dlxm1 = dlog(xm1)
xp1 = 1.0d0 + xx
dlxp1 = dlog(xp1)
dli2a = dilog(xm1)
dli3a = trilog(xm1)
dimz = wgplg(1,1,-xx)
trimz = wgplg(2,1,-xx)
zcom = xm1/xp1
tripco= wgplg(2,1,zcom)
trimco= wgplg(2,1,-zcom)
s1mz = wgplg(1,2,xm1)
smz = wgplg(1,2,-xx)
Z = xx
partIV =
. dlnx**3 * ( 4.D0/3.D0*Z )
. + dlnx**2*dlxp1 * ( - 20.D0 - 10.D0*Z - 20.D0/Z )
. + dlnx**2 * ( 13.D0*Z )
. + dlnx*dlxp1**2 * ( 24.D0 + 12.D0*Z + 24.D0/Z )
. + dlnx*dlxp1 * ( - 20.D0 - 20.D0*Z )
. + dlnx*dli2a * ( - 20.D0 + 2.D0*Z )
. + dlnx*DIMZ * ( - 16.D0*Z - 40.D0/Z )
. + dlnx*z2 * ( - 20.D0 - 2.D0*Z )
. + dlnx * ( 20.D0 + 16.D0*Z )
. + dlxp1*DIMZ * ( 48.D0 + 24.D0*Z + 48.D0/Z )
. + dlxp1*z2 * ( 24.D0 + 12.D0*Z + 24.D0/Z )
partIV = partIV
. + dli2a * ( - 10.D0 + 8.D0*Z )
. + DIMZ * ( - 20.D0 - 20.D0*Z )
. + z2 * ( - 10.D0 - 10.D0*Z )
. + dli3a * ( 12.D0 - 6.D0*Z - 8.D0/Z )
. + TRIMZ * ( - 40.D0 + 12.D0*Z + 40.D0/Z )
. + S1MZ * ( - 16.D0 - 8.D0*Z - 16.D0/Z )
. + SMZ * ( 48.D0 + 24.D0*Z + 48.D0/Z )
. + z3 * ( - 36.D0 + 6.D0*Z + 24.D0/Z )
. + 40.D0 - 40.D0*Z
partIV = cf * partIV
CQQIV = partIV
return
end
double precision function CQQIA(xx,lfh,lfr)
implicit none
double precision partIA
double precision xx,lfh,lfr,z
double precision xm1,xp1,dlnx,dlxm1,dlxp1
double precision dli2a,dli3a,zcom,s1mz,dimz,trimz,smz
double precision trimco,tripco
double precision part1,part2,scale
double precision z2,z3,z4
double precision ca,cf,nf
double precision dilog,trilog,wgplg
double precision cax,cv,cd,sw
external dilog,trilog,wgplg
sw = 0.2228972225239183d0
cax = 1.d0
cv = -1.d0 + 4.d0/3.d0*sw
cd = 1.d0 + cv**2
ca=3.d0
cf=4.d0/ca
nf=5d0
z2 = 1.6449340668482264365d0
z3 = 1.2020569031595942854d0
z4 = 1.0823232337111381915d0
xm1 = 1.d0 - xx
dlnx = dlog(xx)
dlxm1 = dlog(xm1)
xp1 = 1.0d0 + xx
dlxp1 = dlog(xp1)
dli2a = dilog(xm1)
dli3a = trilog(xm1)
dimz = wgplg(1,1,-xx)
trimz = wgplg(2,1,-xx)
zcom = xm1/xp1
tripco= wgplg(2,1,zcom)
trimco= wgplg(2,1,-zcom)
s1mz = wgplg(1,2,xm1)
smz = wgplg(1,2,-xx)
Z = xx
c--- CQQIA
partIA =
. dlnx**3 * ( - 4.D0/3.D0*Z )
. + dlnx**2*dlxp1 * ( 20.D0 + 10.D0*Z )
. + dlnx**2 * ( - Z )
. + dlnx*dlxp1**2 * ( - 24.D0 - 12.D0*Z )
. + dlnx*dlxp1 * ( 4.D0 + 4.D0*Z )
. + dlnx*dli2a * ( 4.D0 + 6.D0*Z )
. + dlnx*DIMZ * ( 32.D0 )
. + dlnx*z2 * ( 4.D0 + 10.D0*Z ) - dlnx * 4.D0
. + dlxp1*DIMZ * ( - 48.D0 - 24.D0*Z )
. + dlxp1*z2 * ( - 24.D0 - 12.D0*Z )
. + dli2a * ( 2.D0 )
. + DIMZ * ( 4.D0 + 4.D0*Z )
. + z2 * ( 2.D0 + 2.D0*Z )
. + dli3a * ( 4.D0 - 2.D0*Z )
. + TRIMZ * ( - 24.D0 + 20.D0*Z )
. + S1MZ * ( 16.D0 + 8.D0*Z )
. + SMZ * ( - 48.D0 - 24.D0*Z )
. + z3 * ( - 12.D0 + 18.D0*Z ) - 8.D0 + 8.D0*Z
partIA = cf * partIA
CQQIA = partIA
return
end
double precision function CQQ(xx,lfh,lfr)
implicit none
double precision xx,lfh,lfr,z
double precision xm1,xp1,dlnx,dlxm1,dlxp1
double precision dli2a,dli3a,zcom,s1mz,dimz,trimz,smz
double precision trimco,tripco
double precision part1,part2,scale
double precision z2,z3,z4
double precision ca,cf,nf
double precision dilog,trilog,wgplg
double precision cax,cv,cd,sw
external dilog,trilog,wgplg
sw = 0.2228972225239183d0
cax = 1.d0
cv = -1.d0 + 4.d0/3.d0*sw
cd = 1.d0 + cv**2
ca=3.d0
cf=4.d0/ca
nf=5d0
z2 = 1.6449340668482264365d0
z3 = 1.2020569031595942854d0
z4 = 1.0823232337111381915d0
xm1 = 1.d0 - xx
dlnx = dlog(xx)
dlxm1 = dlog(xm1)
xp1 = 1.0d0 + xx
dlxp1 = dlog(xp1)
dli2a = dilog(xm1)
dli3a = trilog(xm1)
dimz = wgplg(1,1,-xx)
trimz = wgplg(2,1,-xx)
zcom = xm1/xp1
tripco= wgplg(2,1,zcom)
trimco= wgplg(2,1,-zcom)
s1mz = wgplg(1,2,xm1)
smz = wgplg(1,2,-xx)
Z = xx
part1=-8.0d0*dli3a + 12.0d0*s1mz + 8.0d0*dlxm1*dli2a +
. 2.0d0*dlnx*dli2a - 4.0d0*dlnx*z2 + 1.5d0*dlnx**3
. -4.0d0*dlnx**2*dlxm1 + 4.0d0*dlnx*dlxm1**2
part1= xp1 * part1
part2=dli2a * ( 13.0d0 + 8.0d0/3.0d0*z**2 + 5.0d0*z +
. 16.0d0/3.0d0/z )
. + z2 * ( - 2.0d0 + 8.0d0/3.0d0*z**2 + 2.0d0*z
. - 8.0d0/3.0d0/z )
. + dlxm1**2 * ( 2.0d0 - 8.0d0/3.0d0*z**2 - 2.0d0*z +
. 8.0d0/3.0d0/z )
. + dlnx*dlxm1 * ( 6.0d0 + 8.0d0*z**2 + 12.0d0*z )
. + dlnx**2 * ( - 5.0d0/4.0d0 - 10.0d0/3.0d0*z**2
. - 25.0d0/4.0d0*z )
. + dlxm1 * ( - 26.0d0/3.0d0 - 44.0d0/9.0d0*z**2 +
. 26.0d0/3.0d0*z + 44.0d0/9.0d0/z )
. + dlnx * ( 115.0d0/6.0d0 + 10.0d0/9.0d0*z**2
. - 8.0d0/3.0d0*z )
. + 593.0d0/36.0d0 + 703.0d0/108.0d0*z**2
. - 433.0d0/18.0d0*z + 29.0d0/27.0d0/z
scale =lfh**2 * ( dlnx * ( 1.0d0 + z ) +
. 0.5d0 - 2.0d0/3.0d0*z**2 - 0.5d0*z +
. 2.0d0/3.0d0/z ) +
. lfh * ( dli2a * ( 4.0d0 + 4.0d0*z ) +
. dlnx**2 * ( - 2.0d0 - 2.0d0*z ) +
. dlnx*dlxm1 * ( 4.0d0 + 4.0d0*z ) +
. dlnx * ( 3.0d0 + 4.0d0*z**2 + 6.0d0*z ) +
. dlxm1 * ( 2.0d0 - 8.0d0/3.0d0*z**2
. - 2.0d0*z + 8.0d0/3.0d0/z )
. -13.0d0/3.0d0 - 22.0d0/9.0d0*z**2 +
. 13.0d0/3.0d0*z + 22.0d0/9.0d0/z )
CQQ = ( ca**2 - 1.0d0 )/ ca *( part1 + part2 + scale )
return
end
DOUBLE PRECISION FUNCTION CQBBB(xx,lfh,lfr)
implicit none
double precision xx,lfh,lfr,z
double precision logz,s1mz,dimz,trimz,smz
double precision log1pz,di1mz
double precision trimco,tripco
double precision part1,part2,scale
double precision z2,z3,z4
double precision ca,cf,nf
double precision dilog,trilog,wgplg
external dilog,trilog,wgplg
ca=3.d0
cf=4.d0/ca
nf=5d0
z2 = 1.6449340668482264365d0
z3 = 1.2020569031595942854d0
z4 = 1.0823232337111381915d0
Z = xx
LOGZ=DLOG(Z)
LOG1PZ=DLOG(1.0D0+Z)
DIMZ=WGPLG(1,1,-Z)
CQBBB =CF * (
. LOGZ**2 * ( 4.D0/3.D0 + 4.D0/3.D0*Z**2 + 8.D0/3.D0*Z )
. + LOGZ*LOG1PZ * ( - 16.D0/3.D0 - 16.D0/3.D0*Z**2
. - 32.D0/3.D0*Z )
. + LOGZ * ( 4.D0 + 4.D0*Z**2 + 16.D0/3.D0*Z )
. + DIMZ * ( - 16.D0/3.D0 - 16.D0/3.D0*Z**2 - 32.D0/3.D0*Z )
. + Z2 * ( - 8.D0/3.D0 - 8.D0/3.D0*Z**2 - 16.D0/3.D0*Z )
. + 20.D0/3.D0 - 20.D0/3.D0*Z**2 )
RETURN
END
DOUBLE PRECISION FUNCTION CQBBCD(xx,lfh,lfr)
implicit none
double precision xx,lfh,lfr,z
double precision logz,s1mz,dimz,trimz,smz
double precision log1pz,di1mz
double precision trimco,tripco
double precision part1,part2,scale
double precision z2,z3,z4
double precision ca,cf,nf
double precision dilog,trilog,wgplg
external dilog,trilog,wgplg
ca=3.d0
cf=4.d0/ca
nf=5d0
z2 = 1.6449340668482264365d0
z3 = 1.2020569031595942854d0
z4 = 1.0823232337111381915d0
Z = xx
LOGZ=DLOG(Z)
LOG1PZ=DLOG(1.0D0+Z)
DIMZ=WGPLG(1,1,-Z)
DI1MZ=WGPLG(1,1,1.0D0-Z)
TRIMZ=WGPLG(2,1,-Z)
SMZ=WGPLG(1,2,-Z)
S1MZ=WGPLG(1,2,1.0D0-Z)
CQBBCD =CF * ( CF - CA/2 ) * (
. LOGZ**3 * ( 4.D0/3.D0 + 4.D0/3.D0*Z**2 + 16.D0/3.D0*Z )
. + LOGZ**2*LOG1PZ * ( 20.D0 + 20.D0*Z**2 + 40.D0*Z )
. + LOGZ**2 * ( 12.D0 - 30.D0*Z**2 - 16.D0*Z )
. + LOGZ*LOG1PZ**2 * ( - 24.D0 - 24.D0*Z**2 - 48.D0*Z )
. + LOGZ*LOG1PZ * ( 24.D0 + 24.D0*Z**2 + 48.D0*Z )
. + LOGZ*DI1MZ * ( 16.D0 + 16.D0*Z**2 + 48.D0*Z )
. + LOGZ*DIMZ * ( 24.D0 + 24.D0*Z**2 + 48.D0*Z )
. + LOGZ*Z2 * ( 8.D0 + 8.D0*Z**2 + 16.D0*Z )
. + LOGZ * ( 36.D0 + 44.D0*Z )
. + LOG1PZ*DIMZ * ( - 48.D0 - 48.D0*Z**2 - 96.D0*Z )
. + LOG1PZ*Z2 * ( - 24.D0 - 24.D0*Z**2 - 48.D0*Z )
. + DI1MZ * ( 36.D0 - 36.D0*Z**2 )
. + DIMZ * ( 24.D0 + 24.D0*Z**2 + 48.D0*Z )
. + Z2 * ( 12.D0 + 12.D0*Z**2 + 24.D0*Z )
. + TRIMZ * ( - 8.D0 - 8.D0*Z**2 - 16.D0*Z )
. + S1MZ * ( 32.D0 + 32.D0*Z**2 + 96.D0*Z )
. + SMZ * ( - 48.D0 - 48.D0*Z**2 - 96.D0*Z )
. + 54.D0 - 26.D0*Z**2 - 28.D0*Z )
RETURN
END
double precision function CQBAB(xx,lfh,lfr)
implicit none
double precision partBA
double precision xx,lfh,lfr,z
double precision xm1,xp1,dlnx,dlxm1,dlxp1
double precision dli2a,dli3a,zcom,s1mz,dimz,trimz,smz
double precision trimco,tripco
double precision part1,part2,scale
double precision z2,z3,z4
double precision ca,cf,nf
double precision dilog,trilog,wgplg
double precision cax,cv,cd,sw
external dilog,trilog,wgplg
ca=3.d0
cf=4.d0/ca
nf=5d0
z2 = 1.6449340668482264365d0
z3 = 1.2020569031595942854d0
z4 = 1.0823232337111381915d0
xm1 = 1.d0 - xx
dlnx = dlog(xx)
dlxm1 = dlog(xm1)
xp1 = 1.0d0 + xx
dlxp1 = dlog(xp1)
dli2a = dilog(xm1)
dli3a = trilog(xm1)
dimz = wgplg(1,1,-xx)
trimz = wgplg(2,1,-xx)
zcom = xm1/xp1
tripco= wgplg(2,1,zcom)
trimco= wgplg(2,1,-zcom)
s1mz = wgplg(1,2,xm1)
smz = wgplg(1,2,-xx)
Z = xx
partBA = dlnx * ( - 8.d0 + 8.d0*xx + 16.d0/xm1 )
. + 24.d0 - 8.d0*xx
partBA =CF * partBA
CQBAB = partBA
return
end
double precision function deltabba(xx,lfh,lfr)
implicit none
- double precision xx,lfh,lfr,z,lfac,lfren
+ double precision xx,lfh,lfr,z,lfac,lfren,pi
double precision CQQIV,CQQ,CQQIA,CQBBB
double precision CQBBCD,CQQCRF,CQQCRI
double precision cax,cv,cd,sw
external CQQIV,CQQ,CQQIA,CQBBB
external CQBBCD,CQQCRF,CQQCRI
sw = 0.2228972225239183d0
cax = 1.d0
cv = -1.d0 + 4.d0/3.d0*sw
cd = 1.d0 + cv**2
c lfh here is actually -lfh,while lfr = -lfh+lfr
lfren = -lfh+lfr
lfac = -lfh
-
+ pi = 3.1415926535897932385d0
+
deltabba = xx*(cd*CQBBB(xx,lfac,lfren)
? + 2.d0*(cd*(CQQ(xx,lfac,lfren)
? + CQQCRF(xx,lfac,lfren)
? + CQQCRI(xx,lfac,lfren))
? +cax**2*CQQIA(xx,lfac,lfren)
- ? +cv**2 *CQQIV(xx,lfac,lfren)))
+ ? +cv**2 *CQQIV(xx,lfac,lfren)))/(4.d0*pi)**2
return
end
c.......................................................................
double precision function deltabbf(xx,lfh,lfr)
implicit none
double precision xx,lfh,lfr,CQQ
external CQQ
deltabbf = 0.d0
return
end
*******************************************
double precision function deltabqa(xx,lfh,lfr)
implicit none
- double precision xx,lfh,lfr,lfac,lfren
+ double precision xx,lfh,lfr,lfac,lfren,pi
double precision CQQ
external CQQ
c lfh here is actually -lfh,while lfr = -lfh+lfr
lfren = -lfh+lfr
lfac = -lfh
- deltabqa = xx*(CQQ(xx,lfac,lfren) )
+ deltabqa = xx*(CQQ(xx,lfac,lfren) )/(4.d0*pi)**2
end
c.......................................................................
double precision function deltabqf(xx,lfh,lfr)
implicit none
double precision xx,lfh,lfr
deltabqf = 0.d0
end
*******************************************
double precision function deltaqqba(xx,lfh,lfr)
implicit none
- double precision xx,lfh,lfr,lfac,lfren
+ double precision xx,lfh,lfr,lfac,lfren,pi
double precision CQBBB
external CQBBB
c lfh here is actually -lfh,while lfr = -lfh+lfr
lfren = -lfh+lfr
lfac = -lfh
-
- deltaqqba = xx*CQBBB(xx,lfac,lfren)
+ pi = 3.1415926535897932385d0
+
+ deltaqqba = xx*CQBBB(xx,lfac,lfren)/(4.0d0*pi)**2
return
end
c.......................................................................
double precision function deltaqqbf(xx,lfh,lfr)
implicit none
double precision xx,lfh,lfr
deltaqqbf = 0.d0
return
end
c......................................................................
c......................................................................
c................... FUNCTIONS ........................................
c......................................................................
c......................................................................
double precision function dilog(xx)
c..
c.. Dilogarithm: dilog(x) = Li_2(x)
c..
double precision xx
double precision wgplg
dilog = wgplg(1,1,xx)
end
c.. ------------------------------------------------------------
C-}}}
C-{{{ function trilog:
c.. ------------------------------------------------------------
double precision function trilog(xx)
c..
c.. Trilogarithm: trilog(x) = Li_3(x)
c..
double precision xx
double precision wgplg
trilog = wgplg(2,1,xx)
end
c.. ------------------------------------------------------------
C-}}}
C-{{{ FUNCTION WGPLG:
double precision FUNCTION WGPLG(N,P,X)
c..
c.. Nielson's Generalized Polylogarithm
c..
c.. taken from zwprod.f (R. Hamberg, T. Matsuura and W.L. van Neerven)
c.. (originally from CERNLIB?)
c..
INTEGER P,P1,NC(10),INDEX(31)
DOUBLE PRECISION FCT(0:4),SGN(0:4),U(0:4),S1(4,4),C(4,4)
DOUBLE PRECISION A(0:30,10)
DOUBLE PRECISION X,X1,H,ALFA,R,Q,C1,C2,B0,B1,B2,ZERO,HALF
DOUBLE PRECISION V(0:5),SK,SM
DATA FCT /1.0D0,1.0D0,2.0D0,6.0D0,24.0D0/
DATA SGN /1.0D0,-1.0D0,1.0D0,-1.0D0,1.0D0/
DATA ZERO /0.0D0/, HALF /0.5D0/
DATA C1 /1.33333 33333 333D0/, C2 /0.33333 33333 3333D0/
DATA S1(1,1) /1.64493 40668 482D0/
DATA S1(1,2) /1.20205 69031 596D0/
DATA S1(1,3) /1.08232 32337 111D0/
DATA S1(1,4) /1.03692 77551 434D0/
DATA S1(2,1) /1.20205 69031 596D0/
DATA S1(2,2) /2.70580 80842 778D-1/
DATA S1(2,3) /9.65511 59989 444D-2/
DATA S1(3,1) /1.08232 32337 111D0/
DATA S1(3,2) /9.65511 59989 444D-2/
DATA S1(4,1) /1.03692 77551 434D0/
DATA C(1,1) / 1.64493 40668 482D0/
DATA C(1,2) / 1.20205 69031 596D0/
DATA C(1,3) / 1.08232 32337 111D0/
DATA C(1,4) / 1.03692 77551 434D0/
DATA C(2,1) / 0.00000 00000 000D0/
DATA C(2,2) /-1.89406 56589 945D0/
DATA C(2,3) /-3.01423 21054 407D0/
DATA C(3,1) / 1.89406 56589 945D0/
DATA C(3,2) / 3.01423 21054 407D0/
DATA C(4,1) / 0.00000 00000 000D0/
DATA INDEX /1,2,3,4,6*0,5,6,7,7*0,8,9,8*0,10/
DATA NC /24,26,28,30,22,24,26,19,22,17/
DATA A( 0,1) / .96753 21504 3498D0/
DATA A( 1,1) / .16607 30329 2785D0/
DATA A( 2,1) / .02487 93229 2423D0/
DATA A( 3,1) / .00468 63619 5945D0/
DATA A( 4,1) / .00100 16274 9616D0/
DATA A( 5,1) / .00023 20021 9609D0/
DATA A( 6,1) / .00005 68178 2272D0/
DATA A( 7,1) / .00001 44963 0056D0/
DATA A( 8,1) / .00000 38163 2946D0/
DATA A( 9,1) / .00000 10299 0426D0/
DATA A(10,1) / .00000 02835 7538D0/
DATA A(11,1) / .00000 00793 8705D0/
DATA A(12,1) / .00000 00225 3670D0/
DATA A(13,1) / .00000 00064 7434D0/
DATA A(14,1) / .00000 00018 7912D0/
DATA A(15,1) / .00000 00005 5029D0/
DATA A(16,1) / .00000 00001 6242D0/
DATA A(17,1) / .00000 00000 4827D0/
DATA A(18,1) / .00000 00000 1444D0/
DATA A(19,1) / .00000 00000 0434D0/
DATA A(20,1) / .00000 00000 0131D0/
DATA A(21,1) / .00000 00000 0040D0/
DATA A(22,1) / .00000 00000 0012D0/
DATA A(23,1) / .00000 00000 0004D0/
DATA A(24,1) / .00000 00000 0001D0/
DATA A( 0,2) / .95180 88912 7832D0/
DATA A( 1,2) / .43131 13184 6532D0/
DATA A( 2,2) / .10002 25071 4905D0/
DATA A( 3,2) / .02442 41559 5220D0/
DATA A( 4,2) / .00622 51246 3724D0/
DATA A( 5,2) / .00164 07883 1235D0/
DATA A( 6,2) / .00044 40792 0265D0/
DATA A( 7,2) / .00012 27749 4168D0/
DATA A( 8,2) / .00003 45398 1284D0/
DATA A( 9,2) / .00000 98586 9565D0/
DATA A(10,2) / .00000 28485 6995D0/
DATA A(11,2) / .00000 08317 0847D0/
DATA A(12,2) / .00000 02450 3950D0/
DATA A(13,2) / .00000 00727 6496D0/
DATA A(14,2) / .00000 00217 5802D0/
DATA A(15,2) / .00000 00065 4616D0/
DATA A(16,2) / .00000 00019 8033D0/
DATA A(17,2) / .00000 00006 0204D0/
DATA A(18,2) / .00000 00001 8385D0/
DATA A(19,2) / .00000 00000 5637D0/
DATA A(20,2) / .00000 00000 1735D0/
DATA A(21,2) / .00000 00000 0536D0/
DATA A(22,2) / .00000 00000 0166D0/
DATA A(23,2) / .00000 00000 0052D0/
DATA A(24,2) / .00000 00000 0016D0/
DATA A(25,2) / .00000 00000 0005D0/
DATA A(26,2) / .00000 00000 0002D0/
DATA A( 0,3) / .98161 02799 1365D0/
DATA A( 1,3) / .72926 80632 0726D0/
DATA A( 2,3) / .22774 71490 9321D0/
DATA A( 3,3) / .06809 08329 6197D0/
DATA A( 4,3) / .02013 70118 3064D0/
DATA A( 5,3) / .00595 47848 0197D0/
DATA A( 6,3) / .00176 76901 3959D0/
DATA A( 7,3) / .00052 74821 8502D0/
DATA A( 8,3) / .00015 82746 1460D0/
DATA A( 9,3) / .00004 77492 2076D0/
DATA A(10,3) / .00001 44792 0408D0/
DATA A(11,3) / .00000 44115 4886D0/
DATA A(12,3) / .00000 13500 3870D0/
DATA A(13,3) / .00000 04148 1779D0/
DATA A(14,3) / .00000 01279 3307D0/
DATA A(15,3) / .00000 00395 9070D0/
DATA A(16,3) / .00000 00122 9055D0/
DATA A(17,3) / .00000 00038 2658D0/
DATA A(18,3) / .00000 00011 9459D0/
DATA A(19,3) / .00000 00003 7386D0/
DATA A(20,3) / .00000 00001 1727D0/
DATA A(21,3) / .00000 00000 3687D0/
DATA A(22,3) / .00000 00000 1161D0/
DATA A(23,3) / .00000 00000 0366D0/
DATA A(24,3) / .00000 00000 0116D0/
DATA A(25,3) / .00000 00000 0037D0/
DATA A(26,3) / .00000 00000 0012D0/
DATA A(27,3) / .00000 00000 0004D0/
DATA A(28,3) / .00000 00000 0001D0/
DATA A( 0,4) /1.06405 21184 614 D0/
DATA A( 1,4) /1.06917 20744 981 D0/
DATA A( 2,4) / .41527 19325 1768D0/
DATA A( 3,4) / .14610 33293 6222D0/
DATA A( 4,4) / .04904 73264 8784D0/
DATA A( 5,4) / .01606 34086 0396D0/
DATA A( 6,4) / .00518 88935 0790D0/
DATA A( 7,4) / .00166 29871 7324D0/
DATA A( 8,4) / .00053 05827 9969D0/
DATA A( 9,4) / .00016 88702 9251D0/
DATA A(10,4) / .00005 36832 8059D0/
DATA A(11,4) / .00001 70592 3313D0/
DATA A(12,4) / .00000 54217 4374D0/
DATA A(13,4) / .00000 17239 4082D0/
DATA A(14,4) / .00000 05485 3275D0/
DATA A(15,4) / .00000 01746 7795D0/
DATA A(16,4) / .00000 00556 7550D0/
DATA A(17,4) / .00000 00177 6234D0/
DATA A(18,4) / .00000 00056 7224D0/
DATA A(19,4) / .00000 00018 1313D0/
DATA A(20,4) / .00000 00005 8012D0/
DATA A(21,4) / .00000 00001 8579D0/
DATA A(22,4) / .00000 00000 5955D0/
DATA A(23,4) / .00000 00000 1911D0/
DATA A(24,4) / .00000 00000 0614D0/
DATA A(25,4) / .00000 00000 0197D0/
DATA A(26,4) / .00000 00000 0063D0/
DATA A(27,4) / .00000 00000 0020D0/
DATA A(28,4) / .00000 00000 0007D0/
DATA A(29,4) / .00000 00000 0002D0/
DATA A(30,4) / .00000 00000 0001D0/
DATA A( 0,5) / .97920 86066 9175D0/
DATA A( 1,5) / .08518 81314 8683D0/
DATA A( 2,5) / .00855 98522 2013D0/
DATA A( 3,5) / .00121 17721 4413D0/
DATA A( 4,5) / .00020 72276 8531D0/
DATA A( 5,5) / .00003 99695 8691D0/
DATA A( 6,5) / .00000 83806 4065D0/
DATA A( 7,5) / .00000 18684 8945D0/
DATA A( 8,5) / .00000 04366 6087D0/
DATA A( 9,5) / .00000 01059 1733D0/
DATA A(10,5) / .00000 00264 7892D0/
DATA A(11,5) / .00000 00067 8700D0/
DATA A(12,5) / .00000 00017 7654D0/
DATA A(13,5) / .00000 00004 7342D0/
DATA A(14,5) / .00000 00001 2812D0/
DATA A(15,5) / .00000 00000 3514D0/
DATA A(16,5) / .00000 00000 0975D0/
DATA A(17,5) / .00000 00000 0274D0/
DATA A(18,5) / .00000 00000 0077D0/
DATA A(19,5) / .00000 00000 0022D0/
DATA A(20,5) / .00000 00000 0006D0/
DATA A(21,5) / .00000 00000 0002D0/
DATA A(22,5) / .00000 00000 0001D0/
DATA A( 0,6) / .95021 85196 3952D0/
DATA A( 1,6) / .29052 52916 1433D0/
DATA A( 2,6) / .05081 77406 1716D0/
DATA A( 3,6) / .00995 54376 7280D0/
DATA A( 4,6) / .00211 73389 5031D0/
DATA A( 5,6) / .00047 85947 0550D0/
DATA A( 6,6) / .00011 33432 1308D0/
DATA A( 7,6) / .00002 78473 3104D0/
DATA A( 8,6) / .00000 70478 8108D0/
DATA A( 9,6) / .00000 18278 8740D0/
DATA A(10,6) / .00000 04838 7492D0/
DATA A(11,6) / .00000 01303 3842D0/
DATA A(12,6) / .00000 00356 3769D0/
DATA A(13,6) / .00000 00098 7174D0/
DATA A(14,6) / .00000 00027 6586D0/
DATA A(15,6) / .00000 00007 8279D0/
DATA A(16,6) / .00000 00002 2354D0/
DATA A(17,6) / .00000 00000 6435D0/
DATA A(18,6) / .00000 00000 1866D0/
DATA A(19,6) / .00000 00000 0545D0/
DATA A(20,6) / .00000 00000 0160D0/
DATA A(21,6) / .00000 00000 0047D0/
DATA A(22,6) / .00000 00000 0014D0/
DATA A(23,6) / .00000 00000 0004D0/
DATA A(24,6) / .00000 00000 0001D0/
DATA A( 0,7) / .95064 03218 6777D0/
DATA A( 1,7) / .54138 28546 5171D0/
DATA A( 2,7) / .13649 97959 0321D0/
DATA A( 3,7) / .03417 94232 8207D0/
DATA A( 4,7) / .00869 02788 3583D0/
DATA A( 5,7) / .00225 28408 4155D0/
DATA A( 6,7) / .00059 51608 9806D0/
DATA A( 7,7) / .00015 99561 7766D0/
DATA A( 8,7) / .00004 36521 3096D0/
DATA A( 9,7) / .00001 20747 4688D0/
DATA A(10,7) / .00000 33801 8176D0/
DATA A(11,7) / .00000 09563 2476D0/
DATA A(12,7) / .00000 02731 3129D0/
DATA A(13,7) / .00000 00786 6968D0/
DATA A(14,7) / .00000 00228 3195D0/
DATA A(15,7) / .00000 00066 7205D0/
DATA A(16,7) / .00000 00019 6191D0/
DATA A(17,7) / .00000 00005 8018D0/
DATA A(18,7) / .00000 00001 7246D0/
DATA A(19,7) / .00000 00000 5151D0/
DATA A(20,7) / .00000 00000 1545D0/
DATA A(21,7) / .00000 00000 0465D0/
DATA A(22,7) / .00000 00000 0141D0/
DATA A(23,7) / .00000 00000 0043D0/
DATA A(24,7) / .00000 00000 0013D0/
DATA A(25,7) / .00000 00000 0004D0/
DATA A(26,7) / .00000 00000 0001D0/
DATA A( 0,8) / .98800 01167 2229D0/
DATA A( 1,8) / .04364 06760 9601D0/
DATA A( 2,8) / .00295 09117 8278D0/
DATA A( 3,8) / .00031 47780 9720D0/
DATA A( 4,8) / .00004 31484 6029D0/
DATA A( 5,8) / .00000 69381 8230D0/
DATA A( 6,8) / .00000 12464 0350D0/
DATA A( 7,8) / .00000 02429 3628D0/
DATA A( 8,8) / .00000 00504 0827D0/
DATA A( 9,8) / .00000 00109 9075D0/
DATA A(10,8) / .00000 00024 9467D0/
DATA A(11,8) / .00000 00005 8540D0/
DATA A(12,8) / .00000 00001 4127D0/
DATA A(13,8) / .00000 00000 3492D0/
DATA A(14,8) / .00000 00000 0881D0/
DATA A(15,8) / .00000 00000 0226D0/
DATA A(16,8) / .00000 00000 0059D0/
DATA A(17,8) / .00000 00000 0016D0/
DATA A(18,8) / .00000 00000 0004D0/
DATA A(19,8) / .00000 00000 0001D0/
DATA A( 0,9) / .95768 50654 6350D0/
DATA A( 1,9) / .19725 24967 9534D0/
DATA A( 2,9) / .02603 37031 3918D0/
DATA A( 3,9) / .00409 38216 8261D0/
DATA A( 4,9) / .00072 68170 7110D0/
DATA A( 5,9) / .00014 09187 9261D0/
DATA A( 6,9) / .00002 92045 8914D0/
DATA A( 7,9) / .00000 63763 1144D0/
DATA A( 8,9) / .00000 14516 7850D0/
DATA A( 9,9) / .00000 03420 5281D0/
DATA A(10,9) / .00000 00829 4302D0/
DATA A(11,9) / .00000 00206 0784D0/
DATA A(12,9) / .00000 00052 2823D0/
DATA A(13,9) / .00000 00013 5066D0/
DATA A(14,9) / .00000 00003 5451D0/
DATA A(15,9) / .00000 00000 9436D0/
DATA A(16,9) / .00000 00000 2543D0/
DATA A(17,9) / .00000 00000 0693D0/
DATA A(18,9) / .00000 00000 0191D0/
DATA A(19,9) / .00000 00000 0053D0/
DATA A(20,9) / .00000 00000 0015D0/
DATA A(21,9) / .00000 00000 0004D0/
DATA A(22,9) / .00000 00000 0001D0/
DATA A( 0,10) / .99343 65167 1347D0/
DATA A( 1,10) / .02225 77012 6826D0/
DATA A( 2,10) / .00101 47557 4703D0/
DATA A( 3,10) / .00008 17515 6250D0/
DATA A( 4,10) / .00000 89997 3547D0/
DATA A( 5,10) / .00000 12082 3987D0/
DATA A( 6,10) / .00000 01861 6913D0/
DATA A( 7,10) / .00000 00317 4723D0/
DATA A( 8,10) / .00000 00058 5215D0/
DATA A( 9,10) / .00000 00011 4739D0/
DATA A(10,10) / .00000 00002 3652D0/
DATA A(11,10) / .00000 00000 5082D0/
DATA A(12,10) / .00000 00000 1131D0/
DATA A(13,10) / .00000 00000 0259D0/
DATA A(14,10) / .00000 00000 0061D0/
DATA A(15,10) / .00000 00000 0015D0/
DATA A(16,10) / .00000 00000 0004D0/
DATA A(17,10) / .00000 00000 0001D0/
IF(N .LT. 1 .OR. N .GT. 4 .OR. P .LT. 1 .OR. P .GT. 4 .OR.
1 N+P .GT. 5) THEN
WGPLG=ZERO
PRINT 1000, N,P
RETURN
END IF
IF(X .EQ. SGN(0)) THEN
WGPLG=S1(N,P)
RETURN
END IF
IF(X .GT. FCT(2) .OR. X .LT. SGN(1)) THEN
X1=SGN(0)/X
H=C1*X1+C2
ALFA=H+H
V(0)=SGN(0)
V(1)=LOG(DCMPLX(-X,ZERO))
DO 33 L = 2,N+P
33 V(L)=V(1)*V(L-1)/L
SK=ZERO
DO 34 K = 0,P-1
P1=P-K
R=X1**P1/(FCT(P1)*FCT(N-1))
SM=ZERO
DO 35 M = 0,K
N1=N+K-M
L=INDEX(10*N1+P1-10)
B1=ZERO
B2=ZERO
DO 31 I = NC(L),0,-1
B0=A(I,L)+ALFA*B1-B2
B2=B1
31 B1=B0
Q=(FCT(N1-1)/FCT(K-M))*(B0-H*B2)*R/P1**N1
35 SM=SM+V(M)*Q
34 SK=SK+SGN(K)*SM
SM=ZERO
DO 36 M = 0,N-1
36 SM=SM+V(M)*C(N-M,P)
WGPLG=SGN(N)*SK+SGN(P)*(SM+V(N+P))
RETURN
END IF
IF(X .GT. HALF) THEN
X1=SGN(0)-X
H=C1*X1+C2
ALFA=H+H
V(0)=SGN(0)
U(0)=SGN(0)
V(1)=LOG(DCMPLX(X1,ZERO))
U(1)=LOG(X)
DO 23 L = 2,P
23 V(L)=V(1)*V(L-1)/L
DO 26 L = 2,N
26 U(L)=U(1)*U(L-1)/L
SK=ZERO
DO 24 K = 0,N-1
P1=N-K
R=X1**P1/FCT(P1)
SM=ZERO
DO 25 M = 0,P-1
N1=P-M
L=INDEX(10*N1+P1-10)
B1=ZERO
B2=ZERO
DO 12 I = NC(L),0,-1
B0=A(I,L)+ALFA*B1-B2
B2=B1
12 B1=B0
Q=SGN(M)*(B0-H*B2)*R/P1**N1
25 SM=SM+V(M)*Q
24 SK=SK+U(K)*(S1(P1,P)-SM)
WGPLG=SK+SGN(P)*U(N)*V(P)
RETURN
END IF
L=INDEX(10*N+P-10)
H=C1*X+C2
ALFA=H+H
B1=ZERO
B2=ZERO
DO 11 I = NC(L),0,-1
B0=A(I,L)+ALFA*B1-B2
B2=B1
11 B1=B0
WGPLG=(B0-H*B2)*X**P/(FCT(P)*P**N)
RETURN
1000 FORMAT(/' ***** CERN SUBROUTINE WGPLG ... ILLEGAL VALUES',
1 ' N = ',I3,' P = ',I3)
END
C-}}}
c..
c.. Coefficient of the delta-function at NNLO.
c..
double precision function delta2(lfr,lfh)
implicit none
double precision lfh,lfr,lfac,lfren
- double precision z2,z3
+ double precision z2,z3,pi
double precision consCF,consCA,consNF
double precision cf,ca,nf
cf = 4.d0/3.d0
ca = 3.d0
nf = 5.d0
z2 = 1.6449340668482264365d0
z3 = 1.2020569031595942854d0
lfren = -lfh+lfr
- lfac = -lfh
-
+ lfac = -lfh
+
+ pi = 3.1415926535897932385d0
+
consCF = 8.0d0/5.0d0 * z2**2 - 60.0d0 * z3 - 70.0d0 * z2 +
. 511.0d0/4.0d0
consCF = consCF + (18.0d0 - 32.0d0 * z2 ) * lfac**2 +
. ( -93.0d0 + 24.0d0 * z2 + 176.0d0 * z3 ) * lfac
consCA = -1535.0d0/12.0d0 + 592.0d0/9.0d0 * z2 + 28.0d0 * z3
. -12.0d0/5.0d0 * z2**2
consCA = consCA + 11.0d0*lfac**2 - 22.0d0*lfac*lfren +
. ( 17.0d0/3.0d0 + 88.0d0/3.0d0 * z2 - 24.0d0 * z3)*lfac +
. ( 176.0d0/3.0d0 - 88.0d0/3.0d0 * z2 )*lfren
consNF = 8.0d0 * z3 - 112.0d0/9.0d0 * z2 + 127.0d0/6.0d0
consNF = consNF - 2.0d0*lfac**2 + 4.0D0*lfac*lfren
. - ( 2.0d0/3.0d0 + 16.0d0/3.0d0 * z2 ) * lfac +
. (-32.0d0/3.0d0 + 16.0d0/3.0d0 * z2 ) * lfren
delta2 = cf**2 * consCF + cf * ca * consCA + nf*cf * consNF
-
+ delta2 = delta2/(4.d0*pi)**2
return
end
c..
c.. The plus-distributions at NNLO.
c..
double precision function dterms2(xt,lfh,lfr)
implicit none
- double precision cf,nf,ca
+ double precision cf,nf,ca,pi
double precision partCF,partCA,partNF
double precision xt,lfh,lfr,lfac,lfren
double precision log1mz,dd0,dd1,dd2,dd3
double precision z2,z3,z4,z
z=xt
cf = 4.d0/3.d0
nf = 5.d0
ca = 3.d0
z2 = 1.6449340668482264365d0
z3 = 1.2020569031595942854d0
z4 = 1.0823232337111381915d0
lfren = -lfh+lfr
lfac = -lfh
-
+ pi = 3.1415926535897932385d0
log1mz = dlog(1.0d0-z)
dd0 = 1.0d0/(1.0d0-z)
dd1 = dd0 * log1mz
dd2 = dd0 * log1mz**2
dd3 = dd0 * log1mz**3
partCF = 128.0d0 * dd3 - (128.0d0*z2 + 256.0d0 )*dd1 +
. 256.0d0 * z3 * dd0
partCF = partCF + 192.0d0 * lfac * dd2 +
. ( 64.0d0 * lfac + 96.0d0 ) * lfac * dd1 +
. ( 48.0D0 * lfac - 128.0d0 - 64.0d0*z2 ) * lfac * dd0
partCF = cf**2 * partCF
c--- CA
partCA = -176.0d0/3.0d0 * dd2 + ( 1072.0d0/9.0d0 - 32.0d0*z2)*dd1+
. ( 56.0d0*z3 + 176.0d0/3.0d0*z2 -1616.0d0/27.0d0 )*dd0
partCA = partCA - 176.0d0/3.0d0*lfren*dd1 +
. ( 44.0d0/3.0d0*lfac**2 - 88.0d0/3.0d0*lfac*lfren +
. ( 536.0d0/9.0d0 - 16.0d0*z2 )*lfac ) * dd0
partCA = ca*cf * partCA
c -- NF
partNF = 32.0d0/3.0d0 * dd2 - 160.0d0/9.0d0 * dd1 +
. ( 224.0d0/27.0d0 - 32.0d0/3.0d0*z2 )*dd0
partNF = partNF + 32.0d0/3.0d0*lfren * dd1 +
. ( -8.0d0/3.0d0*lfac**2 + 16.0d0/3.0d0*lfac*lfren
. -80.0d0/9.0d0*lfac ) * dd0
partNF = nf*cf * partNF
dterms2 = partCF + partCA + partNF
-
+ dterms2 = dterms2/(4.d0*pi)**2
return
end
c.......................................................................
c............ Unchanged from bbh .......................................
c.......................................................................
double precision function dtsub2(lfh,lfr,tauh)
c..
c.. Contributions arising from the fact that the integrals over
c.. plus-distributions do not run from 0 to 1, but from z to 1.
c..
implicit double precision (a-z)
dtsub2 = 0.
end
C........................Running mass (Marius) ....................
double precision function runmass(mass0,api0,apif,nf,nloop)
c..
c.. evaluates the running of the MS-bar quark mass
c.. by expanding the equation
c..
c.. m(mu) = m(mu0) * exp( \int_a0^af dx gammam(x)/x/beta(x) )
c..
c.. in terms of alpha_s. The results agree with RunDec.m.
c..
c..
c.. Input:
c.. ------
c.. mass0 : m(mu0)
c.. api0 : alpha_s(mu0)/pi
c.. apif : alpha_s(muf)/pi
c.. nf : number of flavors
c.. nloop : order of calculation (nloop=1..4)
c..
c.. Output:
c.. -------
c.. massout: m(muf)
c..
implicit real*8 (a-h,o-z)
real*8 mass0,massfun
real*8 nf
external massfun
parameter(accmass=1.d-6)
data z3/1.2020569031595942853997/,
& z5/1.0369277551433699263/,
& pi/3.1415926535897932381/
beta0 = (33 - 2*nf)/12.d0
beta1 = (102 - (38*nf)/3.d0)/16.d0
beta2 = (2857/2.d0 - (5033*nf)/18.d0 + (325*nf**2)/54.d0)/64.d0
beta3 = (149753/6.d0 + (1093*nf**3)/729.d0 + 3564*z3 + nf**2
& *(50065/162.d0 + (6472*z3)/81.d0) - nf*(1078361/162.d0 +
& (6508*z3)/27.d0))/256.d0
gamma0 = 1.d0
gamma1 = (67.33333333333333d0 - (20*nf)/9.d0)/16.d0
gamma2 = (1249.d0 - (140*nf**2)/81.d0 + 2*nf*(-20.59259259259259d0
& - 48*z3) +(8*nf*(-46 + 48*z3))/9.d0)/64.d0
gamma3 = (28413.91975308642d0 + (135680*z3)/27.d0 + nf**3*(-1
& .3662551440329218d0 + (64*z3)/27.d0) + nf**2*(21
& .57201646090535d0 - (16*Pi**4)/27.d0 + (800*z3)/9.d0) - 8800
& *z5 + nf*(-3397.1481481481483d0 + (88*Pi**4)/9.d0 - (34192
& *z3)/9.d0 + (18400*z5)/9.d0))/256.d0
bb1 = beta1/beta0
bb2 = beta2/beta0
bb3 = beta3/beta0
cc0 = gamma0/beta0
cc1 = gamma1/beta0
cc2 = gamma2/beta0
cc3 = gamma3/beta0
cfunc1 = 1.d0
cfunc2 = cc1 - bb1*cc0
cfunc3 = 1/2.d0*((cc1-bb1*cc0)**2 + cc2 - bb1*cc1 + bb1**2*cc0 -
& bb2*cc0)
cfunc4 = (1/6*(cc1 - bb1*cc0)**3 + 1/2*(cc1 - bb1*cc0)*(cc2 - bb1
& *cc1 + bb1**2*cc0 - bb2*cc0) + 1/3*(cc3 - bb1*cc2 + bb1**2
& *cc1 - bb2*cc1 - bb1**3*cc0 + 2*bb1*bb2*cc0 - bb3*cc0))
if (nloop.lt.4) then
cfunc4 = 0.d0
if (nloop.lt.3) then
cfunc3 = 0.d0
if (nloop.lt.2) then
cfunc2 = 0.d0
if (nloop.lt.1) then
cfunc1 = 0.d0
endif
endif
endif
endif
cfuncmu0 = cfunc1 + cfunc2*api0 + cfunc3*api0**2 + cfunc4*api0**3
cfuncmuf = cfunc1 + cfunc2*apif + cfunc3*apif**2 + cfunc4*apif**3
runmass = mass0*(apif/api0)**cc0*cfuncmuf/cfuncmu0
end
DOUBLE PRECISION FUNCTION CQBCAR(xx,lfh,lfr)
implicit none
double precision xx,lfh,lfr,z,help,zmin
double precision logz,log1mz,di1mz,s1mz,dimz,trimz,smz
double precision log1pz,tri1mz
double precision trimco,tripco
double precision part1,part2,scale
double precision z2,z3,z4
double precision ca,cf,nf
double precision dilog,trilog,wgplg
external dilog,trilog,wgplg
ca=3.d0
cf=4.d0/ca
nf=5d0
z2 = 1.6449340668482264365d0
z3 = 1.2020569031595942854d0
z4 = 1.0823232337111381915d0
Z = xx
C.
LOGZ=DLOG(Z)
ZMIN=1.0D0-Z
LOG1MZ=DLOG(ZMIN)
DI1MZ=WGPLG(1,1,ZMIN)
TRI1MZ=WGPLG(2,1,ZMIN)
S1MZ=WGPLG(1,2,ZMIN)
HELP=( 1.0D0 + Z ) *
. ( 20.0D0*S1MZ - 28.0D0*z3 +
. 16.0D0*LOG1MZ*DI1MZ - 8.0D0*LOGZ*DI1MZ +
. 16.0D0*LOG1MZ*z2 - 8.0D0*LOGZ*z2 +
. 88.0D0/3.0D0*LOG1MZ**2 ) +
. (-8.0D0*S1MZ - 24.0D0*TRI1MZ
. -16.0D0*LOG1MZ*DI1MZ + 16.0D0*LOGZ*DI1MZ +
. 16.0D0*LOGZ*z2 + 8.0D0/3.0D0*DI1MZ +
. 280.0D0/3.0D0*LOG1MZ*LOGZ - 29.0D0*LOGZ**2
. -208.0D0/3.0D0*LOGZ )/ZMIN
. -DI1MZ*( 32.0D0/3.0D0 + 56.0D0/3.0D0*Z )
. -z2*( 76.0D0/3.0D0 + 100.0D0/3.0D0*Z )
. -LOG1MZ*LOGZ*( 176.0D0/3.0D0 + 176.0D0/3.0D0*Z ) +
. LOGZ**2*( 55.0D0/3.0D0 + 55.0D0/3.0D0*Z )
. -LOG1MZ*( 152.0D0/9.0D0 + 956.0D0/9.0D0*Z ) +
. LOGZ*( 52.0D0/3.0D0 + 218.0D0/3.0D0*Z )
. -446.0D0/27.0D0 + 2278.0D0/27.0D0*Z
HELP = HELP +
. lfh**2*( - 22.D0/3.D0 - 22.D0/3.D0*Z ) +
. lfh*lfr * ( 44.D0/3.D0 + 44.D0/3.D0*Z ) +
. lfh *( LOGZ * ( - 44.D0/3.D0 - 44.D0/3.D0*Z +
. 52.D0/3.D0/ZMIN )
. + DI1MZ * ( 8.D0 + 8.D0*Z - 16.D0/ZMIN )
. + z2 * ( 8.D0 + 8.D0*Z )
. - 76/9.D0 - 496.D0/9.D0*Z ) +
. lfr *( LOGZ * ( - 44.D0/3.D0 - 44.D0/3.D0*Z +
. 88.D0/3.D0/ZMIN )
. + LOG1MZ * ( 88.D0/3.D0 + 88.D0/3.D0*Z ) )
CQBCAR=CA*CF*HELP
RETURN
END
DOUBLE PRECISION FUNCTION CQBNF(xx,lfh,lfr)
implicit none
double precision xx,lfh,lfr,z,help
double precision logz,log1mz,s1mz,dimz,trimz,smz
double precision d0,d1,d2
double precision log1pz,di1mz
double precision trimco,tripco
double precision part1,part2,scale
double precision z2,z3,z4
double precision ca,cf,nf
double precision dilog,trilog,wgplg
external dilog,trilog,wgplg
ca=3.d0
cf=4.d0/ca
nf=5d0
z2 = 1.6449340668482264365d0
z3 = 1.2020569031595942854d0
z4 = 1.0823232337111381915d0
Z = xx
LOG1MZ=DLOG(1.0D0-Z)
D0=1.0D0/(1.0D0-Z)
D1=D0*LOG1MZ
D2=D1*LOG1MZ
HELP=32.0D0/3.0D0*D2 - 160.0D0/9.0D0*D1 +
. ( 224.0D0/27.0D0 - 32.0D0/3.0D0*z2 )*D0
HELP=HELP + 32.0D0/3.0D0*lfr*D1 +
. ( -8.0D0/3.0D0*lfh**2 + 16.0D0/3.0D0*lfh*lfr
. -80.0D0/9.0D0*lfh )*D0
CQBNF=NF*CF * HELP
RETURN
END
DOUBLE PRECISION FUNCTION CQBNFR(xx,lfh,lfr)
implicit none
double precision xx,lfh,lfr,z,zmin,help
double precision logz,s1mz,dimz,trimz,smz,log1mz
double precision log1pz,di1mz
double precision trimco,tripco
double precision part1,part2,scale
double precision z2,z3,z4
double precision ca,cf,nf
double precision dilog,trilog,wgplg
external dilog,trilog,wgplg
ca=3.d0
cf=4.d0/ca
nf=5d0
z2 = 1.6449340668482264365d0
z3 = 1.2020569031595942854d0
z4 = 1.0823232337111381915d0
Z = xx
LOGZ=DLOG(Z)
ZMIN=1.0D0-Z
LOG1MZ=DLOG(ZMIN)
DI1MZ=WGPLG(1,1,ZMIN)
HELP=( 1.0D0 + Z ) * ( 16.0D0/3.0D0*z2 + 8.0D0/3.0D0*DI1MZ -
. 10.0D0/3.0D0*LOGZ**2 + 32.0D0/3.0D0* LOG1MZ*LOGZ -
. 16.0D0/3.0D0*LOG1MZ**2 ) +
. ( -64.0D0/3.0D0*LOG1MZ*LOGZ - 8.0D0/3.0D0*DI1MZ +
. 8.0D0*LOGZ**2 + 40.0D0/3.0D0*LOGZ ) / ZMIN +
. LOG1MZ * ( - 16.0D0/9.0D0 + 176.0D0/9.0D0*Z ) +
. LOGZ * ( - 4.0D0/3.0D0 - 44.0D0/3.0D0*Z ) +
. 188.0D0/27.0D0 - 412.0D0/27.0D0*Z
HELP = HELP +
. lfh**2*( 4.D0/3.D0 + 4.D0/3.D0*Z ) +
. lfh*lfr * ( - 8.D0/3.D0 - 8.D0/3.D0*Z ) +
. lfh *( LOGZ * ( 8.D0/3.D0 + 8.D0/3.D0*Z
. - 16.D0/3.D0/ZMIN )
. - 8.D0/9.D0 + 88.D0/9.D0*Z ) +
. lfr *( LOGZ * ( 8.D0/3.D0 + 8.D0/3.D0*Z
. - 16.D0/3.D0/ZMIN )
. + LOG1MZ * ( - 16.D0/3.D0 - 16.D0/3.D0*Z ) )
CQBNFR=NF*CF*HELP
RETURN
END
DOUBLE PRECISION FUNCTION CQBCFR(xx,lfh,lfr)
implicit none
double precision xx,lfh,lfr,z,zmin,help
double precision logz,s1mz,dimz,trimz,smz,log1mz
double precision log1pz,di1mz,tri1mz
double precision trimco,tripco
double precision part1,part2,scale
double precision z2,z3,z4
double precision ca,cf,nf
double precision dilog,trilog,wgplg
external dilog,trilog,wgplg
ca=3.d0
cf=4.d0/ca
nf=5d0
z2 = 1.6449340668482264365d0
z3 = 1.2020569031595942854d0
z4 = 1.0823232337111381915d0
Z = xx
LOGZ=DLOG(Z)
ZMIN=1.0D0-Z
LOG1MZ=DLOG(ZMIN)
DI1MZ=WGPLG(1,1,ZMIN)
TRI1MZ=WGPLG(2,1,ZMIN)
S1MZ=WGPLG(1,2,ZMIN)
HELP=( 1.0D0 + Z ) *
. ( 48.0D0*S1MZ - 32.0D0*TRI1MZ - 128.0D0*z3 +
. 24.0D0*LOG1MZ*DI1MZ + 24.0D0*LOGZ*DI1MZ +
. 64.0D0*LOG1MZ*z2 -96.0D0*LOGZ*z2
. -64.0D0*LOG1MZ**3 + 156.0D0*LOG1MZ**2*LOGZ
. -96.0D0*LOG1MZ*LOGZ**2 + 50.0D0/3.0D0*LOGZ**3 ) +
. (-64.0D0*S1MZ - 16.0D0*TRI1MZ +
. 48.0D0*LOG1MZ*DI1MZ - 48.0D0*LOGZ*DI1MZ +
. 128.0D0*LOGZ*z2
. -248.0D0*LOG1MZ**2*LOGZ + 144.0D0*LOG1MZ*LOGZ**2
. -24.0D0*LOGZ**3 + 112.0D0*LOGZ )/ZMIN +
. DI1MZ*( -24.0D0 - 32.0D0*Z ) + z2*( 64.0D0 - 64.0D0*Z ) +
. LOG1MZ**2*( - 64.0D0 + 64.0D0*Z ) +
. LOG1MZ*LOGZ*( 64.0D0 - 112.0D0*Z ) +
. LOGZ**2*( - 8.0D0 + 24.0D0*Z ) +
. LOG1MZ*( 256.0D0 + 12.0D0*Z ) +
. LOGZ*( - 104.0D0 + 48.0D0*Z ) - 72.0D0 + 48.0D0*Z
HELP = HELP +
. LFH**2*( LOGZ * ( 24.D0 + 24.D0*Z - 32.D0/ZMIN )
. + LOG1MZ * ( - 32.D0 - 32.D0*Z )
. - 40.D0 - 8.D0*Z ) +
. LFH *( LOGZ**2 * ( - 36.D0 - 36.D0*Z + 48.D0/ZMIN )
. + LOGZ*LOG1MZ * ( 144.D0 + 144.D0*Z - 224.D0/ZMIN )
. + LOGZ * ( 56.D0 - 24.D0*Z - 48.D0/ZMIN )
. + LOG1MZ**2 * ( - 96.D0 - 96.D0*Z )
. + LOG1MZ * ( - 112.D0 + 16.D0*Z )
. + DI1MZ * ( 16.D0 + 16.D0*Z + 32.D0/ZMIN )
. + z2 * ( 32.D0 + 32.D0*Z )
. + 120.D0 + 16.D0*Z )
CQBCFR=CF*CF*HELP
RETURN
END
DOUBLE PRECISION FUNCTION CQBACD(xx,lfh,lfr)
implicit none
double precision xx,lfh,lfr,z,zmin,help
double precision logz,s1mz,dimz,trimz,smz,log1mz
double precision log1pz,di1mz,tri1mz
double precision trimco,tripco
double precision part1,part2,scale
double precision z2,z3,z4
double precision ca,cf,nf
double precision delacd,frterm
double precision dilog,trilog,wgplg
external dilog,trilog,wgplg
ca=3.d0
cf=4.d0/ca
nf=5d0
z2 = 1.6449340668482264365d0
z3 = 1.2020569031595942854d0
z4 = 1.0823232337111381915d0
Z = xx
LOGZ=DLOG(Z)
ZMIN=1.0D0-Z
LOG1MZ=DLOG(ZMIN)
DI1MZ=WGPLG(1,1,ZMIN)
TRI1MZ=WGPLG(2,1,ZMIN)
S1MZ=WGPLG(1,2,ZMIN)
DELACD=CF * ( CF - CA/2 ) * (
. ( LOGZ**3 * 16.D0/3.D0 - LOGZ**2*LOG1MZ * 16.D0
. + LOGZ**2 * 15.D0 - LOGZ*LOG1MZ * 24.D0
. - LOGZ*DI1MZ * 24.D0 + LOGZ * 24.D0
. - LOG1MZ*DI1MZ * 32.D0 - DI1MZ * 12.D0
. + TRI1MZ * 32.D0 - S1MZ * 72.D0 )/ZMIN
. + LOGZ**3 * ( - 2.D0 - 2.D0*Z )
. + LOGZ**2*LOG1MZ * ( 8.D0 + 8.D0*Z )
. + LOGZ**2 * ( 4.D0 + 4.D0*Z )
. + LOGZ*LOG1MZ * ( - 16.D0 - 16.D0*Z )
. + LOGZ*DI1MZ * ( 16.D0 + 16.D0*Z )
. + LOGZ * ( 32.D0 - 30.D0*Z )
. + LOG1MZ*DI1MZ * ( 16.D0 + 16.D0*Z )
. + LOG1MZ * ( - 64.D0 + 56.D0*Z )
. + DI1MZ * ( - 20.D0 - 20.D0*Z )
. + TRI1MZ * ( - 24.D0 - 24.D0*Z )
. + S1MZ * ( 36.D0 + 36.D0*Z )
. + 94.D0 - 78.D0*Z )
FRTERM=CF * ( CF - CA/2 ) * lfh * (
. LOGZ**2 * ( 4.D0 + 4.D0*Z - 8.D0/ZMIN )
. + LOGZ * ( - 8.D0 - 8.D0*Z - 12.D0/ZMIN )
. + DI1MZ * ( 8.D0 + 8.D0*Z - 16.D0/ZMIN )
. - 32.D0 + 28.D0*Z )
CQBACD=DELACD+FRTERM
RETURN
END
double precision function CQQCRI(xx,lfh,lfr)
c--- CQQCRI
implicit none
double precision partCRI
double precision xx,lfh,lfr,z
double precision xm1,xp1,dlnx,dlxm1,dlxp1
double precision dli2a,dli3a,zcom,s1mz,dimz,trimz,smz
double precision trimco,tripco
double precision part1,part2,scale
double precision z2,z3,z4
double precision ca,cf,nf
double precision dilog,trilog,wgplg
double precision cax,cv,cd,sw
external dilog,trilog,wgplg
sw = 0.2228972225239183d0
cax = 1.d0
cv = -1.d0 + 4.d0/3.d0*sw
cd = 1.d0 + cv**2
ca=3.d0
cf=4.d0/ca
nf=5d0
z2 = 1.6449340668482264365d0
z3 = 1.2020569031595942854d0
z4 = 1.0823232337111381915d0
xm1 = 1.d0 - xx
dlnx = dlog(xx)
dlxm1 = dlog(xm1)
xp1 = 1.0d0 + xx
dlxp1 = dlog(xp1)
dli2a = dilog(xm1)
dli3a = trilog(xm1)
dimz = wgplg(1,1,-xx)
trimz = wgplg(2,1,-xx)
zcom = xm1/xp1
tripco= wgplg(2,1,zcom)
trimco= wgplg(2,1,-zcom)
s1mz = wgplg(1,2,xm1)
smz = wgplg(1,2,-xx)
Z = xx
partCRI = CF * ( CF - CA/2 ) * (
. dlnx**3 * ( - 4.D0/3.D0 - 4.D0/3.D0*Z**2 + 8.D0/3.D0*Z )
. + dlnx**2 * ( - 6.D0 - 6.D0*Z**2 + 12.D0*Z )
. + dlnx*dli2a * ( - 8.D0 - 8.D0*Z**2 + 16.D0*Z )
. + dlnx * ( - 14.D0 + 12.D0*Z )
. + dli2a * ( - 12.D0 - 12.D0*Z**2 + 24.D0*Z )
. + dli3a * ( 8.D0 + 8.D0*Z**2 - 16.D0*Z )
. + S1MZ * ( - 8.D0 - 8.D0*Z**2 + 16.D0*Z )
. - 15.D0 - 13.D0*Z**2 + 28.D0*Z )
CQQCRI = partCRI
return
end
double precision function CQQCRF(xx,lfh,lfr)
implicit none
double precision partCRF
double precision xx,lfh,lfr,z
double precision xm1,xp1,dlnx,dlxm1,dlxp1
double precision dli2a,dli3a,zcom,s1mz,dimz,trimz,smz
double precision trimco,tripco
double precision part1,part2,scale
double precision z2,z3,z4
double precision ca,cf,nf
double precision dilog,trilog,wgplg
double precision cax,cv,cd,sw
external dilog,trilog,wgplg
sw = 0.2228972225239183d0
cax = 1.d0
cv = -1.d0 + 4.d0/3.d0*sw
cd = 1.d0 + cv**2
ca=3.d0
cf=4.d0/ca
nf=5d0
z2 = 1.6449340668482264365d0
z3 = 1.2020569031595942854d0
z4 = 1.0823232337111381915d0
xm1 = 1.d0 - xx
dlnx = dlog(xx)
dlxm1 = dlog(xm1)
xp1 = 1.0d0 + xx
dlxp1 = dlog(xp1)
dli2a = dilog(xm1)
dli3a = trilog(xm1)
dimz = wgplg(1,1,-xx)
trimz = wgplg(2,1,-xx)
zcom = xm1/xp1
tripco= wgplg(2,1,zcom)
trimco= wgplg(2,1,-zcom)
s1mz = wgplg(1,2,xm1)
smz = wgplg(1,2,-xx)
Z = xx
partCRF =
. dlnx**3 * ( 2.D0 - 2.D0*Z - 16.D0/3.D0/xp1 )
. + dlnx**2*dlxm1 * ( - 8.D0 + 8.D0*Z + 16.D0/xp1 )
. + dlnx**2*dlxp1 * ( - 24.D0 + 24.D0*Z + 56.D0/xp1 )
. + dlnx**2 * ( - 4.D0 - 12.D0*Z )
. + dlnx*dlxm1*dlxp1 * ( 32.D0 - 32.D0*Z - 64.D0/xp1 )
. + dlnx*dlxm1 * ( 16.D0 + 16.D0*Z )
. + dlnx*dlxp1**2 * ( - 16.D0/xp1 )
. + dlnx*dlxp1 * ( 8.D0 + 8.D0*Z )
. + dlnx*dli2a * ( - 24.D0 + 24.D0*Z + 48.D0/xp1 )
. + dlnx*DIMZ * ( - 32.D0 + 32.D0*Z + 64.D0/xp1 )
. + dlnx*z2 * ( - 8.D0 + 8.D0*Z + 24.D0/xp1 )
. + dlnx * ( - 18.D0 + 14.D0*Z )
. + dlxm1*DIMZ * ( 32.D0 - 32.D0*Z - 64.D0/xp1 )
. + dlxm1*z2 * ( 16.D0 - 16.D0*Z - 32.D0/xp1 )
. + dlxm1 * ( 32.D0 - 32.D0*Z )
. + dlxp1*DIMZ * ( - 32.D0/xp1 )
. + dlxp1*z2 * ( - 16.D0/xp1 )
partCRF = partCRF
. + dli2a * ( 24.D0 + 8.D0*Z )
. + DIMZ * ( 8.D0 + 8.D0*Z )
. + z2 * ( 4.D0 + 4.D0*Z )
. + dli3a * ( 32.D0 - 32.D0*Z - 64.D0/xp1 )
. + TRIMZ * ( 16.D0 - 16.D0*Z - 16.D0/xp1 )
. + TRIMCO * ( 32.D0 - 32.D0*Z - 64.D0/xp1 )
. + TRIPCO * ( - 32.D0 + 32.D0*Z + 64.D0/xp1 )
. + S1MZ * ( - 32.D0 + 32.D0*Z + 64.D0/xp1 )
. + SMZ * ( - 32.D0/xp1 )
. + z3 * ( 12.D0 - 12.D0*Z - 8.D0/xp1 )
. - 34.D0 + 34.D0*Z
partCRF = partCRF + lfh * (
. dlnx**2 * ( - 4.D0 + 4.D0*Z + 8.D0/xp1 )
. + dlnx*dlxp1 * ( 16.D0 - 16.D0*Z - 32.D0/xp1 )
. + dlnx * ( 8.D0 + 8.D0*Z )
. + DIMZ * ( 16.D0 - 16.D0*Z - 32.D0/xp1 )
. + z2 * ( 8.D0 - 8.D0*Z - 16.D0/xp1 )
. + 16.D0 - 16.D0*Z )
partCRF = cf * ( cf - ca/2 ) * partCRF
CQQCRF = partCRF
return
end
Index: branches/bbz/src/Massless_lim.C
===================================================================
--- branches/bbz/src/Massless_lim.C (revision 37)
+++ branches/bbz/src/Massless_lim.C (revision 38)
@@ -1,81 +1,80 @@
#include "Massless_lim.H"
#include "Massless_Kernels.H"
#include <iostream>
#include <vector>
#include <cmath>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_vegas.h>
using namespace LHAPDF;
using namespace std;
-Massless_Limit::Massless_Limit(const char *pdfname, const Scales scale) :
- i_nev(scale.i_evs),
+Massless_Limit::Massless_Limit(const char *pdfname, const Scales scale) :
m2r(pow(scale.mu_R,2)),m2f(pow(scale.mu_F,2)),
- m2h(pow(scale.m_H,2)),S(pow(scale.ss,2)), m2b(pow(scale.mbmb,2))
+ m2h(pow(scale.m_H,2)),m2b(pow(scale.mbmb,2)),S(pow(scale.ss,2)),i_nev(scale.i_evs)
{
m_mode = scale.i_fonll;
m_amin = scale.m_amin;
p_lumi=new LUMI(pdfname,scale.i_member,scale.mu_F,scale.mbmb,false);
m_cu = m_cd = 1.;
m_proc = 1; // bbh
// if(scale.m_proc_name == "bbh" || scale.m_proc_name != "bbz"){
// if(p_lumi->Init_Lumi()){
// m_inte=m_err=0;
// lfr = log(m2f/m2r);
// lfh = log(m2f/m2h);
// tauh = m2h/S;
// double mbr, cv, qmb, vev;
// qmb = 4.18;
// mbr = p_lumi->mb(qmb,sqrt(m2r),4);
// cv = 3.8937966e8;
// vev = 246.221;
// cout << "bbH ==== \n Parameters: \n MB("<< sqrt(m2r)<< " ): " << mbr << "\n"
// << "tauh("<< sqrt(m2h) << "): " << tauh << "\n"
// << endl;
// s0 = M_PI/6. * cv/m2h * pow(mbr/vev,2);
// m_cu = m_cd = 1.;
// }
// }
// else if(scale.m_proc_name == "bbz"){
if(p_lumi->Init_Lumi()){
m_inte=m_err=0;
lfr = log(m2f/m2r);
lfh = log(m2f/m2h);
tauh = m2h/S;
double mbr, cv, qmb, vev;
qmb = 4.18;
-
+ vev = 246.221;
mbr = p_lumi->mb(qmb,sqrt(m2r),4);
cv = 3.8937966e8;
vev = 246.221;
cout << "bbZ ==== \n Parameters: \n MB("<< sqrt(m2r)<< " ): " << mbr << "\n"
- << "tauh("<< sqrt(m2h) << "): " << tauh << "\n"
+ << "tauh("<< sqrt(m2h) << "): " << tauh/vev << "\n"
<< endl;
double GF = 1.1663787e-5;
-
+
double mw = 80.385, mz = sqrt(m2h);
double sw = 1. - sqr(mw/mz);
m_cu = 1 + sqr(-8./3.*sw + 1.);
m_cd = 1 + sqr(4./3.*sw - 1.);
m_proc = 1;
s0 = cv*M_PI*sqrt(2.)*GF/12.*m_cd;
// cout <<scientific << setprecision(15)<<sw << endl;
}
// }
}
Massless_Limit::~Massless_Limit()
{
if(p_lumi) delete p_lumi;
}
Index: branches/bbz/src/Makefile
===================================================================
--- branches/bbz/src/Makefile (revision 37)
+++ branches/bbz/src/Makefile (revision 38)
@@ -1,39 +1,43 @@
CXXFLAGS = $(shell lhapdf-config --cflags)
CXXFLAGS += $(shell gsl-config --cflags)
LDFLAGS = $(shell lhapdf-config --libs)
LDFLAGS += $(shell gsl-config --libs)
-LDFLAGS += -L. -lbbh -lgfortran
-PYFLAGS = -fPIC -shared
-PYFLAGS += ${CXXFLAGS}
+LDFLAGS += -L. -lbbh -L$(HOME)/work/local/lib -lgfortran #-L/usr/lib64 -lgfortran
+
+PYPREFIX=$(shell python3-config --prefix)
+
+PYFLAGS += $(shell python3-config --cflags)
+PYFLAGS += -L${PYPREFIX}/lib -lpython3.6m -ldl
PYFLAGS += $(shell python3-config --cflags --ldflags)
+
PYFLAGS += ${LDFLAGS}
INSTALLPATH=plots_HXSWG
CXX = g++
libbbh.a: bits.o # bbh.o bbhnnloexact.o
ar r $@ $^
bits.o: bbznnlo/bits.f # bbh.o bbhnnloexact.o: bbhnnlo/bbh.f # bbhnnlo/bbhnnloexact.f
gfortran -c -g -fPIC $^
%.o: %.C
$(CXX) -fPIC -c -g $< ${CXXFLAGS}
test: 4F.o Calculate.o 4F0_LO.o 4F0_NLO.o 5F.o Massless_lim.o PDF.o
$(CXX) -g $^ ${LDFLAGS} -o $@
bbHFONLL_wrap.cxx: bbHFONLL.i
swig -c++ -python $<
python: 4F.C Calculate.C 4F0_LO.C 4F0_NLO.C 5F.C Massless_lim.C PDF.C bbHFONLL_wrap.cxx
- $(CXX) -g -fPIC -o _bbHFONLL.so $^ ${PYFLAGS}
+ $(CXX) -g -fPIC -o _bbHFONLL.so $^ ${PYFLAGS} -Wl,-rpath,${PYPREFIX}/lib
install: _bbHFONLL.so bbHFONLL.py
cp $^ ${INSTALLPATH}/.
clean:
rm -rf *.o test *~ ../lib/libbbh.a *_wrap.cxx ../lib/_*.so ../lib/bbHFONLL.py *.dSYM libbbh.a _*.so bbHFONLL.py
Index: branches/bbz/src/Calculate.C
===================================================================
--- branches/bbz/src/Calculate.C (revision 37)
+++ branches/bbz/src/Calculate.C (revision 38)
@@ -1,390 +1,427 @@
#include "Massless_lim.H"
#include "Massless_Kernels.H"
#include "PDF.H"
#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_vegas.h>
using namespace std;
Massless_Limit* mlim_instance;
-
double wrap_LO(double x[], size_t dim, void* p){
return mlim_instance->dsigmaLO(x, dim, p);
}
double wrap_NLO(double x[], size_t dim, void* p){
if(dim == 1) return mlim_instance->dsigmaNLO(x, dim, p, 0);
if(dim == 2)
return
mlim_instance->dsigmaNLO(x, dim, p, 0)
- mlim_instance->dsigmaNLO(x, dim, p, 1);
else return 0.;
}
double wrap_NNLO(double x[], size_t dim, void* p){
if(dim == 1) return mlim_instance->dsigmaNNLO(x, dim, p, 0);
if(dim == 2)
return
mlim_instance->dsigmaNNLO(x, dim, p, 0)
- mlim_instance->dsigmaNNLO(x, dim, p, 1);
else return 0.;
}
double xsec5F_LO(const Scales scal)
{
double inte(0.), w2(0.);
size_t l=1;
gsl_monte_function F;
double xmin[l], xmax[l];
for(size_t i(0); i<l; ++i){
xmin[i] = 0.; xmax[i] =1.;
}
F.f = &(wrap_LO);
F.dim = l;
gsl_monte_vegas_state* state = gsl_monte_vegas_alloc(l);
gsl_monte_vegas_init(state);
gsl_rng * rng = gsl_rng_alloc(gsl_rng_taus);
gsl_monte_vegas_integrate(&F,xmin,xmax,l,scal.i_evs,rng,state,&inte,&w2);
gsl_monte_vegas_integrate(&F,xmin,xmax,l,scal.i_evs,rng,state,&inte,&w2);
gsl_rng_free(rng);
gsl_monte_vegas_free(state);
// cout <<" Delta-LO ... : " << inte << " ± " << w2 << std::endl;
return inte;
}
double xsec5F_NLO(const Scales scal)
{
double inte_1(0.), w2_1(0.);
double inte_2(0.), w2_2(0.);
size_t l1(1), l2(2);
gsl_monte_function Fd, Fp;
double xmin1[l1], xmax1[l1];
double xmin2[l2], xmax2[l2];
xmin1[0]=0.; xmax1[0]=1.;
for(size_t i(0); i<l2; ++i){
xmin2[i] = 0.; xmax2[i] =1.;
}
Fd.dim = l1; // delta-terms
Fp.dim = l2; // other terms
Fd.f = &(wrap_NLO);
Fp.f = &(wrap_NLO);
gsl_monte_vegas_state* state1 = gsl_monte_vegas_alloc(l1);
gsl_monte_vegas_init(state1);
gsl_rng * rng1 = gsl_rng_alloc(gsl_rng_taus);
gsl_monte_vegas_integrate(&Fd,xmin1,xmax1,l1,
scal.i_evs,rng1,state1,
&inte_1,&w2_1);
gsl_monte_vegas_integrate(&Fd,xmin1,xmax1,l1,
scal.i_evs,rng1,state1,
&inte_1,&w2_1);
gsl_rng_free(rng1);
gsl_monte_vegas_free(state1);
// cout <<" Delta-NLO ... : " << inte_1 << " ± " << w2_1 << std::endl;
gsl_monte_vegas_state* state2 = gsl_monte_vegas_alloc(l2);
gsl_monte_vegas_init(state2);
gsl_rng * rng2 = gsl_rng_alloc(gsl_rng_taus);
gsl_monte_vegas_integrate(&Fp,xmin2,xmax2,l2,
scal.i_evs,rng2,state2,
&inte_2,&w2_2);
gsl_monte_vegas_integrate(&Fp,xmin2,xmax2,l2,
scal.i_evs,rng2,state2,
&inte_2,&w2_2);
gsl_rng_free(rng2);
gsl_monte_vegas_free(state2);
// cout <<" Rest ........ : " << inte_2 << " ± " << w2_2 << std::endl;
// cout << " ----------------- " << endl;
// cout << inte_1+inte_2 << " ± " << w2_1 + w2_2 << endl;
return inte_1+inte_2;
}
double xsec5F_NNLO(const Scales scal)
{
double inte_1(0.), w2_1(0.);
double inte_2(0.), w2_2(0.);
size_t l1(1), l2(2);
gsl_monte_function Fd, Fp;
double xmin1[l1], xmax1[l1];
double xmin2[l2], xmax2[l2];
xmin1[0]=0.; xmax1[0]=1.;
for(size_t i(0); i<l2; ++i){
xmin2[i] = 0.; xmax2[i] =1.;
}
Fd.dim = l1; // delta-terms
Fp.dim = l2; // other terms
Fd.f = &(wrap_NNLO);
Fp.f = &(wrap_NNLO);
gsl_monte_vegas_state* state1 = gsl_monte_vegas_alloc(l1);
gsl_monte_vegas_init(state1);
gsl_rng * rng1 = gsl_rng_alloc(gsl_rng_taus);
gsl_monte_vegas_integrate(&Fd,xmin1,xmax1,l1,
scal.i_evs,rng1,state1,
&inte_1,&w2_1);
gsl_monte_vegas_integrate(&Fd,xmin1,xmax1,l1,
scal.i_evs,rng1,state1,
&inte_1,&w2_1);
gsl_rng_free(rng1);
gsl_monte_vegas_free(state1);
// cout <<" Delta-NNLO ... : " << inte_1 << " ± " << w2_1 << std::endl;
gsl_monte_vegas_state* state2 = gsl_monte_vegas_alloc(l2);
gsl_monte_vegas_init(state2);
gsl_rng * rng2 = gsl_rng_alloc(gsl_rng_taus);
gsl_monte_vegas_integrate(&Fp,xmin2,xmax2,l2,
scal.i_evs,rng2,state2,
&inte_2,&w2_2);
gsl_monte_vegas_integrate(&Fp,xmin2,xmax2,l2,
scal.i_evs,rng2,state2,
&inte_2,&w2_2);
gsl_rng_free(rng2);
gsl_monte_vegas_free(state2);
// cout <<" Rest ........ : " << inte_2 << " ± " << w2_2 << std::endl;
// cout << " ----------------- " << endl;
return inte_1 + inte_2;
}
double xsec5F(const char *argv, const Scales scal)
{
mlim_instance = new Massless_Limit(argv, scal);
double xs(0.);
if ( scal.i_or == 0 ){
// cout << " ----------- Computing LO 5F XS --------- " << endl;
xs = xsec5F_LO(scal);
}
if ( scal.i_or == 1 ){
// cout << " ----------- Computing NLO 5F XS --------- " << endl;
xs = xsec5F_LO(scal) + xsec5F_NLO(scal);
// the nlo-res is slowly convergent
// cout << "NLO xs = " << xs << endl;
}
if ( scal.i_or == 2 ){
// cout << " ----------- Computing NNLO 5F XS --------- " << endl;
xs = xsec5F_LO(scal) + xsec5F_NLO(scal) + xsec5F_NNLO(scal);
// cout << "NNLO xs = " << xs << endl;
}
if (scal.i_or > 2 || scal.i_or < 0 ) xs = 0;
delete mlim_instance;
return xs;
}
double wrap_0LO(double x[], size_t dim, void* p){
if(dim == 3) return mlim_instance->dsigma0LO(x, dim, p);
else return 0.;
}
double wrap_0NLO(double x[], size_t dim, void* p){
if(dim == 4) return mlim_instance->dsigma0NLO(x, dim, p);
else return 0.;
}
double xsec4F0_LO(const Scales scal)
{
double inte(0.), w2(0.);
size_t l=3;
gsl_monte_function F;
double xmin[l], xmax[l];
for(size_t i(0); i<l; ++i){
xmin[i] = 0.; xmax[i] =1.;
}
F.f = &(wrap_0LO);
F.dim = l;
gsl_monte_vegas_state* state = gsl_monte_vegas_alloc(l);
gsl_monte_vegas_init(state);
gsl_rng * rng = gsl_rng_alloc(gsl_rng_taus);
gsl_monte_vegas_integrate(&F,xmin,xmax,l,scal.i_evs,rng,state,&inte,&w2);
gsl_monte_vegas_integrate(&F,xmin,xmax,l,scal.i_evs,rng,state,&inte,&w2);
gsl_rng_free(rng);
gsl_monte_vegas_free(state);
// cout << " XS 4F O(as^2) .... " << inte << " ± " << w2 << std::endl;
//cout << " ----------------- " << endl;
- if (scal.i_or < 0){
- // std::cout << " WARNING M-Limit Only Calculation!!!! " << std::endl;
- size_t l1(2);
- double inte_1(0.), w2_1(0.);
- gsl_monte_function F1;
- double xmin1[l], xmax1[l];
- for(size_t i(0); i<l1; ++i){
- xmin1[i] = 0.; xmax1[i] =1.;
- }
-
- F1.f = &(wrap_0LO);
- F1.dim = l1;
- gsl_monte_vegas_state* state1 = gsl_monte_vegas_alloc(l1);
- gsl_monte_vegas_init(state1);
- gsl_rng * rng1 = gsl_rng_alloc(gsl_rng_taus);
- gsl_monte_vegas_integrate(&F1,xmin1,xmax1,l1,scal.i_evs,rng1,state1,&inte_1,&w2_1);
- gsl_monte_vegas_integrate(&F1,xmin1,xmax1,l1,scal.i_evs,rng1,state1,&inte_1,&w2_1);
-
- gsl_rng_free(rng1);
- gsl_monte_vegas_free(state1);
- // cout << " gg + qq : " << inte_1 << endl;
- inte+=inte_1;
+ // if (scal.i_or < 0){
+ // // std::cout << " WARNING M-Limit Only Calculation!!!! " << std::endl;
+ // size_t l1(2);
+ // double inte_1(0.), w2_1(0.);
+ // gsl_monte_function F1;
+ // double xmin1[l], xmax1[l];
+ // for(size_t i(0); i<l1; ++i){
+ // xmin1[i] = 0.; xmax1[i] =1.;
+ // }
+
+ // F1.f = &(wrap_0LO);
+ // F1.dim = l1;
+ // gsl_monte_vegas_state* state1 = gsl_monte_vegas_alloc(l1);
+ // gsl_monte_vegas_init(state1);
+ // gsl_rng * rng1 = gsl_rng_alloc(gsl_rng_taus);
+ // gsl_monte_vegas_integrate(&F1,xmin1,xmax1,l1,scal.i_evs,rng1,state1,&inte_1,&w2_1);
+ // gsl_monte_vegas_integrate(&F1,xmin1,xmax1,l1,scal.i_evs,rng1,state1,&inte_1,&w2_1);
+
+ // gsl_rng_free(rng1);
+ // gsl_monte_vegas_free(state1);
+ // // cout << " gg + qq : " << inte_1 << endl;
+ // inte+=inte_1;
- }
+ // }
return inte;
}
double xsec4F0_NLO(const Scales scal)
{
double inte(0.), w2(0.);
size_t l(4);
gsl_monte_function F;
double xmin[l], xmax[l];
- int Evs= scal.i_evs;//*100;
+
+ int Evs= scal.i_evs;
for(size_t i(0); i<l; ++i){
xmin[i] = 0.; xmax[i] =1.;
}
F.f = &(wrap_0NLO);
F.dim=l;
gsl_monte_vegas_state* state = gsl_monte_vegas_alloc(l);
gsl_monte_vegas_init(state);
gsl_rng * rng = gsl_rng_alloc(gsl_rng_taus);
- gsl_monte_vegas_integrate(&F,xmin,xmax,l,
- Evs,rng,state,
- &inte,&w2);
- gsl_monte_vegas_integrate(&F,xmin,xmax,l,
- Evs,rng,state,
- &inte,&w2);
+ gsl_monte_vegas_params vegas_par;
+ gsl_monte_vegas_params_get(state,&vegas_par);
+ vegas_par.iterations = 20;
+ gsl_monte_vegas_params_set(state,&vegas_par);
+ gsl_monte_vegas_integrate(&F,xmin,xmax,l,Evs,rng,state,&inte,&w2);
+
+ size_t rounds(0);
+ cout << "iteration : " << rounds << "\n"
+ << " % err : " << std::fabs(w2/inte)*100
+ << " , chi2 : " << gsl_monte_vegas_chisq(state) << endl;
+
+ while (fabs (gsl_monte_vegas_chisq(state) - 1.0) > 0.5){
+ vegas_par.iterations += 10;
+ if(rounds <3) vegas_par.stage = rounds+1;
+ gsl_monte_vegas_params_set(state,&vegas_par);
+ //apply this mod to the xs4F0 and include muF in that!
+ gsl_monte_vegas_integrate(&F,xmin,xmax,l,
+ size_t(Evs*(1+rounds)),rng,state,
+ &inte,&w2);
+ rounds++;
+ cout << "iteration : " << rounds << "\n"
+ << " % err : " << std::fabs(w2/inte)*100
+ << " , chi2 : " << gsl_monte_vegas_chisq(state) << endl;
+ }
+ cout << " success ! vegas converged after " << rounds << " iterations !" << endl;
+ cout << " % err : " << std::fabs(w2/inte)*100
+ << " , chi2 : " << gsl_monte_vegas_chisq(state) << endl;
gsl_rng_free(rng);
gsl_monte_vegas_free(state);
// cout << " XS 4F O(as^3) .... " << inte << " ± " <<w2 << std::endl;
// cout << " ----------------- " << endl;
return inte;
}
double wrap_muFTerms(double x[], size_t dim, void* p)
{
if(dim == 4) return mlim_instance->muFct(x, dim, p);
else return 0.;
}
double muFTerm(const Scales scal)
{
double inte(0.), w2(0.);
size_t l(4);
gsl_monte_function F;
double xmin[l], xmax[l];
int Evs= scal.i_evs;
for(size_t i(0); i<l; ++i){
xmin[i] = 0.; xmax[i] =1.;
}
F.f = &(wrap_muFTerms);
F.dim=l;
gsl_monte_vegas_state* state = gsl_monte_vegas_alloc(l);
gsl_monte_vegas_init(state);
+ gsl_monte_vegas_params vegas_par;
+ gsl_monte_vegas_params_get(state,&vegas_par);
+ vegas_par.iterations = 20;
+ gsl_monte_vegas_params_set(state,&vegas_par);
+
gsl_rng * rng = gsl_rng_alloc(gsl_rng_taus);
gsl_monte_vegas_integrate(&F,xmin,xmax,l,
- Evs*10,rng,state,
- &inte,&w2);
- gsl_monte_vegas_integrate(&F,xmin,xmax,l,
- Evs*10,rng,state,
+ Evs,rng,state,
&inte,&w2);
+ size_t rounds(0);
+ while (fabs (gsl_monte_vegas_chisq(state) - 1.0) > 0.5){
+ vegas_par.iterations += 10;
+ if(rounds <3) vegas_par.stage = rounds+1;
+ gsl_monte_vegas_params_set(state,&vegas_par);
+ //apply this mod to the xs4F0 and include muF in that!
+ gsl_monte_vegas_integrate(&F,xmin,xmax,l,
+ size_t(Evs*(1+rounds)),rng,state,
+ &inte,&w2);
+ rounds++;
+ cout << "iteration : " << rounds << "\n"
+ << " % err : " << std::fabs(w2/inte)*100
+ << " , chi2 : " << gsl_monte_vegas_chisq(state) << endl;
+ }
+ cout << " success ! vegas converged after " << rounds << " iterations !" << endl;
+ cout << " % err : " << std::fabs(w2/inte)*100
+ << " , chi2 : " << gsl_monte_vegas_chisq(state) << endl;
gsl_rng_free(rng);
gsl_monte_vegas_free(state);
- // cout << " muF-Terms O(as^3) .... " << inte << " ± " <<w2 << std::endl;
- // cout << " ----------------- " << endl;
-
return inte;
}
double xsec4F0(const char *argv, const Scales scal)
{
double xs(0.), xs40(0.);
mlim_instance = new Massless_Limit(argv, scal);
if( scal.i_fonll == 1) {
// cout << " ----------- Computing 4F0 XS O(as^2) --------- " << endl;
xs40 = xsec4F0_LO(scal);
xs = xs40;
// cout << "Massless-lim O(as^2) xs = .... " << xs << endl;
}
if( scal.i_fonll == 2) {
// cout << " ----------- Computing 4F0 XS O(as^3) --------- " << endl;
double m2r = pow(scal.mu_R,2);
double m2f = pow(scal.mu_F,2);
double m2h = pow(scal.m_H,2);
double b0 = (33.-10.)/12./M_PI;
cout << sqrt(m2r) << " , " << sqrt(m2h) << " , " << sqrt(m2f) << endl;
double asr = mlim_instance->masr(m2r);
double asf = mlim_instance->masr(m2f);
double A = xsec4F0_LO(scal);
double B = xsec4F0_NLO(scal);
double C = 0.;
if(m2f!=m2h)
C = muFTerm(scal);
xs40 = A +(B + 2.*asr*b0*A*log(m2r/m2h)
-asf*C*log(m2f/m2h)/2./M_PI);
- cout << " sLO : " << A << " , sNLO : " << B << " , smuR : " << 2.*asr*b0*A*log(m2r/m2h)
+ cout << " sLO : " << A << " , sNLO : " << B << " , smuR : " << 2.*asr*b0*A*log(m2r/m2h)
<< " , smuF : " << asf*C*log(m2f/m2h)/2./M_PI << endl;
xs = xs40;
// cout << " -------------------------------------------\n" ;
// cout << " ------ Massless lim O(as^3) xs = ... " << xs<< " ------"<<endl;
// cout << " -------------------------------------------\n" ;
}
if (scal.i_fonll > 2 || scal.i_fonll < 0) xs40=0;
xs = xs40;
delete mlim_instance;
return xs;
}
double xsecDiff(const char *argv, const Scales scal)
{
double xs(0.), xs5(0.), xs40(0.);
xs5 = xsec5F(argv,scal); xs40 = xsec4F0(argv,scal);
// cout << " ############################################### " << endl
// << " XS 5F : \t " << xs5 << endl
// << " XS 4F : \t " << xs40 << endl
// << " ############################################### " << endl;
xs = xs5 - xs40;
return xs;
}
int main()
{
//order 5f, order fonll, mb_pole, muR, muF, mV, eCM, evs, member, amin, process
- Scales sc(-1,2,4.92, 33., 91.188, 91.188, 13000., 100000,0,0.);
+ Scales sc(-1,2,4.92, 33., 4., 91.188, 13000., 100000,0,1.e-12);
// cout << xsec5F("NNPDF31_nnlo_as_0118",sc) << endl;
xsec4F0("NNPDF31_nnlo_as_0118",sc);
return 0;
}
Index: branches/bbz/src/Massless_lim.H
===================================================================
--- branches/bbz/src/Massless_lim.H (revision 37)
+++ branches/bbz/src/Massless_lim.H (revision 38)
@@ -1,121 +1,122 @@
#ifndef MASSLESS_LIM_H
#define MASSLESS_LIM_H
#include <vector>
#include <cmath>
#include <string>
#include "PDF.H"
class Scales {
public:
- double mu_R, mu_F, m_H, ss;
- double mbmb, m_vev, m_amin;
- int i_evs, i_or, i_fonll, i_member;
+ int i_or, i_fonll;
+ double mbmb, mu_R, mu_F, m_H, ss;
+ double m_vev, m_amin;
+ int i_evs, i_member;
public:
Scales(const int order, const int ofonll,
const double mf, const double mur, const double muF,
const double mH, const double S, const int evs,
const int PDFmem, const double amin)://, char *proc_name):
i_or(order),i_fonll(ofonll),mbmb(mf),
mu_R(mur),mu_F(muF),m_H(mH),ss(S),i_evs(evs),
i_member(PDFmem)
// m_proc_name(proc_name)
{
m_amin = amin;
}
};
class Massless_Limit {
private:
LUMI * p_lumi;
double lfr, lfh, tauh,m_inte,m_err;
double m2r,m2f,m2h,m2b,S;
double s0,m_amin,m_cu,m_cd;
int i_nev, m_mode,m_proc;
//m_mode = 0 ---> 5F or Massless_lim
//m_mode = 1 ---> Difference Term
public:
Massless_Limit(const char *pdfname, const Scales scale);
//const int mode);
~Massless_Limit();
double dsigmaLO(double x[], size_t dim, void* p);
double dsigmaNLO(double x[], size_t dim,
void* p, const size_t plus);
double dsigmaNNLO(double x[], size_t dim,
void* p, const size_t plus);
double dsigma0LO(double x[], size_t dim, void* p);
double dsigma0NLO(double x[], size_t dim, void* p);
double dsigmagg_lo(double y,double z);
double dsigmaqqb_lo(double y, double z);
double Agg_2_2(double y, double z);
double Agg_2_1(double y, double z);
double Agg_2_0(double y);
double dsigmagg_nlo(double y,double z, double t);
double dsigmaqqb_nlo(double y, double z);
double dsigmagq_nlo(double y, double z);
double dsigmaqg_nlo(double y, double z);
double dsigmagg_st(double y,double z, double kp);
double dsigmaqq_st(double y,double z);
double dsigmagq_st(double y,double z, double kp);
double muFct(double z[], size_t dim, void* p);
double slogg(double y,double z, double r);
double sloqq(double y, double z);
double bb(double y, double z, double t);
double bg(double y, double z);
double bb_LO(double y, double z);
double bb_NLO(double y, double z, double t);
double bg_LO(double y, double z);
double bg_NLO(double y, double z);
double bbs(double y, double z);
double bgs(double y, double z);
double bqs(double y, double z);
double fun(double y, double z);
double L_qq(const double x, const double y);
double L_gq(const double x, const double y);
inline double masr(const double Q2) { return 2.*M_PI*p_lumi->ass(Q2); }
};
double wrap_0LO(double x[], size_t dim, void* p);
double wrap_0NLO(double x[], size_t dim, void* p);
double wrap_muFTerms(double x[], size_t dim, void* p);
double wrap_LO(double x[], size_t dim, void* p);
double wrap_NLO(double x[], size_t dim, void* p);
double wrap_NNLO(double x[], size_t dim, void* p);
double xsec5F_LO(const Scales scal);
double xsec5F_NLO(const Scales scal);
double xsec5F_NNLO(const Scales scal);
double xsec4F0_LO(const Scales scal);
double xsec4F0_NLO(const Scales scal);
double muFTerm(const Scales scal);
double xsec5F(const char *argv, const Scales scal);
double xsec4F0(const char *argv, const Scales scal);
double xsecDiff(const char *argv, const Scales scal);
#endif
Index: branches/bbz/src/4F0_NLO.C
===================================================================
--- branches/bbz/src/4F0_NLO.C (revision 37)
+++ branches/bbz/src/4F0_NLO.C (revision 38)
@@ -1,390 +1,392 @@
#include "Massless_lim.H"
#include "Massless_Kernels.H"
#include "PDF.H"
#include <cmath>
#include <iostream>
using namespace std;
double Massless_Limit::dsigma0NLO(double z[], size_t dim, void *p)
{
double res(0.);
// first round of variables
double tau,xi,ximax,jac(1.),x,y;
// lumis and MEs
double lgg,lqg,lgq,lqq,
sgg,sgq,sqg, sqq;
//second variables
double zz, tt;
double asr(2.*M_PI*Massless_Limit::p_lumi->ass(m2r));
//cout << asr << " , " <<sqrt(m2r) << endl;
double acut = m_amin; //0.7e-2;
double MAX = 1. - acut;
- if(z[0] < MAX && z[1] < MAX && z[2] < MAX && z[3] < MAX){
- // tau = tauh + (1. - tauh)*z[0];
- // jac *= (1. - tauh);
- tau = tauh*pow(MAX/tauh,z[0]);
- jac *= - log(tauh);
-
- ximax = -1./2.*log(tau);
- xi = -ximax + 2.*ximax*z[1];
- jac *= 2.*ximax;
-
- x = sqrt(tau)*exp(xi);
- y = sqrt(tau)*exp(-xi);
-
- lgg = Massless_Limit::p_lumi->Lumi(5,5,x,tauh/tau,1.,1.,10);
- lgq = Massless_Limit::L_gq(x, tauh/tau);
- lqg = Massless_Limit::L_gq(tauh/tau, x);
- lqq = Massless_Limit::L_qq(x,tauh/tau);
-
- // zz = y + (1. - y)*z[2];
- // jac *= (1. - y);
- zz = y*pow(1/y,z[2]);
- jac *= - log(y);
+ // tau = tauh + (1. - tauh)*z[0];
+ // jac *= (1. - tauh);
+ tau = tauh*pow(MAX/tauh,z[0]);
+ jac *= log(MAX/tauh)*tau;
- tt = y/zz + (MAX - y/zz)*z[3];
- jac *= (MAX - y/zz);
+ ximax = -1./2.*log(tau);
+ xi = -ximax + 2.*ximax*z[1];
+ jac *= 2.*ximax;
+
+ x = sqrt(tau)*exp(xi);
+ y = sqrt(tau)*exp(-xi);
+
+ lgg = Massless_Limit::p_lumi->Lumi(5,5,x,tauh/tau,1.,1.,10);
+ lgq = Massless_Limit::L_gq(x, tauh/tau);
+ lqg = Massless_Limit::L_gq(tauh/tau, x);
+ lqq = Massless_Limit::L_qq(x,tauh/tau);
+
+ // zz = y + (1. - y)*z[2];
+ // jac *= (1. - y);
+ zz = y*pow(MAX/y,z[2]);
+ jac *= log(MAX/y)*zz;
- sgg = dsigmagg_nlo(y,zz,tt);
- sgq = dsigmagq_nlo(y,zz)/(1. - y/zz);
+ // tt = y/zz + (MAX - y/zz)*z[3];
+ // jac *= (MAX - y/zz);
- sqg = dsigmaqg_nlo(y,zz)/(1. - y/zz);
+ tt = y/zz*pow(MAX*zz/y,z[3]);
+ jac *= log(MAX*zz/y)*tt;
+
+ sgg = dsigmagg_nlo(y,zz,tt);
- sqq = 0.;
+ sgq = dsigmagq_nlo(y,zz)/(1. - y/zz);
+
+ sqg = dsigmaqg_nlo(y,zz)/(1. - y/zz);
- res = pow(asr,3)*jac*s0*zz*(lgg*sgg + lqq*sqq
- + lqg*sqg + lgq*sgq);// /tau;
- }
+ sqq = 0.;
+
+ //cout << "zz, tt, MAX, "<< << endl;
+ res = pow(asr,3)*jac*s0*(lgg*sgg + lqq*sqq
+ + lqg*sqg + lgq*sgq)/tau;
+
return res;
}
double Massless_Limit::dsigmaqqb_nlo(double y, double z)
{
return 0.;
}
double Massless_Limit::dsigmagg_nlo(double y, double z, double t)
{
double res(0.);
res = 2.*bb(y,z,t)+ 4.*bg(y,z)/(1. - y/z);
return res;
}
double Massless_Limit::bb(double y, double z, double t){
return bb_LO(y,z)/(1.-y/z) + bb_NLO(y,z,t);
}
double Massless_Limit::bb_LO(double y, double z)
{
double res(0.);
double L(log(m2h/m2b)), agb, pqg;
pqg = L*Massless_Limit::p_lumi->Pqg(y/z)/(2.*M_PI);
agb = (L*L*Massless_Limit::p_lumi->agb(z,2,2)
+ L*Massless_Limit::p_lumi->agb(z,2,1)
+ Massless_Limit::p_lumi->agb(z,2,0))/pow(2.*M_PI,2);
res = 2.*y*pqg*agb/z;
return res;
}
double Massless_Limit::bb_NLO(double y, double z, double t)
{
double res(0.);
double sc(0.);
double L(log(m2h/m2b)), pqg1, pqg2, pqg20;
double k=y/z;
double sbb_d, sbb_p, sbb_sub, sbb_rest;
- sbb_d = delta1_(&sc)/M_PI/4.; sbb_p = dterms1_(&t,&sc)/M_PI/4.;
- sbb_sub = - dtsub1_(&sc,&k)/M_PI/4.; sbb_rest = sbbq1_(&t,&sc)/M_PI/4.;
+ sbb_d = delta1_(&sc); sbb_p = dterms1_(&t,&sc);
+ sbb_sub = - dtsub1_(&sc,&k); sbb_rest = sbbq1_(&t,&sc);
pqg1 = Massless_Limit::p_lumi->Pqg(z)/(2.*M_PI)*L;
pqg20 = Massless_Limit::p_lumi->Pqg(y/z)/(2.*M_PI)*L;
pqg2 = Massless_Limit::p_lumi->Pqg(y/z/t)/(2.*M_PI)*L;
// delta bit:
res = sbb_d*pqg20/(1. - y/z);
// plus bit:
- double x = y/z;
- double fac = 0;
+ // double x = y/z;
+ // double fac = 0;
// // I only want the coeffs of the plus as the integral is done anal!
// sbb_p *= (1.-t)/log(1.-t);
// // here goes the integral of pqg(y/zt)/t^2 times the plus terms.
// // prefactor of the splitting function, note that the division for 4pi is included in
// // sbb_p
// fac = sbb_p*L/2./M_PI;
// fac*= (2. - 8.*x + (6. + x*(-9. + 11.*x))*log(x) + 6.*(1. + (-1. + x)*x)*dilog(x) +
// log(1. - x)*(-9. + 15.*x + 5.*pow(x,-1) - 11.*pow(x,2)) +
// pow(M_PI,2)*(-1. + x - pow(x,2)) + 6.*pow(x,2))/12.;
// res += fac/(1.-y/z);
// /*
// alternative, purely numerical implementation of the previous formula:
// sbb_p = dterms1_(&t,&sc)/M_PI/4.;
- res += sbb_p*(pqg2/t/t-pqg20);
+ res += sbb_p*(pqg2/t/t - pqg20);
// */
// plus rest:
res += sbb_sub*pqg20/(1.-y/z);
// rest:
res += pqg2*sbb_rest/t/t;
res*=y*pqg1/z;
return res;
}
double Massless_Limit::bg(double y, double z)
{
return bg_LO(y,z) + bg_NLO(y,z);
}
double Massless_Limit::bg_LO(double y, double z)
{
double res(0.);
double L(log(m2h/m2b)), agb;
double sbg, sc(0.);
// double k=y/z;
agb = (Massless_Limit::p_lumi->agb(y/z,2,2)*L*L
+ Massless_Limit::p_lumi->agb(y/z,2,1)*L
+ Massless_Limit::p_lumi->agb(y/z,2,0))/pow(2.*M_PI,2);
- sbg = sbg1_(&z,&sc)/M_PI/4.;
+ sbg = sbg1_(&z,&sc);
res = y*sbg*agb/z/z;
return res;
}
double Massless_Limit::bg_NLO(double y, double z)
{
double res(0.);
double L(log(m2h/m2b)), pqg;
double sbg, sc(0.), nf(5);
// double k=y/z;
pqg = Massless_Limit::p_lumi->Pqg(y/z)*L/(2.*M_PI);
sbg = (deltabga_(&z,&sc,&sc) +
- nf*deltabgf_(&z,&sc,&sc))/pow(4.*M_PI,2);
+ nf*deltabgf_(&z,&sc,&sc));
res = y*pqg*sbg/z/z;
return res;
}
double Massless_Limit::dsigmagq_nlo(double y, double z)
{
double res(0.);
res = bbs(y,z) + 2.*bgs(y,z) + 2.*bqs(y,z);
return res;
}
double Massless_Limit::dsigmaqg_nlo(double y, double z)
{
return Massless_Limit::dsigmagq_nlo(y,z);
}
double Massless_Limit::bbs(double y, double z)
{
double res(0.);
double L(log(m2h/m2b)), asb, pqg;
pqg = Massless_Limit::p_lumi->Pqg(y/z)/(2.*M_PI)*L;
asb = (Massless_Limit::p_lumi->asb(z,2,2)*L*L
+ Massless_Limit::p_lumi->asb(z,2,1)*L
+ Massless_Limit::p_lumi->asb(z,2,0))/pow(2.*M_PI,2);
res = 2.*y*pqg*asb/z;
return res;
}
double Massless_Limit::bgs(double y, double z)
{
double res(0.);
double L(log(m2h/m2b)), asb;
double sbg, sc(0.);
- // double k=y/z;
asb = (Massless_Limit::p_lumi->asb(y/z,2,2)*L*L
+ Massless_Limit::p_lumi->asb(y/z,2,1)*L
+ Massless_Limit::p_lumi->asb(y/z,2,0))/pow(2.*M_PI,2);
- sbg = sbg1_(&z,&sc)/M_PI/4.;
+ sbg = sbg1_(&z,&sc);
res = y*sbg*asb/z/z;
return res;
}
double Massless_Limit::bqs(double y, double z)
{
double res(0.);
double L(log(m2h/m2b)), pqg;
double sbq, sc(0.), nf(5);
- // double k=y/z;
pqg = Massless_Limit::p_lumi->Pqg(y/z)/(2.*M_PI)*L;
sbq = (deltaqqba_(&z,&sc,&sc) +
- nf*deltaqqbf_(&z,&sc,&sc))/pow(4.*M_PI,2);
+ nf*deltaqqbf_(&z,&sc,&sc));
res = y*pqg*sbq/z/z;
return res;
}
double Massless_Limit::muFct(double z[], size_t dim, void *p)
{
double res(0.);
// first round of variables
double tau,xi,ximax,jac(1.),x,y;
// lumis and MEs
double lgg,lqg,lgq,lqq,
sgg,sgq,sqg, sqq;
//second variables
double zz, tt;
tau = tauh + (1. - tauh)*z[0];
jac *= (1. - tauh);
ximax = -1./2.*log(tau);
xi = -ximax + 2.*ximax*z[1];
jac *= 2.*ximax;
x = sqrt(tau)*exp(xi);
y = sqrt(tau)*exp(-xi);
lgg = Massless_Limit::p_lumi->Lumi(5,5,x,tauh/tau,1.,1.,10);
lgq = Massless_Limit::L_gq(x,tauh/tau);
lqg = Massless_Limit::L_gq(tauh/tau,x);
lqq = Massless_Limit::L_qq(x,tauh/tau);
zz = y + (1. - y)*z[2];
jac *= (1. - y);
tt = z[3];
sgg = dsigmagg_st(y,zz,tt);
sqg = dsigmagq_st(y,zz,tt);
sgq = sqg;
sqq = dsigmaqq_st(y,zz);
res = jac*s0*///
(lgg*sgg + lqq*sqq + lqg*sqg + lgq*sgq)/tau;
return res;
}
double Massless_Limit::slogg(double y, double z, double r)
{
double res(0.);
double t = y/z + (1. - y/z)*r;
double jac = (1. - y/z);
res = jac*Massless_Limit::dsigmagg_lo(y/z,y/z/t)/t/t;
return res;
}
double Massless_Limit::sloqq(double y, double z)
{
double res(0.);
res = (1. - y/z)*Massless_Limit::dsigmaqqb_lo(y/z,1.);
return res;
}
/// Sub term
double Massless_Limit::dsigmagg_st(double y,double z, double kp)
{
double ds(0.);
double t = (1.-y/z)*kp + y/z,
k = (1.-y)*kp + y;
double j = (1.-y/z);
double ppc(p_lumi->Pgg_pt(z)),ppc0(p_lumi->Pgg_pt(1.));
double pgg(p_lumi->Pgg(z)), d(p_lumi->Pgg_d());
ds = 2.*j*
( pgg*dsigmagg_lo(y/z,t) + (d + ppc0*log(1.-y))*dsigmagg_lo(y,z)/(1.-y/z) +
(ppc*dsigmagg_lo(y/z,t) - (1.-y)*(ppc0*dsigmagg_lo(y,k))/(1.-y/z))/(1.-z) );
return ds;
}
double Massless_Limit::dsigmaqq_st(double y, double z)
{
double ds(0.);
double h1(p_lumi->Pqq_pt(1.)), h(p_lumi->Pqq_pt(z)), sqq(dsigmaqqb_lo(y/z,1.)),
sqq1(dsigmaqqb_lo(y,1.)), d(p_lumi->Pqq_d());
ds = 2.*((d + h1*log(1.-y))*sqq1/(1.-y) +
(h*sqq - h1*sqq1)/(1.-z) );
return ds;
}
double Massless_Limit::dsigmagq_st(double y,double z, double kp)
{
double t = (1.-y/z)*kp + y/z;
double ds(0.),
dsgg(dsigmagg_lo(y/z,t)), dsqq(dsigmaqqb_lo(y/z,1.)),
nf = 5.;
ds = (1.-y/z)*(p_lumi->Pgq(z)*dsgg
+nf*p_lumi->Pqg(z)*dsqq);
return ds;
}
double Massless_Limit::fun(double y, double z){
double res = pow(y - 1.*z,-1)*pow(z,-4)*
(y*dilog(1. - 1.*y*pow(z,-1))*pow(z,2)*
(1.333*y*z - 2.67*pow(y,2) + 1.333*pow(z,2)) +
y*(11.85*z*pow(y,3) - 0.1852*pow(y,4) -
27.9*pow(y,2)*pow(z,2) + 22.6*y*pow(z,3) +
log(y*pow(z,-1))*
((0.667*y - 0.667*z)*z*log(y*pow(z,-1))*
pow(y,2) + 5.*z*pow(y,3) -
0.889*pow(y,4) -
13.44*pow(y,2)*pow(z,2) +
12.*y*pow(z,3) - 2.67*pow(z,4)) -
6.32*pow(z,4)) +
log(1. - 1.*y*pow(z,-1))*
(-10.89*z*pow(y,4) + 1.778*pow(y,5) +
23.8*pow(y,3)*pow(z,2) -
23.8*pow(y,2)*pow(z,3) +
y*z*log(y*pow(z,-1))*
(2.67*z*pow(y,2) - 2.67*pow(y,3) -
2.67*y*pow(z,2) + 2.67*pow(z,3)) +
10.89*y*pow(z,4) - 1.778*pow(z,5)) +
y*z*(-5.33*z*pow(y,2) + 2.67*pow(y,3) +
4.*y*pow(z,2) - 1.333*pow(z,3))*
pow(log(1. - 1.*y*pow(z,-1)),2));
return res;
}
// this was valid only for bbh
// fac =L*(4*pow(x,-1)*(2*(2 + 3*x*(-2 + 5*x))*log(1 - x) + 6*(1 - 3*x)*x*log(x) +
// 6*x*(1 + 2*(-1 + x)*x)*dilog(x) +
// x*(4*(-1 + x)*(-1 + 3*x) + (-1 - 2*(-1 + x)*x)*pow(M_PI,2)) -
// 44*atanh(1 - 2*x)*pow(x,3)))/9./M_PI/2./M_PI;
Index: branches/bbz/src/PDF.C
===================================================================
--- branches/bbz/src/PDF.C (revision 37)
+++ branches/bbz/src/PDF.C (revision 38)
@@ -1,389 +1,379 @@
#include "LHAPDF/LHAPDF.h"
#include "PDF.H"
#include "Massless_Kernels.H"
#include <cmath>
#define ZETA2 1.6449340668482264365
#define ZETA3 1.2020569031595942854
using namespace LHAPDF;
using namespace std;
LUMI::LUMI(const char* setname, const int setmem,
const double muF, const double mbext,
const bool ismlim) :
m_pdf(mkPDF(setname,setmem)),m_as(0.),Q(muF),m_flag(false), m_lim(ismlim)
{
if(m_pdf) m_flag=true;
m_as = m_pdf->alphasQ(Q)/(2.*M_PI);
mm2b = mbext*mbext;
m2b = m_pdf->alphaS().quarkMass(5)*m_pdf->alphaS().quarkMass(5);
}
LUMI::~LUMI()
{
if(m_pdf) delete m_pdf;
}
double LUMI::mb(const double mf, const double mu, const int nloop)
{
- double beta0, beta1,beta2,beta3, gamma0,gamma1,gamma2,gamma3;
- double A1,as,asmf,l2, nf(5.);
+ double as,asmf,nf(5.);
double mfrun(0.);
double mu0=mf;
- double Z3 = ZETA3, Z5 = 1.0369277551433699263;
-
- double cfunc1(0.), cfunc2(0.), cfunc3(0.), cfunc4(0.);
- double cfuncmu0(1.), cfuncmuf(1.);
-
-
as = m_pdf->alphasQ(mu)/M_PI;
asmf = m_pdf->alphasQ(mu0)/M_PI;
mfrun = runmass_(&mf,&asmf,&as,&nf,&nloop);
-
-
- return mfrun;
+ return mfrun;
}
double LUMI::Lumi(const int i, const int j,
const double x, const double y,
const double z1, const double z2,
const int order)
{
double Q2 = Q*Q;
// cout << " muF = " << Q << "\n";
if(m_lim == true) m2b = mm2b;
if(order==10){
return m_pdf->xfxQ2(i-5,x,Q2)*m_pdf->xfxQ2(j-5,y,Q2);
}
if(order==1){
if( (i == 10 || i == 0) && (j!=10 && j!=0)){
// double btildex = m_as*log(Q2/m2b)*Pqg(z1)*m_pdf->xfxQ2(0,x/z1,Q2);
double btildex = log(Q2/m2b)*Pqg(z1)*m_pdf->xfxQ2(0,x/z1,Q2);
return btildex*m_pdf->xfxQ2(j-5,y,Q2);
}
if( (i!=10 && i!=0) && (j == 10 || j == 0)){
// double btildey = m_as*log(Q2/m2b)*Pqg(z2)*m_pdf->xfxQ2(0,y/z2,Q2);
double btildey = log(Q2/m2b)*Pqg(z2)*m_pdf->xfxQ2(0,y/z2,Q2);
return m_pdf->xfxQ2(i-5,x,Q2)*btildey;
}
if( (i == 10 || i == 0) && (j == 10 || j == 0) ){
// double btildex = m_as*log(Q2/m2b)*Pqg(z1)*m_pdf->xfxQ2(0,x/z1,Q2);
// double btildey = m_as*log(Q2/m2b)*Pqg(z2)*m_pdf->xfxQ2(0,y/z2,Q2);
double btildex = log(Q2/m2b)*Pqg(z1)*m_pdf->xfxQ2(0,x/z1,Q2);
double btildey = log(Q2/m2b)*Pqg(z2)*m_pdf->xfxQ2(0,y/z2,Q2);
return btildex*btildey;
}
}
if(order==2){
double L = log(Q2/m2b);
if( (i == 10 || i == 0) && (j!=10 && j!=0 )){
// double bt2x = m_as*m_as*(L*L*phi(x,z1,Q2,2)+L*phi(x,z1,Q2,1));
double bt2x = L*L*phi(x,z1,Q2,2) + L*phi(x,z1,Q2,1);
if(sqrt(Q2) != 125.){
bt2x += phi(x,z1,Q2,0);
}
return bt2x*m_pdf->xfxQ2(j-5,y,Q2);
}
if( (i!=10 && i!=0) && (j == 10 || j == 0)){
// double bt2y = m_as*m_as*(L*L*phi(y,z2,Q2,2)+L*phi(y,z2,Q2,1));
double bt2y = L*L*phi(y,z2,Q2,2)+L*phi(y,z2,Q2,1);
if(sqrt(Q2) != 125.){
bt2y += phi(y,z2,Q2,0);
}
return m_pdf->xfxQ2(i-5,x,Q2)*bt2y;
}
if( (i == 10 || i == 0) && (j == 10 || j == 0) ){
// double bt1x = m_as*L*Pqg(z1)*m_pdf->xfxQ2(0,x/z1,Q2);
// double bt1y = m_as*L*Pqg(z2)*m_pdf->xfxQ2(0,y/z2,Q2);
// double bt2x = m_as*m_as*(L*L*phi(x,z1,Q2,2)+L*phi(x,z1,Q2,1));
// double bt2y = m_as*m_as*(L*L*phi(y,z2,Q2,2)+L*phi(y,z2,Q2,1));
double bt1x = L*Pqg(z1)*m_pdf->xfxQ2(0,x/z1,Q2);
double bt1y = L*Pqg(z2)*m_pdf->xfxQ2(0,y/z2,Q2);
double bt2x = L*L*phi(x,z1,Q2,2) + L*phi(x,z1,Q2,1);
double bt2y = L*L*phi(y,z2,Q2,2) + L*phi(y,z2,Q2,1);
if(sqrt(Q2) != 125.){
bt2x += phi(x,z1,Q2,0);
bt2y += phi(y,z2,Q2,0);
}
return bt1x*bt2y + bt2x*bt1y;
}
}
- else return 0.;
+ return 0.;
}
double polevl(double x, double* coef, int N )
{
double ans;
int i;
double *p;
p = coef;
ans = *p++;
i = N;
do
ans = ans * x + *p++;
while( --i );
return ans;
}
double dilog(double x)
{
static double cof_A[8] = {
4.65128586073990045278E-5,
7.31589045238094711071E-3,
1.33847639578309018650E-1,
8.79691311754530315341E-1,
2.71149851196553469920E0,
4.25697156008121755724E0,
3.29771340985225106936E0,
1.00000000000000000126E0,
};
static double cof_B[8] = {
6.90990488912553276999E-4,
2.54043763932544379113E-2,
2.82974860602568089943E-1,
1.41172597751831069617E0,
3.63800533345137075418E0,
5.03278880143316990390E0,
3.54771340985225096217E0,
9.99999999999999998740E-1,
};
if( x >1. ) {
return -dilog(1./x)+M_PI*M_PI/3.-0.5*sqr(log(x));
}
x = 1.-x;
double w, y, z;
int flag;
if( x == 1.0 )
return( 0.0 );
if( x == 0.0 )
return( M_PI*M_PI/6.0 );
flag = 0;
if( x > 2.0 )
{
x = 1.0/x;
flag |= 2;
}
if( x > 1.5 )
{
w = (1.0/x) - 1.0;
flag |= 2;
}
else if( x < 0.5 )
{
w = -x;
flag |= 1;
}
else
w = x - 1.0;
y = -w * polevl( w, cof_A, 7) / polevl( w, cof_B, 7 );
if( flag & 1 )
y = (M_PI * M_PI)/6.0 - log(x) * log(1.0-x) - y;
if( flag & 2 )
{
z = log(x);
y = -0.5 * z * z - y;
}
return y;
}
double trilog(double x)
{
double z = x;
return trilog_(&z);
}
double LUMI::Pqg(const double x)
{
return 0.5*(x*x+pow((1.0-x),2));
}
double LUMI::Pqq_pt(const double x)
{
return 4./3.*(1.+x*x);
}
double LUMI::Pqq_plus(const double x)
{
return 1./(1.-x);
}
double LUMI::Pqq_psub(const double x)
{
return -log(1. - x);
}
double LUMI::Pqq_d()
{
return 2.;
}
double LUMI::Pgq(const double x)
{
return 4./3.*(1.+(1.-x)*(1.-x))/x;
}
double LUMI::Pgg_pt(const double x)
{
return 6.*x;
}
double LUMI::Pgg_plus(const double x)
{
return Pqq_plus(x);
}
double LUMI::Pgg_psub(const double x)
{
return -log(1. - x);
}
double LUMI::Pgg_d()
{
return (33.-10.)/6.;
}
double LUMI::Pgg(const double x)
{
return 6.*((1.-x)/x + x*(1.-x));
}
double LUMI::phi(const double x, const double z,
const double Q2, const int o)
{
double sigma(0.);
for(size_t i(0);i<4;++i){
sigma+=m_pdf->xfxQ2(i+1,x/z,Q2)+m_pdf->xfxQ2(-(i+1),x/z,Q2);
}
return agb(z,2,o)*m_pdf->xfxQ2(0,x/z,Q2) + asb(z,2,o)*sigma;
}
double LUMI::agb(const double z, const int fo, const int lo)
{
double CF(4./3.), CA(3.), TR(0.5), lz(log(z)),lIz(log(1.-z)),
z2(z*z);
double res(0.), cftr(0.), catr(0.), tr2(0.);
if(fo==1 && lo==1) res=Pqg(z);
if(fo==2 && lo==0) {
double Iz(1. - z);
int n(1),p(2);
cftr = 2.*Pqg(z)*(8.*ZETA3 + 4./3.*pow(lIz,3) - 8.*lIz*dilog(Iz)
+8.*ZETA2*lz - 4.*lz*pow(lIz,2) + 2./3.*pow(lz,3)
- 8.*lz*dilog(1.-z) + 8.*trilog(1.-z) - 24.*wgplg_(&n,&p,&Iz))
+ z2*(-16.*ZETA2*lz + 4./3.*pow(lz,3) + 16.*lz*dilog(Iz)+32*wgplg_(&n,&p,&Iz))
- (4. + 96.*z - 64.*z2)*dilog(Iz) - (4.-48.*z + 40.*z2)*ZETA2 - (8.+48.*z - 24.*z2)*lz*lIz
+ (4.+8.*z- 12.*z2)*pow(lIz,2) - (1.+12.*z-20.*z2)*pow(lz,2)
- (52.*z - 48.*z2)*lIz - (16.+18.*z+48.*z2)*lz
+ 26. - 82.*z + 80.*z2;
double mz = -z;
catr = 2.*Pqg(z)*(-4./3.*pow(lIz,3) + 8.*lIz*dilog(Iz) - 8.*trilog(Iz))
+(1. + 2.*z + 2.*z2)*(-8.*ZETA2*log(1.+z) - 16.*log(1.+z)*dilog(-z)
-8.*lz*pow(log(1.+z),2) + 4.*pow(lz,2)*log(1.+z) + 8.*lz*dilog(-z)
-8.*trilog(-z) - 16.*wgplg_(&n,&p,&mz))
+ (16. + 64.*z)*(2.*wgplg_(&n,&p,&Iz)+lz*dilog(Iz))
-(4./3. + 8./3.*z)*pow(lz,3)
+ (8.-32.*z + 16.*z2)*ZETA3
-(16.+64.*z)*ZETA2*lz
+ (16.*z + 16.*z2)*(dilog(-z) +
lz*log(1.+z))
+ (32./3./z + 12. + 64.*z - 272./3.*z2)*dilog(Iz)
- (12. + 48.*z - 260./3.*z2 + 32./3./z)*ZETA2
- 4.*z2*lz*lIz
- (2.+8.*z- 10.*z2)*pow(lIz,2)
+ (2. + 8.*z + 46./3.*z2)*pow(lz,2)
+ (4.+16.*z - 16.*z2)*lIz
- (56./3. + 172./3.*z + 1600./9.*z2)*lz
- 448./27./z - 4./3. - 628./3.*z + 6352./27.*z2;
res = (CF*TR*cftr + CA*TR*catr)/8.; // 8 is = 4 wrt Maria & 2 wrt to Buza...
}
if(fo==2 && lo==1) {
cftr = 8.*Pqg(z)*(2.*lz*lIz - lIz*lIz + 2.*ZETA2)
-(2. - 4.*z + 8.*z2)*lz*lz - 16.*z*(1.-z)*lIz
-(6. - 8.*z + 16.*z2)*lz - 28. + 58.*z - 40.*z2;
catr = (8. + 16.*z + 16.*z2)*(dilog(-z)+log(z)*log(1+z))
+ 8.*Pqg(z)*lIz*lIz
+ (4. + 8.*z)*lz*lz + 16.*z*ZETA2 + 16.*z*(1.-z)*lIz
- (4. + 32.*z + 176.*z2/3.)*lz
- 80./(9.*z) + 8. - 100.*z + 872.*z2/9.;
res = - (CF*TR*cftr + CA*TR*catr)/4.;
}
if(fo==2 && lo==2){
cftr = 8.*Pqg(z)*lIz - (2. - 4.*z + 8.*z2)*lz - (1. - 4.*z);
catr = 8.*Pqg(z)*lIz + (4. + 16.*z)*lz
+ 8./(3.*z) + 2. + 16.*z - 62.*z2/3.;
tr2 = -8.*(z2+(1.-z)*(1.-z))/3.;
double tr(0.);
tr = 8.*Pqg(z)/3.; // 8./3.*Pqg(z);
res = (CF*TR*cftr - CA*TR*catr + TR*TR*tr2 + TR*tr)/4.;
}
return res;
}
double LUMI::asb(const double z, const int fo, const int lo)
{
double CF(4./3.), TR(0.5), lz(log(z)), z2(z*z);
- double res(0.), cftr(0.), lIz(log(1.-z));
+ double res(0.), cftr(0.);
if(fo==2 && lo==0) {
double a1(0.), a2(0.), a3(0.), a4(0.), a5(0.),
a6(0.);
double Iz(1. - z);
int n(1),p(2);
a1 = 32.*wgplg_(&n,&p,&Iz) + 16.*lz*dilog(Iz) - 16.*ZETA2*lz
- 4./3.*pow(lz,3);
a2 = 32./3./z + 8. -8.*z - 32./3.*z2;
a3 = -a2;
a4 = 2. + 10.*z+ 16./3.*z2;
a5 = -(56./3.+88./3.*z + 448./9.*z2);
a6 = -448./27./z - 4./3. - 124./3.*z + 1600./27.*z2;
res =CF*TR*((1.+z)*a1 + dilog(Iz)*a2 + ZETA2*a3 + pow(lz,2)*a4
+ lz*a5 + a6);
res /= 8.;
// res = 0.;
}
if(fo==2 && lo==2) {
cftr = - 4.*(1. + z)*lz - 8./(3.*z) - 2. + 2.*z + 8.*z2/3. ;
res = (CF*TR*cftr)/4.;
}
if(fo==2 && lo==1){
cftr =- 4.*(1. + z)*lz*lz + (4. + 20.*z + 32.*z2/3.)*lz
+ 80./(9.*z) - 8. + 24.*z - 224.*z2/9.;
res = (CF*TR*cftr)/4.;
}
return res;
}
-
Index: branches/bbz/src/5F.C
===================================================================
--- branches/bbz/src/5F.C (revision 37)
+++ branches/bbz/src/5F.C (revision 38)
@@ -1,227 +1,225 @@
////
// Calculate the massless scheme at LO, NLO and NNLO
///
#include "Massless_lim.H"
#include "Massless_Kernels.H"
#include <iostream>
#include <vector>
#include <cmath>
using namespace LHAPDF;
using namespace std;
double Massless_Limit::dsigmaLO(double x[], size_t dim, void* p)
{
double dsg(0.);
double tau, jac(1.);
tau = tauh + (1. - tauh)*x[0];
jac *= (1. - tauh);
dsg = 2.*jac*p_lumi->Lumi(10,0,tau,tauh/tau,1.,1.,10)/tau;
return s0*dsg;
}
double Massless_Limit::dsigmaNLO(double z[], size_t dim,
void* p, const size_t plus)
// here plus means should I do the subtraction? yes=1,no=0
{
- double as = (m_proc)?p_lumi->ass(m2r)/2.:
- 2.*p_lumi->ass(m2r);
+ double as = 2.*M_PI*p_lumi->ass(m2r);
double dsg(0.);
// delta contr
if(dim==1){
double tau, jac(1.);
double sbb(0.);
tau = tauh + (1. - tauh)*z[0];
jac *= (1. - tauh);
if(m_proc==0)
sbb = delta1_(&lfr);
else if(m_proc==1) // bb->Z
sbb = delta1_(&lfh);;
dsg = 2.*as*
sbb*jac*p_lumi->Lumi(10,0,tau,tauh/tau,1.,1.,10)/tau;
return s0*dsg;
}
// main thing
if(dim==2 && plus == 0){
double tau, jac(1.), xi, x, y;
double sbbp(0.), sbb(0.), sbg(0.);
tau = tauh + (1. - tauh)*z[0];
jac *= (1. - tauh);
double ximax = - 0.5*log(tau);
xi = -ximax + 2.*ximax*z[1];
jac *= 2.*ximax;
// standard variables x, y ;
x = sqrt(tau)*exp(xi);
y = sqrt(tau)*exp(-xi);
sbbp = dterms1_(&y, &lfh);
sbb = sbbq1_(&y, &lfh);
sbg = sbg1_(&y,&lfh); /////////
dsg =
as*jac*(
2.*p_lumi->Lumi(10,0,x,tauh/tau,1.,1.,10)*
(sbbp*y + sbb)
- 2.*p_lumi->Lumi(10,0,x,tauh/x,1.,1.,10)*(y*sbbp)
+ (p_lumi->Lumi(10,5,x,tauh/tau,1.,1.,10) +
p_lumi->Lumi(5,10,x,tauh/tau,1.,1.,10) +
p_lumi->Lumi(0,5,x,tauh/tau,1.,1.,10) +
p_lumi->Lumi(5,0,x,tauh/tau,1.,1.,10))*(sbg)
)/tau;
return s0*dsg;
}
if(dim==2 && plus == 1){
double tau, jac(1.), xi, x, y;
double sbbp(0.);
tau = tauh + (1. - tauh)*z[0];
jac *= (1. - tauh);
xi = (tauh/tau)*z[1];
jac *= (tauh/tau);
x = tau; y = xi;
sbbp = dterms1_(&y, &lfh);
dsg = 2.*jac*as*sbbp*p_lumi->Lumi(10,0,x,tauh/x,1.,1.,10)/x;
return s0*dsg;
}
else if (dim > 2){
cout << " No more than 2 dim for dsigmaNLO! " << std::endl;
return 0.;
}
return s0*dsg;
}
double Massless_Limit::dsigmaNNLO(double z[], size_t dim,
void* p, const size_t plus)
// here plus means should I do the subtraction? yes=1,no=0
{
//cout << " MODE:: "<< m_mode << std::endl;
- double as = (m_proc)?p_lumi->ass(m2r)/2.:
- 2.*p_lumi->ass(m2r);
+ double as = 2.*M_PI*p_lumi->ass(m2r);
double a2s = pow(as,2);
double dsg(0.), nf(5.);
// delta contr
if(dim==1){
- double tau, jac(1.); double y(1.);
- double sbba(0.), sbbf(0.);
+ double tau, jac(1.);
+ double sbba(0.);
tau = tauh + (1. - tauh)*z[0];
jac *= (1. - tauh);
sbba = delta2_(&lfr, &lfh);
dsg = 2.*a2s*sbba *
jac*p_lumi->Lumi(10,0,tau,tauh/tau,1.,1.,10)/tau;
return s0*dsg;
}
// main thing
if(dim==2){
if (plus == 0) {
double tau, jac(1.), xi, x, y;
double lbb, lbb0, lbg, lbq, lgg,lqq, lbs;
double sbbp(0.), sbba(0.), sbbf(0.),
sbg(0.), sbq(0.), sgg(0.), sbs(0.), sqq(0.);
tau = tauh + (1. - tauh)*z[0];
jac *= (1. - tauh);
double ximax = - 0.5*log(tau);
xi = -ximax + 2.*ximax*z[1];
jac *= 2.*ximax;
// standard variables x, y ;
x = sqrt(tau)*exp(xi);
y = sqrt(tau)*exp(-xi);
//soft-terms
sbbp = dterms2_(&y,&lfh,&lfr);
// finite-terms
//b b~
sbba = deltabbqa_(&y,&lfh,&lfr);
sbbf = deltabbqf_(&y,&lfh,&lfr);
// // b g
sbg = deltabga_(&y,&lfh,&lfr) +
nf* deltabgf_(&y,&lfh,&lfr);
// // g g
sgg = deltagga_(&y,&lfh,&lfr);
// // b b
sbs = deltabba_(&y,&lfh,&lfr);
// // b q
sbq = deltabqa_(&y,&lfh,&lfr);
// // q q~
sqq = deltaqqba_(&y,&lfh,&lfr);
// Luminosities
lbb = 2.*p_lumi->Lumi(10,0,x,tauh/tau,1.,1.,10);
lbb0 = 2.*p_lumi->Lumi(10,0,x,tauh/x,1.,1.,10);
lbg = p_lumi->Lumi(10,5,x,tauh/tau,1.,1.,10) +
p_lumi->Lumi(5,10,x,tauh/tau,1.,1.,10) +
p_lumi->Lumi(0,5,x,tauh/tau,1.,1.,10) +
p_lumi->Lumi(5,0,x,tauh/tau,1.,1.,10);
lbq = 0.;
lbq += 2.*( p_lumi->Lumi(10,4,x,tauh/tau,1.,1.,10) // d
+ p_lumi->Lumi(10,6,x,tauh/tau,1.,1.,10));
lbq += 2.*( p_lumi->Lumi(4,10,x,tauh/tau,1.,1.,10)
+ p_lumi->Lumi(6,10,x,tauh/tau,1.,1.,10));
lbq += 2.*( p_lumi->Lumi(10,3,x,tauh/tau,1.,1.,10) // u
+ p_lumi->Lumi(10,7,x,tauh/tau,1.,1.,10));
lbq += 2.*( p_lumi->Lumi(3,10,x,tauh/tau,1.,1.,10)
+ p_lumi->Lumi(7,10,x,tauh/tau,1.,1.,10));
lbq += 2.*( p_lumi->Lumi(10,2,x,tauh/tau,1.,1.,10) // s
+ p_lumi->Lumi(10,8,x,tauh/tau,1.,1.,10));
lbq += 2.*( p_lumi->Lumi(3,10,x,tauh/tau,1.,1.,10)
+ p_lumi->Lumi(8,10,x,tauh/tau,1.,1.,10));
// (c = c~)
lbq += 4.*( p_lumi->Lumi(10,9,x,tauh/tau,1.,1.,10) );
lbq += 4.*( p_lumi->Lumi(9,10,x,tauh/tau,1.,1.,10) );
lgg = p_lumi->Lumi(5,5,x,tauh/tau,1.,1.,10);
lqq = L_qq(x,tauh/tau);
lbs = 2.*p_lumi->Lumi(10,10,x,tauh/tau,1.,1.,10);
dsg =
a2s*jac*( lbb * ( sbbp*y + sbba + nf * sbbf)
- lbb0 * sbbp *y
+ lbg * sbg
+ lgg * sgg
+ lbs * sbs
+ lbq * sbq
+ lqq * sqq
)/tau;
return s0*dsg;
}
if(plus == 1){
double tau, jac(1.), xi, x, y;
double sbbp(0.);
tau = tauh + (1. - tauh)*z[0];
jac *= (1. - tauh);
xi = (tauh/tau)*z[1];
jac *= (tauh/tau);
x = tau; y = xi;
sbbp = dterms2_(&y,&lfh,&lfr);
dsg = 2.*jac*a2s*sbbp*p_lumi->Lumi(10,0,x,tauh/x,1.,1.,10)/x;
return s0*dsg;
}
}
else if (dim > 2){
cout << " No more than 2 dim for dsigmaNLO! " << std::endl;
return 0.;
}
return s0*dsg;
}
Index: branches/bbz/src/PDF.H
===================================================================
--- branches/bbz/src/PDF.H (revision 37)
+++ branches/bbz/src/PDF.H (revision 38)
@@ -1,58 +1,57 @@
#ifndef PDF_H
#define PDF_H
#include "LHAPDF/LHAPDF.h"
#include <cmath>
#include <vector>
class LUMI {
private :
LHAPDF::PDF *m_pdf;
- // std::vector< std::vector<double> > m_lumi;
double m_as,Q,m2b;
double mm2b;
bool m_flag, m_lim;
public :
LUMI(const char* setname, const int setmem, const double muF,
const double mbext, const bool ismlim);
~LUMI();
double Lumi(const int i, const int j,
const double x, const double y,
const double z1, const double z2,
const int order);
inline double alphaS() { return m_as; }
inline double ass(const double Q2) { return m_pdf->alphasQ(sqrt(Q2))/(2.*M_PI); }
double mb(const double mf, const double Q,const int nloop);
inline double m2bpdf() { return m2b; }
inline double m2bext() {return mm2b; }
double Pqg(const double x);
double Pgq(const double x);
double Pgg(const double x); // finite bits
double Pgg_d(); // delta term P_gg
double Pgg_plus(const double x); // plus distribution
double Pgg_pt(const double x); // plus coefficient
double Pgg_psub(const double x); // plus reminder
double Pqq_plus(const double x); // same as above
double Pqq_pt(const double x);
double Pqq_psub(const double x);
double Pqq_d();
double phi(const double x, const double z,
const double Q2, const int o);
double agb(const double z, const int fo, const int lo);
double asb(const double z, const int fo, const int lo);
void SetmuF(const double muf) {Q=muf;}
bool Init_Lumi() { return m_flag;}
};
#endif
double dilog(double x);
double trilog(double x);

File Metadata

Mime Type
text/x-diff
Expires
Thu, Apr 24, 6:36 AM (1 d, 19 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4887770
Default Alt Text
(122 KB)

Event Timeline