StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
main51.cc
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 }