StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
main52.cc
1 // main52.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2014 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 // Studies of hadron-level and parton-level minimum-bias quantities,
7 // comparing internal default PDF with one from LHAPDF.
8 // Major differences indicate need for major retuning, e.g. pT0Ref.
9 
10 // Access time information.
11 #include <ctime>
12 
13 #include "Pythia8/Pythia.h"
14 
15 using namespace Pythia8;
16 
17 int main() {
18 
19  // Machine: 1 = Tevatron, 2 = LHC. Statistics.
20  int machine = 1;
21  int nEvent = 10000;
22 
23  // Select new PDF set; LHAPDF file name conventions.
24  //string pdfSet = "cteq5l.LHgrid";
25  //string pdfSet = "cteq61.LHpdf";
26  //string pdfSet = "cteq61.LHgrid";
27  //string pdfSet = "MRST2004nlo.LHgrid";
28  string pdfSet = "MRST2001lo.LHgrid";
29 
30  // Histograms for hadron-level quantities.
31  double nMax = (machine == 1) ? 199.5 : 399.5;
32  Hist nChargedOld("n_charged old PDF", 100, -0.5, nMax);
33  Hist nChargedNew("n_charged new PDF", 100, -0.5, nMax);
34  Hist nChargedRat("n_charged new/old PDF", 100, -0.5, nMax);
35  Hist ySpecOld("y charged distribution old PDF", 100, -10., 10.);
36  Hist ySpecNew("y charged distribution new PDF", 100, -10., 10.);
37  Hist ySpecRat("y charged distribution new/old PDF", 100, -10., 10.);
38  Hist pTSpecOld("pT charged distribution old PDF", 100, 0., 20.);
39  Hist pTSpecNew("pT charged distribution new PDF", 100, 0., 20.);
40  Hist pTSpecRat("pT charged distribution new/old PDF", 100, 0., 20.);
41  Hist avgPTnChOld("<pT>(n_charged) old PDF", 100, -0.5, nMax);
42  Hist avgPTnChNew("<pT>(n_charged) new PDF", 100, -0.5, nMax);
43  Hist avgPTnChRat("<pT>(n_charged) new/old PDF", 100, -0.5, nMax);
44 
45  // Histograms for parton-level quantities.
46  Hist xDistOld("MPI log(x) distribution old PDF", 80, -8., 0.);
47  Hist xDistNew("MPI log(x) distribution new PDF", 80, -8., 0.);
48  Hist xDistRat("MPI log(x) distribution new/old PDF", 80, -8., 0.);
49  Hist pTDistOld("MPI pT (=Q) distribution old PDF", 100, 0., 20.);
50  Hist pTDistNew("MPI pT (=Q) distribution new PDF", 100, 0., 20.);
51  Hist pTDistRat("MPI pT (=Q) distribution new/old PDF", 100, 0., 20.);
52 
53  // Loop over one default run and one with new PDF.
54  for (int iRun = 0; iRun < 2; ++iRun) {
55 
56  // Get starting time in seconds.
57  time_t tBegin = time(0);
58 
59  // Generator.
60  Pythia pythia;
61  Event& event = pythia.event;
62 
63  // Generate minimum-bias events, with or without double diffraction.
64  pythia.readString("SoftQCD:nonDiffractive = on");
65  //pythia.readString("SoftQCD:doubleDiffractive = on");
66 
67  // Generate QCD jet events, above some threshold.
68  //pythia.readString("HardQCD:all = on");
69  //pythia.readString("PhaseSpace:pTHatMin = 50.");
70 
71  // In second run pick new PDF set.
72  if (iRun == 1) {
73  pythia.readString("PDF:useLHAPDF = on");
74  pythia.readString("PDF:LHAPDFset = " + pdfSet);
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  pythia.readString("PDF:extrapolateLHAPDF = on");
79 
80  // Need to change pT0Ref depending on choice of PDF.
81  // One possibility: retune to same <n_charged>.
82  //pythia.readString("MultipartonInteractions:pT0Ref = 2.17");
83  }
84 
85  // Tevatron/LHC initialization.
86  double eCM = (machine == 1) ? 1960. : 7000.;
87  pythia.settings.parm("Beams:eCM", eCM);
88  if (machine == 1) pythia.readString("Beams:idB = -2212");
89  pythia.init();
90 
91  // Begin event loop.
92  for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
93 
94  // Generate events. Skip if error.
95  if (!pythia.next()) continue;
96 
97  // Statistics on multiplicity and pT.
98  int nCh = 0;
99  double pTsum = 0.;
100  for (int i = 0; i < event.size(); ++i)
101  if (event[i].isFinal() && event[i].isCharged()) {
102  ++nCh;
103  pTsum += event[i].pT();
104 
105  // Fill histograms for charged y and pT spectra.
106  if (iRun == 0) {
107  ySpecOld.fill( event[i].y() );
108  pTSpecOld.fill( event[i].pT() );
109  } else {
110  ySpecNew.fill( event[i].y() );
111  pTSpecNew.fill( event[i].pT() );
112  }
113  }
114 
115  // Fill histograms for summed quantities.
116  if (iRun == 0) {
117  nChargedOld.fill( nCh );
118  avgPTnChOld.fill( nCh, pTsum / max(1, nCh) );
119  } else {
120  nChargedNew.fill( nCh );
121  avgPTnChNew.fill( nCh, pTsum / max(1, nCh) );
122  }
123 
124  // Loop through event record and fill x of all incoming partons.
125  for (int i = 1; i < event.size(); ++i)
126  if (event[i].status() == -21 || event[i].status() == -31) {
127  double x = 2. * event[i].e() / eCM;
128  if (iRun == 0) xDistOld.fill( log10(x) );
129  else xDistNew.fill( log10(x) );
130  }
131 
132  // Loop through multiparton interactions list and fill pT of all MPI's.
133  for (int i = 0; i < pythia.info.nMPI(); ++i) {
134  double pT = pythia.info.pTMPI(i);
135  if (iRun == 0) pTDistOld.fill( pT );
136  else pTDistNew.fill( pT );
137  }
138 
139  // End of event loop.
140  }
141 
142  // Statistics.
143  pythia.readString("Stat:showPartonLevel = on");
144  pythia.stat();
145 
146  // Get finishing time in seconds. Print used time.
147  time_t tEnd = time(0);
148  cout << "\n This subrun took " << tEnd - tBegin << " seconds \n" << endl;
149 
150  // End of loop over two runs.
151  }
152 
153  // Form <pT>(n_charged) ratios.
154  avgPTnChOld /= nChargedOld;
155  avgPTnChNew /= nChargedNew;
156 
157  // Take ratios of new to old distributions.
158  nChargedRat = nChargedNew / nChargedOld;
159  ySpecRat = ySpecNew / ySpecOld;
160  pTSpecRat = pTSpecNew / pTSpecOld;
161  avgPTnChRat = avgPTnChNew / avgPTnChOld;
162  xDistRat = xDistNew / xDistOld;
163  pTDistRat = pTDistNew / pTDistOld;
164 
165  // Print histograms.
166  cout << nChargedOld << nChargedNew << nChargedRat
167  << ySpecOld << ySpecNew << ySpecRat
168  << pTSpecOld << pTSpecNew << pTSpecRat
169  << avgPTnChOld << avgPTnChNew << avgPTnChRat
170  << xDistOld << xDistNew << xDistRat
171  << pTDistOld << pTDistNew << pTDistRat;
172 
173  // Second part of study, as simple extra check:
174  // Begin fill shape of effective PDF at typical MPI Q2 = 10 scale:
175  // F_effective(x) = (9/4) x*g(x) + Sum_i (x*q_i(x) + x*qbar_i(x)).
176  double Q2 = 10.;
177  PDF* oldPDF = new CTEQ5L(2212);
178  PDF* newPDF = new LHAPDF(2212, pdfSet, 0);
179 
180  // Histograms.
181  Hist effFlinOld("F_effective( x, Q2 = 10) old", 100 , 0., 1.);
182  Hist effFlinNew("F_effective( x, Q2 = 10) new", 100 , 0., 1.);
183  Hist effFlinRat("F_effective( x, Q2 = 10) new/old", 100 , 0., 1.);
184  Hist effFlogOld("F_effective( x, Q2 = 10) old", 80 , -8., 0.);
185  Hist effFlogNew("F_effective( x, Q2 = 10) new", 80 , -8., 0.);
186  Hist effFlogRat("F_effective( x, Q2 = 10) new/old", 80 , -8., 0.);
187 
188  // Loop over x values, in a linear scale.
189  for (int iX = 0; iX < 99; ++iX) {
190  double x = 0.005 + 0.01 * iX;
191 
192  // Evaluate old summed PDF.
193  double oldSum = 2.25 * oldPDF->xf( 21, x, Q2);
194  for (int i = 1; i < 6; ++i)
195  oldSum += oldPDF->xf( i, x, Q2) + oldPDF->xf( -i, x, Q2);
196  effFlinOld.fill ( x, oldSum );
197 
198  // Evaluate new summed PDF.
199  double newSum = 2.25 * newPDF->xf( 21, x, Q2);
200  for (int i = 1; i < 6; ++i)
201  newSum += newPDF->xf( i, x, Q2) + newPDF->xf( -i, x, Q2);
202  effFlinNew.fill ( x, newSum );
203 
204  //End loop over x values, in a linear scale.
205  }
206 
207  // Loop over x values, in a logarithmic scale
208  for (int iX = 0; iX < 80; ++iX) {
209  double xLog = -(0.1 * iX + 0.05);
210  double x = pow( 10., xLog);
211 
212  // Evaluate old summed PDF.
213  double oldSum = 2.25 * oldPDF->xf( 21, x, Q2);
214  for (int i = 1; i < 6; ++i)
215  oldSum += oldPDF->xf( i, x, Q2) + oldPDF->xf( -i, x, Q2);
216  effFlogOld.fill ( xLog, oldSum );
217 
218  // Evaluate new summed PDF.
219  double newSum = 2.25 * newPDF->xf( 21, x, Q2);
220  for (int i = 1; i < 6; ++i)
221  newSum += newPDF->xf( i, x, Q2) + newPDF->xf( -i, x, Q2);
222  effFlogNew.fill ( xLog, newSum );
223 
224  //End loop over x values, in a logarithmic scale.
225  }
226 
227  // Take ratios of new to old distributions.
228  effFlinRat = effFlinNew / effFlinOld;
229  effFlogRat = effFlogNew / effFlogOld;
230 
231  // Print histograms.
232  cout << effFlinOld << effFlinNew << effFlinRat
233  << effFlogOld << effFlogNew << effFlogRat;
234 
235  // Done.
236  delete oldPDF;
237  delete newPDF;
238  return 0;
239 }