Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

can be used to produce results in x-space #11

Open
wants to merge 5 commits into
base: main
Choose a base branch
from
Open

can be used to produce results in x-space #11

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
bc93697
save
RoyStegeman Jun 25, 2021
1df5c5a
save
RoyStegeman Jun 26, 2021
fa94d6a
save
RoyStegeman Jun 27, 2021
89c74af
fix x-space agreement with HE
RoyStegeman Jun 28, 2021
42e1e96
Committing clang-format changes
RoyStegeman Jul 1, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 4 additions & 4 deletions higgsfodptN.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -145,18 +145,18 @@ int main(int argc, char* argv[]) {
if (sectype == "hadronic") {
results = higgshard.higgsdpt(pt,y1,y2);
} else if (sectype == "partonic") {
results = higgspart.partonichiggsdpt(pt,nn);
results = higgspart.partonichiggsdpt(pt,std::pow(10,nn));
} else {
std::cout << "Error in entry sectype!" << std::endl;
exit(EXIT_FAILURE);
}

// Generate some output logs & write to output file
printf("N=%e: dHdpt = %e +- %e. \n", nn, results[0], results[1]);
printf("N=%e: dHdpt = %e +- %e. %e\n", std::pow(10,nn), results[0], results[1], results[2]);
output_file.setf(std::ios_base::scientific);
output_file << nn << std::setw(space)
output_file << std::pow(10,nn) << std::setw(space)
<< results[0] << std::setw(space)
<< results[1] << "\n";
<< "\n";
output_file.flush();

nn += Nbin;
Expand Down
4 changes: 2 additions & 2 deletions runcards/Higgs-FO-as-N.yaml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,8 +11,8 @@ order : 1

# N bins, now x-bins
Nmin : -12
Nmax : -3
Nbin : 0.1
Nmax : -1
Nbin : 0.5

# relevant scales in GeV
mh : 125 # mass of colorless final state
Expand Down
109 changes: 54 additions & 55 deletions src/higgsptpartonic.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -692,7 +692,7 @@ double HiggsDpTpartonic::deltapartonic(double pt, double nn, double zz) {
///////////////////////////////////////////////////////////////
double xx = nn;
double result = 0;
zz=1; // unused integration variable
zz = 1; // unused integration variable

// This function calculates the terms of the cross section proportional to
// delta(Q²)
Expand Down Expand Up @@ -740,22 +740,21 @@ double HiggsDpTpartonic::deltapartonic(double pt, double nn, double zz) {
// part of the xsec that is also calculated in this function.
case (0): // gg-channel
{
// dominant contribution in small-x are the terms multiplied by uu.
// specifically the logs (not Li2) dominate uu
// dominant contribution in small-x are the terms multiplied by uu.
// specifically the logs (not Li2) dominate uu
double uu =
0.5 * std::pow(log(uh / th), 2) + std::pow(M_PI, 2) / 3. -
log(sh / MH2) * log(-th / MH2) - log(sh / MH2) * log(-uh / MH2) -
log(-uh / MH2) * log(-th / MH2) + std::pow(log(MH2 / sh), 2) +
std::pow(log(MH2 / (MH2 - th)), 2) +
std::pow(log(MH2 / (MH2 - uh)), 2) + 2. * Li2(1. - MH2 / sh) +
2. * Li2(MH2 / (MH2 - th)) + 2. * Li2(MH2 / (MH2 - uh));
uu =
0.5 * std::pow(log(uh / th), 2) + std::pow(M_PI, 2) / 3. -
log(sh / MH2) * log(-th / MH2) - log(sh / MH2) * log(-uh / MH2) -
log(-uh / MH2) * log(-th / MH2) + std::pow(log(MH2 / sh), 2) +
std::pow(log(MH2 / (MH2 - th)), 2) +
std::pow(log(MH2 / (MH2 - uh)), 2) + 2. * Li2(1. - MH2 / sh) +
2. * Li2(MH2 / (MH2 - th)) + 2. * Li2(MH2 / (MH2 - uh));
uu = 0.5 * std::pow(log(uh / th), 2) + std::pow(M_PI, 2) / 3. -
log(sh / MH2) * log(-th / MH2) - log(sh / MH2) * log(-uh / MH2) -
log(-uh / MH2) * log(-th / MH2) + std::pow(log(MH2 / sh), 2) +
std::pow(log(MH2 / (MH2 - th)), 2) +
std::pow(log(MH2 / (MH2 - uh)), 2) + 2. * Li2(1. - MH2 / sh) +
2. * Li2(MH2 / (MH2 - th)) + 2. * Li2(MH2 / (MH2 - uh));
double de = 1.5 * beta0 * (log(-MUR2 / th) + log(-MUR2 / uh)) +
67. / 6. - 5. / 9. * NF;
result += aass / (2. * M_PI) *
Expand Down Expand Up @@ -860,9 +859,9 @@ double HiggsDpTpartonic::distrpartonic(double pt, double nn, double zz1,
// multiplied by either delta(Q) or delta(Q^2) (both singular //
// and non-Singular) //
////////////////////////////////////////////////////////////////
double nonsingular = 0, a1 = 0, b1 = 0, c1 = 0, a10 = 0, b10 = 0;
double nonsingular = 0, a1 = 0, b1 = 0, c1 = 0, a10 = 0, b10 = 0, d10 = 0;

zz2=1; //unused integration variable
zz2 = 1; // unused integration variable
double qq = zz1; // qq = QQ2/QQ2max is an integration variable used to
// integrate out rapidity
double xx = nn; // xx = Q²/sh
Expand Down Expand Up @@ -891,15 +890,16 @@ double HiggsDpTpartonic::distrpartonic(double pt, double nn, double zz1,
// double a1factor = log(1. - za) / (1. - za);
// double b1factor = 1. / (1. - za);
// Alternatively one can write `a1factor' and `b1factor' in terms of qq:
double a1factor = (log(qq)/qq+log(QQ2max*za/-th)/qq)/QQ2max*(-th/za);
double b1factor = 1./qq/QQ2max*(-th/za);
double a1factor =
(log(qq) / qq + log(QQ2max * za / -th) / qq) / QQ2max * (-th / za);
double b1factor = 1. / qq / QQ2max * (-th / za);

coeff(pt, uh, th, sh, MH2);
REG(pt, uh, th, sh, MH2);

double shnew = za * sh;
double thnew = th;
double uhnew = za*(uh-MH2)+MH2;
double uhnew = za * (uh - MH2) + MH2;

switch (CHANNEL) {
// a1:: (log(1-za)/(1-za))+ terms
Expand Down Expand Up @@ -988,14 +988,20 @@ double HiggsDpTpartonic::distrpartonic(double pt, double nn, double zz1,
za = -th / (QQ2 - th);
double jac10 = QQ2max / sqrt(std::pow(sh + MH2 - QQ2, 2) - 4. * sh * mt2);

double a10factor =
double apolefactor =
(log(qq) / qq + log(QQ2max * za / -th) / qq) / QQ2max * (-th / za);
double b10factor = 1. / qq / QQ2max * (-th / za);
double adeltafactor =
0.5 * std::pow(log(-QQ2max / th), 2) / QQ2max * (-th / za);
double a10factor = apolefactor - adeltafactor;
double bdeltafactor = log(-QQ2max / th) / QQ2max * (-th / za);
double bpolefactor = 1. / qq / QQ2max * (-th / za);
double b10factor = bpolefactor - bdeltafactor;

double d10factor = -th / za / QQ2max;

shnew = za * sh;
thnew = th;
uhnew = za*(uh-MH2)+MH2;

uhnew = za * (uh - MH2) + MH2;

coeff(pt, uh, th, sh, MH2);

Expand All @@ -1008,6 +1014,8 @@ double HiggsDpTpartonic::distrpartonic(double pt, double nn, double zz1,
(1. / th * pgg(za) * log(-MUF2 * za / th) *
gg0(shnew, thnew, uhnew, MH2) +
za / th * big1 * log((QQ2 + pt * pt) * za / (-th)) + za / th * big2);
d10 += 1. / th * beta0 * log(-MUF2 * za / th) *
gg0(shnew, thnew, uhnew, MH2);
} break;
case (1): // gq-channel
{
Expand Down Expand Up @@ -1039,21 +1047,13 @@ double HiggsDpTpartonic::distrpartonic(double pt, double nn, double zz1,
} break;
}

double adeltafactor =
0.5 * std::pow(log(-QQ2max / th), 2) / QQ2max * (-th / za);
double bdeltafactor = log(-QQ2max / th) / QQ2max * (-th / za);

double bfinal = b1 * jac1 * b1factor - b10 * jac10 * b10factor +
b10 * jac10 * bdeltafactor;
double afinal = a1 * jac1 * a1factor - a10 * jac10 * a10factor +
a10 * jac10 * adeltafactor;
double bfinal = b1 * jac1 * b1factor - b10 * jac10 * b10factor;
double afinal = a1 * jac1 * a1factor - a10 * jac10 * a10factor;
double cfinal = c1 * jac1;
double dfinal = d10 * jac10 * d10factor;
double nonsingularfinal = nonsingular * jac1;



// nonsingular = 0;
double result = afinal + bfinal + cfinal + nonsingularfinal;
double result = afinal + bfinal + cfinal + dfinal + nonsingularfinal;

// dsigma/pt² to dsigma/dpt
result *= 2. * pt;
Expand All @@ -1070,17 +1070,16 @@ double HiggsDpTpartonic::distrpartonic(double pt, double nn, double zz1,

//==============================================================================================//


double HiggsDpTpartonic::distrpartoniccross(double pt, double nn, double zz1,
double zz2) {
double zz2) {
////////////////////////////////////////////////////////////////
// This function computes the remaining terms that are not //
// multiplied by either delta(Q) or delta(Q^2) (both singular //
// and non-Singular) //
////////////////////////////////////////////////////////////////
double nonsingular = 0, a1 = 0, b1 = 0, c1 = 0, a10 = 0, b10 = 0;
double nonsingular = 0, a1 = 0, b1 = 0, c1 = 0, a10 = 0, b10 = 0, d10 = 0;

zz2=1; //unused integration variable
zz2 = 1; // unused integration variable
double qq = zz1; // qq = QQ2/QQ2max is an integration variable used to
// integrate out rapidity
double xx = nn; // xx = Q²/sh
Expand Down Expand Up @@ -1109,15 +1108,16 @@ double HiggsDpTpartonic::distrpartoniccross(double pt, double nn, double zz1,
// double a1factor = log(1. - za) / (1. - za);
// double b1factor = 1. / (1. - za);
// Alternatively one can write `a1factor' and `b1factor' in terms of qq:
double a1factor = (log(qq)/qq+log(QQ2max*za/-th)/qq)/QQ2max*(-th/za);
double b1factor = 1./qq/QQ2max*(-th/za);
double a1factor =
(log(qq) / qq + log(QQ2max * za / -th) / qq) / QQ2max * (-th / za);
double b1factor = 1. / qq / QQ2max * (-th / za);

coeff(pt, uh, th, sh, MH2);
REG(pt, uh, th, sh, MH2);

double shnew = za * sh;
double thnew = th;
double uhnew = za*(uh-MH2)+MH2;
double uhnew = za * (uh - MH2) + MH2;

switch (CHANNEL) {
// a1:: (log(1-za)/(1-za))+ terms
Expand Down Expand Up @@ -1206,14 +1206,20 @@ double HiggsDpTpartonic::distrpartoniccross(double pt, double nn, double zz1,
za = -th / (QQ2 - th);
double jac10 = QQ2max / sqrt(std::pow(sh + MH2 - QQ2, 2) - 4. * sh * mt2);

double a10factor =
double apolefactor =
(log(qq) / qq + log(QQ2max * za / -th) / qq) / QQ2max * (-th / za);
double b10factor = 1. / qq / QQ2max * (-th / za);
double adeltafactor =
0.5 * std::pow(log(-QQ2max / th), 2) / QQ2max * (-th / za);
double a10factor = apolefactor - adeltafactor;
double bdeltafactor = log(-QQ2max / th) / QQ2max * (-th / za);
double bpolefactor = 1. / qq / QQ2max * (-th / za);
double b10factor = bpolefactor - bdeltafactor;

double d10factor = -th / za / QQ2max;

shnew = za * sh;
thnew = th;
uhnew = za*(uh-MH2)+MH2;

uhnew = za * (uh - MH2) + MH2;

coeff(pt, uh, th, sh, MH2);

Expand All @@ -1226,6 +1232,8 @@ double HiggsDpTpartonic::distrpartoniccross(double pt, double nn, double zz1,
(1. / th * pgg(za) * log(-MUF2 * za / th) *
gg0(shnew, thnew, uhnew, MH2) +
za / th * big1 * log((QQ2 + pt * pt) * za / (-th)) + za / th * big2);
d10 += 1. / th * beta0 * log(-MUF2 * za / th) *
gg0(shnew, thnew, uhnew, MH2);
} break;
case (1): // gq-channel
{
Expand Down Expand Up @@ -1257,21 +1265,13 @@ double HiggsDpTpartonic::distrpartoniccross(double pt, double nn, double zz1,
} break;
}

double adeltafactor =
0.5 * std::pow(log(-QQ2max / th), 2) / QQ2max * (-th / za);
double bdeltafactor = log(-QQ2max / th) / QQ2max * (-th / za);

double bfinal = b1 * jac1 * b1factor - b10 * jac10 * b10factor +
b10 * jac10 * bdeltafactor;
double afinal = a1 * jac1 * a1factor - a10 * jac10 * a10factor +
a10 * jac10 * adeltafactor;
double bfinal = b1 * jac1 * b1factor - b10 * jac10 * b10factor;
double afinal = a1 * jac1 * a1factor - a10 * jac10 * a10factor;
double cfinal = c1 * jac1;
double dfinal = d10 * jac10 * d10factor;
double nonsingularfinal = nonsingular * jac1;



// nonsingular = 0;
double result = afinal + bfinal + cfinal + nonsingularfinal;
double result = afinal + bfinal + cfinal + dfinal + nonsingularfinal;

// dsigma/pt² to dsigma/dpt
result *= 2. * pt;
Expand All @@ -1287,4 +1287,3 @@ double HiggsDpTpartonic::distrpartoniccross(double pt, double nn, double zz1,
}

//==============================================================================================//

27 changes: 14 additions & 13 deletions src/integration.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,18 +50,19 @@ int integration::CubaIntegrator(
static const int nvec = 1;
// static const int spin = -1;
static const int maxpass = 5;
static const int mineval = 0;
static const int nnew = 1000;
static const int nstart = 1000;
static const int nbatch = 1000;
static const int nincrease = 500;
static const int mineval = 1000;
static const int nnew = 1e5;
static const int nmin = 1;
static const int nstart = 1e6;
static const int nbatch = 1e4;
static const int nincrease = 5e5;
static const int ldxgiven = ndim;
static const int maxeval = 100000;
static const int maxeval = 1e9;
static const int verbose = verbos;
static const double border = 0.1;
static const double maxchisq = 10.;
static const double epsabs = 1e-15;
static const double flatness = 25.;
static const double epsabs = 1e-100;
static const double flatness = 1;
static const double mindeviation = 0.25;
const char *statefile = NULL;

Expand All @@ -82,7 +83,7 @@ int integration::CubaIntegrator(
break;
case 1:
Suave(ndim, ncomp, func, par, nvec, prec, epsabs, verbose, seed, mineval,
maxeval, nnew, 5, flatness, statefile, NULL, &nregions, &neval,
maxeval, nnew, nmin, flatness, statefile, NULL, &nregions, &neval,
fail, ris, err, proba);
break;
case 2:
Expand Down Expand Up @@ -119,7 +120,7 @@ double integration::IntegrateDeltaPartonic(int method,
double(Func)(double, void *),
void *pp, double *error) {
int fail;
double prec = 1e-8;
double prec = 1e-3;
double *res = NULL;
double *err = NULL;
double *prb = NULL;
Expand Down Expand Up @@ -156,7 +157,7 @@ double integration::IntegrateDistrPartonic(int method,
double(Func)(double, double, void *),
void *pp, double *error) {
int fail;
double prec = 1e-8;
double prec = 1e-3;
double *res = NULL;
double *err = NULL;
double *prb = NULL;
Expand Down Expand Up @@ -198,7 +199,7 @@ double integration::IntegrateDeltaHadronic(int method,
double(Func)(double, double, void *),
void *pp, double *error) {
int fail;
double prec = 1e-8;
double prec = 1e-3;
double *res = NULL;
double *err = NULL;
double *prb = NULL;
Expand Down Expand Up @@ -241,7 +242,7 @@ double integration::IntegrateDistrHadronic(int method,
void *),
void *pp, double *error) {
int fail;
double prec = 1e-8;
double prec = 1e-3;
double *res = NULL;
double *err = NULL;
double *prb = NULL;
Expand Down
2 changes: 1 addition & 1 deletion src/partonic.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -88,7 +88,7 @@ std::vector<double> CrossHiggs::partonichiggsdpt(double pt, double nn) {
// y2: upper boundary of the rapidity yh //
//////////////////////////////////////////////////

int method = 0;
int method = 1;
double deltres, delterr;
double distres = 0, disterr = 0;
std::vector<double> results;
Expand Down