// -*- C++ -*- // // This file is part of LHAPDF // Copyright (C) 2012-2021 The LHAPDF collaboration (see AUTHORS for details) // #include "LHAPDF/Info.h" #include "LHAPDF/Config.h" #include "LHAPDF/PDFSet.h" #include "LHAPDF/PDFInfo.h" #include "LHAPDF/PDF.h" #include "LHAPDF/GridPDF.h" #include "LHAPDF/BilinearInterpolator.h" #include "LHAPDF/BicubicInterpolator.h" #include "LHAPDF/LogBilinearInterpolator.h" #include "LHAPDF/LogBicubicInterpolator.h" #include "LHAPDF/ErrExtrapolator.h" #include "LHAPDF/NearestPointExtrapolator.h" #include "LHAPDF/ContinuationExtrapolator.h" #include "LHAPDF/AlphaS.h" namespace LHAPDF { Info& getConfig() { return Config::get(); } PDFSet& getPDFSet(const string& setname) { static map _sets; map::iterator it = _sets.find(setname); if (it != _sets.end()) return it->second; _sets[setname] = PDFSet(setname); return _sets[setname]; } PDFInfo* mkPDFInfo(const std::string& setname, int member) { //return new Info(findpdfmempath(setname, member)); return new PDFInfo(setname, member); } PDFInfo* mkPDFInfo(int lhaid) { const pair setname_memid = lookupPDF(lhaid); return mkPDFInfo(setname_memid.first, setname_memid.second); } PDF* mkPDF(const string& setname, int member) { // Find the member data file and ensure that it exists const string searchpath = findpdfmempath(setname, member); if (searchpath.empty()) { const int setsize = getPDFSet(setname).size(); if (member > setsize-1) throw UserError("PDF " + setname + "/" + to_str(member) + " is out of the member range of set " + setname); throw UserError("Can't find a valid PDF " + setname + "/" + to_str(member)); } // First create an Info object to work out what format of PDF this is: Info info(searchpath); const string fmt = info.get_entry("Format"); // Then use the format information to call the appropriate concrete PDF constructor: if (fmt == "lhagrid1") return new GridPDF(setname, member); /// @todo Throw a deprecation error if format version is too old or new throw FactoryError("No LHAPDF factory defined for format type '" + fmt + "'"); } PDF* mkPDF(const string& setname_nmem) { const pair idpair = lookupPDF(setname_nmem); return mkPDF(idpair.first, idpair.second); } PDF* mkPDF(int lhaid) { const pair setname_nmem = lookupPDF(lhaid); return mkPDF(setname_nmem.first, setname_nmem.second); } void mkPDFs(const string& setname, vector& pdfs) { getPDFSet(setname).mkPDFs(pdfs); } vector mkPDFs(const string& setname) { return getPDFSet(setname).mkPDFs(); } Interpolator* mkInterpolator(const string& name) { // Convert name to lower case for comparisons const string iname = to_lower(name); if (iname == "linear") return new BilinearInterpolator(); else if (iname == "cubic") return new BicubicInterpolator(); else if (iname == "log") return new LogBilinearInterpolator(); else if (iname == "logcubic") return new LogBicubicInterpolator(); else throw FactoryError("Undeclared interpolator requested: " + name); } Extrapolator* mkExtrapolator(const string& name) { // Convert name to lower case for comparisons const string xname = to_lower(name); if (xname == "nearest") return new NearestPointExtrapolator(); else if (xname == "error") return new ErrExtrapolator(); else if (xname == "continuation") return new ContinuationExtrapolator(); else { throw FactoryError("Undeclared extrapolator requested: " + name); } } AlphaS* mkAlphaS(const Info& info) { AlphaS* as = 0; const string itype = to_lower(info.get_entry("AlphaS_Type")); if (itype == "analytic") as = new AlphaS_Analytic(); else if (itype == "ode") as = new AlphaS_ODE(); else if (itype == "ipol") as = new AlphaS_Ipol(); else throw FactoryError("Undeclared AlphaS requested: " + itype); // Configure the QCD params on this AlphaS if (info.has_key("AlphaS_OrderQCD")) as->setOrderQCD(info.get_entry_as("AlphaS_OrderQCD")); /// @todo Fall back to generic OrderQCD? (Default value is 4 at the moment --KN) /// Thresholds are where the flavor transition happens in the ODE/analytic solvers, AlphaS_MQ /// should be given as \bar{MQ} (\bar{MQ}) if it is to be used for thresholds away from the mass. /// Since you need the heavy quark mass to calculate the decoupling in the ODE solver /// when threshold =/= mass, I've implemented this as AlphaS_Threshold* -> Threshold*, AlphaS_M* -> M* if (info.has_key("AlphaS_ThresholdDown") && info.has_key("AlphaS_ThresholdUp") && info.has_key("AlphaS_ThresholdStrange") && info.has_key("AlphaS_ThresholdCharm") && info.has_key("AlphaS_ThresholdBottom") && info.has_key("AlphaS_ThresholdTop")) { as->setQuarkThreshold(1, info.get_entry_as("AlphaS_ThresholdDown")); as->setQuarkThreshold(2, info.get_entry_as("AlphaS_ThresholdUp")); as->setQuarkThreshold(3, info.get_entry_as("AlphaS_ThresholdStrange")); as->setQuarkThreshold(4, info.get_entry_as("AlphaS_ThresholdCharm")); as->setQuarkThreshold(5, info.get_entry_as("AlphaS_ThresholdBottom")); as->setQuarkThreshold(6, info.get_entry_as("AlphaS_ThresholdTop")); } else if (info.has_key("ThresholdDown") && info.has_key("ThresholdUp") && info.has_key("ThresholdStrange") && info.has_key("ThresholdCharm") && info.has_key("ThresholdBottom") && info.has_key("ThresholdTop")) { as->setQuarkThreshold(1, info.get_entry_as("ThresholdDown")); as->setQuarkThreshold(2, info.get_entry_as("ThresholdUp")); as->setQuarkThreshold(3, info.get_entry_as("ThresholdStrange")); as->setQuarkThreshold(4, info.get_entry_as("ThresholdCharm")); as->setQuarkThreshold(5, info.get_entry_as("ThresholdBottom")); as->setQuarkThreshold(6, info.get_entry_as("ThresholdTop")); } if (info.has_key("AlphaS_MDown") && info.has_key("AlphaS_MUp") && info.has_key("AlphaS_MStrange") && info.has_key("AlphaS_MCharm") && info.has_key("AlphaS_MBottom") && info.has_key("AlphaS_MTop")) { as->setQuarkMass(1, info.get_entry_as("AlphaS_MDown")); as->setQuarkMass(2, info.get_entry_as("AlphaS_MUp")); as->setQuarkMass(3, info.get_entry_as("AlphaS_MStrange")); as->setQuarkMass(4, info.get_entry_as("AlphaS_MCharm")); as->setQuarkMass(5, info.get_entry_as("AlphaS_MBottom")); as->setQuarkMass(6, info.get_entry_as("AlphaS_MTop")); } // This falls back to lhapdf.conf so should in theory never throw the MetadataError else if (info.has_key("MDown") && info.has_key("MUp") && info.has_key("MStrange") && info.has_key("MCharm") && info.has_key("MBottom") && info.has_key("MTop")) { as->setQuarkMass(1, info.get_entry_as("MDown")); as->setQuarkMass(2, info.get_entry_as("MUp")); as->setQuarkMass(3, info.get_entry_as("MStrange")); as->setQuarkMass(4, info.get_entry_as("MCharm")); as->setQuarkMass(5, info.get_entry_as("MBottom")); as->setQuarkMass(6, info.get_entry_as("MTop")); } else { throw MetadataError("All quark masses required (either as AlphaS_MQ or MQ) for AlphaS."); } const string fscheme = to_lower(info.get_entry("AlphaS_FlavorScheme", info.get_entry("FlavorScheme", "variable"))); // default is VFNS const int nflavs = info.get_entry_as("AlphaS_NumFlavors", info.get_entry_as("NumFlavors", 5)); // default is 5 flavour evolution if (fscheme == "fixed") as->setFlavorScheme(AlphaS::FIXED, nflavs); else if (fscheme == "variable") as->setFlavorScheme(AlphaS::VARIABLE, nflavs); else as->setFlavorScheme(AlphaS::VARIABLE, 5); // default fallback mode // Required parameter settings for each calculation mode if (as->type() == "ode") { /// @todo Handle FFNS / VFNS if ( (!info.has_key("AlphaS_MZ") || !info.has_key("MZ")) || (!info.has_key("AlphaS_MassReference") || !info.has_key("AlphaS_Reference")) ) throw MetadataError("Requested ODE AlphaS but there is no reference point given: define either AlphaS_MZ and MZ, or AlphaS_MassReference and AlphaS_Reference. The latter is given preference if both are defined."); if (info.has_key("AlphaS_MZ")) as->setAlphaSMZ(info.get_entry_as("AlphaS_MZ")); if (info.has_key("MZ"))as->setMZ(info.get_entry_as("MZ")); if (info.has_key("AlphaS_Reference"))as->setAlphaSReference(info.get_entry_as("AlphaS_Reference")); if (info.has_key("AlphaS_MassReference"))as->setMassReference(info.get_entry_as("AlphaS_MassReference")); if (info.has_key("AlphaS_Qs")) { AlphaS_ODE* as_o = dynamic_cast(as); if (info.has_key("AlphaS_Qs")) as_o->setQValues( info.get_entry_as< vector >("AlphaS_Qs")); } } else if (as->type() == "analytic") { /// @todo Handle FFNS / VFNS if (!info.has_key("AlphaS_Lambda5") && !info.has_key("AlphaS_Lambda4") && !info.has_key("AlphaS_Lambda3") ) throw MetadataError("Requested analytic AlphaS but the required parameters are not defined."); if (info.has_key("AlphaS_Lambda3")) as->setLambda(3, info.get_entry_as("AlphaS_Lambda3")); if (info.has_key("AlphaS_Lambda4")) as->setLambda(4, info.get_entry_as("AlphaS_Lambda4")); if (info.has_key("AlphaS_Lambda5")) as->setLambda(5, info.get_entry_as("AlphaS_Lambda5")); } else if (as->type() == "ipol") { if (!info.has_key("AlphaS_Qs") || !info.has_key("AlphaS_Vals") ) throw MetadataError("Requested ipol AlphaS but the required parameters are not defined."); AlphaS_Ipol* as_i = dynamic_cast(as); if (info.has_key("AlphaS_Qs")) as_i->setQValues( info.get_entry_as< vector >("AlphaS_Qs")); if (info.has_key("AlphaS_Vals")) as_i->setAlphaSValues( info.get_entry_as< vector >("AlphaS_Vals")); } return as; } AlphaS* mkAlphaS(const std::string& setname) { return mkAlphaS(getPDFSet(setname)); } AlphaS* mkAlphaS(const std::string& setname, int member) { unique_ptr info( mkPDFInfo(setname, member) ); return mkAlphaS(*info); } AlphaS* mkAlphaS(int lhaid) { unique_ptr info( mkPDFInfo(lhaid) ); return mkAlphaS(*info); } /// @todo Merge to mkAlphaS(0)? AlphaS* mkBareAlphaS(const std::string& type) { AlphaS* as = 0; const string itype = to_lower(type); if (itype == "analytic") as = new AlphaS_Analytic(); else if (itype == "ode") as = new AlphaS_ODE(); else if (itype == "ipol") as = new AlphaS_Ipol(); else throw FactoryError("Undeclared AlphaS requested: " + itype); return as; } }