PYTHIA  8.312
LHAPDF5.h
1 // LHAPDF5.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2024 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 
6 // This file contains the LHAPDF5 PDF plugin class.
7 
8 #ifndef Pythia8_LHAPDF5_H
9 #define Pythia8_LHAPDF5_H
10 
11 #include "Pythia8/PartonDistributions.h"
12 #include "Pythia8/Plugins.h"
13 
14 namespace Pythia8 {
15 
16 //==========================================================================
17 
18 // Interfaces to to make the C++ calls similar to f77.
19 
20 //--------------------------------------------------------------------------
21 
22 // Declare the LHAPDF5 f77 subroutines that are needed.
23 
24 extern "C" {
25 
26  extern void initpdfsetm_(int&, const char*, int);
27 
28  extern void initpdfsetbynamem_(int&, const char*, int);
29 
30  extern void initpdfm_(int&, int&);
31 
32  extern void evolvepdfm_(int&, double&, double&, double*);
33 
34  extern void evolvepdfpm_(int&, double&, double&, double&, double&, double*);
35 
36  extern void evolvepdfphotonm_(int&, double&, double&, double*, double&);
37 
38  extern void setlhaparm_(const char*, int);
39 
40  extern void getxminm_(int &, int &, double &);
41 
42  extern void getxmaxm_(int &, int &, double &);
43 
44  extern void getq2minm_(int &, int &, double &);
45 
46  extern void getq2maxm_(int &, int &, double &);
47 
48 }
49 
50 //--------------------------------------------------------------------------
51 
52 // Map the f77 routines to C++.
53 
54 namespace LHAPDF5Interface {
55 
56  // Initialize set with full pathname, allowing multiple sets.
57  void initPDFsetM( int& nSet, string name) {
58  const char* cName = name.c_str(); int lenName = name.length();
59  initpdfsetm_( nSet, cName, lenName);
60  }
61 
62  // Initialize set with simple name, allowing multiple sets.
63  void initPDFsetByNameM( int& nSet, string name) {
64  const char* cName = name.c_str(); int lenName = name.length();
65  initpdfsetbynamem_( nSet, cName, lenName);
66  }
67 
68  // Initialize member of set.
69  void initPDFM(int& nSet, int member) {
70  initpdfm_(nSet, member);
71  }
72 
73  // Evaluate x f_i(x, Q).
74  void evolvePDFM( int& nSet, double x, double Q, double* xfArray) {
75  evolvepdfm_( nSet, x, Q, xfArray);
76  }
77 
78  // Evaluate x f_i(x, Q) for photon beams.
79  void evolvePDFpM( int& nSet, double x, double Q, double P2, double IP2,
80  double* xfArray) {
81  evolvepdfpm_( nSet, x, Q, P2, IP2, xfArray);
82  }
83 
84  // Evaluate x f_i(x, Q) including photon
85  void evolvePDFPHOTONM(int& nSet, double x, double Q, double* xfArray,
86  double& xPhoton) {
87  evolvepdfphotonm_( nSet, x, Q, xfArray, xPhoton);
88  }
89 
90  // Extrapolate PDF set beyond boundaries, or freeze them there.
91  void setPDFparm(string name) {
92  const char* cName = name.c_str(); int lenName = name.length();
93  setlhaparm_( cName, lenName);
94  }
95 
96  // Simple structure to hold LHAPDF set information.
97  struct LHAPDFInfo {
98  string name;
99  int member;
100  bool photon;
101  };
102 
103  // Global tracking of opened PDF sets.
104  map<int, LHAPDFInfo> initializedSets;
105 
106  // Method to find the nSet number corresponding to a name and member.
107  // Returns -1 if no such LHAPDF5 set has been initialized.
108  int findNSet(string setName, int member) {
109  for (map<int, LHAPDFInfo>::const_iterator i = initializedSets.begin();
110  i != initializedSets.end(); ++i) {
111  int iSet = i->first;
112  string iName = i->second.name;
113  int iMember = i->second.member;
114  if (iName == setName && iMember == member) return iSet;
115  }
116  return -1;
117  }
118 
119  // Method to return the lowest non-occupied nSet number.
120  int freeNSet() {
121  for (int iSet = 1; iSet <= int(initializedSets.size()); ++iSet) {
122  if (initializedSets.find(iSet) == initializedSets.end())
123  return iSet;
124  }
125  return initializedSets.size() + 1;
126  }
127 
128 }
129 
130 //==========================================================================
131 
132 // Plugin interface to the LHAPDF5 library.
133 
134 //==========================================================================
135 
136 // Provide plugin interface to the LHAPDF5 library of parton densities.
137 
138 class LHAPDF5 : public PDF {
139 
140 public:
141 
142  // Constructor.
143  LHAPDF5(Pythia*, Settings* settingsPtr, Logger*) :
144  PDF(), doExtraPol(false) {
145  if (settingsPtr == nullptr) return;
146  sSymmetric(settingsPtr->flag("LHAPDF:sSymmetric"));
147  cSymmetric(settingsPtr->flag("LHAPDF:cSymmetric"));
148  bSymmetric(settingsPtr->flag("LHAPDF:bSymmetric"));
149  }
150 
151  // Initialization of PDF set.
152  bool init(int idBeamIn, string setName, int member, Logger* loggerPtr)
153  override;
154 
155  // Allow extrapolation beyond boundaries. This is optional.
156  void setExtrapolate(bool extrapol);
157 
158 private:
159 
160  // Update all PDF values.
161  void xfUpdate(int , double x, double Q2);
162 
163  // Current set and pdf values.
164  bool doExtraPol;
165  int nSet;
166  double xfArray[13];
167  bool hasPhoton, isPhoton;
168  double xPhoton;
169 
170 };
171 
172 //--------------------------------------------------------------------------
173 
174 // Initialize a parton density function from LHAPDF5.
175 
176 bool LHAPDF5::init(int idBeamIn, string setName, int member,
177  Logger* loggerPtr) {
178  idBeam = idBeamIn;
179  idBeamAbs = abs(idBeamIn);
180  isPhoton = idBeamIn == 22;
181  nSet = LHAPDF5Interface::findNSet(setName, member);
182  if (nSet == -1) nSet = LHAPDF5Interface::freeNSet();
183 
184  // If already initialized then need not do anything further.
185  LHAPDF5Interface::LHAPDFInfo initializedInfo =
187  string initializedSetName = initializedInfo.name;
188  int initializedMember = initializedInfo.member;
189  hasPhoton = initializedInfo.photon;
190  if (setName == initializedSetName && member == initializedMember)
191  return true;
192 
193  // Initialize set. If first character is '/' then assume that name
194  // is given with path, else not.
195  if (setName[0] == '/') LHAPDF5Interface::initPDFsetM( nSet, setName);
196  else LHAPDF5Interface::initPDFsetByNameM( nSet, setName);
197  isSet = (nSet >= 0);
198 
199  // Initialize member.
200  LHAPDF5Interface::initPDFM(nSet, member);
201 
202  // Do not collect statistics on under/overflow to save time and space.
203  LHAPDF5Interface::setPDFparm( "NOSTAT" );
204  LHAPDF5Interface::setPDFparm( "LOWKEY" );
205 
206  // Check if photon PDF available (has_photon does not work properly).
207  xPhoton = 0;
208  LHAPDF5Interface::evolvePDFPHOTONM(nSet, 0.01, 1, xfArray, xPhoton);
209  hasPhoton = xPhoton != 0;
210 
211  // Save values to avoid unnecessary reinitializations.
212  initializedInfo.name = setName;
213  initializedInfo.member = member;
214  initializedInfo.photon = hasPhoton;
215  if (nSet > 0) LHAPDF5Interface::initializedSets[nSet] = initializedInfo;
216  return true;
217 
218 }
219 
220 //--------------------------------------------------------------------------
221 
222 // Allow optional extrapolation beyond boundaries.
223 
224 void LHAPDF5::setExtrapolate(bool extrapol) {
225 
226  doExtraPol = extrapol;
227  LHAPDF5Interface::setPDFparm( (extrapol) ? "EXTRAPOLATE" : "18" );
228 
229 }
230 
231 //--------------------------------------------------------------------------
232 
233 // Give the parton distribution function set from LHAPDF5.
234 
235 void LHAPDF5::xfUpdate(int, double x, double Q2) {
236 
237  // Freeze at boundary value if PDF is evaluated outside the fit region.
238  int member = LHAPDF5Interface::initializedSets[nSet].member;
239  double xMin, xMax, q2Min, q2Max;
240  getxminm_( nSet, member, xMin);
241  getxmaxm_( nSet, member, xMax);
242  getq2minm_(nSet, member, q2Min);
243  getq2maxm_(nSet, member, q2Max);
244  if (x < xMin && !doExtraPol) x = xMin;
245  if (x > xMax) x = xMax;
246  if (Q2 < q2Min) Q2 = q2Min;
247  if (Q2 > q2Max) Q2 = q2Max;
248  double Q = sqrt( max( 0., Q2));
249 
250  // Use special call if photon included in proton.
251  if (hasPhoton) {
252  LHAPDF5Interface::evolvePDFPHOTONM( nSet, x, Q, xfArray, xPhoton);
253  }
254 
255  // Use special call with photon beams. No virtualities implemented yet.
256  else if (isPhoton) {
257  LHAPDF5Interface::evolvePDFpM( nSet, x, Q, 0., 0., xfArray);
258  }
259 
260  // Else use default LHAPDF5 call.
261  else {
262  LHAPDF5Interface::evolvePDFM( nSet, x, Q, xfArray);
263  xPhoton=0.0;
264  }
265 
266  // Update values.
267  xg = xfArray[6];
268  xd = xfArray[7];
269  xdbar = xfArray[5];
270  xu = xfArray[8];
271  xubar = xfArray[4];
272  xs = xfArray[9];
273  xc = xfArray[10];
274  xb = xfArray[11];
275  xsbar = sSymmetricSave ? xs : xfArray[3];
276  xcbar = cSymmetricSave ? xc : xfArray[2];
277  xbbar = bSymmetricSave ? xb : xfArray[1];
278  xgamma = xPhoton;
279 
280  // idSav = 9 to indicate that all flavours reset.
281  idSav = 9;
282 
283 }
284 
285 //--------------------------------------------------------------------------
286 
287 // Declare the plugin.
288 
289 PYTHIA8_PLUGIN_CLASS(PDF, LHAPDF5, false, false, false)
290 PYTHIA8_PLUGIN_VERSIONS(PYTHIA_VERSION_INTEGER)
291 
292 //==========================================================================
293 
294 } // end namespace Pythia8
295 
296 #endif // end Pythia8_LHAPDF5_H
bool flag(string keyIn)
Give back current value, with check that key exists.
Definition: Settings.cc:1478
Plugin interface to the LHAPDF5 library.
Definition: LHAPDF5.h:138
void initPDFsetByNameM(int &nSet, string name)
Initialize set with simple name, allowing multiple sets.
Definition: LHAPDF5.h:63
void initPDFsetM(int &nSet, string name)
Initialize set with full pathname, allowing multiple sets.
Definition: LHAPDF5.h:57
map< int, LHAPDFInfo > initializedSets
Global tracking of opened PDF sets.
Definition: LHAPDF5.h:104
void initPDFM(int &nSet, int member)
Initialize member of set.
Definition: LHAPDF5.h:69
void evolvePDFPHOTONM(int &nSet, double x, double Q, double *xfArray, double &xPhoton)
Evaluate x f_i(x, Q) including photon.
Definition: LHAPDF5.h:85
bool init(int idBeamIn, string setName, int member, Logger *loggerPtr) override
Initialization of PDF set.
Definition: LHAPDF5.h:176
Definition: Logger.h:23
int freeNSet()
Method to return the lowest non-occupied nSet number.
Definition: LHAPDF5.h:120
void setExtrapolate(bool extrapol)
Allow extrapolation beyond boundaries. This is optional.
Definition: LHAPDF5.h:224
void evolvePDFM(int &nSet, double x, double Q, double *xfArray)
Evaluate x f_i(x, Q).
Definition: LHAPDF5.h:74
void setPDFparm(string name)
Extrapolate PDF set beyond boundaries, or freeze them there.
Definition: LHAPDF5.h:91
void initpdfsetm_(int &, const char *, int)
Interfaces to to make the C++ calls similar to f77.
Base class for parton distribution functions.
Definition: PartonDistributions.h:49
Header for classes to set beam momentum and interaction vertex spread.
Definition: Analysis.h:20
The Pythia class contains the top-level routines to generate an event.
Definition: Pythia.h:71
Simple structure to hold LHAPDF set information.
Definition: LHAPDF5.h:97
LHAPDF5(Pythia *, Settings *settingsPtr, Logger *)
Constructor.
Definition: LHAPDF5.h:143
int findNSet(string setName, int member)
Definition: LHAPDF5.h:108
void evolvePDFpM(int &nSet, double x, double Q, double P2, double IP2, double *xfArray)
Evaluate x f_i(x, Q) for photon beams.
Definition: LHAPDF5.h:79
Definition: Settings.h:195