Update to pythi8.170
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8170 / examples / main52.cc
1 // main52.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 // 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 "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:minBias = 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   return 0;
237 }