source: HiSusy/trunk/Pythia8/pythia8170/examples/main04.cc @ 1

Last change on this file since 1 was 1, checked in by zerwas, 11 years ago

first import of structure, PYTHIA8 and DELPHES
File size: 8.2 KB
Line 
1// main04.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// This is a simple test program.
7// It illustrates how to generate and study "total cross section" processes,
8// i.e. elastic, single and double diffractive, and the "minimum-bias" rest.
9// All input is specified in the main06.cmnd file.
10// Note that the "total" cross section does NOT include
11// the Coulomb contribution to elastic scattering, as switched on here.
12
13#include "Pythia.h"
14
15using namespace Pythia8; 
16
17//==========================================================================
18
19int main() {
20
21  // Generator. Shorthand for the event.
22  Pythia pythia;
23  Event& event = pythia.event;
24
25  // Read in commands from external file.
26  pythia.readFile("main04.cmnd");   
27
28  // Extract settings to be used in the main program.
29  int    nEvent    = pythia.mode("Main:numberOfEvents");
30  int    nAbort    = pythia.mode("Main:timesAllowErrors");
31 
32  // Initialize.
33  pythia.init();
34
35  // Book histograms: multiplicities and mean transverse momenta.
36  Hist yChg("rapidity of charged particles; all",      100, -10., 10.);
37  Hist nChg("number of charged particles; all",        100, -0.5, 799.5);
38  Hist nChgSD("number of charged particles; single diffraction", 
39                                                       100, -0.5, 799.5);
40  Hist nChgDD("number of charged particles, double diffractive", 
41                                                       100, -0.5, 799.5);
42  Hist nChgCD("number of charged particles, central diffractive", 
43                                                       100, -0.5, 799.5);
44  Hist nChgND("number of charged particles, non-diffractive", 
45                                                       100, -0.5, 799.5);
46  Hist pTnChg("<pt>(n_charged) all",                   100, -0.5, 799.5);
47  Hist pTnChgSD("<pt>(n_charged) single diffraction",  100, -0.5, 799.5);
48  Hist pTnChgDD("<pt>(n_charged) double diffraction",  100, -0.5, 799.5);
49  Hist pTnChgCD("<pt>(n_charged) central diffraction", 100, -0.5, 799.5);
50  Hist pTnChgND("<pt>(n_charged) non-diffractive   ",  100, -0.5, 799.5);
51
52  // Book histograms: ditto as function of separate subsystem mass.
53  Hist mLogInel("log10(mass), by diffractive system",  100, 0., 5.);
54  Hist nChgmLog("<n_charged>(log10(mass))",            100, 0., 5.);
55  Hist pTmLog("<pT>_charged>(log10(mass))",            100, 0., 5.);
56
57  // Book histograms: elastic/diffractive.
58  Hist tSpecEl("elastic |t| spectrum",              100, 0., 1.);
59  Hist tSpecElLog("elastic log10(|t|) spectrum",    100, -5., 0.);
60  Hist tSpecSD("single diffractive |t| spectrum",   100, 0., 2.); 
61  Hist tSpecDD("double diffractive |t| spectrum",   100, 0., 5.); 
62  Hist tSpecCD("central diffractive |t| spectrum",  100, 0., 5.); 
63  Hist mSpec("diffractive mass spectrum",           100, 0., 100.); 
64  Hist mLogSpec("log10(diffractive mass spectrum)", 100, 0., 4.); 
65
66  // Book histograms: inelastic nondiffractive "minbias".
67  double pTmax = 20.;
68  double bMax  = 4.;
69  Hist pTspec("total pT_hard spectrum",             100, 0., pTmax); 
70  Hist pTspecND("nondiffractive pT_hard spectrum",  100, 0., pTmax); 
71  Hist bSpec("b impact parameter spectrum",         100, 0., bMax);
72  Hist enhanceSpec("b enhancement spectrum",        100, 0., 10.);
73  Hist number("number of interactions",             100, -0.5, 99.5);
74  Hist pTb1("pT spectrum for b < 0.5",              100, 0., pTmax); 
75  Hist pTb2("pT spectrum for 0.5 < b < 1",          100, 0., pTmax); 
76  Hist pTb3("pT spectrum for 1 < b < 1.5",          100, 0., pTmax); 
77  Hist pTb4("pT spectrum for 1.5 < b",              100, 0., pTmax); 
78  Hist bpT1("b spectrum for pT < 2",                100, 0., bMax);
79  Hist bpT2("b spectrum for 2 < pT < 5",            100, 0., bMax);
80  Hist bpT3("b spectrum for 5 < pT < 15",           100, 0., bMax);
81  Hist bpT4("b spectrum for 15 < pT",               100, 0., bMax);
82 
83  // Begin event loop.
84  int iAbort = 0; 
85  for (int iEvent = 0; iEvent < nEvent; ++iEvent) {
86
87    // Generate events. Quit if too many failures.
88    if (!pythia.next()) {
89      if (++iAbort < nAbort) continue;
90      cout << " Event generation aborted prematurely, owing to error!\n"; 
91      break;
92    }
93
94    // Extract event classification.
95    int code = pythia.info.code();
96   
97    // Charged multiplicity and mean pT: all and by event class.
98    int nch = 0;
99    double pTsum = 0.; 
100    for (int i = 1; i < event.size(); ++i)
101    if (event[i].isFinal() && event[i].isCharged()) {
102      yChg.fill( event[i].y() );
103      ++nch; 
104      pTsum += event[i].pT();
105    }
106    nChg.fill( nch );
107    if (nch > 0) pTnChg.fill( nch, pTsum/nch);
108    if (code == 103 || code == 104) {
109      nChgSD.fill( nch );
110      if (nch > 0) pTnChgSD.fill( nch, pTsum/nch);
111    } else if (code == 105) {
112      nChgDD.fill( nch );
113      if (nch > 0) pTnChgDD.fill( nch, pTsum/nch);
114    } else if (code == 106) {
115      nChgCD.fill( nch );
116      if (nch > 0) pTnChgCD.fill( nch, pTsum/nch);
117    } else if (code == 101) {
118      nChgND.fill( nch );
119      if (nch > 0) pTnChgND.fill( nch, pTsum/nch);
120      double mLog = log10( event[0].m() ); 
121      mLogInel.fill( mLog ); 
122      nChgmLog.fill( mLog, nch ); 
123      if (nch > 0) pTmLog.fill( mLog, pTsum / nch ); 
124    }
125
126    // Charged multiplicity and mean pT: per diffractive system.
127    for (int iDiff = 0; iDiff < 3; ++iDiff) 
128    if ( (iDiff == 0 && pythia.info.isDiffractiveA()) 
129      || (iDiff == 1 && pythia.info.isDiffractiveB()) 
130      || (iDiff == 2 && pythia.info.isDiffractiveC()) ) {
131      int ndiff = 0;
132      double pTdiff = 0.;
133      int nDoc = (iDiff < 2) ? 4 : 5;   
134      for (int i = nDoc + 1; i < event.size(); ++i) 
135      if (event[i].isFinal() && event[i].isCharged()) {
136        // Trace back final particle to see which system it comes from.
137        int k = i;
138        do k = event[k].mother1(); 
139        while (k > nDoc);
140        if (k == iDiff + 3) {
141          ++ndiff;
142          pTdiff += event[i].pT();
143        }
144      } 
145      // Study diffractive mass spectrum.
146      double mDiff = event[iDiff+3].m();
147      double mLog  = log10( mDiff); 
148      mLogInel.fill( mLog ); 
149      nChgmLog.fill( mLog, ndiff ); 
150      if (ndiff > 0) pTmLog.fill( mLog, pTdiff / ndiff ); 
151      mSpec.fill( mDiff ); 
152      mLogSpec.fill( mLog );
153    }
154
155    // Study pT spectrum of all hard collisions, no distinction.
156    double pT = pythia.info.pTHat();
157    pTspec.fill( pT );
158
159    // Study t distribution of elastic/diffractive events.
160    if (code > 101) {
161      double tAbs = abs(pythia.info.tHat());
162      if (code == 102) {
163        tSpecEl.fill(tAbs);
164        tSpecElLog.fill(log10(tAbs));
165      }
166      else if (code == 103 || code == 104) tSpecSD.fill(tAbs);
167      else if (code == 105) tSpecDD.fill(tAbs);
168      else if (code == 106) {
169        double t1Abs = abs( (event[3].p() - event[1].p()).m2Calc() );
170        double t2Abs = abs( (event[4].p() - event[2].p()).m2Calc() );
171        tSpecCD.fill(t1Abs);
172        tSpecCD.fill(t2Abs);
173      }
174
175    // Study nondiffractive inelastic events in (pT, b) space.
176    } else {
177      double b = pythia.info.bMPI();
178      double enhance = pythia.info.enhanceMPI();
179      int nMPI = pythia.info.nMPI();
180      pTspecND.fill( pT );
181      bSpec.fill( b );
182      enhanceSpec.fill( enhance );
183      number.fill( nMPI );
184      if (b < 0.5) pTb1.fill( pT );
185      else if (b < 1.0) pTb2.fill( pT );
186      else if (b < 1.5) pTb3.fill( pT );
187      else pTb4.fill( pT );
188      if (pT < 2.) bpT1.fill( b );
189      else if (pT < 5.) bpT2.fill( b );
190      else if (pT < 15.) bpT3.fill( b );
191      else bpT4.fill( b );
192    }
193
194  // End of event loop.
195  }
196
197  // Final statistics and histograms.
198  pythia.stat();
199  pTnChg   /= nChg;
200  pTnChgSD /= nChgSD;
201  pTnChgDD /= nChgDD;
202  pTnChgCD /= nChgCD;
203  pTnChgND /= nChgND;
204  nChgmLog /= mLogInel;
205  pTmLog   /= mLogInel;
206  cout << yChg << nChg << nChgSD << nChgDD << nChgCD << nChgND
207       << pTnChg << pTnChgSD << pTnChgDD << pTnChgCD << pTnChgND
208       << mLogInel << nChgmLog << pTmLog
209       << tSpecEl << tSpecElLog << tSpecSD << tSpecDD << tSpecCD
210       << mSpec << mLogSpec
211       << pTspec << pTspecND << bSpec << enhanceSpec << number
212       << pTb1 << pTb2 << pTb3 << pTb4 << bpT1 << bpT2 << bpT3 << bpT4;
213
214  // Done.
215  return 0;
216}
Note: See TracBrowser for help on using the repository browser.