[1] | 1 | // main14.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 | // Comparison with some PYTHIA 6.413 cross sections process by process. |
---|
| 7 | // Several processes have been left out to keep reasonable execution time. |
---|
| 8 | // Some processes are not handled absolutely identically, so minor |
---|
| 9 | // systematic differences may occur in addition to the statistical ones. |
---|
| 10 | // (For some MSSM Higgs processes 6.413 has been modified to use |
---|
| 11 | // running quark masses in loops, like 8.1, to allow proper comparison.) |
---|
| 12 | // Subruns 0 - 5 : QCD jets |
---|
| 13 | // 6 - 10 : prompt photons. |
---|
| 14 | // 11 - 12 : t-channel gamma/Z/W exchange. |
---|
| 15 | // 13 - 23 : gamma*/Z^0/W^+-, singly, in pairs or with parton |
---|
| 16 | // 24 - 25 : onia. |
---|
| 17 | // 26 - 30 : top. |
---|
| 18 | // 31 - 40 : Standard Model Higgs. |
---|
| 19 | // 41 - 45 : MSSM Higgses (trivial couplings). |
---|
| 20 | // 46 - 47 : Z' and W' |
---|
| 21 | // 48 - 51 : Left-right-symmetric scenario. |
---|
| 22 | // 52 - 52 : Leptoquark. |
---|
| 23 | // 53 - 55 : Excited fermions (compositeness). |
---|
| 24 | // 56 - 56 : excited Graviton (RS extra dimensions). |
---|
| 25 | |
---|
| 26 | #include "Pythia.h" |
---|
| 27 | |
---|
| 28 | using namespace Pythia8; |
---|
| 29 | |
---|
| 30 | int main() { |
---|
| 31 | |
---|
| 32 | // First and last process to test: can run from 0 through 40. |
---|
| 33 | int iFirst = 0; |
---|
| 34 | int iLast = 56; |
---|
| 35 | |
---|
| 36 | // Statistics. Pythia6 run was with 10000, so no point to use more. |
---|
| 37 | int nEvent = 10000; |
---|
| 38 | |
---|
| 39 | // Normally one subprocess word per subrun, but exceptions exist. |
---|
| 40 | int nSub[100]; |
---|
| 41 | for (int i = 0; i < 100; ++i) nSub[i] = 1; |
---|
| 42 | nSub[1] = 3; |
---|
| 43 | nSub[5] = 3; |
---|
| 44 | // Starting positions in subprocess words list, recursively defined. |
---|
| 45 | int iBeg[101] = { 0 }; |
---|
| 46 | for (int i = 0; i < 100; ++i) iBeg[i + 1] = iBeg[i] + nSub[i]; |
---|
| 47 | // List of subprocess words. |
---|
| 48 | string processes[61] = { "HardQCD:gg2gg", "HardQCD:gg2qqbar", |
---|
| 49 | "HardQCD:gg2ccbar", "HardQCD:gg2bbbar","HardQCD:qg2qg" , |
---|
| 50 | "HardQCD:qq2qq", "HardQCD:qqbar2gg", "HardQCD:qqbar2qqbarNew", |
---|
| 51 | "HardQCD:qqbar2ccbar", "HardQCD:qqbar2bbbar", "PromptPhoton:qg2qgamma", |
---|
| 52 | "PromptPhoton:qqbar2ggamma", "PromptPhoton:gg2ggamma", |
---|
| 53 | "PromptPhoton:ffbar2gammagamma", "PromptPhoton:gg2gammagamma", |
---|
| 54 | "WeakBosonExchange:ff2ff(t:gmZ)", "WeakBosonExchange:ff2ff(t:W)", |
---|
| 55 | "WeakSingleBoson:ffbar2gmZ", "WeakSingleBoson:ffbar2W", |
---|
| 56 | "WeakDoubleBoson:ffbar2gmZgmZ", "WeakDoubleBoson:ffbar2ZW", |
---|
| 57 | "WeakDoubleBoson:ffbar2WW", "WeakBosonAndParton:qqbar2gmZg", |
---|
| 58 | "WeakBosonAndParton:qg2gmZq", "WeakBosonAndParton:ffbar2gmZgm", |
---|
| 59 | "WeakBosonAndParton:qqbar2Wg", "WeakBosonAndParton:qg2Wq", |
---|
| 60 | "WeakBosonAndParton:ffbar2Wgm", "Charmonium:all", "Bottomonium:all", |
---|
| 61 | "Top:gg2ttbar", "Top:qqbar2ttbar", "Top:qq2tq(t:W)", |
---|
| 62 | "Top:ffbar2ttbar(s:gmZ)", "Top:ffbar2tqbar(s:W)", |
---|
| 63 | "HiggsSM:ffbar2H", "HiggsSM:gg2H", "HiggsSM:ffbar2HZ", |
---|
| 64 | "HiggsSM:ffbar2HW", "HiggsSM:ff2Hff(t:ZZ)", "HiggsSM:ff2Hff(t:WW)", |
---|
| 65 | "HiggsSM:qg2Hq", "HiggsSM:gg2Hg(l:t)", "HiggsSM:qg2Hq(l:t)", |
---|
| 66 | "HiggsSM:qqbar2Hg(l:t)", "HiggsBSM:allH1", "HiggsBSM:allH2", |
---|
| 67 | "HiggsBSM:allA3", "HiggsBSM:allH+-", "HiggsBSM:allHpair", |
---|
| 68 | "NewGaugeBoson:ffbar2gmZZprime", "NewGaugeBoson:ffbar2Wprime", |
---|
| 69 | "LeftRightSymmmetry:ffbar2ZR", "LeftRightSymmmetry:ffbar2WR", |
---|
| 70 | "LeftRightSymmmetry:ffbar2HLHL", "LeftRightSymmmetry:ffbar2HRHR", |
---|
| 71 | "LeptoQuark:all", "ExcitedFermion:dg2dStar", |
---|
| 72 | "ExcitedFermion:qq2dStarq", "ExcitedFermion:qqbar2eStare", |
---|
| 73 | "ExtraDimensionsG*:all" }; |
---|
| 74 | |
---|
| 75 | // List of cross sections from Pythia6. |
---|
| 76 | double sigma6[57] = { 4.960e-01, 1.627e-02, 2.790e-01, 2.800e-02, |
---|
| 77 | 3.310e-04, 3.653e-04, 1.697e-04, 1.163e-05, 1.065e-07, 8.259e-08, |
---|
| 78 | 8.237e-08, 2.544e-05, 5.321e-06, 5.571e-05, 1.621e-04, 9.039e-09, |
---|
| 79 | 2.247e-08, 5.893e-08, 3.781e-06, 1.078e-05, 4.551e-08, 1.025e-05, |
---|
| 80 | 3.208e-05, 5.435e-08, 1.038e-04, 3.929e-05, 4.155e-07, 6.685e-08, |
---|
| 81 | 1.898e-07, 4.240e-10, 7.142e-09, 1.547e-10, 7.064e-09, 1.316e-10, |
---|
| 82 | 2.332e-10, 5.105e-10, 1.316e-09, 4.462e-11, 5.557e-09, 1.966e-09, |
---|
| 83 | 8.725e-12, 2.450e-08, 5.839e-09, 1.687e-08, 8.950e-11, 4.188e-11, |
---|
| 84 | 1.980e-07, 4.551e-07, 6.005e-09, 1.102e-07, 7.784e-11, 3.488e-11, |
---|
| 85 | 6.006e-08, 3.235e-06, 1.689e-05, 5.986e-07, 3.241e-10 }; |
---|
| 86 | |
---|
| 87 | // Generator. |
---|
| 88 | Pythia pythia; |
---|
| 89 | |
---|
| 90 | // Standard set of masses for comparison with Fortran code. |
---|
| 91 | pythia.readString("5:m0 = 4.2"); |
---|
| 92 | pythia.readString("6:m0 = 175."); |
---|
| 93 | pythia.readString("23:m0 = 91.2"); |
---|
| 94 | pythia.readString("24:m0 = 80."); |
---|
| 95 | |
---|
| 96 | // Same kinematics cuts as Fortran code. |
---|
| 97 | pythia.readString("PhaseSpace:pTHatMin = 20."); |
---|
| 98 | pythia.readString("6:mMin = 20."); |
---|
| 99 | pythia.readString("23:mMin = 20."); |
---|
| 100 | pythia.readString("24:mMin = 20."); |
---|
| 101 | pythia.readString("25:mMin = 20."); |
---|
| 102 | pythia.readString("32:mMin = 400."); |
---|
| 103 | pythia.readString("34:mMin = 400."); |
---|
| 104 | pythia.readString("42:mMin = 50."); |
---|
| 105 | pythia.readString("5000039:mMin = 50."); |
---|
| 106 | |
---|
| 107 | // Also same renormalization and factorization scale. |
---|
| 108 | pythia.readString("SigmaProcess:renormScale2 = 3"); |
---|
| 109 | pythia.readString("SigmaProcess:factorScale2 = 3"); |
---|
| 110 | |
---|
| 111 | // Switch off unnecessary parts. |
---|
| 112 | pythia.readString("PartonLevel:all = off"); |
---|
| 113 | pythia.readString("ProcessLevel:resonanceDecays = off"); |
---|
| 114 | |
---|
| 115 | // No printing of settings, particle data or events. |
---|
| 116 | pythia.readString("Init:showProcesses = off"); |
---|
| 117 | pythia.readString("Init:showChangedSettings = off"); |
---|
| 118 | pythia.readString("Init:showChangedParticleData = off"); |
---|
| 119 | pythia.readString("Next:numberCount = 0"); |
---|
| 120 | pythia.readString("Next:numberShowInfo = 0"); |
---|
| 121 | pythia.readString("Next:numberShowProcess = 0"); |
---|
| 122 | pythia.readString("Next:numberShowEvent = 0"); |
---|
| 123 | |
---|
| 124 | // Debug: show information on cross section maximum and violation. |
---|
| 125 | //pythia.readString("PhaseSpace:showSearch = on"); |
---|
| 126 | //pythia.readString("PhaseSpace:showViolation = on"); |
---|
| 127 | |
---|
| 128 | // Loop over processes. |
---|
| 129 | for (int iProc = iFirst; iProc <= iLast; ++iProc) { |
---|
| 130 | cout << "\n Begin subrun number " << iProc << " : "; |
---|
| 131 | |
---|
| 132 | // Switch off previous process(es) and switch on new one(s). |
---|
| 133 | if (iProc > iFirst) for (int i = iBeg[iProc - 1]; i < iBeg[iProc]; ++i) |
---|
| 134 | pythia.readString( processes[i] + " = off" ); |
---|
| 135 | for (int i = iBeg[iProc]; i < iBeg[iProc + 1]; ++i) { |
---|
| 136 | pythia.readString( processes[i] + " = on" ); |
---|
| 137 | if (i > iBeg[iProc]) cout << " + "; |
---|
| 138 | cout << processes[i]; |
---|
| 139 | } |
---|
| 140 | cout << endl; |
---|
| 141 | |
---|
| 142 | // Switch between SM and MSSM Higgs scenario. |
---|
| 143 | if (iProc <= 40) { |
---|
| 144 | pythia.readString("Higgs:useBSM = off"); |
---|
| 145 | pythia.readString("25:m0 = 200."); |
---|
| 146 | } else { |
---|
| 147 | pythia.readString("Higgs:useBSM = on"); |
---|
| 148 | pythia.readString("25:m0 = 115."); |
---|
| 149 | pythia.readString("35:m0 = 300."); |
---|
| 150 | pythia.readString("36:m0 = 300."); |
---|
| 151 | pythia.readString("37:m0 = 320."); |
---|
| 152 | } |
---|
| 153 | |
---|
| 154 | // Initialize for LHC. |
---|
| 155 | pythia.readString("Beams:eCM = 14000."); |
---|
| 156 | pythia.init(); |
---|
| 157 | |
---|
| 158 | // Debug: show initialized resonance data first time around. |
---|
| 159 | //if (iProc == iFirst) pythia.particleData.listChanged(true); |
---|
| 160 | |
---|
| 161 | // Generate events to get cross section statistics. |
---|
| 162 | for (int iEvent = 0; iEvent < nEvent; ++iEvent) pythia.next(); |
---|
| 163 | |
---|
| 164 | // Show statistics. |
---|
| 165 | //pythia.stat(); |
---|
| 166 | double sigma = pythia.info.sigmaGen(); |
---|
| 167 | cout << " Cross section is " << scientific << setprecision(3) |
---|
| 168 | << sigma << " and in Pythia6 was " << sigma6[iProc] |
---|
| 169 | << ",\n i.e. now is factor >>> " << fixed |
---|
| 170 | << sigma / sigma6[iProc] << " <<< different" <<endl; |
---|
| 171 | |
---|
| 172 | // End of loop over processes. |
---|
| 173 | } |
---|
| 174 | |
---|
| 175 | // Done. |
---|
| 176 | return 0; |
---|
| 177 | } |
---|