Line data Source code
1 : // PartonDistributions.h is a part of the PYTHIA event generator.
2 : // Copyright (C) 2015 Torbjorn Sjostrand.
3 : // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4 : // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 :
6 : // Header file for parton densities.
7 : // PDF: base class.
8 : // LHAPDF: derived class for interface to the LHAPDF library.
9 : // GRV94L: derived class for the GRV 94L parton densities.
10 : // CTEQ5L: derived class for the CTEQ 5L parton densities.
11 : // MSTWpdf: derived class for MRST LO*, LO**, MSTW 2008 LO, NLO.
12 : // CTEQ6pdf: derived class for CTEQ 6L, 6L1, 66, CT09 MC1, MC2, (MCS?).
13 : // ProtonPoint: unresolved proton with equivalent photon spectrum.
14 : // GRVpiL: derived class for the GRV LO pion parton densities.
15 : // PomFix: derived class for Q2-independent Pomeron parton densities.
16 : // PomH1FitAB: derived class for the H1 2006 Fit A and Fit B Pomeron PDFs.
17 : // PomH1Jets: derived class for the H1 2007 Jets Pomeron PDFs.
18 : // Lepton: derived class for parton densities inside a lepton.
19 : // LeptonPoint: derived class for unresolved lepton (mainly dummy).
20 : // NNPDF: derived class for the NNPDF2.3 QCD+QED PDF sets.
21 :
22 : #ifndef Pythia8_PartonDistributions_H
23 : #define Pythia8_PartonDistributions_H
24 :
25 : #include "Pythia8/Basics.h"
26 : #include "Pythia8/Info.h"
27 : #include "Pythia8/ParticleData.h"
28 : #include "Pythia8/PythiaStdlib.h"
29 :
30 : namespace Pythia8 {
31 :
32 : //==========================================================================
33 :
34 : // Base class for parton distribution functions.
35 :
36 : class PDF {
37 :
38 : public:
39 :
40 : // Constructor.
41 0 : PDF(int idBeamIn = 2212) {idBeam = idBeamIn; idBeamAbs = abs(idBeam);
42 0 : setValenceContent(); idSav = 9; xSav = -1.; Q2Sav = -1.;
43 0 : xu = 0.; xd = 0.; xs = 0.; xubar = 0.; xdbar = 0.; xsbar = 0.; xc = 0.;
44 0 : xb = 0.; xg = 0.; xlepton = 0.; xgamma = 0.; xuVal = 0.; xuSea = 0.;
45 0 : xdVal = 0.; xdSea = 0.; isSet = true; isInit = false;}
46 :
47 : // Destructor.
48 0 : virtual ~PDF() {}
49 :
50 : // Confirm that PDF has been set up (important for LHAPDF and H1 Pomeron).
51 0 : virtual bool isSetup() {return isSet;}
52 :
53 : // Dynamic choice of meson valence flavours for pi0, K0S, K0L, Pomeron.
54 : virtual void newValenceContent(int idVal1In, int idVal2In) {
55 0 : idVal1 = idVal1In; idVal2 = idVal2In;}
56 :
57 : // Allow extrapolation beyond boundaries. This is optional.
58 0 : virtual void setExtrapolate(bool) {}
59 :
60 : // Read out parton density
61 : virtual double xf(int id, double x, double Q2);
62 :
63 : // Read out valence and sea part of parton densities.
64 : virtual double xfVal(int id, double x, double Q2);
65 : virtual double xfSea(int id, double x, double Q2);
66 :
67 : // Check whether x and Q2 values fall inside the fit bounds (LHAPDF6 only).
68 0 : virtual bool insideBounds(double, double) {return true;}
69 :
70 : // Access the running alpha_s of a PDF set (LHAPDF6 only).
71 0 : virtual double alphaS(double) { return 1.;}
72 :
73 : // Return quark masses used in the PDF fit (LHAPDF6 only).
74 0 : virtual double mQuarkPDF(int) { return -1.;}
75 :
76 : protected:
77 :
78 : // Allow the LHAPDF class to access these methods.
79 : friend class LHAPDF;
80 :
81 : // Store relevant quantities.
82 : int idBeam, idBeamAbs, idSav, idVal1, idVal2;
83 : double xSav, Q2Sav;
84 : double xu, xd, xs, xubar, xdbar, xsbar, xc, xb, xg, xlepton, xgamma,
85 : xuVal, xuSea, xdVal, xdSea;
86 : bool isSet, isInit;
87 :
88 : // Resolve valence content for assumed meson. Possibly modified later.
89 : void setValenceContent();
90 :
91 : // Update parton densities.
92 : virtual void xfUpdate(int id, double x, double Q2) = 0;
93 :
94 : };
95 :
96 : //==========================================================================
97 :
98 : // Gives the GRV 94L (leading order) parton distribution function set
99 : // in parametrized form. Authors: M. Glueck, E. Reya and A. Vogt.
100 :
101 0 : class GRV94L : public PDF {
102 :
103 : public:
104 :
105 : // Constructor.
106 0 : GRV94L(int idBeamIn = 2212) : PDF(idBeamIn) {}
107 :
108 : private:
109 :
110 : // Update PDF values.
111 : void xfUpdate(int , double x, double Q2);
112 :
113 : // Auxiliary routines used during the updating.
114 : double grvv (double x, double n, double ak, double bk, double a,
115 : double b, double c, double d);
116 : double grvw (double x, double s, double al, double be, double ak,
117 : double bk, double a, double b, double c, double d, double e, double es);
118 : double grvs (double x, double s, double sth, double al, double be,
119 : double ak, double ag, double b, double d, double e, double es);
120 :
121 : };
122 :
123 : //==========================================================================
124 :
125 : // Gives the CTEQ 5L (leading order) parton distribution function set
126 : // in parametrized form. Parametrization by J. Pumplin. Authors: CTEQ.
127 :
128 0 : class CTEQ5L : public PDF {
129 :
130 : public:
131 :
132 : // Constructor.
133 0 : CTEQ5L(int idBeamIn = 2212) : PDF(idBeamIn) {}
134 :
135 : private:
136 :
137 : // Update PDF values.
138 : void xfUpdate(int , double x, double Q2);
139 :
140 : };
141 :
142 : //==========================================================================
143 :
144 : // The MSTWpdf class.
145 : // MRST LO*(*) and MSTW 2008 PDF's, specifically the LO one.
146 : // Original C++ version by Jeppe Andersen.
147 : // Modified by Graeme Watt <watt(at)hep.ucl.ac.uk>.
148 : // Sets available:
149 : // iFit = 1 : MRST LO* (2007).
150 : // iFit = 2 : MRST LO** (2008).
151 : // iFit = 3 : MSTW 2008 LO, central member.
152 : // iFit = 4 : MSTW 2008 NLO, central member. (Warning!)
153 :
154 0 : class MSTWpdf : public PDF {
155 :
156 : public:
157 :
158 : // Constructor.
159 : MSTWpdf(int idBeamIn = 2212, int iFitIn = 1,
160 : string xmlPath = "../share/Pythia8/xmldoc/", Info* infoPtr = 0)
161 0 : : PDF(idBeamIn) {init( iFitIn, xmlPath, infoPtr);}
162 :
163 : private:
164 :
165 : // Constants: could only be changed in the code itself.
166 : static const int np, nx, nq, nqc0, nqb0;
167 : static const double xmin, xmax, qsqmin, qsqmax, xxInit[65], qqInit[49];
168 :
169 : // Data read in from grid file or set at initialization.
170 : int iFit, alphaSorder, alphaSnfmax;
171 : double mCharm, mBottom, alphaSQ0, alphaSMZ, distance, tolerance,
172 : xx[65], qq[49], c[13][64][48][5][5];
173 :
174 : // Initialization of data array.
175 : void init( int iFitIn, string xmlPath, Info* infoPtr);
176 :
177 : // Update PDF values.
178 : void xfUpdate(int , double x, double Q2);
179 :
180 : // Evaluate PDF of one flavour species.
181 : double parton(int flavour,double x,double q);
182 : double parton_interpolate(int flavour,double xxx,double qqq);
183 : double parton_extrapolate(int flavour,double xxx,double qqq);
184 :
185 : // Auxiliary routines for evaluation.
186 : int locate(double xx[],int n,double x);
187 : double polderivative1(double x1, double x2, double x3, double y1,
188 : double y2, double y3);
189 : double polderivative2(double x1, double x2, double x3, double y1,
190 : double y2, double y3);
191 : double polderivative3(double x1, double x2, double x3, double y1,
192 : double y2, double y3);
193 :
194 : };
195 :
196 : //==========================================================================
197 :
198 : // The CTEQ6pdf class.
199 : // Sets available:
200 : // iFit = 1 : CTEQ6L
201 : // iFit = 2 : CTEQ6L1
202 : // iFit = 3 : CTEQ66.00 (NLO, central member)
203 : // iFit = 4 : CT09MC1
204 : // iFit = 5 : CT09MC2
205 : // iFit = 6 : CT09MCS (not yet implemented)
206 :
207 0 : class CTEQ6pdf : public PDF {
208 :
209 : public:
210 :
211 : // Constructor.
212 : CTEQ6pdf(int idBeamIn = 2212, int iFitIn = 1,
213 : string xmlPath = "../share/Pythia8/xmldoc/", Info* infoPtr = 0)
214 0 : : PDF(idBeamIn) {init( iFitIn, xmlPath, infoPtr);}
215 :
216 : private:
217 :
218 : // Constants: could only be changed in the code itself.
219 : static const double EPSILON, XPOWER;
220 :
221 : // Data read in from grid file or set at initialization.
222 : int iFit, order, nQuark, nfMx, mxVal, nX, nT, nG,
223 : iGridX, iGridQ, iGridLX, iGridLQ;
224 : double lambda, mQ[7], qIni, qMax, tv[26], xMin, xv[202], upd[57773],
225 : xvpow[202], xMinEps, xMaxEps, qMinEps, qMaxEps, fVec[5],
226 : tConst[9], xConst[9], xLast, qLast;
227 :
228 : // Initialization of data array.
229 : void init( int iFitIn, string xmlPath, Info* infoPtr);
230 :
231 : // Update PDF values.
232 : void xfUpdate(int id, double x, double Q2);
233 :
234 : // Evaluate PDF of one flavour species.
235 : double parton6(int iParton, double x, double q);
236 :
237 : // Interpolation in grid.
238 : double polint4F(double xgrid[], double fgrid[], double xin);
239 :
240 : };
241 :
242 : //==========================================================================
243 :
244 : // SA Unresolved proton: equivalent photon spectrum from
245 : // V.M. Budnev, I.F. Ginzburg, G.V. Meledin and V.G. Serbo,
246 : // Phys. Rept. 15 (1974/1975) 181.
247 :
248 0 : class ProtonPoint : public PDF {
249 :
250 : public:
251 :
252 : // Constructor.
253 : ProtonPoint(int idBeamIn = 2212, Info* infoPtrIn = 0) :
254 : PDF(idBeamIn), m_infoPtr(infoPtrIn) {}
255 :
256 : private:
257 :
258 : // Stored value for PDF choice.
259 : static const double ALPHAEM, Q2MAX, Q20, A, B, C;
260 :
261 : // Update PDF values.
262 : void xfUpdate(int , double x, double Q2);
263 :
264 : // phi function from Q2 integration.
265 : double phiFunc(double x, double Q);
266 :
267 : // Info and errors
268 : Info* m_infoPtr;
269 :
270 : };
271 :
272 : //==========================================================================
273 :
274 : // Gives the GRV 1992 pi+ (leading order) parton distribution function set
275 : // in parametrized form. Authors: Glueck, Reya and Vogt.
276 :
277 0 : class GRVpiL : public PDF {
278 :
279 : public:
280 :
281 : // Constructor.
282 0 : GRVpiL(int idBeamIn = 221) : PDF(idBeamIn) {}
283 :
284 : private:
285 :
286 : // Update PDF values.
287 : void xfUpdate(int , double x, double Q2);
288 :
289 : };
290 :
291 : //==========================================================================
292 :
293 : // Gives generic Q2-independent Pomeron PDF.
294 :
295 0 : class PomFix : public PDF {
296 :
297 : public:
298 :
299 : // Constructor.
300 : PomFix(int idBeamIn = 990, double PomGluonAIn = 0.,
301 : double PomGluonBIn = 0., double PomQuarkAIn = 0.,
302 : double PomQuarkBIn = 0., double PomQuarkFracIn = 0.,
303 0 : double PomStrangeSuppIn = 0.) : PDF(idBeamIn),
304 0 : PomGluonA(PomGluonAIn), PomGluonB(PomGluonBIn),
305 0 : PomQuarkA(PomQuarkAIn), PomQuarkB(PomQuarkBIn),
306 0 : PomQuarkFrac(PomQuarkFracIn), PomStrangeSupp(PomStrangeSuppIn)
307 0 : {init();}
308 :
309 : private:
310 :
311 : // Stored value for PDF choice.
312 : double PomGluonA, PomGluonB, PomQuarkA, PomQuarkB, PomQuarkFrac,
313 : PomStrangeSupp, normGluon, normQuark;
314 :
315 : // Initialization of some constants.
316 : void init();
317 :
318 : // Update PDF values.
319 : void xfUpdate(int , double x, double);
320 :
321 : };
322 :
323 : //==========================================================================
324 :
325 : // The H1 2006 Fit A and Fit B Pomeron parametrization.
326 : // H1 Collaboration, A. Aktas et al., "Measurement and QCD Analysis of
327 : // the Diffractive Deep-Inelastic Scattering Cross Section at HERA",
328 : // DESY-06-049, Eur. Phys. J. C48 (2006) 715. e-Print: hep-ex/0606004.
329 :
330 0 : class PomH1FitAB : public PDF {
331 :
332 : public:
333 :
334 : // Constructor.
335 : PomH1FitAB(int idBeamIn = 990, int iFit = 1, double rescaleIn = 1.,
336 : string xmlPath = "../share/Pythia8/xmldoc/", Info* infoPtr = 0)
337 0 : : PDF(idBeamIn) {rescale = rescaleIn; init( iFit, xmlPath, infoPtr);}
338 :
339 : private:
340 :
341 : // Limits for grid in x, in Q2, and data in (x, Q2).
342 : int nx, nQ2;
343 : double rescale, xlow, xupp, dx, Q2low, Q2upp, dQ2;
344 : double gluonGrid[100][30];
345 : double quarkGrid[100][30];
346 :
347 : // Initialization of data array.
348 : void init( int iFit, string xmlPath, Info* infoPtr);
349 :
350 : // Update PDF values.
351 : void xfUpdate(int , double x, double );
352 :
353 : };
354 :
355 : //==========================================================================
356 :
357 : // The H1 2007 Jets Pomeron parametrization..
358 : // H1 Collaboration, A. Aktas et al., "Dijet Cross Sections and Parton
359 : // Densities in Diffractive DIS at HERA", DESY-07-115, Aug 2007. 33pp.
360 : // Published in JHEP 0710:042,2007. e-Print: arXiv:0708.3217 [hep-ex]
361 :
362 0 : class PomH1Jets : public PDF {
363 :
364 : public:
365 :
366 : // Constructor.
367 : PomH1Jets(int idBeamIn = 990, double rescaleIn = 1.,
368 : string xmlPath = "../share/Pythia8/xmldoc/", Info* infoPtr = 0)
369 0 : : PDF(idBeamIn) {rescale = rescaleIn; init( xmlPath, infoPtr);}
370 :
371 : private:
372 :
373 : // Arrays for grid in x, in Q2, and data in (x, Q2).
374 : double rescale;
375 : double xGrid[100];
376 : double Q2Grid[88];
377 : double gluonGrid[100][88];
378 : double singletGrid[100][88];
379 : double charmGrid[100][88];
380 :
381 : // Initialization of data array.
382 : void init( string xmlPath, Info* infoPtr);
383 :
384 : // Update PDF values.
385 : void xfUpdate(int id, double x, double );
386 :
387 : };
388 :
389 : //==========================================================================
390 :
391 : // Gives electron (or muon, or tau) parton distribution.
392 :
393 0 : class Lepton : public PDF {
394 :
395 : public:
396 :
397 : // Constructor.
398 0 : Lepton(int idBeamIn = 11) : PDF(idBeamIn) {}
399 :
400 : private:
401 :
402 : // Constants: could only be changed in the code itself.
403 : static const double ALPHAEM, ME, MMU, MTAU;
404 :
405 : // Update PDF values.
406 : void xfUpdate(int id, double x, double Q2);
407 :
408 : // The squared lepton mass, set at initialization.
409 : double m2Lep;
410 :
411 : };
412 :
413 : //==========================================================================
414 :
415 : // Gives electron (or other lepton) parton distribution when unresolved.
416 :
417 0 : class LeptonPoint : public PDF {
418 :
419 : public:
420 :
421 : // Constructor.
422 0 : LeptonPoint(int idBeamIn = 11) : PDF(idBeamIn) {}
423 :
424 : private:
425 :
426 : // Update PDF values in trivial way.
427 0 : void xfUpdate(int , double , double ) {xlepton = 1; xgamma = 0.;}
428 :
429 : };
430 :
431 : //==========================================================================
432 :
433 : // Gives neutrino parton distribution when unresolved (only choice for now).
434 : // Note factor of 2 since only lefthanded implies no spin averaging.
435 :
436 0 : class NeutrinoPoint : public PDF {
437 :
438 : public:
439 :
440 : // Constructor.
441 0 : NeutrinoPoint(int idBeamIn = 12) : PDF(idBeamIn) {}
442 :
443 : private:
444 :
445 : // Update PDF values, with spin factor of 2.
446 0 : void xfUpdate(int , double , double ) {xlepton = 2; xgamma = 0.;}
447 :
448 : };
449 :
450 : //==========================================================================
451 :
452 : // The NNPDF class.
453 : // Sets available:
454 : // Leading order QCD+QED Proton PDF sets
455 : // iFit = 1 : NNPDF2.3 QCD+QED LO, alphas(MZ) = 0.130
456 : // iFit = 2 : NNPDF2.3 QCD+QED LO, alphas(MZ) = 0.119
457 : // (Next-to-)Next-to-Leading order QCD+QED Proton PDF sets
458 : // iFit = 3 : NNPDF2.3 QCD+QED NLO, alphas(MZ) = 0.119
459 : // iFit = 4 : NNPDF2.3 QCD+QED NNLO, alphas(MZ) = 0.119
460 : // Code provided by Juan Rojo and Stefano Carrazza.
461 :
462 : class NNPDF : public PDF {
463 :
464 : public:
465 :
466 : // Constructor.
467 : NNPDF(int idBeamIn = 2212, int iFitIn = 1,
468 : string xmlPath = "../share/Pythia8/xmldoc/", Info* infoPtr = 0)
469 0 : : PDF(idBeamIn), fPDFGrid(NULL), fXGrid(NULL), fLogXGrid(NULL),
470 0 : fQ2Grid(NULL), fLogQ2Grid(NULL), fRes(NULL) {
471 0 : init( iFitIn, xmlPath, infoPtr); };
472 :
473 : // Destructor.
474 0 : ~NNPDF() {
475 0 : if (fPDFGrid) {
476 0 : for (int i = 0; i < fNFL; i++) {
477 0 : for (int j = 0; j < fNX; j++)
478 0 : if (fPDFGrid[i][j]) delete[] fPDFGrid[i][j];
479 0 : if (fPDFGrid[i]) delete[] fPDFGrid[i];
480 : }
481 0 : delete[] fPDFGrid;
482 : }
483 0 : if (fXGrid) delete[] fXGrid;
484 0 : if (fLogXGrid) delete[] fLogXGrid;
485 0 : if (fQ2Grid) delete[] fQ2Grid;
486 0 : if (fLogQ2Grid) delete[] fLogQ2Grid;
487 0 : if (fRes) delete[] fRes;
488 0 : };
489 :
490 : private:
491 :
492 : // Constants: could only be changed in the code itself.
493 : static const double fXMINGRID;
494 :
495 : // Number of flavors (including photon) and interpolation parameters.
496 : static const int fNFL = 14;
497 : static const int fM = 4;
498 : static const int fN = 2;
499 :
500 : // Variables to be set during code initialization.
501 : int iFit, fNX, fNQ2;
502 : double ***fPDFGrid;
503 : double *fXGrid;
504 : double *fLogXGrid;
505 : double *fQ2Grid;
506 : double *fLogQ2Grid;
507 : double *fRes;
508 :
509 : // Initialization of data array.
510 : void init( int iFitIn, string xmlPath, Info* infoPtr);
511 :
512 : // Update PDF values.
513 : void xfUpdate(int id, double x, double Q2);
514 :
515 : // Interpolation in the grid for a given PDF flavour.
516 : void xfxevolve(double x, double Q2);
517 :
518 : // 1D and 2D polynomial interpolation.
519 : void polint(double xa[], double ya[], int n, double x,
520 : double& y, double& dy);
521 : void polin2(double x1a[], double x2a[], double ya[][fN],
522 : double x1, double x2, double& y, double& dy);
523 :
524 : };
525 :
526 : //==========================================================================
527 :
528 : // LHAPDF plugin interface class.
529 :
530 : class LHAPDF : public PDF {
531 :
532 : public:
533 :
534 : // Constructor and destructor.
535 : LHAPDF(int idIn, string pSet, Info* infoPtrIn);
536 : ~LHAPDF();
537 :
538 : // Confirm that PDF has been set up.
539 0 : bool isSetup() {if (pdfPtr) return pdfPtr->isSetup(); return false;}
540 :
541 : // Dynamic choice of meson valence flavours for pi0, K0S, K0L, Pomeron.
542 : void newValenceContent(int idVal1In, int idVal2In) {
543 0 : if (pdfPtr) pdfPtr->newValenceContent(idVal1In, idVal2In);}
544 :
545 : // Allow extrapolation beyond boundaries.
546 : void setExtrapolate(bool extrapolate) {
547 0 : if (pdfPtr) pdfPtr->setExtrapolate(extrapolate);}
548 :
549 : // Read out parton density
550 : double xf(int id, double x, double Q2) {
551 0 : if (pdfPtr) return pdfPtr->xf(id, x, Q2); else return 0;}
552 :
553 : // Read out valence and sea part of parton densities.
554 : double xfVal(int id, double x, double Q2) {
555 0 : if (pdfPtr) return pdfPtr->xfVal(id, x, Q2); else return 0;}
556 : double xfSea(int id, double x, double Q2) {
557 0 : if (pdfPtr) return pdfPtr->xfSea(id, x, Q2); else return 0;}
558 :
559 : // Check whether x and Q2 values fall inside the fit bounds (LHAPDF6 only).
560 : bool insideBounds(double x, double Q2) {
561 0 : if(pdfPtr) return pdfPtr->insideBounds(x, Q2); else return true;}
562 :
563 : // Access the running alpha_s of a PDF set (LHAPDF6 only).
564 : double alphaS(double Q2) {
565 0 : if(pdfPtr) return pdfPtr->alphaS(Q2); else return 1.;}
566 :
567 : // Return quark masses used in the PDF fit (LHAPDF6 only).
568 : double mQuarkPDF(int idIn) {
569 0 : if(pdfPtr) return pdfPtr->mQuarkPDF(idIn); else return -1.;}
570 :
571 : private:
572 :
573 : // Resolve valence content for assumed meson.
574 : void setValenceContent() {if (pdfPtr) pdfPtr->setValenceContent();}
575 :
576 : // Update parton densities.
577 : void xfUpdate(int id, double x, double Q2) {
578 0 : if (pdfPtr) pdfPtr->xfUpdate(id, x, Q2);}
579 :
580 : // Typedefs of the hooks used to access the plugin.
581 : typedef PDF* NewLHAPDF(int, string, int, Info*);
582 : typedef void DeleteLHAPDF(PDF*);
583 : typedef void (*Symbol)();
584 :
585 : // Acccess a plugin library symbol.
586 : Symbol symbol(string symName);
587 :
588 : // The loaded LHAPDF object, info pointer, and plugin library name.
589 : PDF *pdfPtr;
590 : Info *infoPtr;
591 : string libName;
592 :
593 : };
594 :
595 : //==========================================================================
596 :
597 : } // end namespace Pythia8
598 :
599 : #endif // Pythia8_PartonDistributions_H
|
