-
Notifications
You must be signed in to change notification settings - Fork 0
/
dSigmadptN.cpp
155 lines (133 loc) · 4.72 KB
/
dSigmadptN.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
/*
* =====================================================================================
*
* Filename: main.cpp
*
* Description: Main file that computes either the full hadronic or the
* Mellin partonic cross section.
*
* Version: 1.0
* Created: 17/02/2021 23:16:58
* Revision: none
* Compiler: gcc
*
* Author: Tanjona R. Rabemananjara, Roy Stegeman
* Organization: N3PDF
*
* =====================================================================================
*/
#include <stdio.h>
#include <stdlib.h>
#include <exception>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <string>
#include <vector>
#include "higgs-fo/hadronic.h"
#include "higgs-fo/higgspt.h"
#include "higgs-fo/higgsptpartonic.h"
#include "higgs-fo/luminosity.h"
#include "higgs-fo/params.h"
#include "higgs-fo/partonic.h"
#include "include/CombinedRes.h"
#include "include/SmallptExp.h"
#include "include/ThresExp.h"
#include "include/HighEnergyExp.h"
#include "yaml-cpp/yaml.h"
// Exception for wrong inputs
struct err_message : public std::exception {
const char *what() const throw() { return "Wrong Parameters!!"; }
};
int main(int argc, char *argv[]) {
LHAPDF::setVerbosity(0);
YAML::Node node = YAML::LoadFile(argv[1]);
int inorm = node["inorm"].as<int>();
int order = node["order"].as<int>();
int _nf = node["nf"].as<int>();
int channel = node["channel"].as<int>();
int scheme = node["scheme"].as<int>();
double _mh = node["mh"].as<double>();
double _mur = node["mur"].as<double>();
double _muf = node["muf"].as<double>();
double _sroot = node["sroot"].as<double>();
long double Nmin = node["Nmin"].as<long double>();
long double Nmax = node["Nmax"].as<long double>();
long double Nbin = node["Nbin"].as<long double>();
long double pt = node["pt"].as<long double>();
std::string pdfname = node["pdfname"].as<std::string>();
std::string filename = node["outfile"].as<std::string>();
std::string ord_fixod[2] = {"_LO", "_NLO"};
std::string par_chanl[5] = {"_gg_channel", "_gq_channel", "_qq_channel",
"_qqb_channel", "_all_channels"};
std::string matsch[3] = {"_smallpt_asN.dat", "_threshold_asN.dat",
"_combined_asN.dat"};
try {
if (order < 0 || order > 1) throw err_message();
if (channel < 0 || channel > 5) throw err_message();
filename = ord_fixod[order];
filename = par_chanl[channel];
filename = matsch[scheme];
} catch (err_message &err) {
std::cout << err.what() << std::endl;
exit(EXIT_FAILURE);
}
// Factors for Born cross-section
double factor;
double gf = 1.16637e-5; // Fermi Constant
double gevpb = 3.8937966e8; // GeV to pb
if (inorm == 1) {
std::cout << "ERROR in inorm!" << std::endl;
exit(EXIT_FAILURE); // TODO: complete implementation!
} else if (inorm == 0) {
factor = gf / 288. / M_PI / sqrt(2.); // large-top mass limit
} else {
std::cout << "ERROR in inorm!" << std::endl;
exit(EXIT_FAILURE);
}
// Initialize PDF and extract alpha_s
LHAPDF::initPDFSetByName(pdfname);
double _as = LHAPDF::alphasPDF(_mur);
double _sigma0 = factor * gevpb * std::pow(_as, 2);
// Define parameters
PhysParams physparam;
physparam.nc = 3;
physparam.nf = _nf;
physparam.mh = _mh;
physparam.mur = _mur;
physparam.muf = _muf;
physparam.alphas = _as;
physparam.sroot = _sroot;
physparam.sigma0 = _sigma0;
// Init. combined resummation class
CombinedRes combres(order, channel, pdfname, &physparam);
// Construct output fie
std::ofstream output_file(filename);
output_file << "# PDF set name : " << pdfname << "\n"
<< "# Matching scheme : " << scheme << "\n"
<< "# Fixed Order : " << order << "\n"
<< "# Partonic channel : " << channel << "\n"
<< "# Center of M.E. (GeV) : " << _sroot << "\n"
<< "# Higgs mass (GeV) : " << _mh << "\n"
<< "# Renorm. scale (GeV) : " << _mur << "\n"
<< "# Fact. scale (GeV) : " << _muf << "\n";
const int space = 16;
output_file << "# [N values]" << std::setw(space) << "[dHpt (pb)]"
<< "\n";
long double nn = Nmin;
std::complex<long double> results;
while (nn <= Nmax) {
std::complex<long double> Ncmpx(nn, 0.);
results = combres.CombinedResExpr(Ncmpx, pt, scheme);
// Generate some output logs & write to output file
printf("N=%Le: dHdpt = %Le %Le II. \n", nn, results.real(),
results.imag());
output_file.setf(std::ios_base::scientific);
output_file << nn << std::setw(space) << results.real() << "\n";
output_file.flush();
nn = Nbin;
}
output_file.close();
return 0;
}