LHAPDF5.h

// LHAPDF5.h is a part of the PYTHIA event generator.
// Copyright (C) 2024 Torbjorn Sjostrand.
// PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
// Please respect the MCnet Guidelines, see GUIDELINES for details.

// This file contains the LHAPDF5 PDF plugin class.

#ifndef Pythia8_LHAPDF5_H
#define Pythia8_LHAPDF5_H

#include "Pythia8/PartonDistributions.h"
#include "Pythia8/Plugins.h"

namespace Pythia8 {

//==========================================================================

// Interfaces to to make the C++ calls similar to f77.

//--------------------------------------------------------------------------

// Declare the LHAPDF5 f77 subroutines that are needed.

extern "C" {

  extern void initpdfsetm_(int&, const char*, int);

  extern void initpdfsetbynamem_(int&, const char*, int);

  extern void initpdfm_(int&, int&);

  extern void evolvepdfm_(int&, double&, double&, double*);

  extern void evolvepdfpm_(int&, double&, double&, double&, double&, double*);

  extern void evolvepdfphotonm_(int&, double&, double&, double*, double&);

  extern void setlhaparm_(const char*, int);

  extern void getxminm_(int &, int &, double &);

  extern void getxmaxm_(int &, int &, double &);

  extern void getq2minm_(int &, int &, double &);

  extern void getq2maxm_(int &, int &, double &);

}

//--------------------------------------------------------------------------

// Map the f77 routines to C++.

namespace LHAPDF5Interface {

  // Initialize set with full pathname, allowing multiple sets.
  void initPDFsetM( int& nSet, string name) {
    const char* cName = name.c_str(); int lenName = name.length();
    initpdfsetm_( nSet, cName, lenName);
  }

  // Initialize set with simple name, allowing multiple sets.
  void initPDFsetByNameM( int& nSet, string name) {
    const char* cName = name.c_str(); int lenName = name.length();
    initpdfsetbynamem_( nSet, cName, lenName);
  }

  // Initialize member of set.
  void initPDFM(int& nSet, int member) {
    initpdfm_(nSet, member);
  }

  // Evaluate x f_i(x, Q).
  void evolvePDFM( int& nSet, double x, double Q, double* xfArray) {
    evolvepdfm_( nSet, x, Q, xfArray);
  }

  // Evaluate x f_i(x, Q) for photon beams.
  void evolvePDFpM( int& nSet, double x, double Q, double P2, double IP2,
    double* xfArray) {
    evolvepdfpm_( nSet, x, Q, P2, IP2, xfArray);
  }

  // Evaluate x f_i(x, Q) including photon
  void evolvePDFPHOTONM(int& nSet, double x, double Q, double* xfArray,
    double& xPhoton) {
    evolvepdfphotonm_( nSet, x, Q, xfArray, xPhoton);
  }

  // Extrapolate PDF set beyond boundaries, or freeze them there.
  void setPDFparm(string name) {
    const char* cName = name.c_str(); int lenName = name.length();
    setlhaparm_( cName, lenName);
  }

  // Simple structure to hold LHAPDF set information.
  struct LHAPDFInfo {
    string name;
    int member;
    bool photon;
  };

  // Global tracking of opened PDF sets.
  map<int, LHAPDFInfo> initializedSets;

  // Method to find the nSet number corresponding to a name and member.
  // Returns -1 if no such LHAPDF5 set has been initialized.
  int findNSet(string setName, int member) {
    for (map<int, LHAPDFInfo>::const_iterator i = initializedSets.begin();
      i != initializedSets.end(); ++i) {
      int    iSet    = i->first;
      string iName   = i->second.name;
      int    iMember = i->second.member;
      if (iName == setName && iMember == member) return iSet;
    }
    return -1;
  }

  // Method to return the lowest non-occupied nSet number.
  int freeNSet() {
    for (int iSet = 1; iSet <= int(initializedSets.size()); ++iSet) {
      if (initializedSets.find(iSet) == initializedSets.end())
        return iSet;
    }
    return initializedSets.size() + 1;
  }

}

//==========================================================================

// Plugin interface to the LHAPDF5 library.

//==========================================================================

// Provide plugin interface to the LHAPDF5 library of parton densities.

class LHAPDF5 : public PDF {

public:

  // Constructor.
  LHAPDF5(Pythia*, Settings* settingsPtr, Logger*) :
    PDF(), doExtraPol(false) {
    if (settingsPtr == nullptr) return;
    sSymmetric(settingsPtr->flag("LHAPDF:sSymmetric"));
    cSymmetric(settingsPtr->flag("LHAPDF:cSymmetric"));
    bSymmetric(settingsPtr->flag("LHAPDF:bSymmetric"));
  }

  // Initialization of PDF set.
  bool init(int idBeamIn, string setName, int member, Logger* loggerPtr)
    override;

  // Allow extrapolation beyond boundaries. This is optional.
  void setExtrapolate(bool extrapol);

private:

