StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
main14.cc
1 // main14.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2014 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 "Pythia8/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  // With default option Higgs:clipWings = on need to reset mass range.
153  pythia.readString("25:mMin = 50.");
154  pythia.readString("25:mMax = 0.");
155  }
156 
157  // Initialize for LHC.
158  pythia.readString("Beams:eCM = 14000.");
159  pythia.init();
160 
161  // Debug: show initialized resonance data first time around.
162  //if (iProc == iFirst) pythia.particleData.listChanged(true);
163 
164  // Generate events to get cross section statistics.
165  for (int iEvent = 0; iEvent < nEvent; ++iEvent) pythia.next();
166 
167  // Show statistics.
168  //pythia.stat();
169  double sigma = pythia.info.sigmaGen();
170  cout << " Cross section is " << scientific << setprecision(3)
171  << sigma << " and in Pythia6 was " << sigma6[iProc]
172  << ",\n i.e. now is factor >>> " << fixed
173  << sigma / sigma6[iProc] << " <<< different" <<endl;
174 
175  // End of loop over processes.
176  }
177 
178  // Done.
179  return 0;
180 }