1 | // main51.cc is a part of the PYTHIA event generator. |
---|
2 | // Copyright (C) 2012 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 | // Test of LHAPDF interface and whether PDF's behave sensibly. |
---|
7 | |
---|
8 | #include "Pythia.h" |
---|
9 | |
---|
10 | using namespace Pythia8; |
---|
11 | |
---|
12 | //========================================================================== |
---|
13 | |
---|
14 | // Integration to check momentum sum rule. |
---|
15 | |
---|
16 | double integrate(PDF* nowPDF, double Q2) { |
---|
17 | |
---|
18 | // Number of points, x ranges and initial values. |
---|
19 | int nLin = 980; |
---|
20 | int nLog = 1000; |
---|
21 | double xLin = 0.02; |
---|
22 | double xLog = 1e-8; |
---|
23 | double dxLin = (1. - xLin) / nLin; |
---|
24 | double dxLog = log(xLin / xLog) / nLog; |
---|
25 | double sum = 0.; |
---|
26 | double x, sumNow; |
---|
27 | |
---|
28 | // Integration at large x in linear steps. |
---|
29 | for (int iLin = 0; iLin < nLin; ++iLin) { |
---|
30 | x = xLin + (iLin + 0.5) * dxLin; |
---|
31 | sumNow = nowPDF->xf( 21, x, Q2); |
---|
32 | for (int i = 1; i < 6; ++i) |
---|
33 | sumNow += nowPDF->xf( i, x, Q2) + nowPDF->xf( -i, x, Q2); |
---|
34 | sum += dxLin * sumNow; |
---|
35 | } |
---|
36 | |
---|
37 | // Integration at small x in logarithmic steps. |
---|
38 | for (int iLog = 0; iLog < nLog; ++iLog) { |
---|
39 | x = xLog * pow( xLin / xLog, (iLog + 0.5) / nLog ); |
---|
40 | sumNow = nowPDF->xf( 21, x, Q2); |
---|
41 | for (int i = 1; i < 6; ++i) |
---|
42 | sumNow += nowPDF->xf( i, x, Q2) + nowPDF->xf( -i, x, Q2); |
---|
43 | sum += dxLog * x * sumNow; |
---|
44 | } |
---|
45 | |
---|
46 | // Done. |
---|
47 | return sum; |
---|
48 | |
---|
49 | } |
---|
50 | |
---|
51 | //========================================================================== |
---|
52 | |
---|
53 | int main() { |
---|
54 | |
---|
55 | // The Pythia class itself is not used, but some facilities that com along. |
---|
56 | Pythia pythia; |
---|
57 | |
---|
58 | // Chosen new PDF set; LHAPDF file name conventions. |
---|
59 | //string pdfSet = "cteq5l.LHgrid"; |
---|
60 | //string pdfSet = "cteq61.LHpdf"; |
---|
61 | //string pdfSet = "cteq61.LHgrid"; |
---|
62 | //string pdfSet = "MRST2004nlo.LHgrid"; |
---|
63 | //string pdfSet = "MRST2001lo.LHgrid"; |
---|
64 | string pdfSet = "MRST2007lomod.LHgrid"; |
---|
65 | |
---|
66 | // Pointers to old default and new tryout PDF sets. |
---|
67 | PDF* oldPDF = new CTEQ5L(2212); |
---|
68 | PDF* newPDF = new LHAPDF(2212, pdfSet, 0); |
---|
69 | |
---|
70 | // Alternative: compare two Pomeron PDF's. Boost second by factor 2. |
---|
71 | //PDF* oldPDF = new PomFix( 990, -0.2, 2.5, 0., 3., 0.4, 0.5); |
---|
72 | //PDF* newPDF = new PomH1Jets( 990, 2.); |
---|
73 | //PDF* oldPDF = new PomH1FitAB( 990, 2); |
---|
74 | //PDF* newPDF = new PomH1FitAB( 990, 3); |
---|
75 | |
---|
76 | // Allow extrapolation of PDF's beyond x and Q2 boundaries, at own risk. |
---|
77 | // Default behaviour is to freeze PDF's at boundaries. |
---|
78 | newPDF->setExtrapolate(true); |
---|
79 | |
---|
80 | // Histogram F(x, Q2) = (9/4) x*g(x, Q2) + sum_{i = q, qbar} x*f_i(x, Q2) |
---|
81 | // for range 10^{-8} < x < 1 logarithmic in x and for Q2 = 4 and 100. |
---|
82 | Hist oldF4("F( x, Q2 = 4) old", 80 , -8., 0.); |
---|
83 | Hist newF4("F( x, Q2 = 4) new", 80 , -8., 0.); |
---|
84 | Hist ratF4("F( x, Q2 = 4) new/old", 80 , -8., 0.); |
---|
85 | Hist oldF100("F( x, Q2 = 100) old", 80 , -8., 0.); |
---|
86 | Hist newF100("F( x, Q2 = 100) new", 80 , -8., 0.); |
---|
87 | Hist ratF100("F( x, Q2 = 100) new/old", 80 , -8., 0.); |
---|
88 | |
---|
89 | // Loop over the two Q2 values. |
---|
90 | for (int iQ = 0; iQ < 2; ++iQ) { |
---|
91 | double Q2 = (iQ == 0) ? 4. : 100; |
---|
92 | |
---|
93 | // Loop over x values, in a logarithmic scale |
---|
94 | for (int iX = 0; iX < 80; ++iX) { |
---|
95 | double xLog = -(0.1 * iX + 0.05); |
---|
96 | double x = pow( 10., xLog); |
---|
97 | |
---|
98 | // Evaluate old summed PDF. |
---|
99 | double oldSum = 2.25 * oldPDF->xf( 21, x, Q2); |
---|
100 | for (int i = 1; i < 6; ++i) |
---|
101 | oldSum += oldPDF->xf( i, x, Q2) + oldPDF->xf( -i, x, Q2); |
---|
102 | if (iQ == 0) oldF4.fill ( xLog, oldSum ); |
---|
103 | else oldF100.fill ( xLog, oldSum ); |
---|
104 | |
---|
105 | // Evaluate new summed PDF. |
---|
106 | double newSum = 2.25 * newPDF->xf( 21, x, Q2); |
---|
107 | for (int i = 1; i < 6; ++i) |
---|
108 | newSum += newPDF->xf( i, x, Q2) + newPDF->xf( -i, x, Q2); |
---|
109 | if (iQ == 0) newF4.fill ( xLog, newSum ); |
---|
110 | else newF100.fill ( xLog, newSum ); |
---|
111 | |
---|
112 | //End loops over x and Q2 values. |
---|
113 | } |
---|
114 | } |
---|
115 | |
---|
116 | // Show F(x, Q2) and their ratio new/old. |
---|
117 | ratF4 = newF4 / oldF4; |
---|
118 | ratF100 = newF100 / oldF100; |
---|
119 | cout << oldF4 << newF4 << ratF4 << oldF100 << newF100 << ratF100; |
---|
120 | |
---|
121 | // Histogram momentum sum as a function of Q2 (or rather log10(Q2)). |
---|
122 | Hist oldXSum("momentum sum(log10(Q2)) old", 100, -2., 8.); |
---|
123 | Hist newXSum("momentum sum(log10(Q2)) new", 100, -2., 8.); |
---|
124 | |
---|
125 | // Loop over Q2 values. |
---|
126 | for (int iQ = 0; iQ < 100; ++iQ) { |
---|
127 | double log10Q2 = -2.0 + 0.1 * iQ + 0.05; |
---|
128 | double Q2 = pow( 10., log10Q2); |
---|
129 | |
---|
130 | // Evaluate old and new momentum sums. |
---|
131 | double oldSum = integrate( oldPDF, Q2); |
---|
132 | oldXSum.fill( log10Q2, oldSum); |
---|
133 | double newSum = integrate( newPDF, Q2); |
---|
134 | newXSum.fill( log10Q2, newSum); |
---|
135 | } |
---|
136 | |
---|
137 | // Show momentum sum as a function of Q2. |
---|
138 | cout << oldXSum << newXSum; |
---|
139 | |
---|
140 | // Done. |
---|
141 | return 0; |
---|
142 | } |
---|