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