[1] | 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 | |
---|
| 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 | } |
---|