// main22.cc is a part of the PYTHIA event generator. // Copyright (C) 2008 Peter Skands, Torbjorn Sjostrand. // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details. // Please respect the MCnet Guidelines, see GUIDELINES for details. // This is a simple test program. // It illustrates how to run SUSY processes in Pythia8. // All input is specified in the main22.cmnd file. #include "Pythia.h" using namespace Pythia8; int main() { // Generator. Shorthand for the event and the (static) Settings. Pythia pythia; Event& event = pythia.event; Settings& settings = pythia.settings; // Read in commands from external file. pythia.readFile("main22.cmnd"); // Extract settings to be used in the main program. int nEvent = settings.mode("Main:numberOfEvents"); int nList = settings.mode("Main:numberToList"); int nShow = settings.mode("Main:timesToShow"); bool showChangedSettings = settings.flag("Main:showChangedSettings"); bool showAllSettings = settings.flag("Main:showAllSettings"); bool showChangedParticleData = settings.flag("Main:showChangedParticleData"); bool showAllParticleData = settings.flag("Main:showAllParticleData"); double eCM = settings.parm("Beams:eCM"); // Initialize. Beam parameters set in .cmnd file. pythia.init(); // List changed data. if (showChangedSettings) settings.listChanged(); if (showAllSettings) settings.listAll(); // List particle data. if (showChangedParticleData) ParticleDataTable::listChanged(); if (showAllParticleData) ParticleDataTable::listAll(); // Histograms. double epTol = 1e-6 * eCM; Hist epCons("deviation from energy-momentum conservation",100,0.,epTol); Hist nFinal("final particle multiplicity",100,-0.5,799.5); Hist dnparticledy("dn/dy for particles",100,-10.,10.); Hist pxsum1("px1sum", 50000, -10., 10.); Hist pysum1("py1sum", 50000, -10., 10.); Hist pzsum1("pz1sum", 50000, -10., 10.); Hist ET1("missing pt1", 50000, 0., 100.); Hist pxsum2("px2sum", 50000, -100., 100.); Hist pysum2("py2sum", 50000, -100., 100.); Hist pzsum2("pz2sum", 50000, -100., 100.); Hist ET2("missing pt2", 50000, 0., 100.); Hist pxsum3("px3sum", 50000, -100., 100.); Hist pysum3("py3sum", 50000, -100., 100.); Hist pzsum3("pz3sum", 50000, -100., 100.); Hist ET3("missing pt3", 50000, 0., 100.); // Begin event loop. int nPace = max(1,nEvent/nShow); for (int iEvent = 0; iEvent < nEvent; ++iEvent) { if (iEvent%nPace == 0) cout << " Now begin event " << iEvent << endl; // Generate events. Quit if failure. if (!pythia.next()) { cout << " Event generation aborted prematurely, owing to error!\n"; break; } // List first few events, both hard process and complete events. if (iEvent < nList) { pythia.process.list(); event.list(); } double pxs1=0., pys1=0., pzs1=0.;//, es1=0.; double pxs2=0., pys2=0., pzs2=0.;//, es2=0.; double pxs3=0., pys3=0., pzs3=0.;//, es3=0.; // Loop over final particles in the event. int nFin = 0; Vec4 pSum; for (int i = 0; i < event.size(); ++i){ if (event[i].isFinal()) { nFin++; pSum += event[i].p(); dnparticledy.fill(event[i].y()); } if (pythia.event[i].isFinal()) { pxs1 += pythia.event[i].px(); pys1 += pythia.event[i].py(); pzs1 += pythia.event[i].pz(); // es1 += sqrt(pythia.event[i].px()*pythia.event[i].px()+pythia.event[i].py()*pythia.event[i].py()); } if (pythia.event[i].isFinal() && abs(pythia.event[i].eta()) < 5) { pxs2 += pythia.event[i].px(); pys2 += pythia.event[i].py(); pzs2 += pythia.event[i].pz(); // es2 += sqrt(pythia.event[i].px()*pythia.event[i].px()+pythia.event[i].py()*pythia.event[i].py()); } if (pythia.event[i].isFinal() && abs(pythia.event[i].eta()) < 5 && pythia.event[i].isVisible()) { pxs3 += pythia.event[i].px(); pys3 += pythia.event[i].py(); pzs3 += pythia.event[i].pz(); // es3 += sqrt(pythia.event[i].px()*pythia.event[i].px()+pythia.event[i].py()*pythia.event[i].py()); } } // Check and print event with too big energy-momentum deviation. nFinal.fill(nFin); double epDev = abs(pSum.e() - eCM) + abs(pSum.px()) + abs(pSum.py()) + abs(pSum.pz()); epCons.fill(epDev); if (epDev > epTol) { cout << " Warning! Event with epDev = " << scientific << setprecision(4) << epDev << " now listed:"; event.list(); } pxsum1.fill(pxs1); pysum1.fill(pys1); pzsum1.fill(pzs1); ET1.fill(sqrt(pxs1*pxs1+pys1*pys1)); pxsum2.fill(pxs2); pysum2.fill(pys2); pzsum2.fill(pzs2); ET2.fill(sqrt(pxs2*pxs2+pys2*pys2)); pxsum3.fill(pxs3); pysum3.fill(pys3); pzsum3.fill(pzs3); ET3.fill(sqrt(pxs3*pxs3+pys3*pys3)); // End of event loop. } // Final statistics and histogram output. pythia.statistics(); cout << epCons << nFinal << dnparticledy << pxsum1 << pysum1 << pzsum1 << ET1 << pxsum2 << pysum2 << pzsum2 << ET2 << pxsum3 << pysum3 << pzsum3 << ET3; ET3.table(cout); return 0; } /* Pythia8107 with main22.spc | q qbar -> ~chi_10 ~chi_10 1201 | 195282 18237 18237 | 1.050e-11 4.361e-14 | | q qbar -> ~chi_10 ~chi_20 1202 | 20837 1806 1806 | 1.036e-12 1.466e-14 | | q qbar -> ~chi_10 ~chi_30 1203 | 14590 3534 3534 | 2.030e-12 1.970e-14 | | q qbar -> ~chi_10 ~chi_40 1204 | 4549 789 789 | 4.422e-13 9.488e-15 | | q qbar -> ~chi_20 ~chi_20 1205 | 249758 32059 32059 | 1.821e-11 5.781e-14 | | q qbar -> ~chi_20 ~chi_30 1206 | 42509 10542 10542 | 6.040e-12 3.347e-14 | | q qbar -> ~chi_20 ~chi_40 1207 | 19282 2964 2964 | 1.705e-12 1.809e-14 | | q qbar -> ~chi_30 ~chi_30 1208 | 14 1 1 | 1.512e-15 1.512e-15 | | q qbar -> ~chi_30 ~chi_40 1209 | 126415 29986 29986 | 1.707e-11 5.545e-14 | | q qbar -> ~chi_40 ~chi_40 1210 | 498 82 82 | 4.403e-14 2.890e-15 | Pythia 6.4.16 with same spectrum and neutralinos forced stable. ! 216 f + fbar -> ~chi1 + ~chi1 I 177200 1875296 I 1.030E-11 I ! 217 f + fbar -> ~chi2 + ~chi2 I 324265 2542319 I 1.882E-11 I ! 218 f + fbar -> ~chi3 + ~chi3 I 38 163 I 2.196E-15 I ! 219 f + fbar -> ~chi4 + ~chi4 I 789 4886 I 4.536E-14 I ! 220 f + fbar -> ~chi1 + ~chi2 I 17541 207626 I 1.017E-12 I ! 221 f + fbar -> ~chi1 + ~chi3 I 35838 146111 I 2.064E-12 I ! 222 f + fbar -> ~chi1 + ~chi4 I 7859 44297 I 4.502E-13 I ! 223 f + fbar -> ~chi2 + ~chi3 I 105194 426610 I 6.115E-12 I ! 224 f + fbar -> ~chi2 + ~chi4 I 30338 198133 I 1.751E-12 I ! 225 f + fbar -> ~chi3 + ~chi4 I 300938 1266278 I 1.746E-11 I */