Page MenuHomeHEPForge
Files F7879077

No OneTemporary
Actions

View Options

Index: nnlo-bridge/trunk/src/fastnlo_utils.cxx
===================================================================
--- nnlo-bridge/trunk/src/fastnlo_utils.cxx (revision 1626)
+++ nnlo-bridge/trunk/src/fastnlo_utils.cxx (revision 1627)
@@ -1,939 +1,953 @@
//
// @file fastnlo_utils.cxx
//
//
// @author D. Britzger, K. Rabbertz
//
//
//
#include "amconfig.h"
#ifdef HAVE_FASTNLO
#include <iostream>
#include <fstream>
#include <algorithm> //c++98
#include <utility> //c++11
#include <sys/stat.h> //c++98
/// fastnlo headers
#include "fastnlo_utils.h"
//#include "fastnlotk/fastNLOCreate.h"
#include "nnlo_common.h"
#include "nnlo_utils.h"
//-----------------------------------------------------------------------------
//
// --- --- --- fnloUtils --- --- --
//
//-----------------------------------------------------------------------------
extern "C" void evolvepdf_(const double& x, const double& Q, double* xf);
extern "C" double alphaspdf_(const double& Q);
//
// --- fastNLO specific functions
//
// _____________________________________________________________________ //
fastNLO::GeneratorConstants fnloUtils::GetNNLOJET_GenConsts() {
//!< Get NNLOJET specific generator constants for this interface
fastNLO::GeneratorConstants gc;
gc.Name = "NNLOJET"; //!< Name and version of generator
gc.References.push_back("JHEP 1607 (2016) 133 [arXiv:1605.04295]"); //!< References for generator
gc.References.push_back("Phys. Rev. Lett. 117, 042001 (2016) [arXiv:1606.03991]"); //!< References for generator
gc.UnitsOfCoefficients = 15;
return gc;
}
// _____________________________________________________________________ //
//fastNLO::ProcessConstants fnloUtils::GetNNLOJET_ProcConstsPP(int LOord, const std::vector<std::vector<int> >& partLiCos) {
fastNLO::ProcessConstants fnloUtils::GetNNLOJET_ProcConstsPP() {
//!< Get default process constants for pp processes
// --- external input from NNLOJET
std::vector<std::vector<int> > partLiCos = fnloUtils::sp.GetPartonCombinations();
if ( partLiCos.empty() ) {
std::cerr<<"[fnloUtils::GetNNLOJET_ProcConstsPP] Error. No parton combinations found. "
<< "Please read config file into Subprocess-class: fnloUtils::sp .Exiting."<<std::endl;
exit(5);
}
fastNLO::ProcessConstants pc;
pc.LeadingOrder = nnlo::GetPowLO(); // Order in alpha_s of leading order process
pc.NPDF = 2; // No. of PDFs involved
pc.NPDFDim = 2; // 2: full-matrix (then we don't need 'asymmetric processes
int nprocs = partLiCos.size();
pc.IPDFdef1 = 3; // Flag 1 to define PDF linear combinations
pc.IPDFdef2 = 0; // Flag 2 to define PDF linear combinations
pc.IPDFdef3LO = nprocs; // No. of subprocesses
pc.IPDFdef3NLO = nprocs;
pc.IPDFdef3NNLO = nprocs;
pc.NSubProcessesLO = nprocs; //No. of LO subprocesses
pc.NSubProcessesNLO = nprocs;
pc.NSubProcessesNNLO= nprocs;
// parton combinations in two different formats
pc.PDFCoeffLO = VectorToPairs(partLiCos);//!< PDF Linear combinations
pc.PDFCoeffNLO = pc.PDFCoeffLO;
pc.PDFCoeffNNLO = pc.PDFCoeffLO;
pc.PDFLiCoInLO = partLiCos; //!< PDF Linear combinations
pc.PDFLiCoInNLO = pc.PDFLiCoInLO;
pc.PDFLiCoInNNLO = pc.PDFLiCoInLO;
//AsymmetricProcesses = ; // empty
pc.Name = "Name of process not set and only specified in NNLOJET runcard.";
pc.References.push_back("Please see NNLOJET manual for more reference."); // references for processes
return pc;
}
// _____________________________________________________________________ //
fastNLO::ProcessConstants fnloUtils::GetNNLOJET_ProcConstsDIS() {
//!< Get default process constants for pp processes
// --- external input from NNLOJET
//int LOord = inppar_.njets - 1; // DIS: minus 1
std::cout<< "\n\n LOord = "<<nnlo::GetPowLO()<<std::endl<<std::endl;
std::vector<std::vector<int> > partLiCos = fnloUtils::sp.GetPartonCombinations();
if ( partLiCos.empty() ) {
std::cerr<<"[fnloUtils::GetNNLOJET_ProcConstsDIS] Error. No parton combinations found. "
<< "Please read config file into Subprocess-class: fnloUtils::sp .Exiting."<<std::endl;
exit(5);
}
fastNLO::ProcessConstants pc;
pc.LeadingOrder = nnlo::GetPowLO(); // Order in alpha_s of leading order process
pc.NPDF = 1; // No. of PDFs involved
pc.NPDFDim = 0; // 2: full-matrix (then we don't need 'asymmetric processes
int nprocs = partLiCos.size();
pc.IPDFdef1 = 2; // Flag 1 to define PDF linear combinations
pc.IPDFdef2 = 0; // Flag 2 to define PDF linear combinations
pc.IPDFdef3LO = nprocs; // No. of subprocesses
pc.IPDFdef3NLO = nprocs;
pc.IPDFdef3NNLO = nprocs;
pc.NSubProcessesLO = nprocs; //No. of LO subprocesses
pc.NSubProcessesNLO = nprocs;
pc.NSubProcessesNNLO= nprocs;
// parton combinations in two different formats
pc.PDFCoeffLO = VectorToPairs(partLiCos);//!< PDF Linear combinations
pc.PDFCoeffNLO = pc.PDFCoeffLO;
pc.PDFCoeffNNLO = pc.PDFCoeffLO;
pc.PDFLiCoInLO = partLiCos; //!< PDF Linear combinations
pc.PDFLiCoInNLO = pc.PDFLiCoInLO;
pc.PDFLiCoInNNLO = pc.PDFLiCoInLO;
//AsymmetricProcesses = ; // empty
pc.Name = "Name of process not set and only specified in NNLOJET runcard.";
pc.References.push_back("Please see NNLOJET manual for more reference."); // references for processes
return pc;
}
// _____________________________________________________________________ //
fastNLO::ScenarioConstants fnloUtils::GetNNLOJET_ScenConsts() {
//!< Get default process constants for pp processes
fastNLO::ScenarioConstants sc;
sc.ScenarioName = "fnlXYZ"; // no white spaces here!
sc.ScenarioDescription.push_back("test for nnlojet"); //< Description of the scenario
sc.PublicationUnits = 12; //< we are usually working in units of [pb]
// Dimensionality of binning
// also decides if SingleDifferentialBinning or DoubleDifferentialBinning is used
// 1: single-differential,
// 2: double-differential;
sc.DifferentialDimension = 1;
// Labels (symbol and unit) for the measurement dimensions (from outer to inner "loop"),
// e.g. "|y|" and "p_T [GeV]".
// This may also help to define the observables to be calculated in an automatized way!
sc.DimensionLabels.push_back("varName");
// Specify:
// 0 : the cross section is NOT differential,
// i.e. there are two bin borders (but NO division (normalization) by bin width);
// 1 : the cross section is point-wise differential, i.e. only one point is given;
// 2 : the cross section is bin-wise differential, i.e. there are two bin borders and division by bin width
sc.DimensionIsDifferential.push_back(2);
sc.CalculateBinSize = true; //< Calculate bin width from lower and upper bin boundaries
sc.BinSizeFactor = 1; //< Possibility to provide additional normalization factor, e.g. of 2 for bins in |y|
//sc.BinSize = ; //< If 'CalculateBinSize' is 'false' provide table with bin widths for normalization.
sc.ScaleDescriptionScale1 = "scale1"; //< "<pT_1,2>_[GeV]" # This defines the scale to be used (Note: The 1st scale should always be in units of [GeV]!)
sc.ScaleDescriptionScale2 = "scale2"; //< "pT_max_[GeV]" # Specify 2nd scale name and unit (ONLY for flexible-scale tables)
//DB: binning remains empty here and has to be set later.
// Observable binning Use either 'SingleDifferentialBinning' or
// 'DoubleDifferentialBinning' or 'TripleDifferentialBinning' in
// accordance with 'DifferentialDimension' above
// sc.SingleDifferentialBinning;
// sc.DoubleDifferentialBinning; //< Observable binning
// sc.TripleDifferentialBinning; //< Observable binning
sc.CenterOfMassEnergy = inputpar_.roots; //todo //< Center-of-mass energy in GeV. LHC Next Run II: 13000
sc.PDF1 = 2212; //< PDF of 1st hadron (following PDG convention: proton 2212).
sc.PDF2 = 2212; //< PDF of 2nd hadron (following PDG convention: proton 2212).
sc.OutputFilename = "test.tab";//< Filename of fastNLO output table
sc.OutputPrecision = 8; //< Number of decimal digits to store in output table (def.=8).
// Create table fully flexible in mu_f
// larger size, and requires scale independent weights during creation
// or table with fixed number of mu_f scale factors, def.=false.
sc.FlexibleScaleTable = false;
// Factorization scale variations
// only needed for fixed-scale tables,
// List of scale factors must include factor '1',
// Scale factors will be ordered according to fastNLO convention:
// (1, min, ... , max). Defaults: {0.5, 1, 2}
sc.ScaleVariationFactors.push_back(1);
sc.ReadBinningFromSteering = true; // Specify if binning is read from fScenConst or from warmup
sc.ApplyPDFReweighting = true; // Apply reweighting of pdfs for an optimized interpolation, def.=true.
// For warmup-run! Set limits for scale nodes to bin borders, if possible
sc.CheckScaleLimitsAgainstBins = true;
// Choose fastNLO interpolation kernels and distance measures
sc.X_Kernel = "Catmull";
sc.X_DistanceMeasure = "sqrtlog10";
sc.X_NNodes = 18;
sc.X_NoOfNodesPerMagnitude = false;
sc.Mu1_Kernel = "Lagrange";
sc.Mu1_DistanceMeasure = "log10";//"loglog025";
sc.Mu1_NNodes = 12;
sc.Mu2_Kernel = "OneNode"; //"Lagrange";//< Lagrange # Scale2 not used for fixed-scale tables
sc.Mu2_DistanceMeasure = "linear";//"loglog025";//< "loglog025"
sc.Mu2_NNodes = 1;//< 6
return sc;
}
// _____________________________________________________________________ //
void fnloUtils::SetNNLOJETDefaultBinning(fastNLO::ScenarioConstants& sc, const int& nbins, const double& lo, const double& hi) {
//!< set binning as used by NNLOJET
sc.SingleDifferentialBinning.resize(nbins+1);
for (int i=0; i<=nbins; i++) sc.SingleDifferentialBinning[i] = (&lo)[i];
std::cout<<"DefaultBinning: "<<lo;
for (int i=1; i<=nbins; i++) std::cout<<"\t"<<sc.SingleDifferentialBinning[i];
std::cout<<std::endl;
return ;
// double binwidth = (hi-lo)/nbins;
// sc.SingleDifferentialBinning[0] = lo;
// std::cout<<"DefaultBinning: "<<lo;
// for (int i=1; i<=nbins; i++) {
// sc.SingleDifferentialBinning[i] = sc.SingleDifferentialBinning[i-1]+binwidth;
// std::cout<<"\t"<<sc.SingleDifferentialBinning[i];
// }
// std::cout<<std::endl;
}
// // _____________________________________________________________________ //
// std::string fnloUtils::GetProcessName() {
// //!< get name of the process from NNLOJET
// std::string process(process_.sproc);
// process.erase(std::remove(process.begin(), process.end(), ' '),process.end());
// return process;
// }
// _____________________________________________________________________ //
std::string fnloUtils::GetContribName(int nloops, int ne) {
//!<
//!< Get name of contributions
//!<
/*
* LO = #looops=0 #extra-emissions=0
* V = #looops=1 #extra-emissions=0
* R = #looops=0 #extra-emissions=1
* VV = #looops=2 #extra-emissions=0
* RV = #looops=1 #extra-emissions=1
* RR = #looops=0 #extra-emissions=2
*/
if ( nloops==-1 && ne==-1 ) {
std::cout<<"Works only for having a CURRENT event !!!"<<std::endl;
nloops = ndimcurrent_.nv;
static const int& nord = nnlo::GetOrder();
// int njet = inppar_.njets;
// if ( nnlo::IsDIS() ) njet -= 1;
ne = nord - nloops;
}
if ( nloops == 0 && ne == 0) return "LO";
else if ( nloops == 1 && ne == 0) return "V";
else if ( nloops == 0 && ne == 1) return "R";
else if ( nloops == 2 && ne == 0) return "VV";
else if ( nloops == 1 && ne == 1) return "RV";
else if ( nloops == 0 && ne == 2) return "RR";
else return "Contrib Name NOT IDENTIFIED: "+std::to_string(nnlo::GetOrder());
}
// // _____________________________________________________________________ //
// int fnloUtils::GetOrder() {
// static const int& nord = ndimcurrent_.norder ;
// static const int& nv = ndimcurrent_.nv ;
// return nord + nv;
// }
// _____________________________________________________________________ //
std::string fnloUtils::GetOrderName( std::set<std::pair<int,int> >* nlne ) {
//!< get order of calculation as 'string'
if ( nlne == NULL) {
int nord = nnlo::GetOrder();
//static const int& nord = ndimcurrent_.norder;
if ( nord == 0 ) return "LO";
else if ( nord == 1 ) return "NLO";
else if ( nord == 2 ) return "NNLO";
}
- else if ( fIDs.size()==1 ) {
+ else if ( fnloUtils::fIDs.size()==1 ) {
// only one channel was calculated
// take: 'chZYZ'
char buf[3];
- sprintf(buf, "%03d", *fIDs.begin());
+ sprintf(buf, "%03d", *fnloUtils::fIDs.begin());
return "ch"+std::string(buf);
}
else {
std::set<std::string> names;
- for ( auto i : *nlne ) names.insert(GetContribName(i.first,i.second));
+ for ( auto i : *nlne ) names.insert(fnloUtils::GetContribName(i.first,i.second));
if ( names.count("R") && names.count("V") ) {
names.erase("R");
names.erase("V");
names.insert("NLO");
}
if ( names.count("RR") && names.count("VV") && names.count("RV") ) {
names.erase("RR");
names.erase("RV");
names.erase("VV");
names.insert("NNLO");
}
std::string ret;
for ( auto i : names )
if ( ret.empty() ) ret = i;
else ret += "+"+i;
return ret;
}
return ""; //?
}
// _____________________________________________________________________ //
std::string fnloUtils::GetConfigDir() {
//!< Get directory where config files are stored
// CONFIGDIR is set by NNLOJET
//
std::string configdir(getenv( "CONFIGDIR" ));
if ( configdir == "" ) {
std::cout << "[fnloUtils::GetConfigDir()] env $CONFIGDIR could not be found, "
<< "using $PWD/../driver/process instead" << std::endl;
configdir = std::string(getenv("PWD"))+"/../driver/process";
//std::exit(3);
}
// --- test directory
struct stat myStat;
if ( !((stat(configdir.c_str(), &myStat) == 0) && (((myStat.st_mode) & S_IFMT) == S_IFDIR)) ) {
std::cout << "[fnloUtils::GetConfigDir()] Warning. Could not find a config directory" << std::endl;
// another guess
configdir = std::string(getenv("PWD"))+"/driver/process";
if ( !((stat(configdir.c_str(), &myStat) == 0) && (((myStat.st_mode) & S_IFMT) == S_IFDIR)) ) {
std::cerr << "[fnloUtils::GetConfigDir()] Error. Could really not find a config directory" << std::endl;
exit(2);
}
}
std::cout<<"[fnloUtils::GetConfigDir()] Using config directory: "<<configdir<<std::endl;
return configdir;
}
// _____________________________________________________________________ //
std::string fnloUtils::GetConfigFile() {
//!< get name of the process
// "ZJ-LO.config"
std::string fname = nnlo::GetProcessName();
fname += "-" + fnloUtils::GetOrderName() +".config";
fname = fnloUtils::GetConfigDir()+"/"+nnlo::GetProcessName()+"/"+fname;
std::cout<<"[fnloUtils::GetConfigFile()] Using config file: "<<fname<<std::endl;
return fname;
}
// _____________________________________________________________________ //
std::vector<std::string> fnloUtils::GetConfigFiles() {
//!< get name of the process
// "ZJ-LO.config"
std::vector<std::string> fnames;// = fnloUtils::GetProcessName();
std::string pref = fnloUtils::GetConfigDir()+"/"+nnlo::GetProcessName()+"/"+ nnlo::GetProcessName();
if ( nnlo::GetOrder() == 0 ) {
fnames.push_back( pref+"-LO.config");
}
else if ( nnlo::GetOrder() == 1 ){
fnames.push_back( pref+"-R.config");
fnames.push_back( pref+"-V.config");
}
else if ( nnlo::GetOrder() == 2 ){
fnames.push_back( pref+"-RR.config");
fnames.push_back( pref+"-RV.config");
fnames.push_back( pref+"-VV.config");
}
std::cout<<"[fnloUtils::GetConfigFile()] Found order "<<nnlo::GetOrder()<<". Using config files: "<<std::endl;
for ( auto f : fnames ) std::cout<<"\t"<<f<<std::endl;
return fnames;
}
// _____________________________________________________________________ //
void fnloUtils::RememberContribution() {
//!< fill the variable fNlNe
//!< and remember the calculated contributions
static const int& iproc = currentprocess_.iproc;
int nloops = channel_.iloops[iproc-1];
int nreal = ndimcurrent_.norder;
// static const int& nloops = ndimcurrent_.nv;
// int nord = nnlo::GetOrder();//ndimcurrent_.norder;
// // int njet = inppar_.njets;
// // if ( nnlo::IsDIS() ) njet -= 1;
// // int ne = nord - nloops;
// //std::cout<<"nl="<<nloops<<"\tne="<<ne<<"\t\tnord="<<nord<<"\tIsDIS="<<nnlo::IsDIS()<<"\tnjet="<<njet<<std::endl;
fnloUtils::fNlNe.insert({ nloops,nreal });
fnloUtils::fIDs.insert(iproc);
}
// _____________________________________________________________________ //
double fnloUtils::CalculateReferenceWeight(double wgt, double x1, double x2) {
//!< calculate the reference weight
//!< and compare it later with the ref-weight from NNLOJET
//!<
//!< Set value: fnloUtils::_myfillwgt
//!< Need access to evolvepdf_()
//!< needs either strong_.as or alphaspdf_()
//!<
double PDF = 0;
double muf2 = muf2_hook(1);
const std::vector<std::vector<int > >& pc = fnloUtils::sp.GetInitialValues();
if ( nnlo::IsDIS() ) {
double xfx[13] = {0};
evolvepdf_(x1,sqrt(muf2),xfx);
for ( unsigned int ip = 0 ; ip < pc[currentprocess_.iproc].size() ; ip+=2 ) {
PDF += xfx[pc[currentprocess_.iproc][ip]+6];
}
}
else {
double xfx1[13] = {0};
double xfx2[13] = {0};
evolvepdf_(x1,sqrt(muf2),xfx1);
evolvepdf_(x2,sqrt(muf2),xfx2);
for ( unsigned int ip = 0 ; ip < pc[currentprocess_.iproc].size() ; ip+=2 ) {
PDF += xfx1[pc[currentprocess_.iproc][ip]+6]*xfx2[pc[currentprocess_.iproc][ip+1]+6];
}
}
static const double& as = strong_.as; //alphaspdf_(sqrt(mur2));
fnloUtils::_myfillwgt = wgt * PDF * pow( as/(2*M_PI) , nnlo::GetNPow() ) ;
return fnloUtils::_myfillwgt;
}
// _____________________________________________________________________ //
void fnloUtils::InitFastNLO( const int& id ) {
using namespace std;
// static const int& iproc = currentprocess_.iproc;
// int nloops = channel_.iloops[iproc-1]; /// loop order
// int ipars = channel_.ipars[iproc-1]; /// number of additional particles + 1
// std::cout<<" nloops "<<nloops<<std::endl;
// std::cout<<" ipars "<<ipars<<std::endl;
// std::cout<<" ndim.norder "<<ndimcurrent_.norder<<std::endl;
// std::cout<<" inppar.njets "<<inppar_.njets<<std::endl;
// std::cout<<" order_.ieorder "<<order_.ieorder<<std::endl;
// std::cout<<" GetNPos "<<nnlo::GetNPow()<<std::endl;
// std::cout<<" GetOrder() "<<nnlo::GetOrder()<<std::endl;
// std::cout<<std::endl;
// --- welcome
std::cout << std::endl;
std::cout << "[nnlo::init_fastnlo()] fastNLO initialisation" << std::endl;
if ( fnloUtils::ftable[id] != NULL ) return;
std::cout << "[nnlo::init_fastnlo()] NNLOJET provides process: " << nnlo::GetProcessName() << std::endl;
const std::string& gridname = fnloUtils::fInits[id].gridname;
const int& nbins = fnloUtils::fInits[id].nbins;
const double& lo = *(fnloUtils::fInits[id].lo);
const double& hi = *(fnloUtils::fInits[id].hi);
cout<<"fnloUtils::GetProcess(): "<<nnlo::GetProcessName()<<endl;
// --- read config file with parton combinations
if ( fnloUtils::sp.GetConfigFile().empty() )
fnloUtils::sp.ReadConfigFiles(fnloUtils::GetConfigFiles());
// --- get NNLOJET specific generator constants
fastNLO::GeneratorConstants gc = fnloUtils::GetNNLOJET_GenConsts();
// --- build process specific 'processConstants'
fastNLO::ProcessConstants pc = nnlo::IsDIS() ? //fnloUtils::GetProcess() == "DIS" ) ?
fnloUtils::GetNNLOJET_ProcConstsDIS() :
fnloUtils::GetNNLOJET_ProcConstsPP() ;
// --- build scenario specific 'scenarioConstants'
fastNLO::ScenarioConstants sc = fnloUtils::GetNNLOJET_ScenConsts();
fnloUtils::SetNNLOJETDefaultBinning(sc,nbins,lo,hi);
// --- temporary special cases.
if ( nnlo::IsDIS() ){
sc.FlexibleScaleTable = true;
sc.CheckScaleLimitsAgainstBins = false;
}
// ---- initialize fastNLOCreate
std::cout<<"\n INPUT TABLENAME: "<<gridname<<std::endl;
std::string tablename = gridname.substr(0,gridname.find_last_of("."));
- //std::string runid = runid_.srunid;
+ std::string runid = runid_.srunid;
+ std::string jname = jobname_.sjobname;
+ jname.erase(std::remove(jname.begin(), jname.end(), ' '),jname.end());
std::string fname = outfile_.fname;
fname.erase(std::remove(fname.begin(), fname.end(), ' '),fname.end());
+
//tablename = nnlo::GetProcessName()+".EXM."+runid_.srunid+"."+tablename;
- tablename = fname+tablename;
- std::string wrmfile = tablename+".wrm";
+ //tablename = fname+tablename;
+ tablename = nnlo::GetProcessName()+"."+jname+"."+runid+"."+tablename;
+ //std::string wrmfile = tablename+".wrm";
+ std::string wrmfile = nnlo::GetProcessName()+"."+jname+"."+tablename+".wrm";
+
+
sc.OutputFilename = tablename;//+"."+fnloUtils::GetOrderName()+".tab";
//std::cout<<"wrmfile="<<wrmfile<<"\ttablename="<<tablename<<"\ttabname="<<sc.OutputFilename<<std::endl;
fnloUtils::ftablename[id] = tablename;
fnloUtils::ftable[id] = new fastNLOCreate(wrmfile,gc,pc,sc);
//cout<<"\n\n\n ORDER OF alpha_s(): "<<pc.LeadingOrder+nnlo::GetOrder()<<"\t\t ORDER OF CALCULATION: "<<nnlo::GetOrder()<<"\n\n"<<endl;
cout<<"\n\n\n ORDER OF alpha_s(): "<<nnlo::GetNPow()<<"\t\t ORDER OF CALCULATION: "<<nnlo::GetOrder()<<"\n\n"<<endl;
fnloUtils::ftable[id]->SetOrderOfAlphasOfCalculation(nnlo::GetNPow());
//fnloUtils::ftable[id]->fastNLOTable::SetNcontrib(1); // to be removed later!
std::cout << "[nnlo::init_fastnlo] fastNLO initialized. " << std::endl;
return ;
}
// _____________________________________________________________________ //
namespace fnloUtils {
using namespace std;
// _____________________________________________________________________ //
fnloSuprocesses::fnloSuprocesses(const std::string& config){
//!< constructor reading config file with subprocesses
ReadConfigFile(config);
}
// _____________________________________________________________________ //
void fnloSuprocesses::ReadConfigFile(const std::string& config) {
fConfigFile = config;
this->resize(1000);
fInitialValues.resize(1000);
ifstream cf(config);
if ( !cf.is_open() ){
std::cerr<<"[fnloSubprocesses::ReadConfigFile()] Error. Could not open file: "<<config<<std::endl;
exit(2);
}
int id, np, v1, v2;
int nl = 0;
cf>>id; // first line is '0'
while ( !cf.eof() ) {
cf>>id>>np;
if(cf.fail()) break;
nl++;
//cout<<"nl: "<<nl<<"\tid="<<id<<"\tnp="<<np<<endl;
for ( int i = 0 ; i<np ; i++ ) {
cf>>v1>>v2;
//cout<<"\t"<<v1<<"\t"<<v2<<endl;
fInitialValues[id].push_back(v1);
fInitialValues[id].push_back(v2);
}
}
std::cout<<"[fnloSuprocesses::ReadConfigFile()] found "<<nl<<" subprocess definitions."<<std::endl;
SetUniqueValues();
}
// _____________________________________________________________________ //
bool fnloSuprocesses::CompareVectors(const std::vector<int>& v1, const std::vector<int>& v2 ) const {
//!< compare to vectors if all elements are equal
if ( v1.size() != v2.size() ) return false;
for ( unsigned int i=0 ; i<v1.size() ; i++ ) {
if ( v1[i] != v2[i] ) return false;
}
return true;
}
// _____________________________________________________________________ //
void fnloSuprocesses::SetUniqueValues() {
//!< init member SetUniqueValues and 'mother' vector
for ( unsigned int i = 0 ; i<fInitialValues.size() ; i++ ) {
if ( !fInitialValues[i].empty() ){
// look if values are already present ?
unsigned int umax = fUniqueValues.size();
for ( unsigned int u = 0 ; u <= umax ; u++ ) {
if ( u==fUniqueValues.size() ) fUniqueValues.push_back(fInitialValues[i]);
else if ( CompareVectors(fInitialValues[i],fUniqueValues[u]) ) break;
}
}
}
// sort a bit ?!
// put g-g to 0 if available
vector<int> gg(2);//{0,0}; gg -> X
for ( unsigned int u = 0 ; u < fUniqueValues.size() ; u++ ) {
if ( CompareVectors(gg,fUniqueValues[u])) {
std::swap(fUniqueValues[u],fUniqueValues[0]);
}
}
gg[1]=99;//{0,99}; // DIS eg -> X
for ( unsigned int u = 0 ; u < fUniqueValues.size() ; u++ ) {
if ( CompareVectors(gg,fUniqueValues[u])) {
std::swap(fUniqueValues[u],fUniqueValues[0]);
}
}
// define our map (which is actually a vector)
for ( unsigned int i = 0 ; i<fInitialValues.size() ; i++ ) {
(*this)[i] = -1; // init values
if ( !fInitialValues[i].empty() ){
for ( unsigned int u = 0 ; u < fUniqueValues.size() ; u++ ) {
if ( CompareVectors(fInitialValues[i],fUniqueValues[u]) )
(*this)[i] = u; // define our final map
}
}
}
std::cout<<"[fnloSuprocesses::SetUniqueValues()] found "<<fUniqueValues.size()<<" unique subprocess definitions."<<std::endl;
for ( unsigned int i = 0 ; i<fInitialValues.size() ; i++ ) cout<<i<<"\t"<<fInitialValues[i] <<endl;
for ( unsigned int i = 0 ; i<fUniqueValues.size() ; i++ ) cout<<i<<"\t"<<fUniqueValues[i] <<endl<<endl;
cout<<"[fnloSuprocesses::SetUniqueValues()] Mapping: "<<(*this)<<endl<<endl;
}
// _____________________________________________________________________ //
std::vector<std::vector<int > > fnloSuprocesses::GetPartonCombinations() const {
auto out = fUniqueValues;
int n = 0;
for (auto& i : out ) i.emplace(i.begin(),n++);
//for ( auto& i : PartonCombinationsLO ) cout<<i<<endl;
return out;
}
} // end namespace
// _____________________________________________________________________ //
// Temporary pointer to APPLGRID lumi function for fastNLO
//const lumi_pdf *fnlopdf;
//-----------------------------------------------------------------------------
//
// --- --- --- fastNLO_nnlo --- --- --
//
//-----------------------------------------------------------------------------
// _____________________________________________________________________ //
void nnlo::init_fastnlo( const int& id, const std::string& gridname, const int& nbins, const double& lo, const double& hi ) {
//!<
//!< Remember init variables for later initialization
//!<
fnloUtils::fInits[id] = fnloUtils::fnloInit{int(id), std::string(gridname), int(nbins), &lo, &hi };
fnloUtils::ftable[id] = NULL;
}
// _____________________________________________________________________ //
void nnlo::fill_fastnlo( const int& id, const double& obs, const double* wt ) {
//!<
//!< Fill fastNLO table
//!<
say::SetGlobalVerbosity(say::WARNING);
// { // --- debugging
// double x1 = parfrac_.x1;
// double x2 = parfrac_.x2;
// double Q2 = mur2_hook(1);
// double xfx[13] = {0};
// evolvepdf_(x1,sqrt(Q2),xfx);
// double as = alphaspdf_(sqrt(Q2));
// std::cout<<"id="<<id<<"\tfnloProc="<<fnloUtils::sp[currentprocess_.iproc]
// <<"\tx1="<<x1
// <<"\tx2="<<x2
// <<"\tQ2="<<Q2
// <<"\tobs="<<obs
// <<"\tgluon-x1="<<xfx[6]
// <<"\tas="<<as
// <<"\twt="<<wt[0]<<std::endl;
// }
// --- init fastNLO
if ( fnloUtils::ftable[id]==NULL ) fnloUtils::InitFastNLO(id);
// --- keep track on what we have calculated
fnloUtils::RememberContribution();
/*
//static const int& ieorder = order_.ieorder;
static const int& nloops = ndimcurrent_.nv;
int nord = nnlo::GetOrder();//ndimcurrent_.norder;
int njet = inppar_.njets;
if ( nnlo::IsDIS() ) njet -= 1;
int ne = nord - nloops;
//std::cout<<"nl="<<nloops<<"\tne="<<ne<<"\t\tnord="<<nord<<"\tIsDIS="<<nnlo::IsDIS()<<"\tnjet="<<njet<<std::endl;
fnloUtils::fNlNe.insert({ nloops,ne });
*/
// --- calculate weights and x-values
double __w2 = wt[0];//
double __x1, __x2;
nnlo::GetFillWeights(__w2,__x1,__x2);
// --- set reference weight for cross check with NNLOJET
fnloUtils::CalculateReferenceWeight(__w2,__x1,__x2);
/*
double mur2 = mur2_hook(1);
double muf2 = muf2_hook(1);
const double& as = strong_.as; //alphaspdf_(sqrt(mur2));
double wgt = wt[0];
wgt /= parfrac_.x1;
if ( parfrac_.x2 != 0 )
wgt /= parfrac_.x2;
//wgt*=runinfo_.nshota;
wgt /= pow( as/(2*M_PI) , fnloUtils::ftable[id]->GetTheCoeffTable()->GetNpow() );// we have to divide by (alpha_s/2pi)^n
// do not use fnlo PDF
// {
// double x1 = parfrac_.x1;
// double xfx[13] = {0};
// evolvepdf_(x1,sqrt(mur2),xfx);
// wgt *= xfx[6];
// }
// do not use fnlo as
//wgt *= pow( as/(2*M_PI) , fnloUtils::ftable[id]->GetTheCoeffTable()->GetNpow() );// we have to divide by (alpha_s/2pi)^n
// !* remaining xreg[1|2] mapping
// ./core/sig.f:L578
if ( ndimcurrent_.nv!=0 ) {
wgt /= (1.-parfrac_.x1);
if ( nnlo::IsDIS() )
wgt /= (1.-parfrac_.x2);
}
//std::cout<<"wNew/wOld="<<__w2/wgt<<std::endl;
{
//std::cout<<"as="<<as<<"\tnpow="<<fnloUtils::ftable[id]->GetTheCoeffTable()->GetNpow()<<std::endl;
// double x1 = parfrac_.x1;
// double xfx[13] = {0};
// double xfxReg[13] = {0};
double PDF = 0;
//evolvepdf_(x1,sqrt(mur2),xfx);
//xval = x1_kin / xreg1_kin
//evolvepdf_(x1/xregions_.xreg1,sqrt(mur2),xfxReg);
if ( nnlo::IsDIS() ) {
double xfx[13] = {0};
evolvepdf_(__x1,sqrt(muf2),xfx);
const std::vector<std::vector<int > >& pc = fnloUtils::sp.GetInitialValues();
for ( unsigned int ip = 0 ; ip < pc[currentprocess_.iproc].size() ; ip+=2 ) {
//std::cout<<"\tip="<<ip<<"\tp="<<pc[currentprocess_.iproc][ip]<<"\tx="<<x1<<"\t1./x="<<1./x1<<"\txreg1="<<xregions_.xreg1<<"\t1/xreg="<<1./xregions_.xreg1<<"\tPDFreg/PDF="<<xfxReg[pc[currentprocess_.iproc][ip]+6]/xfx[pc[currentprocess_.iproc][ip]+6]<<std::endl;
PDF += xfx[pc[currentprocess_.iproc][ip]+6];
//PDF += xfx[pc[currentprocess_.iproc][ip]+6];
//PDF += xfxReg[pc[currentprocess_.iproc][ip]+6];
}
// ----------------
// int ix = gridix_.igx;
// double _x1 = x1;
// double _x2 = x2;
// if ( ix==2 ) _x2 = x2/xregions_.xreg2;
// else if ( ix==3 ) _x1 = x1/xregions_.xreg1;
// else if ( ix==4 ) {
// _x1 = x1/xregions_.xreg1;
// _x2 = x2/xregions_.xreg2;
// }
// //if ( jacobian_.jflag ) jac *= 1/((1-x1)*(1-x2)); /// *not* 1/((1-_x1)(1-_x2))
// if ( ix==2 ) jac *= 1./xregions_.xreg2;
// else if ( ix==3 ) jac *= 1./xregions_.xreg1;
// else if ( ix==4 ) jac *= 1./(xregions_.xreg1*xregions_.xreg2);
// // --- shortcuts to NNLOJET parameters
// // static const int& nv = ndimcurrent_.nv;
// // static const int& iproc = currentprocess_.iproc;
// static const int& ix = gridix_.igx;
// // --- variables for calculation
// int order = nnlo::GetOrder();//getorder( iproc );
// int as_order = inppar_.njets + order;
// double wtjac = colourflag_.colflag ? 1./xcoljac_.xcoljac : 1;
// std::cout<<"\twtjac="<<wtjac<<"\tnv="<<ndimcurrent_.nv<<"\tigx="<<gridix_.igx<<"\t(1-x1)="<<(1.-x1)<<"\t1/(1-x1)="<<1./(1.-x1)<<std::endl;
// //PDF += xfx[pc[currentprocess_.iproc][0]+6];
// //fnloUtils::_myfillwgt = wgt * PDF * pow( as/(2*M_PI) , fnloUtils::ftable[id]->GetTheCoeffTable()->GetNpow() ) ; // valid for gluons only
}
else {
double xfx1[13] = {0};
double xfx2[13] = {0};
evolvepdf_(__x1,sqrt(muf2),xfx1);
evolvepdf_(__x2,sqrt(muf2),xfx2);
const std::vector<std::vector<int > >& pc = fnloUtils::sp.GetInitialValues();
for ( unsigned int ip = 0 ; ip < pc[currentprocess_.iproc].size() ; ip+=2 ) {
PDF += xfx1[pc[currentprocess_.iproc][ip]+6]*xfx2[pc[currentprocess_.iproc][ip+1]+6];
}
}
fnloUtils::_myfillwgt = __w2 * PDF * pow( as/(2*M_PI) , nnlo::GetNPow() //fnloUtils::ftable[id]->GetTheCoeffTable()->GetNpow() ) ; // valid for gluons only
//std::cout<<"\tprocid="<<currentprocess_.iproc<<std::endl;
//std::cout<<"myfil="<<fnloUtils::_myfillwgt<<std::endl;
}
*/
// fnloUtils::ftable[id]->fEvent.SetWeight(wgt);
// fnloUtils::ftable[id]->fEvent.SetX1(parfrac_.x1);
// fnloUtils::ftable[id]->fEvent.SetX2(parfrac_.x2);
// --- sanity check
if ( fnloUtils::sp[currentprocess_.iproc]==-1 ) std::cout<<"Error! iproc="<<currentprocess_.iproc<<"\tProcId="<<fnloUtils::sp[currentprocess_.iproc]<<std::endl;
// --- pass variables to fastNLO
double mur2 = mur2_hook(1);
double muf2 = muf2_hook(1);
fnloUtils::ftable[id]->fEvent.SetProcessId(fnloUtils::sp[currentprocess_.iproc]);
fnloUtils::ftable[id]->fEvent.SetWeight(__w2);
fnloUtils::ftable[id]->fEvent.SetX1(__x1);
fnloUtils::ftable[id]->fEvent.SetX2(__x2);
fnloUtils::ftable[id]->fScenario.SetObservable0(obs);
//double muf2 = muf2_hook(1);
if ( nnlo::IsDIS() ) {
fnloUtils::ftable[id]->fScenario.SetObsScale1(sqrt(mur2));
fnloUtils::ftable[id]->fScenario.SetObsScale2(sqrt(muf2));
}
else {
fnloUtils::ftable[id]->fScenario.SetObsScale1(sqrt(muf2));
fnloUtils::ftable[id]->fScenario.SetObsScale2(sqrt(mur2));
}
// --- fill event within fastNLO
fnloUtils::ftable[id]->Fill();
fnloUtils::ftable[id]->fEvent.Reset(); // needed?
return ;
}
// _____________________________________________________________________ //
void nnlo::fill_ref_fastnlo( const int& id, const double& obs, const double* wt ) {
// Do nothing for the moment
using namespace std;
// test weights for real radiation
// if ( wt[0]!=0 ) std::cout<<" wt/ref="<<fnloUtils::_myfillwgt/wt[0]<<"\t\ttest wt="<<wt[0]<<"\twtfill="<<fnloUtils::_myfillwgt<<std::endl;
// if ( wt[0]!=0 ) std::cout<<" "<<"\t\ttest wt="<<wt[0]<<"\twtfill="<<fnloUtils::_myfillwgt<<std::endl;
if ( wt[0]!=0 && fabs(1-fnloUtils::_myfillwgt/wt[0])>1.e-5 ) {
std::cout<<"ERROR !! Your fill weight is not correctly calculated!"<<std::endl;
cout<<" wt/ref="<<fnloUtils::_myfillwgt/wt[0]<<"\t\ttest wt="<<wt[0]<<"\twtfill="<<fnloUtils::_myfillwgt<<std::endl;
cout<<" currentprocess_.iproc: "<<currentprocess_.iproc<<endl;
cout<<" xregions_.xreg1 : "<<xregions_.xreg1<<"\t1./xr1="<<1./xregions_.xreg1<<endl;
cout<<" xregions_.xreg2 : "<<xregions_.xreg2<<"\t1./xr1="<<1./xregions_.xreg2<<endl;
cout<<" parfrac_.x1 : "<<parfrac_.x1<<"\t1./x1="<<1./parfrac_.x1<<"\t(1-x1)="<<1-parfrac_.x1<<"\t1./(1-x1)="<<1./(1-parfrac_.x1)<<endl;
cout<<" parfrac_.x2 : "<<parfrac_.x2<<"\t1./x2="<<1./parfrac_.x2<<"\t(1-x2)="<<1-parfrac_.x2<<"\t1./(1-x2)="<<1./(1-parfrac_.x2)<<endl;
//cout<<" ndim : "<<ndimcurrent_.ndim<<endl;
cout<<" nv : "<<ndimcurrent_.nv<<endl;
cout<<" iproc : "<<currentprocess_.iproc<<endl;
cout<<" igx : "<<gridix_.igx<<endl;
cout<<" colflag : "<<colourflag_.colflag<<endl;
cout<<" xcoljac : "<<xcoljac_.xcoljac<<"\t\t1./xcoljac="<<1./xcoljac_.xcoljac<<endl;
cout<<" jacobian_.jflag : "<<jacobian_.jflag<<endl;
static int ii=0;
if ( ii++ >50 ) {
std::cout <<" ende test weight and 'x' value for real radiation. use gluon channels only! " <<std::endl;
exit(1);
}
}
}
// _____________________________________________________________________ //
void nnlo::term_fastnlo( const int& id, int nweights ) {
//!<
//!< Terminate fastNLO: Write table to disk
//!<
if ( vegasiterweight_.it != runinfo_.itmax2 ) return;
std::cout << "[term_fastnlo()] Write fastNLO table id="<<id<<" with nweights" <<nweights << std::endl;
//! --- Set fastNLO verbosity level
say::SetGlobalVerbosity(say::DEBUG);
// --- set number of event
std::cout<<"[nnlo::term_fastnlo] VegasIterwgt: "<<vegasiterweight_.swgt<<"\tniave="<<vegasiterweight_.iave<<"\tit="<<vegasiterweight_.it<<std::endl;
+
+ if ( fnloUtils::ftable[id]==NULL ) {
+ std::cout<<"[nnlo::term_fastnlo] Warning! No event was passed into fastNLO table! Do not write empty file!"<<std::endl;
+ return;
+ }
//std::string runid = runid_.srunid;
fnloUtils::ftable[id]->MultiplyCoefficientsByConstant(vegasiterweight_.swgt/runinfo_.itmax2);
fnloUtils::ftable[id]->SetNumberOfEvents(vegasiterweight_.swgt); // *runinfo_.nshota
+
//std::string filename = fnloUtils::ftablename[id] +"."+fnloUtils::GetOrderName( &fnloUtils::fNlNe )+"."+runid+".tab";
std::string filename = fnloUtils::ftablename[id] +"."+fnloUtils::GetOrderName( &fnloUtils::fNlNe )+".tab";
fnloUtils::ftable[id]->SetFilename(filename);
-
+
+
// --- write table
if (fnloUtils::ftable[id]->GetIsWarmup())
std::cout << "[nnlo::term()] fastNLO warmup run finished ..." << std::endl;
else
std::cout << "[nnlo::term()] fastNLO run finished, writing table '" << fnloUtils::ftablename[id] << "'." << std::endl;
// fnloUtils::ftable[id]->WriteTable(filename); // buggy. This needs a fix!
fnloUtils::ftable[id]->WriteTable();
return;
}
#endif

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 7:29 PM (1 d, 7 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805831
Default Alt Text
(38 KB)

Event Timeline