Page MenuHomeHEPForge
Files F19252254

No OneTemporary
Actions

Size
18 KB
Referenced Files
None
Subscribers
None

View Options

Index: amcatnlo-bridge/trunk/src/amcbridge_interface.cxx
===================================================================
--- amcatnlo-bridge/trunk/src/amcbridge_interface.cxx (revision 1162)
+++ amcatnlo-bridge/trunk/src/amcbridge_interface.cxx (revision 1163)
@@ -1,545 +1,545 @@
/*
NLO mode of aMCatNLO
*/
#include <cstdlib>
#include <iostream>
#include <cmath>
// applgrid header files
#include "appl_grid/appl_grid.h"
#include "appl_grid/lumi_pdf.h"
// LHAPDF
#include "LHAPDF/LHAPDF.h"
// amcbridge header files
#include "amcbridge_interface.h"
#include "amcatnlo_common.h"
// ROOT header files
#include "TPad.h"
// lhapdf functions
extern "C" double alphasPDF_(const double&);
extern "C" void evolvePDF_(const double&, const double&, double*);
// External hook functions for aMC@NLO control
extern void (*appl_initptr)();
extern void (*appl_fillptr)();
extern void (*appl_fillrefptr)();
extern void (*appl_termptr)();
extern void (*appl_fillrefoutptr)();
extern void (*appl_recoptr)();
// Hook functions
extern "C" bool setup_amcbridge_() {
std::cout << "setup_amcatnlo() setting up hook functions ..." << std::endl;
appl_initptr = amcbridge::init;
appl_fillptr = amcbridge::fill;
appl_fillrefptr = amcbridge::fill_ref;
appl_termptr = amcbridge::term;
appl_fillrefoutptr = amcbridge::fill_ref_out;
appl_recoptr = amcbridge::reco;
return true;
}
namespace amcbridge {
// Initialise the hook functions
bool amcbridge_ready = setup_amcbridge_();
// Parton luminosity
double lhapdf_lumi(double, double, double, int);
// Check if a file exists
bool file_exists(const std::string& s) {
if(std::FILE* testfile=std::fopen(s.c_str(),"r")) {
std::fclose(testfile);
return true;
}
else return false;
}
// Different names for reading and writing grids to avoid overwriting
appl::grid* grid_obsA = 0;
std::string grid_filenameA_in = "grid_obsA_in.root";
std::string grid_filenameA_out = "grid_obsA_out.root";
appl::grid* grid_obsB = 0;
std::string grid_filenameB_in = "grid_obsB_in.root";
std::string grid_filenameB_out = "grid_obsB_out.root";
//
// Initialization
//
void init() {
//std::cout << "C++, entering init() ..." << std::endl;
// leading_order -> Power of alphas of Born events
int const leading_order = appl_common_fixed_.bpower;
//std::cout << "bpower = " << leading_order << std::endl;
// Number of loops in the APPLgrid formalism.
// Here we store four grids, W0, WR, WF and W0, in the notation of arXiv:1110.4738.
// Always book with nloops = 3
int const nloops = 3; // 4 grids
// Define the settings for the interpolation in x and Q2.
// These are common to all the grids computed.
// (Should we implement a dynamical tuning of the following parameters?)
// Number of points for the Q2 interpolation
int NQ2 = 10;
// Max and min value of Q2
double Q2min = 100;
double Q2max = 1000000;
// Order of the polynomial interpolation in Q2
int Q2order = 3;
// Number of points for the x interpolation
int Nx = 50;
// Min and max value of x
double xmin = 2e-7;
double xmax = 1;
// Order of the polynomial interpolation in x
int xorder = 3;
// Set up the APPLgrid PDF luminosities
std::vector<int> pdf_luminosities;
pdf_luminosities.push_back(appl_common_lumi_.nlumi);
// Loop over parton luminosities
for(int ilumi=0; ilumi<appl_common_lumi_.nlumi; ilumi++) {
int nproc = appl_common_lumi_.nproc[ilumi];
pdf_luminosities.push_back(ilumi);
pdf_luminosities.push_back(nproc);
// Loop over parton-parton combinations within each luminosity
for(int iproc=0; iproc<nproc; iproc++) {
pdf_luminosities.push_back(appl_common_lumi_.lumimap[ilumi][iproc][0]);
pdf_luminosities.push_back(appl_common_lumi_.lumimap[ilumi][iproc][1]);
}
}
// Use a name for the PDF combination type ending with .config.
// This will configure from a file with the same format as the
// aMC@NLO initial_states_map.dat file unless a serialised
// vector including the combinations is also passed to the
// constructor, as is the case here.
std::string genpdf = "amcatnlo.config";
new lumi_pdf(genpdf,pdf_luminosities);
// Grid Initialization and definition of the observables.
// The default is an histogram with equally spaced bins.
// Check that the grid files exists. If so read the grid from the file,
// otherwise create a new grid from scratch.
// Observable A
if(file_exists(grid_filenameA_in)) {
std::cout << "Reading existing APPLgrid from file " << grid_filenameA_in << " ..." << std::endl;
// Open the existing grid
grid_obsA = new appl::grid(grid_filenameA_in);
grid_obsA->optimise(true);
}
else {
std::cout << "Booking grid from scratch with name " << grid_filenameA_in << " ..." << std::endl;
// Create a grid
int Nobs = appl_common_histokin_.obsA_nbins;
double obsmin = appl_common_histokin_.obsA_min;
double obsmax = appl_common_histokin_.obsA_max;
grid_obsA = new appl::grid(NQ2, Q2min, Q2max, Q2order,
Nx, xmin, xmax, xorder,
Nobs, obsmin, obsmax,
genpdf, leading_order, nloops);
// Use the reweighting function
grid_obsA->reweight(true);
// The grid is an aMC@NLO type
grid_obsA->amcatnlo();
// Optimise the grid
//grid_obsA->optimise(true);
}
// Observable B
if(file_exists(grid_filenameB_in)) {
std::cout << "Reading existing applgrid from file " << grid_filenameB_in << " ..." << std::endl;
// Open the existing grid
grid_obsB = new appl::grid(grid_filenameB_in);
grid_obsB->optimise(true);
}
else {
std::cout << "Booking grid from scratch with name " << grid_filenameB_in << " ..." << std::endl;
// Create a grid
int Nobs = appl_common_histokin_.obsB_nbins;
double obsmin = appl_common_histokin_.obsB_min;
double obsmax = appl_common_histokin_.obsB_max;
grid_obsB = new appl::grid(NQ2, Q2min, Q2max, Q2order,
Nx, xmin, xmax, xorder,
Nobs, obsmin, obsmax,
genpdf, leading_order, nloops);
// Use the reweighting function
grid_obsB->reweight(true);
// The grid is an aMC@NLO type
grid_obsB->amcatnlo();
// Optimise the grid
//grid_obsB->optimise(true);
}
//std::cout << "C++, ... exiting init()" << std::endl;
}
//
// Fill the APPLgrid grid with the weights
//
void fill() {
//std::cout << "C++, entering fill() ..." << std::endl;
// Check that itype ranges from 1 to 3.
int itype = appl_common_histokin_.itype_histo;
if(itype != 1 && itype != 2 && itype != 3) {
std::cout << "ERROR: Invalid value of itype = " << itype << std::endl;
exit(-10);
}
// aMC@NLO weights. Four grids, ordered as {W0,WR,WF,WB}.
double* W0 = appl_common_weights_.W0;
double* WR = appl_common_weights_.WR;
double* WF = appl_common_weights_.WF;
double* WB = appl_common_weights_.WB;
int ilumi;
int nlumi = appl_common_lumi_.nlumi;
double ttol = 1e-100;
double x1,x2;
double scale2;
double obsA = appl_common_histokin_.obsA_histo;
double obsB = appl_common_histokin_.obsB_histo;
// Weight vector whose size is the total number of subprocesses
std::vector<double> weight(nlumi,0);
// (n+1)-body contribution (corresponding to xsec11 in aMC@NLO)
// It uses only Events (k=0) and the W0 weight.
if(itype == 1) {
// Get Bjorken x's
x1 = appl_common_weights_.x1[0];
x2 = appl_common_weights_.x2[0];
// Energy scale
scale2 = appl_common_weights_.muF2[0];
// Relevant parton luminosity combinations (-1 offset in c++)
ilumi = appl_common_weights_.flavmap[0] - 1;
if(x1 < 0 || x1 > 1 || x2 < 0 || x2 > 1) {
std::cout << "ERROR: Invalid value of x1 and/or x2 = " << x1 << " " << x2 << std::endl;
exit(-10);
}
// Fill grid only if x1 and x1 are non-zero
if(x1 == 0 && x2 == 0) return;
// Fill only if W0 is non zero
if(fabs(W0[0]) < ttol) return;
// Fill the grid with the values of the observables
// W0
weight.at(ilumi) = W0[0];
grid_obsA->fill_grid(x1,x2,scale2,obsA,&weight[0],0);
grid_obsB->fill_grid(x1,x2,scale2,obsB,&weight[0],0);
weight.at(ilumi) = 0;
}
// n-body contribution without Born (corresponding to xsec12 in aMC@NLO)
// It uses all the CounterEvents (k=1,2,3) and the weights W0, WR and WF.
else if(itype == 2) {
for(int k=1; k<=3; k++) {
x1 = appl_common_weights_.x1[k];
x2 = appl_common_weights_.x2[k];
scale2 = appl_common_weights_.muF2[k];
ilumi = appl_common_weights_.flavmap[k] - 1;
if(x1 < 0 || x1 > 1 || x2 < 0 || x2 > 1) {
std::cout << "ERROR: Invalid value of x1 and/or x2 = " << x1 << " " << x2 << std::endl;
exit(-10);
}
if(x1 == 0 && x2 == 0) continue;
if(fabs(W0[k]) < ttol && fabs(WR[k]) < ttol && fabs(WF[k]) < ttol) continue;
// W0
weight.at(ilumi) = W0[k];
grid_obsA->fill_grid(x1,x2,scale2,obsA,&weight[0],0);
grid_obsB->fill_grid(x1,x2,scale2,obsB,&weight[0],0);
weight.at(ilumi) = 0;
// WR
weight.at(ilumi) = WR[k];
grid_obsA->fill_grid(x1,x2,scale2,obsA,&weight[0],1);
grid_obsB->fill_grid(x1,x2,scale2,obsB,&weight[0],1);
weight.at(ilumi) = 0;
// WF
weight.at(ilumi) = WF[k];
grid_obsA->fill_grid(x1,x2,scale2,obsA,&weight[0],2);
grid_obsB->fill_grid(x1,x2,scale2,obsB,&weight[0],2);
weight.at(ilumi) = 0;
}
}
// Born (n-body) contribution (corresponding to xsec20 in aMC@NLO)
// It uses only the soft kinematics (k=1) and the weight WB.
else if(itype == 3) {
x1 = appl_common_weights_.x1[1];
x2 = appl_common_weights_.x2[1];
scale2 = appl_common_weights_.muF2[1];
ilumi = appl_common_weights_.flavmap[1] - 1;
if(x1 < 0 || x1 > 1 || x2 < 0 || x2 > 1) {
std::cout << "ERROR: Invalid value of x1 and/or x2 = " << x1 << " " << x2 << std::endl;
exit(-10);
}
if(x1 == 0 && x2 == 0) return;
if(fabs(WB[1]) < ttol) return;
// WB
weight.at(ilumi) = WB[1];
grid_obsA->fill_grid(x1,x2,scale2,obsA,&weight[0],3);
grid_obsB->fill_grid(x1,x2,scale2,obsB,&weight[0],3);
weight.at(ilumi) = 0;
}
//std::cout << "C++, ... exiting fill()" << std::endl;
}
//
// Fill the APPLgrid reference histogram.
//
void fill_ref() {
//std::cout << "C++, entering fill_ref() ..." << std::endl;
// Event weights to fill the histograms
double www = appl_common_histokin_.www_histo;
// Physical observables
double obsA = appl_common_histokin_.obsA_histo;
double obsB = appl_common_histokin_.obsB_histo;
grid_obsA->getReference()->Fill(obsA,www);
grid_obsB->getReference()->Fill(obsB,www);
//std::cout << "C++, ... exiting fill_ref()" << std::endl;
}
//
// Final normalization of the APPLgrid reference histogram.
//
void fill_ref_out() {
//std::cout << "C++, entering fill_ref_out() ..." << std::endl;
// Normalization factor
double norm = appl_common_histokin_.norm_histo;
// Check normalization value
if(norm <= 0.0 || norm > 1e50) {
std::cout << "ERROR: Invalid value for histogram normalization = " << norm << std::endl;
exit(-10);
}
// APPLgrid does not normalizes the reference histograms to the bin width
// just as in aMC@NLO, so no extra scaling by bin width needed.
grid_obsA->getReference()->Scale(norm);
grid_obsB->getReference()->Scale(norm);
//std::cout << "C++, ... exiting fill_ref_out()" << std::endl;
}
//
// Write out the grids
//
void term() {
//std::cout << "C++, entering term() ..." << std::endl;
// Normalization factor
double norm = appl_common_histokin_.norm_histo;
double n_runs = 1 / norm;
int n_runs_int = int(round(n_runs));
std::cout << "number of runs = " << n_runs_int << " " << n_runs << " " << norm << std::endl;
// Write out all the grids that have been filled.
// Save only if the pointer has been initialized.
if(grid_obsA) {
grid_obsA->run() = n_runs_int;
//grid_obsA->optimise(true);
// Write the grid to file
grid_obsA->Write(grid_filenameA_out);
}
if(grid_obsB) {
grid_obsB->run() = n_runs_int;
//grid_obsB->optimise(true);
// Write the grid to file
grid_obsB->Write(grid_filenameB_out);
}
//std::cout << "C++, ... exiting term()" << std::endl;
}
//
// Test that the original aMC@NLO event weight is properly reproduced
// by the APPLgrid information and LHAPDF.
//
void reco() {
//std::cout << "C++, entering reco() ..." << std::endl;
// Initialization
double xsec11 = 0;
double xsec12 = 0;
double xsec20 = 0;
// Set the Born alphas power
int const bpow = appl_common_fixed_.bpower;
// Bjorken x's
double x1bj = appl_common_weights_.x1[0];
double x2bj = appl_common_weights_.x2[0];
// Renormalization scale
double muR = std::sqrt(appl_common_weights_.muR2[0]);
// Factorization scale
double muF = std::sqrt(appl_common_weights_.muF2[0]);
// Get value of the strong coupling
double g;
if(muR != 0) g = std::sqrt(4 * M_PI * LHAPDF::alphasPDF(muR));
else g = 0;
// Compute the partonic luminosities
double pdflumi;
if(muF != 0 && x1bj != 0 && x2bj != 0) pdflumi = lhapdf_lumi(x1bj,x2bj,muF,0);
else pdflumi = 0;
// Compute the (n+1)-body NLO contribution to the xsec (xsec11)
xsec11 = pdflumi * appl_common_weights_.W0[0] * std::pow(g,2*bpow+2);
// Counter Events.
// Soft (k=1), Collinear (k=2) and Soft-Collinear (k=3)
for(int k=1; k<4; k++) {
double QES = std::sqrt(appl_common_weights_.muQES2[k]);
double muR = std::sqrt(appl_common_weights_.muR2[k]);
double muF = std::sqrt(appl_common_weights_.muF2[k]);
double xlgmuf = 0;
double xlgmur = 0;
if(QES != 0 && muR != 0 && muF != 0) {
xlgmuf = 2 * std::log(muF / QES);
xlgmur = 2 * std::log(muR / QES);
}
-
+ /*
// Check that (for the moment) the scale variation logarithms are
// always set to zero.
// Scale variation is included a posteriori using APPLgrid.
double const eps_scale = 1e-10;
if(fabs(xlgmuf) > eps_scale || fabs(xlgmur) > eps_scale) {
std::cout << "ERROR: scale variations logarithms non-zero" << std::endl;
std::cout << "xlgmuf = " << xlgmuf << std::endl;
std::cout << "xlgmur = " << xlgmur << std::endl;
exit(-10);
}
-
+ */
// Get values of Bjorken's x for the counterevent k
double x1bj = appl_common_weights_.x1[k];
double x2bj = appl_common_weights_.x2[k];
// Get value of the strong coupling
double g;
if(muR != 0) g = std::sqrt(4 * M_PI * LHAPDF::alphasPDF(muR));
else g = 0;
// Compute the partonic luminosities
double pdflumi;
if(muF != 0 && x1bj != 0 && x2bj != 0) pdflumi = lhapdf_lumi(x1bj,x2bj,muF,k);
else pdflumi = 0;
// Compute the n-body NLO contribution to the xsec (xsec12)
xsec12 += pdflumi * ( appl_common_weights_.W0[k] + appl_common_weights_.WF[k] * xlgmuf + appl_common_weights_.WR[k] * xlgmur ) * std::pow(g,2*bpow+2);
// Compute the Born contribution to the xsec (xsec20)
if(k == 1) {
if(bpow > 0) xsec20 += pdflumi * appl_common_weights_.WB[k] * std::pow(g,2*bpow);
else xsec20 += pdflumi * appl_common_weights_.WB[k];
}
}
// Sum up all the contributions
double xsec = xsec20 + xsec11 + xsec12;
// Divide by the vegas weight to compare with the physical cross section
xsec /= appl_common_reco_.vegaswgt;
// Compare with the reference weight
double xsec_ref = appl_common_reco_.event_weight;
double diff;
if(xsec_ref == 0) {
if(xsec == 0) diff = 0;
else {
std::cout << "ERROR: Failed reconstruction of original event weight" << std::endl;
std::cout << "xsec = " << xsec << std::endl;
std::cout << "xsec_ref = " << xsec_ref << std::endl;
exit(-10);
}
}
else diff = fabs( ( xsec - xsec_ref ) / xsec_ref );
// Check the accuracy
if(diff > 1e-5) {
std::cout << "ERROR: Failed reconstruction of original event weight" << std::endl;
std::cout << "xsec = " << xsec << std::endl;
std::cout << "xsec_ref = " << xsec_ref << std::endl;
exit(-10);
}
//std::cout << "C++, ... exiting reco()" << std::endl;
return;
}
//
// PDF luminosity based on LHAPDF and the information on the appl_common.inc block
//
double lhapdf_lumi(double x1bj, double x2bj, double qpdf, int k) {
// Conversion factor from natural units to pb
double const conv = 389379660;
// initialization
double pdflumi = 0;
// Check
if(k<0 || k>3) {
std::cout << "ERROR: Invalid value of k = " << k << std::endl;
exit(-10);
}
// Return zero if the Bjorken's x is equal to zero
if(x1bj == 0 || x2bj == 0 ) return 0;
// Select luminosity
int ilumi = appl_common_weights_.flavmap[k] - 1;
// Check physical range
if(ilumi < 0 || ilumi >= appl_common_lumi_.nlumi) {
std::cout << "ERROR: Invalid value of flavor map" << std::endl;
std::cout << "flavmap[k], nlumi = " << ilumi << " " << appl_common_lumi_.nlumi << std::endl;
exit(-10);
}
// Call the PDFs
for(int iproc=0; iproc<appl_common_lumi_.nproc[ilumi]; iproc++) {
// Get PDF flavors ID
int flav1 = appl_common_lumi_.lumimap[ilumi][iproc][0];
int flav2 = appl_common_lumi_.lumimap[ilumi][iproc][1];
// Check LHAPDF conventions are satisfied
// Assume always working in nf=5 scheme at most
if(flav1 < (-5) || flav2 < (-5) || flav2 > (+5) || flav2 > (+5)) {
std::cout << "ERROR: Invalid value of PDF flavor indices" << std::endl;
std::cout << "flav1, flav2 = " << flav1 << " " << flav2 << std::endl;
exit(-10);
}
// Compute the luminosity
double xpdflh_1 = LHAPDF::xfx(x1bj,qpdf,flav1);
double xpdflh_2 = LHAPDF::xfx(x2bj,qpdf,flav2);
pdflumi += ( xpdflh_1 / x1bj ) * ( xpdflh_2 / x2bj );
}
pdflumi *= conv;
return pdflumi;
}
} // end namespace amcbridge

File Metadata

Mime Type
text/x-diff
Expires
Tue, Sep 30, 6:15 AM (2 h, 38 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6566505
Default Alt Text
(18 KB)

Event Timeline

HEPForge.org · .