1 | // main08.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 methods to emphasize generation at high pT. |
---|
8 | |
---|
9 | #include "Pythia.h" |
---|
10 | |
---|
11 | using namespace Pythia8; |
---|
12 | |
---|
13 | int main() { |
---|
14 | |
---|
15 | // Different modes are illustrated for setting the pT ranges. |
---|
16 | // 1 : Hardcoded in the main program. |
---|
17 | // 2 : Using the Main:subrun keyword in a separate command file. |
---|
18 | // A third method instead biases selection continuously. |
---|
19 | // 3 : Bias high-pT selection by a pT^4 factor. |
---|
20 | // Matching also to low-pT processes is more complicated. |
---|
21 | // 4 : Matching between low- and high-pT. (No diffraction.) |
---|
22 | // 5: As 4, but bias high-pT selection by a pT^4 factor. |
---|
23 | int mode = 5; |
---|
24 | |
---|
25 | // Number of events to generate per bin. |
---|
26 | int nEvent = 100000; |
---|
27 | |
---|
28 | // One does not need complete events to study pThard spectrum only. |
---|
29 | bool completeEvents = false; |
---|
30 | |
---|
31 | // Optionally minimize output (almost) to final results. |
---|
32 | bool smallOutput = true; |
---|
33 | |
---|
34 | // Book histograms. |
---|
35 | int nRange = 100; |
---|
36 | double pTrange = (mode < 4) ? 1000. : 100.; |
---|
37 | Hist pTraw("pTHat distribution, unweighted", nRange, 0., pTrange); |
---|
38 | Hist pTnorm("pTHat distribution, weighted", nRange, 0., pTrange); |
---|
39 | Hist pTpow3("pTHat distribution, pT3*weighted", nRange, 0., pTrange); |
---|
40 | Hist pTpow5("pTHat distribution, pT5*weighted", nRange, 0., pTrange); |
---|
41 | Hist pTnormPart("pTHat distribution, weighted", nRange, 0., pTrange); |
---|
42 | Hist pTpow3Part("pTHat distribution, pT3*weighted", nRange, 0., pTrange); |
---|
43 | Hist pTpow5Part("pTHat distribution, pT5*weighted", nRange, 0., pTrange); |
---|
44 | |
---|
45 | // Generator. |
---|
46 | Pythia pythia; |
---|
47 | |
---|
48 | // Shorthand for some public members of pythia (also static ones). |
---|
49 | Settings& settings = pythia.settings; |
---|
50 | Info& info = pythia.info; |
---|
51 | |
---|
52 | // Optionally limit output to minimal one. |
---|
53 | if (smallOutput) { |
---|
54 | pythia.readString("Init:showProcesses = off"); |
---|
55 | pythia.readString("Init:showMultipartonInteractions = off"); |
---|
56 | pythia.readString("Init:showChangedSettings = off"); |
---|
57 | pythia.readString("Init:showChangedParticleData = off"); |
---|
58 | pythia.readString("Next:numberCount = 1000000000"); |
---|
59 | pythia.readString("Next:numberShowInfo = 0"); |
---|
60 | pythia.readString("Next:numberShowProcess = 0"); |
---|
61 | pythia.readString("Next:numberShowEvent = 0"); |
---|
62 | } |
---|
63 | |
---|
64 | // Number of bins to use. In mode 2 read from main08.cmnd file. |
---|
65 | int nBin = 5; |
---|
66 | if (mode == 2) { |
---|
67 | pythia.readFile("main08.cmnd"); |
---|
68 | nBin = pythia.mode("Main:numberOfSubruns"); |
---|
69 | } |
---|
70 | else if (mode == 3) nBin = 1; |
---|
71 | else if (mode == 4) nBin = 4; |
---|
72 | else if (mode == 5) nBin = 2; |
---|
73 | |
---|
74 | // Mode 1: set up five pT bins - last one open-ended. |
---|
75 | double pTlimit[6] = {100., 150., 250., 400., 600., 0.}; |
---|
76 | |
---|
77 | // Modes 4 & 5: set up pT bins for range [0, 100]. The lowest bin |
---|
78 | // is generated with soft processes, to regularize pT -> 0 blowup. |
---|
79 | // Warning: if pTlimitLow[1] is picked too low there will be a |
---|
80 | // visible discontinuity, since soft processes are generated with |
---|
81 | // dampening and "Sudakov" for pT -> 0, while hard processes are not. |
---|
82 | double pTlimitLow[6] = {0., 20., 40., 70., 100.}; |
---|
83 | double pTlimitTwo[3] = {0., 20., 100.}; |
---|
84 | |
---|
85 | // Loop over number of bins, i.e. number of subruns. |
---|
86 | for (int iBin = 0; iBin < nBin; ++iBin) { |
---|
87 | |
---|
88 | // Normally HardQCD, but in two cases minBias. |
---|
89 | // Need MPI on in minBias to get first interaction, but not else. |
---|
90 | if (mode > 3 && iBin == 0) { |
---|
91 | pythia.readString("HardQCD:all = off"); |
---|
92 | pythia.readString("SoftQCD:minBias = on"); |
---|
93 | if (!completeEvents) { |
---|
94 | pythia.readString("PartonLevel:all = on"); |
---|
95 | pythia.readString("PartonLevel:ISR = off"); |
---|
96 | pythia.readString("PartonLevel:FSR = off"); |
---|
97 | pythia.readString("HadronLevel:all = off"); |
---|
98 | } |
---|
99 | } else { |
---|
100 | pythia.readString("HardQCD:all = on"); |
---|
101 | pythia.readString("SoftQCD:minBias = off"); |
---|
102 | if (!completeEvents) pythia.readString("PartonLevel:all = off"); |
---|
103 | } |
---|
104 | |
---|
105 | // Mode 1: hardcoded here. Use settings.parm for non-string input. |
---|
106 | if (mode == 1) { |
---|
107 | settings.parm("PhaseSpace:pTHatMin", pTlimit[iBin]); |
---|
108 | settings.parm("PhaseSpace:pTHatMax", pTlimit[iBin + 1]); |
---|
109 | } |
---|
110 | |
---|
111 | // Mode 2: subruns stored in the main08.cmnd file. |
---|
112 | else if (mode == 2) pythia.readFile("main08.cmnd", iBin); |
---|
113 | |
---|
114 | // Mode 3: The whole range in one step, but pT-weighted. |
---|
115 | else if (mode == 3) { |
---|
116 | settings.parm("PhaseSpace:pTHatMin", pTlimit[0]); |
---|
117 | settings.parm("PhaseSpace:pTHatMax", 0.); |
---|
118 | pythia.readString("PhaseSpace:bias2Selection = on"); |
---|
119 | pythia.readString("PhaseSpace:bias2SelectionPow = 4."); |
---|
120 | pythia.readString("PhaseSpace:bias2SelectionRef = 100."); |
---|
121 | } |
---|
122 | |
---|
123 | // Mode 4: hardcoded here. Use settings.parm for non-string input. |
---|
124 | else if (mode == 4) { |
---|
125 | settings.parm("PhaseSpace:pTHatMin", pTlimitLow[iBin]); |
---|
126 | settings.parm("PhaseSpace:pTHatMax", pTlimitLow[iBin + 1]); |
---|
127 | } |
---|
128 | |
---|
129 | // Mode 5: hardcoded here. Use settings.parm for non-string input. |
---|
130 | // Hard processes in one step, but pT-weighted. |
---|
131 | else if (mode == 5) { |
---|
132 | settings.parm("PhaseSpace:pTHatMin", pTlimitTwo[iBin]); |
---|
133 | settings.parm("PhaseSpace:pTHatMax", pTlimitTwo[iBin + 1]); |
---|
134 | if (iBin == 1) { |
---|
135 | pythia.readString("PhaseSpace:bias2Selection = on"); |
---|
136 | pythia.readString("PhaseSpace:bias2SelectionPow = 4."); |
---|
137 | pythia.readString("PhaseSpace:bias2SelectionRef = 20."); |
---|
138 | } |
---|
139 | } |
---|
140 | |
---|
141 | // Initialize for LHC at 14 TeV. |
---|
142 | pythia.readString("Beams:eCM = 14000."); |
---|
143 | pythia.init(); |
---|
144 | |
---|
145 | // Reset local histograms (that need to be rescaled before added). |
---|
146 | pTnormPart.null(); |
---|
147 | pTpow3Part.null(); |
---|
148 | pTpow5Part.null(); |
---|
149 | |
---|
150 | // Begin event loop. |
---|
151 | for (int iEvent = 0; iEvent < nEvent; ++iEvent) { |
---|
152 | |
---|
153 | // Generate events. Quit if failure. |
---|
154 | if (!pythia.next()) break; |
---|
155 | |
---|
156 | // Soft events have no upper pT limit. They therefore overlap |
---|
157 | // with hard events, and the overlap must be removed by hand. |
---|
158 | // No overlap for elastic/diffraction, which is only part of soft. |
---|
159 | double pTHat = info.pTHat(); |
---|
160 | if (mode > 3 && iBin == 0 && info.isMinBias() |
---|
161 | && pTHat > pTlimitLow[1]) continue; |
---|
162 | |
---|
163 | // Fill hard scale of event. |
---|
164 | double weight = info.weight(); |
---|
165 | pTraw.fill( pTHat ); |
---|
166 | pTnormPart.fill( pTHat, weight); |
---|
167 | pTpow3Part.fill( pTHat, weight * pow3(pTHat) ); |
---|
168 | pTpow5Part.fill( pTHat, weight * pow5(pTHat) ); |
---|
169 | |
---|
170 | // End of event loop. Statistics. |
---|
171 | } |
---|
172 | if (!smallOutput) pythia.stat(); |
---|
173 | |
---|
174 | // Normalize to cross section for each case, and add to sum. |
---|
175 | double sigmaNorm = (info.sigmaGen() / info.weightSum()) |
---|
176 | * (nRange / pTrange); |
---|
177 | pTnorm += sigmaNorm * pTnormPart; |
---|
178 | pTpow3 += sigmaNorm * pTpow3Part; |
---|
179 | pTpow5 += sigmaNorm * pTpow5Part; |
---|
180 | |
---|
181 | // End of pT-bin loop. |
---|
182 | } |
---|
183 | |
---|
184 | // Output histograms. |
---|
185 | cout << pTraw << pTnorm << pTpow3 << pTpow5; |
---|
186 | |
---|
187 | // Done. |
---|
188 | return 0; |
---|
189 | } |
---|