PYTHIA  8.311
LHAPDF6.h
1 // LHAPDF6.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 LHAPDF6 PDF plugin class.
7 
8 #ifndef Pythia8_LHAPDF6_H
9 #define Pythia8_LHAPDF6_H
10 
11 #include "Pythia8/PartonDistributions.h"
12 #include "Pythia8/Plugins.h"
13 #include "LHAPDF/LHAPDF.h"
14 #include <mutex>
15 
16 namespace Pythia8 {
17 
18 //==========================================================================
19 
20 // Containers for PDF sets.
21 
22 //--------------------------------------------------------------------------
23 
24 // Class to hold a PDF set, its information, and its uncertainty sets.
25 
26 class PdfSets {
27 
28 public:
29 
30  // Constructors.
31  PdfSets() {;}
32  PdfSets(string setName) : info(::LHAPDF::PDFSet(setName)),
33  pdfs(vector< ::LHAPDF::PDF* >(info.size(), 0)) {;}
34 
35  // Access a PDF set.
36  ::LHAPDF::PDF *operator[](unsigned int member) {
37  if (!pdfs[member]) {
38  lock_guard<mutex> lck (mtx);
39  pdfs[member] = info.mkPDF(member);
40  }
41  return pdfs[member];
42  }
43 
44  // Get number of PDF sets.
45  int size() {return pdfs.size();}
46 
47  // PDF sets and info.
48  ::LHAPDF::PDFSet info;
49  vector< ::LHAPDF::PDF* > pdfs;
50  static mutex mtx;
51 
52 };
53 
54 mutex PdfSets::mtx;
55 
56 //==========================================================================
57 
58 // Provide interface to the LHAPDF6 library of parton densities.
59 
60 class LHAPDF6 : public PDF {
61 
62 public:
63 
64  // Constructor.
65  LHAPDF6(Pythia*, Settings* settingsPtr, Logger*) :
66  PDF(), pdf(nullptr), extrapol(false) {
67  if (settingsPtr == nullptr) return;
68  sSymmetric(settingsPtr->flag("LHAPDF:sSymmetric"));
69  cSymmetric(settingsPtr->flag("LHAPDF:cSymmetric"));
70  bSymmetric(settingsPtr->flag("LHAPDF:bSymmetric"));
71  }
72 
73  // Initialization of PDF set.
74  bool init(int idBeamIn, string setName, int member, Logger* loggerPtr)
75  override;
76 
77  // Allow extrapolation beyond boundaries (not implemented).
78  void setExtrapolate(bool extrapolIn) {extrapol = extrapolIn;}
79 
80 private:
81 
82  // The LHAPDF objects.
83  PdfSets pdfs;
84  ::LHAPDF::PDF *pdf;
85  ::LHAPDF::Extrapolator *ext;
86  bool extrapol;
87 
88  // Update parton densities.
89  void xfUpdate(int id, double x, double Q2);
90 
91  // Check whether x and Q2 values fall inside the fit bounds.
92  bool insideBounds(double x, double Q2) {
93  return (x > xMin && x < xMax && Q2 > q2Min && Q2 < q2Max);}
94 
95  // Return the running alpha_s shipped with the LHAPDF set.
96  double alphaS(double Q2) { return pdf->alphasQ2(Q2); }
97 
98  // Return quark masses used in the PDF fit.
99  double muPDFSave, mdPDFSave, mcPDFSave, msPDFSave, mbPDFSave,
100  xMin, xMax, q2Min, q2Max;
101  double mQuarkPDF(int id) {
102  switch(abs(id)){
103  case 1: return mdPDFSave;
104  case 2: return muPDFSave;
105  case 3: return msPDFSave;
106  case 4: return mcPDFSave;
107  case 5: return mbPDFSave;
108  }
109  return -1.;
110  }
111 
112  // Calculate uncertainties using the LHAPDF prescription.
113  void calcPDFEnvelope(int, double, double, int);
114  void calcPDFEnvelope(pair<int,int>, pair<double,double>, double, int);
115  PDFEnvelope pdfEnvelope;
116  PDFEnvelope getPDFEnvelope() {return pdfEnvelope;}
117  static const double PDFMINVALUE;
118 
119  int nMembersSave;
120  int nMembers() { return nMembersSave; }
121 
122 };
123 
124 //--------------------------------------------------------------------------
125 
126 // Constants.
127 
128 const double LHAPDF6::PDFMINVALUE = 1e-10;
129 
130 //--------------------------------------------------------------------------
131 
132 // Initialize a parton density function from LHAPDF6.
133 
134 bool LHAPDF6::init(int idBeamIn, string setName, int member,
135  Logger* loggerPtr) {
136  idBeam = idBeamIn;
137  idBeamAbs = abs(idBeamIn);
138  isSet = false;
139 
140  // Find the PDF set. Note, LHAPDF aborts if the PDF does not exist,
141  // which we avoid with this try/catch statement. Ideally, we would
142  // check with :LHAPDF::lookupLHAPDFID, but this is not thread safe.
143  try {
144  pdfs = PdfSets(setName);
145  } catch (const std::exception &e) {
146  loggerPtr->ERROR_MSG("unknown PDF " + setName);
147  return false;
148  }
149 
150  // Find the PDF member.
151  if (pdfs.size() == 0) {
152  loggerPtr->ERROR_MSG("could not initialize PDF " + setName);
153  return false;
154  } else if (member >= pdfs.size()) {
155  loggerPtr->ERROR_MSG(setName + " does not contain requested member");
156  return false;
157  }
158  pdf = pdfs[member];
159  isSet = true;
160 
161  // Save x and Q2 limits.
162  xMax = pdf->xMax();
163  xMin = pdf->xMin();
164  q2Max = pdf->q2Max();
165  q2Min = pdf->q2Min();
166 
167  // Store quark masses used in PDF fit.
168  muPDFSave = pdf->info().get_entry_as<double>("MUp");
169  mdPDFSave = pdf->info().get_entry_as<double>("MDown");
170  mcPDFSave = pdf->info().get_entry_as<double>("MCharm");
171  msPDFSave = pdf->info().get_entry_as<double>("MStrange");
172  mbPDFSave = pdf->info().get_entry_as<double>("MBottom");
173  nMembersSave = pdf->info().get_entry_as<int>("NumMembers");
174  return true;
175 
176 }
177 
178 //--------------------------------------------------------------------------
179 
180 // Give the parton distribution function set from LHAPDF6.
181 
182 void LHAPDF6::xfUpdate(int, double x, double Q2) {
183  if (!isSet) return;
184 
185  // Freeze at boundary value if PDF is evaluated outside the fit region.
186  if (x < xMin && !extrapol) x = xMin;
187  if (x > xMax) x = xMax;
188  if (Q2 < q2Min) Q2 = q2Min;
189  if (Q2 > q2Max) Q2 = q2Max;
190 
191  // Update values.
192  xg = pdf->xfxQ2(21, x, Q2);
193  xd = pdf->xfxQ2(1, x, Q2);
194  xu = pdf->xfxQ2(2, x, Q2);
195  xdbar = pdf->xfxQ2(-1, x, Q2);
196  xubar = pdf->xfxQ2(-2, x, Q2);
197  xs = pdf->xfxQ2(3, x, Q2);
198  xc = pdf->xfxQ2(4, x, Q2);
199  xb = pdf->xfxQ2(5, x, Q2);
200  xsbar = sSymmetricSave ? xs : pdf->xfxQ2(-3, x, Q2);
201  xcbar = cSymmetricSave ? xc : pdf->xfxQ2(-4, x, Q2);
202  xbbar = bSymmetricSave ? xb : pdf->xfxQ2(-5, x, Q2);
203  xgamma = pdf->xfxQ2(22, x, Q2);
204 
205  // idSav = 9 to indicate that all flavours reset.
206  idSav = 9;
207 
208 }
209 
210 //--------------------------------------------------------------------------
211 
212 // Calculate uncertainties using the LHAPDF prescription.
213 
214 void LHAPDF6::calcPDFEnvelope(int idNow, double xNow, double Q2NowIn,
215  int valSea) {
216  if (!isSet) return;
217 
218  // Freeze at boundary value if PDF is evaluated outside the fit region.
219  double x1 = (xNow < xMin && !extrapol) ? xMin : xNow;
220  if (x1 > xMax) x1 = xMax;
221  double Q2Now = (Q2NowIn < q2Min) ? q2Min : Q2NowIn;
222  if (Q2Now > q2Max) Q2Now = q2Max;
223 
224  // Loop over the members.
225  vector<double> xfCalc(pdfs.size());
226  for(int iMem = 0; iMem < pdfs.size(); ++iMem) {
227  if (valSea==0 || (idNow != 1 && idNow != 2)) {
228  xfCalc[iMem] = pdfs[iMem]->xfxQ2(idNow, x1, Q2Now);
229  } else if (valSea==1 && (idNow == 1 || idNow == 2 )) {
230  xfCalc[iMem] = pdfs[iMem]->xfxQ2(idNow, x1, Q2Now) -
231  pdfs[iMem]->xfxQ2(-idNow, x1, Q2Now);
232  } else if (valSea==2 && (idNow == 1 || idNow == 2 )) {
233  xfCalc[iMem] = pdfs[iMem]->xfxQ2(-idNow, x1, Q2Now);
234  }
235  }
236 
237  // Calculate the uncertainty.
238  ::LHAPDF::PDFUncertainty xfErr = pdfs.info.uncertainty(xfCalc);
239  pdfEnvelope.centralPDF = xfErr.central;
240  pdfEnvelope.errplusPDF = xfErr.errplus;
241  pdfEnvelope.errminusPDF = xfErr.errminus;
242  pdfEnvelope.errsymmPDF = xfErr.errsymm;
243  pdfEnvelope.scalePDF = xfErr.scale;
244 }
245 
246 //--------------------------------------------------------------------------
247 
248 // Calculate uncertainties using the LHAPDF prescription.
249 
250 void LHAPDF6::calcPDFEnvelope(pair<int,int> idNows, pair<double,double> xNows,
251  double Q2NowIn, int valSea) {
252  if (!isSet) return;
253 
254  // Freeze at boundary value if PDF is evaluated outside the fit region.
255  double x1 = (xNows.first < xMin && !extrapol) ? xMin : xNows.first;
256  if (x1 > xMax) x1 = xMax;
257  double x2 = (xNows.second < xMin && !extrapol) ? xMin : xNows.second;
258  if (x2 > xMax) x2 = xMax;
259  double Q2Now = (Q2NowIn < q2Min) ? q2Min : Q2NowIn;
260  if (Q2Now > q2Max) Q2Now = q2Max;
261 
262  // Loop over the members.
263  vector<double> xfCalc(pdfs.size());
264  pdfEnvelope.pdfMemberVars.resize(pdfs.size());
265  for(int iMem = 0; iMem < pdfs.size(); ++iMem) {
266  if (valSea == 0 || (idNows.first != 1 && idNows.first != 2 ) ) {
267  xfCalc[iMem] = pdfs[iMem]->xfxQ2(idNows.first, x1, Q2Now);
268  } else if (valSea == 1 && (idNows.first == 1 || idNows.first == 2)) {
269  xfCalc[iMem] = pdfs[iMem]->xfxQ2(idNows.first, x1, Q2Now)
270  - pdfs[iMem]->xfxQ2(-idNows.first, x1, Q2Now);
271  } else if (valSea == 2 && (idNows.first == 1 || idNows.first == 2 )) {
272  xfCalc[iMem] = pdfs[iMem]->xfxQ2(-idNows.first, x1, Q2Now);
273  }
274  xfCalc[iMem] = max(0.0, xfCalc[iMem]);
275  if (valSea == 0 || (idNows.second != 1 && idNows.second != 2)) {
276  xfCalc[iMem] /= max
277  (PDFMINVALUE, pdfs[iMem]->xfxQ2(idNows.second, x2, Q2Now));
278  } else if (valSea == 1 && (idNows.second == 1 || idNows.second == 2 )) {
279  xfCalc[iMem] /= max
280  (pdfs[iMem]->xfxQ2(idNows.second, x2, Q2Now) - pdfs[iMem]->xfxQ2
281  (-idNows.second, x2, Q2Now), PDFMINVALUE);
282  } else if (valSea == 2 && (idNows.second == 1 || idNows.second == 2 )) {
283  xfCalc[iMem] /= max
284  (pdfs[iMem]->xfxQ2(-idNows.second, x2, Q2Now), PDFMINVALUE);
285  }
286  pdfEnvelope.pdfMemberVars[iMem] = xfCalc[iMem];
287  }
288 
289  // Calculate the uncertainty.
290  ::LHAPDF::PDFUncertainty xfErr = pdfs.info.uncertainty(xfCalc);
291  pdfEnvelope.centralPDF = xfErr.central;
292  pdfEnvelope.errplusPDF = xfErr.errplus;
293  pdfEnvelope.errminusPDF = xfErr.errminus;
294  pdfEnvelope.errsymmPDF = xfErr.errsymm;
295  pdfEnvelope.scalePDF = xfErr.scale;
296 
297 }
298 
299 //--------------------------------------------------------------------------
300 
301 // Declare the plugin.
302 
303 PYTHIA8_PLUGIN_CLASS(PDF, LHAPDF6, false, false, false)
304 PYTHIA8_PLUGIN_VERSIONS(PYTHIA_VERSION_INTEGER)
305 
306 //==========================================================================
307 
308 } // end namespace Pythia8
309 
310 #endif // end Pythia8_LHAPDF6_H
Containers for PDF sets.
Definition: LHAPDF6.h:26
void setExtrapolate(bool extrapolIn)
Allow extrapolation beyond boundaries (not implemented).
Definition: LHAPDF6.h:78
bool flag(string keyIn)
Give back current value, with check that key exists.
Definition: Settings.cc:1473
::LHAPDF::PDFSet info
PDF sets and info.
Definition: LHAPDF6.h:48
LHAPDF6(Pythia *, Settings *settingsPtr, Logger *)
Constructor.
Definition: LHAPDF6.h:65
Provide interface to the LHAPDF6 library of parton densities.
Definition: LHAPDF6.h:60
Definition: Logger.h:23
Error envelope from PDF uncertainty.
Definition: PartonDistributions.h:102
PdfSets()
Constructors.
Definition: LHAPDF6.h:31
::LHAPDF::PDF * operator[](unsigned int member)
Access a PDF set.
Definition: LHAPDF6.h:36
int size()
Get number of PDF sets.
Definition: LHAPDF6.h:45
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
bool init(int idBeamIn, string setName, int member, Logger *loggerPtr) override
Initialization of PDF set.
Definition: LHAPDF6.h:134
Definition: Settings.h:195