include/Pythia8Plugins/LHAPDF6.h

0001 // LHAPDF6.h is a part of the PYTHIA event generator.
0002 // Copyright (C) 2024 Torbjorn Sjostrand.
0003 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
0004 // Please respect the MCnet Guidelines, see GUIDELINES for details.
0005
0006 // This file contains the LHAPDF6 PDF plugin class.
0007
0008 #ifndef Pythia8_LHAPDF6_H
0009 #define Pythia8_LHAPDF6_H
0010
0011 #include "Pythia8/PartonDistributions.h"
0012 #include "Pythia8/Plugins.h"
0013 #include "LHAPDF/LHAPDF.h"
0014 #include <mutex>
0015
0016 namespace Pythia8 {
0017
0018 //==========================================================================
0019
0020 // Containers for PDF sets.
0021
0022 //--------------------------------------------------------------------------
0023
0024 // Class to hold a PDF set, its information, and its uncertainty sets.
0025
0026 class PdfSets {
0027
0028 public:
0029
0030   // Constructors.
0031   PdfSets() {;}
0032   PdfSets(string setName) : info(::LHAPDF::PDFSet(setName)),
0033     pdfs(vector< ::LHAPDF::PDF* >(info.size(), 0)) {;}
0034
0035   // Access a PDF set.
0036   ::LHAPDF::PDF *operator[](unsigned int member) {
0037     if (!pdfs[member]) {
0038       lock_guard<mutex> lck (mtx);
0039       pdfs[member] = info.mkPDF(member);
0040     }
0041     return pdfs[member];
0042   }
0043
0044   // Get number of PDF sets.
0045   int size() {return pdfs.size();}
0046
0047   // PDF sets and info.
0048   ::LHAPDF::PDFSet info;
0049   vector< ::LHAPDF::PDF* > pdfs;
0050   static mutex mtx;
0051
0052 };
0053
0054 mutex PdfSets::mtx;
0055
0056 //==========================================================================
0057
0058 // Provide interface to the LHAPDF6 library of parton densities.
0059
0060 class LHAPDF6 : public PDF {
0061
0062 public:
0063
0064   // Constructor.
0065   LHAPDF6(Pythia*, Settings* settingsPtr, Logger*) :
0066     PDF(), pdf(nullptr), extrapol(false) {
0067     if (settingsPtr == nullptr) return;
0068     sSymmetric(settingsPtr->flag("LHAPDF:sSymmetric"));
0069     cSymmetric(settingsPtr->flag("LHAPDF:cSymmetric"));
0070     bSymmetric(settingsPtr->flag("LHAPDF:bSymmetric"));
0071   }
0072
0073   // Initialization of PDF set.
0074   bool init(int idBeamIn, string setName, int member, Logger* loggerPtr)
0075     override;
0076
0077   // Allow extrapolation beyond boundaries (not implemented).
0078   void setExtrapolate(bool extrapolIn) {extrapol = extrapolIn;}
0079
0080 private:
0081
0082   // The LHAPDF objects.
0083   PdfSets pdfs;
0084   ::LHAPDF::PDF *pdf;
0085   ::LHAPDF::Extrapolator *ext;
0086   bool extrapol;
0087
0088   // Update parton densities.
0089   void xfUpdate(int id, double x, double Q2);
0090
0091   // Check whether x and Q2 values fall inside the fit bounds.
0092   bool insideBounds(double x, double Q2) {
0093     return (x > xMin && x < xMax && Q2 > q2Min && Q2 < q2Max);}
0094
0095   // Return the running alpha_s shipped with the LHAPDF set.
0096   double alphaS(double Q2) { return pdf->alphasQ2(Q2); }
0097
0098   // Return quark masses used in the PDF fit.
0099   double muPDFSave, mdPDFSave, mcPDFSave, msPDFSave, mbPDFSave,
0100          xMin, xMax, q2Min, q2Max;
0101   double mQuarkPDF(int id) {
0102     switch(abs(id)){
0103       case 1: return mdPDFSave;
0104       case 2: return muPDFSave;
0105       case 3: return msPDFSave;
0106       case 4: return mcPDFSave;
0107       case 5: return mbPDFSave;
0108     }
0109     return -1.;
0110  }
0111
0112   // Calculate uncertainties using the LHAPDF prescription.
0113   void calcPDFEnvelope(int, double, double, int);
0114   void calcPDFEnvelope(pair<int,int>, pair<double,double>, double, int);
0115   PDFEnvelope pdfEnvelope;
0116   PDFEnvelope getPDFEnvelope() {return pdfEnvelope;}
0117   static const double PDFMINVALUE;
0118
0119   int nMembersSave;
0120   int nMembers() { return nMembersSave; }
0121
0122 };
0123
0124 //--------------------------------------------------------------------------
0125
0126 // Constants.
0127
0128 const double LHAPDF6::PDFMINVALUE = 1e-10;
0129
0130 //--------------------------------------------------------------------------
0131
0132 // Initialize a parton density function from LHAPDF6.
0133
0134 bool LHAPDF6::init(int idBeamIn, string setName, int member,
0135   Logger* loggerPtr) {
0136   idBeam = idBeamIn;
0137   idBeamAbs = abs(idBeamIn);
0138   isSet = false;
0139
0140   // Find the PDF set. Note, LHAPDF aborts if the PDF does not exist,
0141   // which we avoid with this try/catch statement. Ideally, we would
0142   // check with :LHAPDF::lookupLHAPDFID, but this is not thread safe.
0143   try {
0144     pdfs = PdfSets(setName);
0145   } catch (const std::exception &e) {
0146     loggerPtr->ERROR_MSG("unknown PDF " + setName);
0147     return false;
0148   }
0149
0150   // Find the PDF member.
0151   if (pdfs.size() == 0) {
0152     loggerPtr->ERROR_MSG("could not initialize PDF " + setName);
0153     return false;
0154   } else if (member >= pdfs.size()) {
0155     loggerPtr->ERROR_MSG(setName + " does not contain requested member");
0156     return false;
0157   }
0158   pdf = pdfs[member];
0159   isSet = true;
0160
0161   // Save x and Q2 limits.
0162   xMax  = pdf->xMax();
0163   xMin  = pdf->xMin();
0164   q2Max = pdf->q2Max();
0165   q2Min = pdf->q2Min();
0166
0167   // Store quark masses used in PDF fit.
0168   muPDFSave = pdf->info().get_entry_as<double>("MUp");
0169   mdPDFSave = pdf->info().get_entry_as<double>("MDown");
0170   mcPDFSave = pdf->info().get_entry_as<double>("MCharm");
0171   msPDFSave = pdf->info().get_entry_as<double>("MStrange");
0172   mbPDFSave = pdf->info().get_entry_as<double>("MBottom");
0173   nMembersSave  = pdf->info().get_entry_as<int>("NumMembers");
0174   return true;
0175
0176 }
0177
0178 //--------------------------------------------------------------------------
0179
0180 // Give the parton distribution function set from LHAPDF6.
0181
0182 void LHAPDF6::xfUpdate(int, double x, double Q2) {
0183   if (!isSet) return;
0184
0185   // Freeze at boundary value if PDF is evaluated outside the fit region.
0186   if (x < xMin && !extrapol) x = xMin;
0187   if (x > xMax)    x = xMax;
0188   if (Q2 < q2Min) Q2 = q2Min;
0189   if (Q2 > q2Max) Q2 = q2Max;
0190
0191   // Update values.
0192   xg     = pdf->xfxQ2(21, x, Q2);
0193   xd     = pdf->xfxQ2(1,  x, Q2);
0194   xu     = pdf->xfxQ2(2,  x, Q2);
0195   xdbar  = pdf->xfxQ2(-1, x, Q2);
0196   xubar  = pdf->xfxQ2(-2, x, Q2);
0197   xs     = pdf->xfxQ2(3,  x, Q2);
0198   xc     = pdf->xfxQ2(4,  x, Q2);
0199   xb     = pdf->xfxQ2(5,  x, Q2);
0200   xsbar  = sSymmetricSave ? xs : pdf->xfxQ2(-3, x, Q2);
0201   xcbar  = cSymmetricSave ? xc : pdf->xfxQ2(-4, x, Q2);
0202   xbbar  = bSymmetricSave ? xb : pdf->xfxQ2(-5, x, Q2);
0203   xgamma = pdf->xfxQ2(22, x, Q2);
0204
0205   // idSav = 9 to indicate that all flavours reset.
0206   idSav = 9;
0207
0208 }
0209
0210 //--------------------------------------------------------------------------
0211
0212 // Calculate uncertainties using the LHAPDF prescription.
0213
0214 void LHAPDF6::calcPDFEnvelope(int idNow, double xNow, double Q2NowIn,
0215   int valSea) {
0216   if (!isSet) return;
0217
0218   // Freeze at boundary value if PDF is evaluated outside the fit region.
0219   double x1 = (xNow < xMin && !extrapol) ? xMin : xNow;
0220   if (x1 > xMax) x1 = xMax;
0221   double Q2Now = (Q2NowIn < q2Min) ? q2Min : Q2NowIn;
0222   if (Q2Now > q2Max) Q2Now = q2Max;
0223
0224   // Loop over the members.
0225   vector<double> xfCalc(pdfs.size());
0226   for(int iMem = 0; iMem < pdfs.size(); ++iMem) {
0227     if (valSea==0 || (idNow != 1 && idNow != 2)) {
0228       xfCalc[iMem] = pdfs[iMem]->xfxQ2(idNow, x1, Q2Now);
0229     } else if (valSea==1 && (idNow == 1 || idNow == 2 )) {
0230       xfCalc[iMem] = pdfs[iMem]->xfxQ2(idNow, x1, Q2Now) -
0231         pdfs[iMem]->xfxQ2(-idNow, x1, Q2Now);
0232     } else if (valSea==2 && (idNow == 1 || idNow == 2 )) {
0233       xfCalc[iMem] = pdfs[iMem]->xfxQ2(-idNow, x1, Q2Now);
0234     }
0235   }
0236
0237   // Calculate the uncertainty.
0238   ::LHAPDF::PDFUncertainty xfErr = pdfs.info.uncertainty(xfCalc);
0239   pdfEnvelope.centralPDF = xfErr.central;
0240   pdfEnvelope.errplusPDF = xfErr.errplus;
0241   pdfEnvelope.errminusPDF = xfErr.errminus;
0242   pdfEnvelope.errsymmPDF = xfErr.errsymm;
0243   pdfEnvelope.scalePDF = xfErr.scale;
0244 }
0245
0246 //--------------------------------------------------------------------------
0247
0248 // Calculate uncertainties using the LHAPDF prescription.
0249
0250 void LHAPDF6::calcPDFEnvelope(pair<int,int> idNows, pair<double,double> xNows,
0251   double Q2NowIn, int valSea) {
0252   if (!isSet) return;
0253
0254   // Freeze at boundary value if PDF is evaluated outside the fit region.
0255   double x1 = (xNows.first < xMin && !extrapol) ? xMin : xNows.first;
0256   if (x1 > xMax) x1 = xMax;
0257   double x2 = (xNows.second < xMin && !extrapol) ? xMin : xNows.second;
0258   if (x2 > xMax) x2 = xMax;
0259   double Q2Now = (Q2NowIn < q2Min) ? q2Min : Q2NowIn;
0260   if (Q2Now > q2Max) Q2Now = q2Max;
0261
0262   // Loop over the members.
0263   vector<double> xfCalc(pdfs.size());
0264   pdfEnvelope.pdfMemberVars.resize(pdfs.size());
0265   for(int iMem = 0; iMem < pdfs.size(); ++iMem) {
0266     if        (valSea == 0 || (idNows.first != 1 && idNows.first != 2 ) ) {
0267       xfCalc[iMem] = pdfs[iMem]->xfxQ2(idNows.first, x1, Q2Now);
0268     } else if (valSea == 1 && (idNows.first == 1 || idNows.first == 2)) {
0269       xfCalc[iMem] = pdfs[iMem]->xfxQ2(idNows.first, x1, Q2Now)
0270         - pdfs[iMem]->xfxQ2(-idNows.first, x1, Q2Now);
0271     } else if (valSea == 2 && (idNows.first == 1 || idNows.first == 2 )) {
0272       xfCalc[iMem] = pdfs[iMem]->xfxQ2(-idNows.first, x1, Q2Now);
0273     }
0274     xfCalc[iMem] = max(0.0, xfCalc[iMem]);
0275     if        (valSea == 0 || (idNows.second != 1 && idNows.second != 2)) {
0276       xfCalc[iMem] /= max
0277         (PDFMINVALUE, pdfs[iMem]->xfxQ2(idNows.second, x2, Q2Now));
0278     } else if (valSea == 1 && (idNows.second == 1 || idNows.second == 2 )) {
0279       xfCalc[iMem] /= max
0280         (pdfs[iMem]->xfxQ2(idNows.second, x2, Q2Now) - pdfs[iMem]->xfxQ2
0281          (-idNows.second, x2, Q2Now), PDFMINVALUE);
0282     } else if (valSea == 2 && (idNows.second == 1 || idNows.second == 2 )) {
0283       xfCalc[iMem] /= max
0284         (pdfs[iMem]->xfxQ2(-idNows.second, x2, Q2Now), PDFMINVALUE);
0285     }
0286     pdfEnvelope.pdfMemberVars[iMem] = xfCalc[iMem];
0287   }
0288
0289   // Calculate the uncertainty.
0290   ::LHAPDF::PDFUncertainty xfErr = pdfs.info.uncertainty(xfCalc);
0291   pdfEnvelope.centralPDF = xfErr.central;
0292   pdfEnvelope.errplusPDF = xfErr.errplus;
0293   pdfEnvelope.errminusPDF = xfErr.errminus;
0294   pdfEnvelope.errsymmPDF = xfErr.errsymm;
0295   pdfEnvelope.scalePDF = xfErr.scale;
0296
0297 }
0298
0299 //--------------------------------------------------------------------------
0300
0301 // Declare the plugin.
0302
0303 PYTHIA8_PLUGIN_CLASS(PDF, LHAPDF6, false, false, false)
0304 PYTHIA8_PLUGIN_VERSIONS(PYTHIA_VERSION_INTEGER)
0305
0306 //==========================================================================
0307
0308 } // end namespace Pythia8
0309
0310 #endif // end Pythia8_LHAPDF6_H