StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
LHAPDF6.h
1 // LHAPDF6.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2018 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 "LHAPDF/LHAPDF.h"
13 
14 namespace Pythia8 {
15 
16 //==========================================================================
17 
18 // Global tracking of opened PDF sets.
19 
20 //--------------------------------------------------------------------------
21 
22 namespace LHAPDF6Interface {
23 
24 //--------------------------------------------------------------------------
25 
26 // Class to hold a PDF set, its information, and its uncertainty sets.
27 
28  class PdfSets {
29 
30  public:
31 
32  // Constructors.
33  PdfSets() {;}
34  PdfSets(string setName) : info(::LHAPDF::PDFSet(setName)),
35  pdfs(vector< ::LHAPDF::PDF* >(info.size(), 0)) {;}
36 
37  // Access a PDF set.
38  ::LHAPDF::PDF *operator[](unsigned int member) {
39  if (!pdfs[member]) pdfs[member] = info.mkPDF(member);
40  return pdfs[member];
41  }
42 
43  // Get number of PDF sets.
44  int size() {return pdfs.size();}
45 
46  // PDF sets and info.
47  ::LHAPDF::PDFSet info;
48  vector< ::LHAPDF::PDF* > pdfs;
49 
50  };
51 
52 //--------------------------------------------------------------------------
53 
54 // Class to globally track all open PDF sets.
55 
56  class PdfTracker {
57 
58  public:
59 
60  // Destructor, all PDFs from LHAPDF::mkPDF must be deleted.
61  ~PdfTracker() {
62  for (map<int, PdfSets>::iterator pdf = pdfs.begin();
63  pdf != pdfs.end(); ++pdf)
64  for (int iMem = 0; iMem < (int)pdf->second.size(); ++iMem)
65  if (pdf->second.pdfs[iMem]) delete pdf->second.pdfs[iMem];
66  }
67 
68  // Find and return a requested PDF set.
69  PdfSets *find(string setName) {
70  int id = ::LHAPDF::lookupLHAPDFID(setName, 0);
71  if (id < 0) return 0;
72  else if (pdfs.find(id) == pdfs.end()) pdfs[id] = PdfSets(setName);
73  return &pdfs[id];
74  }
75 
76  private:
77 
78  // Map to hold open PDF set information.
79  map<int, PdfSets> pdfs;
80 
81  };
82 
83 //--------------------------------------------------------------------------
84 
85 // Define opened PDF sets global variable.
86 
87  PdfTracker pdfTracker;
88 
89 }
90 
91 //==========================================================================
92 
93 // Provide interface to the LHAPDF6 library of parton densities.
94 
95 class LHAPDF6 : public PDF {
96 
97 public:
98 
99  // Constructor.
100  LHAPDF6(int idBeamIn, string setName, int member, int, Info* infoPtr)
101  : PDF(idBeamIn), pdf(0), extrapol(false)
102  { init(setName, member, infoPtr); }
103 
104  // Allow extrapolation beyond boundaries (not implemented).
105  void setExtrapolate(bool extrapolIn) {extrapol = extrapolIn;}
106 
107 private:
108 
109  // The LHAPDF objects.
110  LHAPDF6Interface::PdfSets *pdfs;
111  ::LHAPDF::PDF *pdf;
112  ::LHAPDF::Extrapolator *ext;
113  bool extrapol;
114 
115  // Initialization of PDF set.
116  void init(string setName, int member, Info* infoPtr);
117 
118  // Update parton densities.
119  void xfUpdate(int id, double x, double Q2);
120 
121  // Check whether x and Q2 values fall inside the fit bounds.
122  bool insideBounds(double x, double Q2) {
123  return (x > pdf->xMin() && x < pdf->xMax()
124  && Q2 > pdf->q2Min() && Q2 < pdf->q2Max());}
125 
126  // Return the running alpha_s shipped with the LHAPDF set.
127  double alphaS(double Q2) { return pdf->alphasQ2(Q2); }
128 
129  // Return quark masses used in the PDF fit.
130  double muPDFSave, mdPDFSave, mcPDFSave, msPDFSave, mbPDFSave;
131  double mQuarkPDF(int id) {
132  if (abs(id) == 1) return mdPDFSave;
133  if (abs(id) == 2) return muPDFSave;
134  if (abs(id) == 3) return msPDFSave;
135  if (abs(id) == 4) return mcPDFSave;
136  if (abs(id) == 5) return mbPDFSave;
137  return -1.;
138  }
139 
140  // Calculate uncertainties using the LHAPDF prescription.
141  void calcPDFEnvelope(int, double, double, int);
142  void calcPDFEnvelope(pair<int,int>, pair<double,double>, double, int);
143  PDFEnvelope pdfEnvelope;
144  PDFEnvelope getPDFEnvelope() {return pdfEnvelope;}
145  static const double PDFMINVALUE;
146 
147  int nMembersSave;
148  int nMembers() { return nMembersSave; }
149 
150 };
151 
152 //--------------------------------------------------------------------------
153 
154 // Constants.
155 
156 const double LHAPDF6::PDFMINVALUE = 1e-10;
157 
158 //--------------------------------------------------------------------------
159 
160 // Initialize a parton density function from LHAPDF6.
161 
162 void LHAPDF6::init(string setName, int member, Info *info) {
163  isSet = false;
164 
165  // Initialize the LHAPDF sets.
166  pdfs = LHAPDF6Interface::pdfTracker.find(setName);
167  if (!pdfs) {
168  info->errorMsg("Error in LHAPDF6::init: unknown PDF " + setName);
169  return;
170  } else if ((*pdfs).size() == 0) {
171  info->errorMsg("Error in LHAPDF6::init: could not initialize PDF "
172  + setName);
173  return;
174  } else if (member >= (*pdfs).size()) {
175  info->errorMsg("Error in LHAPDF6::init: " + setName
176  + " does not contain requested member");
177  return;
178  }
179  pdf = (*pdfs)[member];
180  isSet = true;
181 
182  // Store quark masses used in PDF fit.
183  muPDFSave = pdf->info().get_entry_as<double>("MUp");
184  mdPDFSave = pdf->info().get_entry_as<double>("MDown");
185  mcPDFSave = pdf->info().get_entry_as<double>("MCharm");
186  msPDFSave = pdf->info().get_entry_as<double>("MStrange");
187  mbPDFSave = pdf->info().get_entry_as<double>("MBottom");
188 
189  nMembersSave = pdf->info().get_entry_as<int>("NumMembers");
190 
191 }
192 
193 //--------------------------------------------------------------------------
194 
195 // Give the parton distribution function set from LHAPDF6.
196 
197 void LHAPDF6::xfUpdate(int, double x, double Q2) {
198 
199  // Freeze at boundary value if PDF is evaluated outside the fit region.
200  if (x < pdf->xMin() && !extrapol) x = pdf->xMin();
201  if (x > pdf->xMax() ) x = pdf->xMax();
202  if (Q2 < pdf->q2Min() ) Q2 = pdf->q2Min();
203  if (Q2 > pdf->q2Max() ) Q2 = pdf->q2Max();
204 
205  // Update values.
206  xg = pdf->xfxQ2(21, x, Q2);
207  xu = pdf->xfxQ2(2, x, Q2);
208  xd = pdf->xfxQ2(1, x, Q2);
209  xs = pdf->xfxQ2(3, x, Q2);
210  xubar = pdf->xfxQ2(-2, x, Q2);
211  xdbar = pdf->xfxQ2(-1, x, Q2);
212  xsbar = pdf->xfxQ2(-3, x, Q2);
213  xc = pdf->xfxQ2(4, x, Q2);
214  xb = pdf->xfxQ2(5, x, Q2);
215  xgamma = pdf->xfxQ2(22, x, Q2);
216 
217  // Subdivision of valence and sea.
218  xuVal = xu - xubar;
219  xuSea = xubar;
220  xdVal = xd - xdbar;
221  xdSea = xdbar;
222 
223  // idSav = 9 to indicate that all flavours reset.
224  idSav = 9;
225 
226 }
227 
228 //--------------------------------------------------------------------------
229 
230 // Calculate uncertainties using the LHAPDF prescription.
231 
232 void LHAPDF6::calcPDFEnvelope(int idNow, double xNow, double Q2NowIn,
233  int valSea) {
234 
235  // Freeze at boundary value if PDF is evaluated outside the fit region.
236  double x1 = (xNow < pdf->xMin() && !extrapol)
237  ? pdf->xMin() : xNow;
238  if (x1 > pdf->xMax() ) x1 = pdf->xMax();
239  double Q2Now = (Q2NowIn < pdf->q2Min() )
240  ? pdf->q2Min() : Q2NowIn;
241  if (Q2Now > pdf->q2Max() ) Q2Now = pdf->q2Max();
242 
243  // Loop over the members.
244  vector<double> xfCalc((*pdfs).size());
245  for(int iMem = 0; iMem < (*pdfs).size(); ++iMem) {
246  if (valSea==0 || (idNow != 1 && idNow != 2)) {
247  xfCalc[iMem] = (*pdfs)[iMem]->xfxQ2(idNow, x1, Q2Now);
248  } else if (valSea==1 && (idNow == 1 || idNow == 2 )) {
249  xfCalc[iMem] = (*pdfs)[iMem]->xfxQ2(idNow, x1, Q2Now) -
250  (*pdfs)[iMem]->xfxQ2(-idNow, x1, Q2Now);
251  } else if (valSea==2 && (idNow == 1 || idNow == 2 )) {
252  xfCalc[iMem] = (*pdfs)[iMem]->xfxQ2(-idNow, x1, Q2Now);
253  }
254  }
255 
256  // Calculate the uncertainty.
257  ::LHAPDF::PDFUncertainty xfErr = (*pdfs).info.uncertainty(xfCalc);
258  pdfEnvelope.centralPDF = xfErr.central;
259  pdfEnvelope.errplusPDF = xfErr.errplus;
260  pdfEnvelope.errminusPDF = xfErr.errminus;
261  pdfEnvelope.errsymmPDF = xfErr.errsymm;
262  pdfEnvelope.scalePDF = xfErr.scale;
263 }
264 
265 //--------------------------------------------------------------------------
266 
267 // Calculate uncertainties using the LHAPDF prescription.
268 
269 void LHAPDF6::calcPDFEnvelope(pair<int,int> idNows, pair<double,double> xNows,
270  double Q2NowIn, int valSea) {
271 
272  // Freeze at boundary value if PDF is evaluated outside the fit region.
273  double x1 = (xNows.first < pdf->xMin() && !extrapol)
274  ? pdf->xMin() : xNows.first;
275  if (x1 > pdf->xMax() ) x1 = pdf->xMax();
276  double x2 = (xNows.second < pdf->xMin() && !extrapol)
277  ? pdf->xMin() : xNows.second;
278  if (x2 > pdf->xMax() ) x2 = pdf->xMax();
279  double Q2Now = (Q2NowIn < pdf->q2Min() )
280  ? pdf->q2Min() : Q2NowIn;
281  if (Q2Now > pdf->q2Max() ) Q2Now = pdf->q2Max();
282 
283  // Loop over the members.
284  vector<double> xfCalc((*pdfs).size());
285  pdfEnvelope.pdfMemberVars.resize((*pdfs).size());
286  for(int iMem = 0; iMem < (*pdfs).size(); ++iMem) {
287  if (valSea == 0 || (idNows.first != 1 && idNows.first != 2 ) ) {
288  xfCalc[iMem] = (*pdfs)[iMem]->xfxQ2(idNows.first, x1, Q2Now);
289  } else if (valSea == 1 && (idNows.first == 1 || idNows.first == 2)) {
290  xfCalc[iMem] = (*pdfs)[iMem]->xfxQ2(idNows.first, x1, Q2Now)
291  - (*pdfs)[iMem]->xfxQ2(-idNows.first, x1, Q2Now);
292  } else if (valSea == 2 && (idNows.first == 1 || idNows.first == 2 )) {
293  xfCalc[iMem] = (*pdfs)[iMem]->xfxQ2(-idNows.first, x1, Q2Now);
294  }
295  xfCalc[iMem] = max(0.0, xfCalc[iMem]);
296  if (valSea == 0 || (idNows.second != 1 && idNows.second != 2)) {
297  xfCalc[iMem] /= max
298  (PDFMINVALUE, (*pdfs)[iMem]->xfxQ2(idNows.second, x2, Q2Now));
299  } else if (valSea == 1 && (idNows.second == 1 || idNows.second == 2 )) {
300  xfCalc[iMem] /= max
301  ((*pdfs)[iMem]->xfxQ2(idNows.second, x2, Q2Now) - (*pdfs)[iMem]->xfxQ2
302  (-idNows.second, x2, Q2Now), PDFMINVALUE);
303  } else if (valSea == 2 && (idNows.second == 1 || idNows.second == 2 )) {
304  xfCalc[iMem] /= max
305  ((*pdfs)[iMem]->xfxQ2(-idNows.second, x2, Q2Now), PDFMINVALUE);
306  }
307  pdfEnvelope.pdfMemberVars[iMem] = xfCalc[iMem];
308  }
309 
310  // Calculate the uncertainty.
311  ::LHAPDF::PDFUncertainty xfErr = (*pdfs).info.uncertainty(xfCalc);
312  pdfEnvelope.centralPDF = xfErr.central;
313  pdfEnvelope.errplusPDF = xfErr.errplus;
314  pdfEnvelope.errminusPDF = xfErr.errminus;
315  pdfEnvelope.errsymmPDF = xfErr.errsymm;
316  pdfEnvelope.scalePDF = xfErr.scale;
317 
318 }
319 
320 //--------------------------------------------------------------------------
321 
322 // Define external handles to the plugin for dynamic loading.
323 
324 extern "C" LHAPDF6* newLHAPDF(int idBeamIn, string setName, int member,
325  Info* infoPtr) {
326  return new LHAPDF6(idBeamIn, setName, member, 1, infoPtr);
327 
328 }
329 
330 extern "C" void deleteLHAPDF(LHAPDF6* pdf) {
331  delete pdf;
332 
333 }
334 
335 //==========================================================================
336 
337 } // end namespace Pythia8
338 
339 #endif // end Pythia8_LHAPDF6_H
Generated by   doxygen 1.8.5