  // Update all PDF values.
  void xfUpdate(int , double x, double Q2);

  // Current set and pdf values.
  bool   doExtraPol;
  int    nSet;
  double xfArray[13];
  bool   hasPhoton, isPhoton;
  double xPhoton;

};

//--------------------------------------------------------------------------

// Initialize a parton density function from LHAPDF5.

bool LHAPDF5::init(int idBeamIn, string setName, int member,
  Logger* loggerPtr) {
  idBeam = idBeamIn;
  idBeamAbs = abs(idBeamIn);
  isPhoton = idBeamIn == 22;
  nSet = LHAPDF5Interface::findNSet(setName, member);
  if (nSet == -1) nSet = LHAPDF5Interface::freeNSet();

  // If already initialized then need not do anything further.
  LHAPDF5Interface::LHAPDFInfo initializedInfo =
    LHAPDF5Interface::initializedSets[nSet];
  string initializedSetName   = initializedInfo.name;
  int    initializedMember    = initializedInfo.member;
  hasPhoton                   = initializedInfo.photon;
  if (setName == initializedSetName && member == initializedMember)
    return true;

  // Initialize set. If first character is '/' then assume that name
  // is given with path, else not.
  if (setName[0] == '/') LHAPDF5Interface::initPDFsetM( nSet, setName);
  else LHAPDF5Interface::initPDFsetByNameM( nSet, setName);
  isSet = (nSet >= 0);

  // Initialize member.
  LHAPDF5Interface::initPDFM(nSet, member);

  // Do not collect statistics on under/overflow to save time and space.
  LHAPDF5Interface::setPDFparm( "NOSTAT" );
  LHAPDF5Interface::setPDFparm( "LOWKEY" );

  // Check if photon PDF available (has_photon does not work properly).
  xPhoton = 0;
  LHAPDF5Interface::evolvePDFPHOTONM(nSet, 0.01, 1, xfArray, xPhoton);
  hasPhoton = xPhoton != 0;

  // Save values to avoid unnecessary reinitializations.
  initializedInfo.name   = setName;
  initializedInfo.member = member;
  initializedInfo.photon = hasPhoton;
  if (nSet > 0) LHAPDF5Interface::initializedSets[nSet] = initializedInfo;
  return true;

}

//--------------------------------------------------------------------------

// Allow optional extrapolation beyond boundaries.

void LHAPDF5::setExtrapolate(bool extrapol) {

  doExtraPol = extrapol;
  LHAPDF5Interface::setPDFparm( (extrapol) ? "EXTRAPOLATE" : "18" );

}

//--------------------------------------------------------------------------

// Give the parton distribution function set from LHAPDF5.

void LHAPDF5::xfUpdate(int, double x, double Q2) {

  // Freeze at boundary value if PDF is evaluated outside the fit region.
  int member = LHAPDF5Interface::initializedSets[nSet].member;
  double xMin, xMax, q2Min, q2Max;
  getxminm_( nSet, member, xMin);
  getxmaxm_( nSet, member, xMax);
  getq2minm_(nSet, member, q2Min);
  getq2maxm_(nSet, member, q2Max);
  if (x < xMin && !doExtraPol) x = xMin;
  if (x > xMax)    x = xMax;
  if (Q2 < q2Min) Q2 = q2Min;
  if (Q2 > q2Max) Q2 = q2Max;
  double Q = sqrt( max( 0., Q2));

  // Use special call if photon included in proton.
  if (hasPhoton) {
    LHAPDF5Interface::evolvePDFPHOTONM( nSet, x, Q, xfArray, xPhoton);
  }

  // Use special call with photon beams. No virtualities implemented yet.
  else if (isPhoton) {
    LHAPDF5Interface::evolvePDFpM( nSet, x, Q, 0., 0., xfArray);
  }

  // Else use default LHAPDF5 call.
  else {
    LHAPDF5Interface::evolvePDFM( nSet, x, Q, xfArray);
    xPhoton=0.0;
  }

  // Update values.
  xg     = xfArray[6];
  xd     = xfArray[7];
  xdbar  = xfArray[5];
  xu     = xfArray[8];
  xubar  = xfArray[4];
  xs     = xfArray[9];
  xc     = xfArray[10];
  xb     = xfArray[11];
  xsbar  = sSymmetricSave ? xs : xfArray[3];
  xcbar  = cSymmetricSave ? xc : xfArray[2];
  xbbar  = bSymmetricSave ? xb : xfArray[1];
  xgamma = xPhoton;

  // idSav = 9 to indicate that all flavours reset.
  idSav = 9;

}

//--------------------------------------------------------------------------

// Declare the plugin.

PYTHIA8_PLUGIN_CLASS(PDF, LHAPDF5, false, false, false)
PYTHIA8_PLUGIN_VERSIONS(PYTHIA_VERSION_INTEGER)

//==========================================================================

} // end namespace Pythia8

#endif // end Pythia8_LHAPDF5_H