#include "ATOOLS/Org/Message.H" #include "ATOOLS/Org/Exception.H" #include "ATOOLS/Org/CXXFLAGS_PACKAGES.H" #include "ATOOLS/Org/CXXFLAGS.H" #include "ATOOLS/Org/MyStrStream.H" #include "ATOOLS/Org/Run_Parameter.H" #include "ATOOLS/Org/Scoped_Settings.H" #include "ATOOLS/Math/Random.H" #include "PDF/Main/PDF_Base.H" #include "ATOOLS/Phys/Flavour.H" #include "LHAPDF/LHAPDF.h" namespace PDF { class LHAPDF_CPP_Interface : public PDF_Base { private: LHAPDF::PDF * p_pdf; int m_smember; std::map m_xfx; std::map m_calculated; double m_x,m_Q2; std::vector m_disallowedflavour; public: LHAPDF_CPP_Interface(const ATOOLS::Flavour,std::string,int); ~LHAPDF_CPP_Interface(); PDF_Base * GetCopy(); void CalculateSpec(const double&,const double&); double AlphaSPDF(const double &); double GetXPDF(const ATOOLS::Flavour&); double GetXPDF(const kf_code&, bool); double GetDefaultAlpha(); double GetDefaultScale(); int GetFlavourScheme(); void SetAlphaSInfo(); void SetPDFMember(); }; } using namespace PDF; using namespace ATOOLS; LHAPDF_CPP_Interface::LHAPDF_CPP_Interface(const ATOOLS::Flavour _bunch, const std::string _set, const int _member) : p_pdf(NULL) { m_set=_set; m_smember=_member; m_type="LHA["+m_set+"]"; Scoped_Settings s{ Settings::GetMainSettings()["LHAPDF"] }; m_bunch = _bunch; static std::set s_init; if (s_init.find(m_set)==s_init.end()) { m_member=abs(m_smember); int lhapdfverb(LHAPDF::verbosity()); LHAPDF::setVerbosity(msg_LevelIsDebugging()?lhapdfverb:0); p_pdf = LHAPDF::mkPDF(m_set,m_smember); LHAPDF::setVerbosity(lhapdfverb); SetAlphaSInfo(); } auto q2lim = s["USE_Q2LIMIT"].SetDefault(1).Get(); // get x,Q2 ranges from PDF m_xmin=p_pdf->xMin(); m_xmax=p_pdf->xMax(); m_q2min=q2lim?p_pdf->q2Min():0.0; m_q2max=q2lim?p_pdf->q2Max():1.0e37; m_nf=m_asinfo.m_nf; // initialise all book-keep arrays etc. std::vector kfcs; kfcs.push_back(kf_d); kfcs.push_back(-kf_d); kfcs.push_back(kf_u); kfcs.push_back(-kf_u); kfcs.push_back(kf_s); kfcs.push_back(-kf_s); kfcs.push_back(kf_c); kfcs.push_back(-kf_c); kfcs.push_back(kf_b); kfcs.push_back(-kf_b); kfcs.push_back(kf_t); kfcs.push_back(-kf_t); kfcs.push_back(kf_e); kfcs.push_back(-kf_e); kfcs.push_back(kf_mu); kfcs.push_back(-kf_mu); kfcs.push_back(kf_tau); kfcs.push_back(-kf_tau); kfcs.push_back(kf_gluon); kfcs.push_back(kf_photon); kfcs.push_back(kf_Z); kfcs.push_back(kf_Wplus); kfcs.push_back(-kf_Wplus); kfcs.push_back(kf_h0); for (int i=0;ihasFlavor(kfcs[i])) { m_partons.insert(Flavour(abs(kfcs[i]),kfcs[i]<0)); m_xfx[kfcs[i]]=0.; m_calculated[kfcs[i]]=false; } if (p_pdf->hasFlavor(kf_d)) { m_partons.insert(Flavour(kf_quark)); if (p_pdf->hasFlavor(kf_gluon)) { m_partons.insert(Flavour(kf_jet)); if (p_pdf->hasFlavor(kf_photon)) { m_partons.insert(Flavour(kf_ewjet)); } } } m_lhef_number = p_pdf->lhapdfID(); m_disallowedflavour = s["DISALLOW_FLAVOUR"].GetVector(); if (m_disallowedflavour.size()) { msg_Info()<gen.AddCitation(1,"LHAPDF6 is published under \\cite{Buckley:2014ana}."); } double LHAPDF_CPP_Interface::GetDefaultAlpha() { return m_asinfo.m_asmz; } int LHAPDF_CPP_Interface::GetFlavourScheme() { int nflav=p_pdf->info().get_entry_as("NumFlavors"); if (p_pdf->info().get_entry_as("FlavorScheme")=="variable") { if (nflav==6){ return -1; } nflav+=10; } return nflav; } double LHAPDF_CPP_Interface::GetDefaultScale() { return m_asinfo.m_mz2; } void LHAPDF_CPP_Interface::SetAlphaSInfo() { if (m_asinfo.m_order>=0) return; // TODO: get alphaS info m_asinfo.m_order=p_pdf->info().get_entry_as("AlphaS_OrderQCD"); m_asinfo.m_nf=p_pdf->info().get_entry_as("NumFlavors",-1); if (m_asinfo.m_nf<0) { Scoped_Settings s{ Settings::GetMainSettings()["LHAPDF"] }; const int nf(s["NUMBER_OF_FLAVOURS"].Get()); msg_Info()<info().get_entry_as("MDown"); else if (i==1) m_asinfo.m_flavs[i].m_mass=m_asinfo.m_flavs[i].m_thres =p_pdf->info().get_entry_as("MUp"); else if (i==2) m_asinfo.m_flavs[i].m_mass=m_asinfo.m_flavs[i].m_thres =p_pdf->info().get_entry_as("MStrange"); else if (i==3) m_asinfo.m_flavs[i].m_mass=m_asinfo.m_flavs[i].m_thres =p_pdf->info().get_entry_as("MCharm"); else if (i==4) m_asinfo.m_flavs[i].m_mass=m_asinfo.m_flavs[i].m_thres =p_pdf->info().get_entry_as("MBottom"); else if (i==5) m_asinfo.m_flavs[i].m_mass=m_asinfo.m_flavs[i].m_thres =p_pdf->info().get_entry_as("MTop"); } m_asinfo.m_asmz=p_pdf->info().get_entry_as("AlphaS_MZ"); m_asinfo.m_mz2=sqr(p_pdf->info().get_entry_as("MZ")); } LHAPDF_CPP_Interface::~LHAPDF_CPP_Interface() { if (p_pdf) { delete p_pdf; p_pdf=NULL; } } PDF_Base * LHAPDF_CPP_Interface::GetCopy() { return new LHAPDF_CPP_Interface(m_bunch,m_set,m_smember); } double LHAPDF_CPP_Interface::AlphaSPDF(const double &scale2) { if (IsBad(scale2) || scale2<0.0) { msg_Error()<alphasQ2(scale2); } void LHAPDF_CPP_Interface::SetPDFMember() { if (m_smember<0) { THROW(not_implemented,"Not implemented yet.") double rn=ran->Get(); m_member=1+Min((int)(rn*abs(m_smember)),-m_smember-1); //p_pdf->initPDF(m_member); } } void LHAPDF_CPP_Interface::CalculateSpec(const double& x,const double& Q2) { for (std::map::iterator it=m_calculated.begin(); it!=m_calculated.end();++it) it->second=false; m_x=x/m_rescale; m_Q2=Q2; } double LHAPDF_CPP_Interface::GetXPDF(const ATOOLS::Flavour& infl) { return GetXPDF(infl.Kfcode(), infl.IsAnti()); } double LHAPDF_CPP_Interface::GetXPDF(const kf_code& kf, bool anti) { if (IsBad(m_x) || IsBad(m_Q2) || m_Q2<0.0) { msg_Error()<xfxQ2(kfc,m_x,m_Q2); m_calculated[kfc]=true; } return m_rescale*m_xfx[kfc]; } DECLARE_PDF_GETTER(LHAPDF_Getter); PDF_Base *LHAPDF_Getter::operator() (const Parameter_Type &args) const { if (!args.m_bunch.IsHadron() && !args.m_bunch.IsPhoton()) return NULL; return new LHAPDF_CPP_Interface(args.m_bunch,args.m_set,args.m_member); } void LHAPDF_Getter::PrintInfo (std::ostream &str,const size_t width) const { str<<"LHAPDF interface"; } std::vector p_get_lhapdf; extern "C" void InitPDFLib() { Scoped_Settings s{ Settings::GetMainSettings()["LHAPDF"] }; if (s["GRID_PATH"].IsSetExplicitly()) LHAPDF::setPaths(s["GRID_PATH"].Get()); const std::vector& sets(LHAPDF::availablePDFSets()); msg_Debugging()<