include/Pythia8Plugins/LHAPDF5.h

0001 // LHAPDF5.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 LHAPDF5 PDF plugin class.
0007
0008 #ifndef Pythia8_LHAPDF5_H
0009 #define Pythia8_LHAPDF5_H
0010
0011 #include "Pythia8/PartonDistributions.h"
0012 #include "Pythia8/Plugins.h"
0013
0014 namespace Pythia8 {
0015
0016 //==========================================================================
0017
0018 // Interfaces to to make the C++ calls similar to f77.
0019
0020 //--------------------------------------------------------------------------
0021
0022 // Declare the LHAPDF5 f77 subroutines that are needed.
0023
0024 extern "C" {
0025
0026   extern void initpdfsetm_(int&, const char*, int);
0027
0028   extern void initpdfsetbynamem_(int&, const char*, int);
0029
0030   extern void initpdfm_(int&, int&);
0031
0032   extern void evolvepdfm_(int&, double&, double&, double*);
0033
0034   extern void evolvepdfpm_(int&, double&, double&, double&, double&, double*);
0035
0036   extern void evolvepdfphotonm_(int&, double&, double&, double*, double&);
0037
0038   extern void setlhaparm_(const char*, int);
0039
0040   extern void getxminm_(int &, int &, double &);
0041
0042   extern void getxmaxm_(int &, int &, double &);
0043
0044   extern void getq2minm_(int &, int &, double &);
0045
0046   extern void getq2maxm_(int &, int &, double &);
0047
0048 }
0049
0050 //--------------------------------------------------------------------------
0051
0052 // Map the f77 routines to C++.
0053
0054 namespace LHAPDF5Interface {
0055
0056   // Initialize set with full pathname, allowing multiple sets.
0057   void initPDFsetM( int& nSet, string name) {
0058     const char* cName = name.c_str(); int lenName = name.length();
0059     initpdfsetm_( nSet, cName, lenName);
0060   }
0061
0062   // Initialize set with simple name, allowing multiple sets.
0063   void initPDFsetByNameM( int& nSet, string name) {
0064     const char* cName = name.c_str(); int lenName = name.length();
0065     initpdfsetbynamem_( nSet, cName, lenName);
0066   }
0067
0068   // Initialize member of set.
0069   void initPDFM(int& nSet, int member) {
0070     initpdfm_(nSet, member);
0071   }
0072
0073   // Evaluate x f_i(x, Q).
0074   void evolvePDFM( int& nSet, double x, double Q, double* xfArray) {
0075     evolvepdfm_( nSet, x, Q, xfArray);
0076   }
0077
0078   // Evaluate x f_i(x, Q) for photon beams.
0079   void evolvePDFpM( int& nSet, double x, double Q, double P2, double IP2,
0080     double* xfArray) {
0081     evolvepdfpm_( nSet, x, Q, P2, IP2, xfArray);
0082   }
0083
0084   // Evaluate x f_i(x, Q) including photon
0085   void evolvePDFPHOTONM(int& nSet, double x, double Q, double* xfArray,
0086     double& xPhoton) {
0087     evolvepdfphotonm_( nSet, x, Q, xfArray, xPhoton);
0088   }
0089
0090   // Extrapolate PDF set beyond boundaries, or freeze them there.
0091   void setPDFparm(string name) {
0092     const char* cName = name.c_str(); int lenName = name.length();
0093     setlhaparm_( cName, lenName);
0094   }
0095
0096   // Simple structure to hold LHAPDF set information.
0097   struct LHAPDFInfo {
0098     string name;
0099     int member;
0100     bool photon;
0101   };
0102
0103   // Global tracking of opened PDF sets.
0104   map<int, LHAPDFInfo> initializedSets;
0105
0106   // Method to find the nSet number corresponding to a name and member.
0107   // Returns -1 if no such LHAPDF5 set has been initialized.
0108   int findNSet(string setName, int member) {
0109     for (map<int, LHAPDFInfo>::const_iterator i = initializedSets.begin();
0110       i != initializedSets.end(); ++i) {
0111       int    iSet    = i->first;
0112       string iName   = i->second.name;
0113       int    iMember = i->second.member;
0114       if (iName == setName && iMember == member) return iSet;
0115     }
0116     return -1;
0117   }
0118
0119   // Method to return the lowest non-occupied nSet number.
0120   int freeNSet() {
0121     for (int iSet = 1; iSet <= int(initializedSets.size()); ++iSet) {
0122       if (initializedSets.find(iSet) == initializedSets.end())
0123         return iSet;
0124     }
0125     return initializedSets.size() + 1;
0126   }
0127
0128 }
0129
0130 //==========================================================================
0131
0132 // Plugin interface to the LHAPDF5 library.
0133
0134 //==========================================================================
0135
0136 // Provide plugin interface to the LHAPDF5 library of parton densities.
0137
0138 class LHAPDF5 : public PDF {
0139
0140 public:
0141
0142   // Constructor.
0143   LHAPDF5(Pythia*, Settings* settingsPtr, Logger*) :
0144     PDF(), doExtraPol(false) {
0145     if (settingsPtr == nullptr) return;
0146     sSymmetric(settingsPtr->flag("LHAPDF:sSymmetric"));
0147     cSymmetric(settingsPtr->flag("LHAPDF:cSymmetric"));
0148     bSymmetric(settingsPtr->flag("LHAPDF:bSymmetric"));
0149   }
0150
0151   // Initialization of PDF set.
0152   bool init(int idBeamIn, string setName, int member, Logger* loggerPtr)
0153     override;
0154
0155   // Allow extrapolation beyond boundaries. This is optional.
0156   void setExtrapolate(bool extrapol);
0157
0158 private:
0159
0160   // Update all PDF values.
0161   void xfUpdate(int , double x, double Q2);
0162
0163   // Current set and pdf values.
0164   bool   doExtraPol;
0165   int    nSet;
0166   double xfArray[13];
0167   bool   hasPhoton, isPhoton;
0168   double xPhoton;
0169
0170 };
0171
0172 //--------------------------------------------------------------------------
0173
0174 // Initialize a parton density function from LHAPDF5.
0175
0176 bool LHAPDF5::init(int idBeamIn, string setName, int member,
0177   Logger* loggerPtr) {
0178   idBeam = idBeamIn;
0179   idBeamAbs = abs(idBeamIn);
0180   isPhoton = idBeamIn == 22;
0181   nSet = LHAPDF5Interface::findNSet(setName, member);
0182   if (nSet == -1) nSet = LHAPDF5Interface::freeNSet();
0183
0184   // If already initialized then need not do anything further.
0185   LHAPDF5Interface::LHAPDFInfo initializedInfo =
0186     LHAPDF5Interface::initializedSets[nSet];
0187   string initializedSetName   = initializedInfo.name;
0188   int    initializedMember    = initializedInfo.member;
0189   hasPhoton                   = initializedInfo.photon;
0190   if (setName == initializedSetName && member == initializedMember)
0191     return true;
0192
0193   // Initialize set. If first character is '/' then assume that name
0194   // is given with path, else not.
0195   if (setName[0] == '/') LHAPDF5Interface::initPDFsetM( nSet, setName);
0196   else LHAPDF5Interface::initPDFsetByNameM( nSet, setName);
0197   isSet = (nSet >= 0);
0198
0199   // Initialize member.
0200   LHAPDF5Interface::initPDFM(nSet, member);
0201
0202   // Do not collect statistics on under/overflow to save time and space.
0203   LHAPDF5Interface::setPDFparm( "NOSTAT" );
0204   LHAPDF5Interface::setPDFparm( "LOWKEY" );
0205
0206   // Check if photon PDF available (has_photon does not work properly).
0207   xPhoton = 0;
0208   LHAPDF5Interface::evolvePDFPHOTONM(nSet, 0.01, 1, xfArray, xPhoton);
0209   hasPhoton = xPhoton != 0;
0210
0211   // Save values to avoid unnecessary reinitializations.
0212   initializedInfo.name   = setName;
0213   initializedInfo.member = member;
0214   initializedInfo.photon = hasPhoton;
0215   if (nSet > 0) LHAPDF5Interface::initializedSets[nSet] = initializedInfo;
0216   return true;
0217
0218 }
0219
0220 //--------------------------------------------------------------------------
0221
0222 // Allow optional extrapolation beyond boundaries.
0223
0224 void LHAPDF5::setExtrapolate(bool extrapol) {
0225
0226   doExtraPol = extrapol;
0227   LHAPDF5Interface::setPDFparm( (extrapol) ? "EXTRAPOLATE" : "18" );
0228
0229 }
0230
0231 //--------------------------------------------------------------------------
0232
0233 // Give the parton distribution function set from LHAPDF5.
0234
0235 void LHAPDF5::xfUpdate(int, double x, double Q2) {
0236
0237   // Freeze at boundary value if PDF is evaluated outside the fit region.
0238   int member = LHAPDF5Interface::initializedSets[nSet].member;
0239   double xMin, xMax, q2Min, q2Max;
0240   getxminm_( nSet, member, xMin);
0241   getxmaxm_( nSet, member, xMax);
0242   getq2minm_(nSet, member, q2Min);
0243   getq2maxm_(nSet, member, q2Max);
0244   if (x < xMin && !doExtraPol) x = xMin;
0245   if (x > xMax)    x = xMax;
0246   if (Q2 < q2Min) Q2 = q2Min;
0247   if (Q2 > q2Max) Q2 = q2Max;
0248   double Q = sqrt( max( 0., Q2));
0249
0250   // Use special call if photon included in proton.
0251   if (hasPhoton) {
0252     LHAPDF5Interface::evolvePDFPHOTONM( nSet, x, Q, xfArray, xPhoton);
0253   }
0254
0255   // Use special call with photon beams. No virtualities implemented yet.
0256   else if (isPhoton) {
0257     LHAPDF5Interface::evolvePDFpM( nSet, x, Q, 0., 0., xfArray);
0258   }
0259
0260   // Else use default LHAPDF5 call.
0261   else {
0262     LHAPDF5Interface::evolvePDFM( nSet, x, Q, xfArray);
0263     xPhoton=0.0;
0264   }
0265
0266   // Update values.
0267   xg     = xfArray[6];
0268   xd     = xfArray[7];
0269   xdbar  = xfArray[5];
0270   xu     = xfArray[8];
0271   xubar  = xfArray[4];
0272   xs     = xfArray[9];
0273   xc     = xfArray[10];
0274   xb     = xfArray[11];
0275   xsbar  = sSymmetricSave ? xs : xfArray[3];
0276   xcbar  = cSymmetricSave ? xc : xfArray[2];
0277   xbbar  = bSymmetricSave ? xb : xfArray[1];
0278   xgamma = xPhoton;
0279
0280   // idSav = 9 to indicate that all flavours reset.
0281   idSav = 9;
0282
0283 }
0284
0285 //--------------------------------------------------------------------------
0286
0287 // Declare the plugin.
0288
0289 PYTHIA8_PLUGIN_CLASS(PDF, LHAPDF5, false, false, false)
0290 PYTHIA8_PLUGIN_VERSIONS(PYTHIA_VERSION_INTEGER)
0291
0292 //==========================================================================
0293
0294 } // end namespace Pythia8
0295
0296 #endif // end Pythia8_LHAPDF5_H