[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8175 / examples / main04.cc

// main04.cc is a part of the PYTHIA event generator.
// Copyright (C) 2013 Torbjorn Sjostrand.
// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
// Please respect the MCnet Guidelines, see GUIDELINES for details.

// This is a simple test program. 
// It illustrates how to generate and study "total cross section" processes,
// i.e. elastic, single and double diffractive, and the "minimum-bias" rest.
// All input is specified in the main06.cmnd file.
// Note that the "total" cross section does NOT include 
// the Coulomb contribution to elastic scattering, as switched on here.

#include "Pythia.h"

using namespace Pythia8; 

//==========================================================================

int main() {

  // Generator. Shorthand for the event.
  Pythia pythia;
  Event& event = pythia.event;

  // Read in commands from external file.
  pythia.readFile("main04.cmnd");    

  // Extract settings to be used in the main program.
  int    nEvent    = pythia.mode("Main:numberOfEvents");
  int    nAbort    = pythia.mode("Main:timesAllowErrors");
 
  // Initialize. 
  pythia.init();

  // Book histograms: multiplicities and mean transverse momenta.
  Hist yChg("rapidity of charged particles; all",      100, -10., 10.);
  Hist nChg("number of charged particles; all",        100, -0.5, 799.5);
  Hist nChgSD("number of charged particles; single diffraction", 
                                                       100, -0.5, 799.5);
  Hist nChgDD("number of charged particles, double diffractive", 
                                                       100, -0.5, 799.5);
  Hist nChgCD("number of charged particles, central diffractive", 
                                                       100, -0.5, 799.5);
  Hist nChgND("number of charged particles, non-diffractive", 
                                                       100, -0.5, 799.5);
  Hist pTnChg("<pt>(n_charged) all",                   100, -0.5, 799.5);
  Hist pTnChgSD("<pt>(n_charged) single diffraction",  100, -0.5, 799.5);
  Hist pTnChgDD("<pt>(n_charged) double diffraction",  100, -0.5, 799.5);
  Hist pTnChgCD("<pt>(n_charged) central diffraction", 100, -0.5, 799.5);
  Hist pTnChgND("<pt>(n_charged) non-diffractive   ",  100, -0.5, 799.5);

  // Book histograms: ditto as function of separate subsystem mass.
  Hist mLogInel("log10(mass), by diffractive system",  100, 0., 5.);
  Hist nChgmLog("<n_charged>(log10(mass))",            100, 0., 5.);
  Hist pTmLog("<pT>_charged>(log10(mass))",            100, 0., 5.);

  // Book histograms: elastic/diffractive.
  Hist tSpecEl("elastic |t| spectrum",              100, 0., 1.);
  Hist tSpecElLog("elastic log10(|t|) spectrum",    100, -5., 0.);
  Hist tSpecSD("single diffractive |t| spectrum",   100, 0., 2.); 
  Hist tSpecDD("double diffractive |t| spectrum",   100, 0., 5.); 
  Hist tSpecCD("central diffractive |t| spectrum",  100, 0., 5.); 
  Hist mSpec("diffractive mass spectrum",           100, 0., 100.); 
  Hist mLogSpec("log10(diffractive mass spectrum)", 100, 0., 4.); 

  // Book histograms: inelastic nondiffractive "minbias".
  double pTmax = 20.;
  double bMax  = 4.;
  Hist pTspec("total pT_hard spectrum",             100, 0., pTmax); 
  Hist pTspecND("nondiffractive pT_hard spectrum",  100, 0., pTmax); 
  Hist bSpec("b impact parameter spectrum",         100, 0., bMax);
  Hist enhanceSpec("b enhancement spectrum",        100, 0., 10.);
  Hist number("number of interactions",             100, -0.5, 99.5);
  Hist pTb1("pT spectrum for b < 0.5",              100, 0., pTmax); 
  Hist pTb2("pT spectrum for 0.5 < b < 1",          100, 0., pTmax); 
  Hist pTb3("pT spectrum for 1 < b < 1.5",          100, 0., pTmax); 
  Hist pTb4("pT spectrum for 1.5 < b",              100, 0., pTmax); 
  Hist bpT1("b spectrum for pT < 2",                100, 0., bMax);
  Hist bpT2("b spectrum for 2 < pT < 5",            100, 0., bMax);
  Hist bpT3("b spectrum for 5 < pT < 15",           100, 0., bMax);
  Hist bpT4("b spectrum for 15 < pT",               100, 0., bMax);
 
  // Begin event loop.
  int iAbort = 0; 
  for (int iEvent = 0; iEvent < nEvent; ++iEvent) {

    // Generate events. Quit if too many failures.
    if (!pythia.next()) {
      if (++iAbort < nAbort) continue;
      cout << " Event generation aborted prematurely, owing to error!\n"; 
      break;
    }

    // Extract event classification.
    int code = pythia.info.code();
    
    // Charged multiplicity and mean pT: all and by event class.
    int nch = 0;
    double pTsum = 0.; 
    for (int i = 1; i < event.size(); ++i)
    if (event[i].isFinal() && event[i].isCharged()) {
      yChg.fill( event[i].y() );
      ++nch; 
      pTsum += event[i].pT();
    }
    nChg.fill( nch );
    if (nch > 0) pTnChg.fill( nch, pTsum/nch);
    if (code == 103 || code == 104) {
      nChgSD.fill( nch );
      if (nch > 0) pTnChgSD.fill( nch, pTsum/nch);
    } else if (code == 105) {
      nChgDD.fill( nch );
      if (nch > 0) pTnChgDD.fill( nch, pTsum/nch);
    } else if (code == 106) {
      nChgCD.fill( nch );
      if (nch > 0) pTnChgCD.fill( nch, pTsum/nch);
    } else if (code == 101) {
      nChgND.fill( nch );
      if (nch > 0) pTnChgND.fill( nch, pTsum/nch);
      double mLog = log10( event[0].m() ); 
      mLogInel.fill( mLog );  
      nChgmLog.fill( mLog, nch );  
      if (nch > 0) pTmLog.fill( mLog, pTsum / nch );  
    }

    // Charged multiplicity and mean pT: per diffractive system.
    for (int iDiff = 0; iDiff < 3; ++iDiff) 
    if ( (iDiff == 0 && pythia.info.isDiffractiveA()) 
      || (iDiff == 1 && pythia.info.isDiffractiveB()) 
      || (iDiff == 2 && pythia.info.isDiffractiveC()) ) {
      int ndiff = 0;
      double pTdiff = 0.;
      int nDoc = (iDiff < 2) ? 4 : 5;   
      for (int i = nDoc + 1; i < event.size(); ++i) 
      if (event[i].isFinal() && event[i].isCharged()) {
        // Trace back final particle to see which system it comes from.
        int k = i;
        do k = event[k].mother1(); 
        while (k > nDoc);
        if (k == iDiff + 3) {
          ++ndiff;
          pTdiff += event[i].pT();
        }
      } 
      // Study diffractive mass spectrum.
      double mDiff = event[iDiff+3].m();
      double mLog  = log10( mDiff);  
      mLogInel.fill( mLog );  
      nChgmLog.fill( mLog, ndiff );  
      if (ndiff > 0) pTmLog.fill( mLog, pTdiff / ndiff );  
      mSpec.fill( mDiff );  
      mLogSpec.fill( mLog );
    }

    // Study pT spectrum of all hard collisions, no distinction.
    double pT = pythia.info.pTHat();
    pTspec.fill( pT );

    // Study t distribution of elastic/diffractive events.
    if (code > 101) {
      double tAbs = abs(pythia.info.tHat());
      if (code == 102) {
        tSpecEl.fill(tAbs);
        tSpecElLog.fill(log10(tAbs));
      }
      else if (code == 103 || code == 104) tSpecSD.fill(tAbs);
      else if (code == 105) tSpecDD.fill(tAbs);
      else if (code == 106) {
        double t1Abs = abs( (event[3].p() - event[1].p()).m2Calc() );
        double t2Abs = abs( (event[4].p() - event[2].p()).m2Calc() );
        tSpecCD.fill(t1Abs);
        tSpecCD.fill(t2Abs);
      }

    // Study nondiffractive inelastic events in (pT, b) space.
    } else {
      double b = pythia.info.bMPI();
      double enhance = pythia.info.enhanceMPI();
      int nMPI = pythia.info.nMPI();
      pTspecND.fill( pT );
      bSpec.fill( b );
      enhanceSpec.fill( enhance );
      number.fill( nMPI );
      if (b < 0.5) pTb1.fill( pT );
      else if (b < 1.0) pTb2.fill( pT );
      else if (b < 1.5) pTb3.fill( pT );
      else pTb4.fill( pT );
      if (pT < 2.) bpT1.fill( b );
      else if (pT < 5.) bpT2.fill( b );
      else if (pT < 15.) bpT3.fill( b );
      else bpT4.fill( b );
    }

  // End of event loop.
  }

  // Final statistics and histograms.
  pythia.stat();
  pTnChg   /= nChg;
  pTnChgSD /= nChgSD;
  pTnChgDD /= nChgDD;
  pTnChgCD /= nChgCD;
  pTnChgND /= nChgND;
  nChgmLog /= mLogInel;
  pTmLog   /= mLogInel;
  cout << yChg << nChg << nChgSD << nChgDD << nChgCD << nChgND
       << pTnChg << pTnChgSD << pTnChgDD << pTnChgCD << pTnChgND
       << mLogInel << nChgmLog << pTmLog
       << tSpecEl << tSpecElLog << tSpecSD << tSpecDD << tSpecCD
       << mSpec << mLogSpec 
       << pTspec << pTspecND << bSpec << enhanceSpec << number
       << pTb1 << pTb2 << pTb3 << pTb4 << bpT1 << bpT2 << bpT3 << bpT4;

  // Done.
  return 0;
}
Commit	Line	Data
c6b60c38	1	// main04.cc is a part of the PYTHIA event generator.
	2	// Copyright (C) 2013 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
	15	using namespace Pythia8;
	16
	17	//==========================================================================
	18
	19	int 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	}