diff --git a/Makefile.am b/Makefile.am index ee288f31b7f283bc64f497415be5576a9676fbcb..c0b7629bab63d885f549eeb3da7af0596d1742e0 100644 --- a/Makefile.am +++ b/Makefile.am @@ -7,7 +7,9 @@ SUBDIRS = minuit/src interfaces/src DY/src DIPOLE/src RT/src EW/src common commo include interfaces/include FastNLO/include FastNLO/include/fastnlotk DiffDIS/include \ DY/include tools/draw/include \ pdf2yaml tools/process \ - reactions/AFB/src \ + pdf2yaml tools/process \ + reactions/HathorSingleTop/src \ + reactions/AFB/src \ reactions/KFactor/src reactions/Fractal_DISNC/src reactions/BaseDISCC/src reactions/Hathor/src reactions/BaseDISNC/src \ reactions/RT_DISNC/src reactions/FFABM_DISNC/src reactions/FFABM_DISCC/src reactions/APPLgrid/src reactions/BaseHVQMNR/src \ reactions/HVQMNR_LHCb_7TeV_beauty/src reactions/HVQMNR_LHCb_7TeV_charm/src reactions/testZMVFNS/src reactions/fastNLO/src/ \ diff --git a/Reactions.txt b/Reactions.txt index de4193dc22e525198698064a19dbebaaf87b44da..e11a6ebcc2f7ed4d549fc6c66ebefd989b7140e8 100644 --- a/Reactions.txt +++ b/Reactions.txt @@ -14,5 +14,6 @@ FFABM_DISCC libffabm_discc_xfitter.so KFactor libkfactor_xfitter.so FONLL_DISNC libfonll_disnc_xfitter.so FONLL_DISCC libfonll_discc_xfitter.so +HathorSingleTop libhathorsingletop_xfitter.so AFB libafb_xfitter.so KMatrix libkmatrix_xfitter.so diff --git a/configure.ac b/configure.ac index 28fece7f484c575e9396f0459201284d6144e4b5..4a3a87f02c293e03589cc542ede11aa94d47b0ed 100644 --- a/configure.ac +++ b/configure.ac @@ -558,6 +558,7 @@ AC_CONFIG_FILES([include/Makefile examples/Makefile python/Makefile xfitter-config + reactions/HathorSingleTop/src/Makefile reactions/AFB/src/Makefile reactions/KFactor/src/Makefile reactions/BaseDISCC/src/Makefile diff --git a/reactions/Hathor/yaml/parameters.yaml b/reactions/Hathor/yaml/parameters.yaml index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0ee32897bd55a7f047b3fcf23abf2adcec58bf73 100644 --- a/reactions/Hathor/yaml/parameters.yaml +++ b/reactions/Hathor/yaml/parameters.yaml @@ -0,0 +1,5 @@ +# Hathor: +MS_MASS: 0 +muR: 1. +muF: 1. +RATIO: 0 diff --git a/reactions/HathorSingleTop/include/HathorPdfxFitter.h b/reactions/HathorSingleTop/include/HathorPdfxFitter.h new file mode 100644 index 0000000000000000000000000000000000000000..20e81c1471538c9f8b36936eab3ccb4ae33cb1fd --- /dev/null +++ b/reactions/HathorSingleTop/include/HathorPdfxFitter.h @@ -0,0 +1,33 @@ +#ifndef Hathor_xFitterPdf_ +#define Hathor_xFitterPdf_ + +#include "HathorPdf.h" +#include "ReactionHathorSingleTop.h" + +class HathorPdfxFitter : public Pdf +{ +protected: + + std::string PDFname; + int imember; + +public: + + HathorPdfxFitter(ReactionHathorSingleTop* ptrReactionTheory); + void GetPdf(double x, double muf, double h[13]); + double GetAlphas(double mu); + void InitMember(int i); + int NumberPdf(void); + std::string GetName(void); + void Foo(); + + // flag which determines if GetPdf() and GetAlphas() methods should call proper routines or return dummy results + // (Hathor uses these at the very early initialisation step, need protection against QCDNUM not initialized) + bool IsValid; + +private: + // pointer to instance inherited from ReactionTheory (allows access to alphas and PDF routines) + ReactionHathorSingleTop* _reactionTheory; +}; + +#endif diff --git a/reactions/HathorSingleTop/include/ReactionHathorSingleTop.h b/reactions/HathorSingleTop/include/ReactionHathorSingleTop.h new file mode 100644 index 0000000000000000000000000000000000000000..64820f61f65854a9a4fb120e275438bc03198b4e --- /dev/null +++ b/reactions/HathorSingleTop/include/ReactionHathorSingleTop.h @@ -0,0 +1,53 @@ +#pragma once + +#include "ReactionTheory.h" + +/** + @class' ReactionHathorSingleTop + + @brief A wrapper class for HathorSingleTop reaction + + Based on the ReactionTheory class. Reads options produces 3d cross section. + + @version 0.1 + @date 2018-07-25 + */ + +// Authors: Laia Parets Peris <laia.parets.peris@desy.de>, Katerina Lipka <katerina.lipka@desy.de> +// transition from pole to MSbar scheme by S. Moch (private communication) + +class HathorSgTopT; +class HathorPdfxFitter; + +class ReactionHathorSingleTop : public ReactionTheory +{ + public: + ReactionHathorSingleTop(); + + ~ReactionHathorSingleTop(); + +// ~ReactionHathorSingleTop(){}; +// ~ReactionHathorSingleTop(const ReactionHathorSingleTop &){}; +// ReactionHathorSingleTop & operator =(const ReactionHathorSingleTop &r){return *(new ReactionHathorSingleTop(r));}; + + public: + virtual string getReactionName() const { return "HathorSingleTop" ;}; + int initAtStart(const string &); + virtual void setDatasetParameters(int dataSetID, map<string,string> pars, map<string,double> dsPars) override; + virtual int compute(int dataSetID, valarray<double> &val, map<string, valarray<double> > &err); + double CalculationRatio(double _mtop, double _mr, double _mf, int _scheme, HathorSgTopT* hathor); + protected: + virtual int parseOptions(){ return 0;}; + + // this is map of key = dataset, value = pointer to Hathor instances, + // one instance per one dataset + std::map<int, HathorSgTopT*> _hathorArray; + + HathorPdfxFitter* _pdf; + int* _rndStore; + int _scheme; + double _mtop; + double _mr; + double _mf; +}; + diff --git a/reactions/HathorSingleTop/src/HathorPdfxFitter.cc b/reactions/HathorSingleTop/src/HathorPdfxFitter.cc new file mode 100644 index 0000000000000000000000000000000000000000..417406af37e493927f7ad0920232a95a13c7ffdd --- /dev/null +++ b/reactions/HathorSingleTop/src/HathorPdfxFitter.cc @@ -0,0 +1,59 @@ +#include "HathorPdfxFitter.h" +#include "ReactionHathorSingleTop.h" + +#include "get_pdfs.h" +#include <iostream> +#include <valarray> + +HathorPdfxFitter::HathorPdfxFitter(ReactionHathorSingleTop *ptrReactionTheory) +{ + PDFname = "xFitterPdf"; + imember = 0; + _reactionTheory = ptrReactionTheory; + IsValid = false; +} + +void HathorPdfxFitter::GetPdf(double x, double muf, double h[13]) +{ + if(!IsValid) + return; + + std::valarray<double> pdf(0.0, 13); + _reactionTheory->xfx(x, muf, &pdf[0]); + for(auto& val : pdf) + val /= x; + + // ordering is: + // 0 1 2 3 4 5 6 7 8 9 10 11 12 + // tb bb cb sb ub db g d u s c b t + + h[AbstractHathor::ATOP] = pdf[0]; + h[AbstractHathor::ABOTTOM] = pdf[1]; + h[AbstractHathor::ACHARM] = pdf[2]; + h[AbstractHathor::ASTRANGE] = pdf[3]; + h[AbstractHathor::AUP] = pdf[4]; + h[AbstractHathor::ADOWN] = pdf[5]; + + h[AbstractHathor::GLUON] = pdf[6]; + + h[AbstractHathor::DOWN] = pdf[7]; + h[AbstractHathor::UP] = pdf[8]; + h[AbstractHathor::STRANGE] = pdf[9]; + h[AbstractHathor::CHARM] = pdf[10]; + h[AbstractHathor::BOTTOM] = pdf[11]; + h[AbstractHathor::TOP ] = pdf[12]; +} + +double HathorPdfxFitter::GetAlphas(double mu) +{ + if(!IsValid) + return 0.0; + + return _reactionTheory->alphaS(mu); +} + +void HathorPdfxFitter::InitMember(int i){ imember=i; } + +int HathorPdfxFitter::NumberPdf(void){ return 1; } + +std::string HathorPdfxFitter::GetName(void){ return(PDFname); } diff --git a/reactions/HathorSingleTop/src/Makefile.am b/reactions/HathorSingleTop/src/Makefile.am new file mode 100644 index 0000000000000000000000000000000000000000..b93f9863c16e3e18c145854596fdc8e54afaa093 --- /dev/null +++ b/reactions/HathorSingleTop/src/Makefile.am @@ -0,0 +1,30 @@ + +# Created by AddReaction.py on 2017-08-07 + +if ENABLE_HATHOR + +# LHAPDF is required by Hathor-2.0 +if ENABLE_LHAPDF + + # HATHOR_ROOT points to the directory containing libHathor.a: + # (1) in Hathor-1.5 this is the same common directory which contains all headers and sources, + # (2) in Hathor-2.0 the library is Hathor-2.0/lib while the headers are in Hathor-2.0/include, + # therefore add both header places below + HATHOR_CXXFLAGS = -I$(HATHOR_ROOT) -I$(HATHOR_ROOT)/../include + + # Hathor-2.0 requires LHAPDF + HATHOR_CXXFLAGS += $(LHAPDF_CPPFLAGS) + + AM_CXXFLAGS = $(HATHOR_CXXFLAGS) -I$(srcdir)/../include -I$(srcdir)/../../../include -I$(srcdir)/../../../interfaces/include -Wall -fPIC -Wno-deprecated -O0 + + lib_LTLIBRARIES = libhathorsingletop_xfitter.la + libhathorsingletop_xfitter_la_SOURCES = ReactionHathorSingleTop.cc HathorPdfxFitter.cc + + # add Hathor library(ies) + libhathorsingletop_xfitter_la_LDFLAGS = $(HATHOR_LIBS) + +endif + +endif + +dist_noinst_HEADERS = ../include ../yaml diff --git a/reactions/HathorSingleTop/src/ReactionHathorSingleTop.cc b/reactions/HathorSingleTop/src/ReactionHathorSingleTop.cc new file mode 100644 index 0000000000000000000000000000000000000000..a2f871b0a6f79e39c43d90b99d80682022413085 --- /dev/null +++ b/reactions/HathorSingleTop/src/ReactionHathorSingleTop.cc @@ -0,0 +1,310 @@ +/* + @file ReactionHathorSingleTop.cc + @date 2018-07-25 + @author AddReaction.py + Created by AddReaction.py on 2018-07-25 +*/ + +#include "ReactionHathorSingleTop.h" +#include "HathorPdfxFitter.h" +#include "Hathor.h" +#include "cstring" + +// the class factories +extern "C" ReactionHathorSingleTop* create() { + return new ReactionHathorSingleTop(); +} + +extern "C" +{ + int rlxd_size(void); + void rlxd_get(int state[]); + void rlxd_reset(int state[]); + void rlxd_init(int level,int seed); +} + +ReactionHathorSingleTop::ReactionHathorSingleTop() +{ + _pdf = NULL; + _rndStore = NULL; + _scheme = -1; + _mtop = -1.0; + _mr = -1.0; + _mf = -1.0; +} + +ReactionHathorSingleTop::~ReactionHathorSingleTop() +{ + if(_rndStore) + delete[] _rndStore; + + // do NOT delete Hathor instances here, because: + // (1) Hathor classes do not have vitual destructors which produces a warning + // (2) Hathor classes do not have destructors at all + + //if(_pdf) + // delete _pdf; + //for(auto item : _hathorArray) + // if(item.second) + // delete item.second; +} + + +// Initialize at the start of the computation +int ReactionHathorSingleTop::initAtStart(const string &s) +{ + // PDFs for Hathor + _pdf = new HathorPdfxFitter(this); + + // random number generator + rlxd_init(1, 1); + int nRnd = rlxd_size(); + //std::cout << " Size of random number array = " << nRnd << "\n"; + _rndStore = new int[nRnd]; + rlxd_get(_rndStore); + + // scheme (perturbative order and pole/MSbar mass treatment) + const string order = GetParamS("Order"); + const int pertubOrder = OrderMap(order); + _scheme = Hathor::LO; + if(pertubOrder > 1) + _scheme = _scheme | Hathor::NLO; + if(pertubOrder > 2) + _scheme = _scheme | Hathor::NNLO; + int msMass = 0; // pole mass by default + if(checkParam("MS_MASS")) + msMass = GetParamI("MS_MASS"); + if(msMass) + _scheme = _scheme | Hathor::MS_MASS; + + // top quark mass + std::string mtopName = "mtp"; + if(!checkParam(mtopName)) + { + std::string str = "F: no top quark mass (\"" + mtopName + "\" parameter) for Hathor"; + hf_errlog_(18081701, str.c_str(), strlen(str.c_str())); + } + _mtop = GetParam("mtp"); + + // renorm. scale + _mr = _mtop; + if(checkParam("muR")) + _mr *= GetParam("muR"); + + // fact. scale + _mf = _mtop; + if(checkParam("muF")) + _mf *= GetParam("muF"); + + std::cout << " Hathor will use:"; + std::cout << " mtop = " << _mtop << "[GeV] "; + std::cout << " renorm. scale = " << _mr << "[GeV] "; + std::cout << " fact. scale = " << _mf << "[GeV]"; + std::cout << std::endl; + return 0; +} + + +void ReactionHathorSingleTop::setDatasetParameters(int dataSetID, map<std::string, std::string> pars, map<std::string, double> dsPars) +{ + + // check if dataset with provided ID already exists + if(_hathorArray.find(dataSetID) != _hathorArray.end()) + { + char str[256]; + sprintf(str, "F: dataset with id = %d already exists", dataSetID); + hf_errlog_(17080701, str, strlen(str)); + } + + // read centre-of-mass energy from provided dataset parameters + // (must be provided) + auto it = pars.find("SqrtS"); + if(it == pars.end()) + { + char str[256]; + sprintf(str, "F: no SqrtS for dataset with id = %d", dataSetID); + hf_errlog_(18081702, str, strlen(str)); + } + double sqrtS = atof(it->second.c_str()); + + // read precision level from provided dataset parameters + // if not specified set to default 2 -> Hathor::MEDIUM + int precisionLevel = Hathor::MEDIUM; + it = pars.find("precisionLevel"); + if(it != pars.end()) + { + precisionLevel = std::pow(10, 2 + atoi(it->second.c_str())); + // check that this setting is allowed + // see in AbstractHathor.h: + // enum ACCURACY { LOW=1000, MEDIUM=10000, HIGH=100000 }; + // and + // precisionLevel = 1 -> Hathor::LOW + // precisionLevel = 2 -> Hathor::MEDIUM + // precisionLevel = 3 -> Hathor::HIGH + if(precisionLevel != Hathor::LOW && precisionLevel != Hathor::MEDIUM && precisionLevel != Hathor::HIGH) + { + char str[256]; + sprintf(str, "F: provided precision level = %d not supported by Hathor", precisionLevel); + hf_errlog_(18081702, str, strlen(str)); + } + } + + // read ppbar from provided dataset parameters + // if not specified assume it is false (pp collisions) + int ppbar = false; + it = pars.find("ppbar"); + if(it != pars.end()) + { + ppbar = atoi(it->second.c_str()); + if(ppbar != 0 && ppbar != 1) + { + char str[256]; + sprintf(str, "F: provided ppbar = %d not recognised (must be 0 or 1)", ppbar); + hf_errlog_(17081103, str, strlen(str)); + } + } + + // read topquark from provided dataset parameters + // if not specified assume it is false (topquark collisions) + int antitopquark = false; + it = pars.find("antitopquark"); + if(it != pars.end()) + { + antitopquark = atoi(it->second.c_str()); + if(antitopquark != 0 && antitopquark != 1) + { + char str[256]; + sprintf(str, "F: provided antitopquark = %d not recognised (must be 0 or 1)", antitopquark); + hf_errlog_(17081103, str, strlen(str)); + } + + } + + // instantiate Hathor + HathorSgTopT* hathor = new HathorSgTopT(*_pdf); + + hathor->sethc2(0.38937911e9); + + // set collision type + if(ppbar) + hathor->setColliderType(Hathor::PPBAR); + else + hathor->setColliderType(Hathor::PP); + + std::cout << "ReactionHathorSingleTop: PP/PPBAR parameter set to " << ppbar << std::endl; + + // set centre-of-mass energy + hathor->setSqrtShad(sqrtS); + std::cout << "ReactionHathorSingleTop: center of mass energy set to " << sqrtS << std::endl; + + // choose TOPQUARK/ANTITOPQUARK + std::cout << "ReactionHathorSingleTop: TOPQUARK/ANTITOPQUARK se to " << antitopquark << std::endl; + if(antitopquark) + { + std::cout << " Antitopquark is set" << std::endl; + hathor->setParticle(SgTop::ANTITOPQUARK); + } + else + { + std::cout << " Topquark is selected" << std::endl; + hathor->setParticle(SgTop::TOPQUARK); + } + + // set scheme + hathor->setScheme(_scheme); + std::cout << "ReactionHathorSingleTop: Setting the scheme" << std::endl; + // set precision level + hathor->setPrecision(precisionLevel); + + // done + hathor->PrintOptions(); + _hathorArray[dataSetID] = hathor; +} + +// Main function to compute results at an iteration +int ReactionHathorSingleTop::compute(int dataSetID, valarray<double> &val, map<string, valarray<double> > &err) +{ + _pdf->IsValid = true; + rlxd_reset(_rndStore); + + HathorSgTopT* hathor = _hathorArray.at(dataSetID); + + int msMass = 0; // pole mass by default + if(checkParam("MS_MASS")) + msMass = GetParamI("MS_MASS"); + if(msMass){ + + double dmtms; + double Lrbar,nfl; + double d1dec,d2dec,aspi; + double const pi = 3.141592653589793; + double const z2 = 1.644934066848226; + double const z3 = 1.202056903159594; + double const ln2= 0.693147180559945; + + double crst; + double valtclo, valtclop, valtclom, valtcnlo, valtcnlop, valtcnlom; + double err1, chi1; + + aspi = hathor->getAlphas(_mr)/(pi); + dmtms = _mtop/100.; + + // decoupling coefficients + nfl = 5.; + + Lrbar = 0.; + d1dec = ( 4./3. + Lrbar ); + d2dec = ( 307./32. + 2.*z2 + 2./3.*z2*ln2 - z3/6. + + 509./72.*Lrbar + 47./24.*pow(Lrbar,2) + - nfl*(71./144. + z2/3. + 13./36.*Lrbar + pow(Lrbar,2)/12.) ); + + + _scheme = Hathor::LO; + hathor->setScheme(_scheme); + + // LO + hathor->getXsection(_mtop,_mr,_mf); + hathor->getResult(0,valtclo,err1,chi1); + + // LO derivatives + hathor->getXsection(_mtop+dmtms,_mr,_mf); + hathor->getResult(0,valtclop,err1,chi1); + + hathor->getXsection(_mtop-dmtms,_mr,_mf); + hathor->getResult(0,valtclom,err1,chi1); + + + _scheme = Hathor::LO | Hathor::NLO; + hathor->setScheme(_scheme) ; + + // NLO + hathor->getXsection(_mtop,_mr,_mf); + hathor->getResult(0,valtcnlo,err1,chi1); + + // NLO derivatives + hathor->getXsection(_mtop+dmtms,_mr,_mf); + hathor->getResult(0,valtcnlop,err1,chi1); + + hathor->getXsection(_mtop-dmtms,_mr,_mf); + hathor->getResult(0,valtcnlom,err1,chi1); + + // add things up + crst = valtcnlo + + aspi* d1dec*_mtop/(2.*dmtms)* (valtcnlop-valtcnlom) + + pow(aspi,2)* d2dec*_mtop/(2.*dmtms)* (valtclop-valtclom) + + pow(aspi*d1dec*_mtop/dmtms,2)/2.* (valtclop-2.*valtclo+valtclom); + + val[0]=crst; + + } + else{ + + hathor->getXsection(_mtop, _mr, _mf); + double dum = 0.0; + val[0] = 0.0; + hathor->getResult(0, val[0], dum); + + } + return 0; +} + diff --git a/reactions/HathorSingleTop/yaml/parameters.yaml b/reactions/HathorSingleTop/yaml/parameters.yaml new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391