1 | // main07.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 | // Illustration how to generate various two-body channels from |
---|
7 | // astroparticle processes, e.g. neutralino annihilation or decay. |
---|
8 | // To this end a "blob" of energy is created with unit cross section, |
---|
9 | // from the fictitious collision of two non-radiating incoming e+e-. |
---|
10 | // In the accompanying main29.cmnd file the decay channels of this |
---|
11 | // blob can be set up. Furthermore, only gamma, e+-, p/pbar and |
---|
12 | // neutrinos are stable, everything else is set to decay. |
---|
13 | // (The "single-particle gun" of main21.cc offers another possible |
---|
14 | // approach to the same problem.) |
---|
15 | |
---|
16 | #include "Pythia.h" |
---|
17 | |
---|
18 | using namespace Pythia8; |
---|
19 | |
---|
20 | //========================================================================== |
---|
21 | |
---|
22 | // A derived class for (e+ e- ->) GenericResonance -> various final states. |
---|
23 | |
---|
24 | class Sigma1GenRes : public Sigma1Process { |
---|
25 | |
---|
26 | public: |
---|
27 | |
---|
28 | // Constructor. |
---|
29 | Sigma1GenRes() {} |
---|
30 | |
---|
31 | // Evaluate sigmaHat(sHat): dummy unit cross section. |
---|
32 | virtual double sigmaHat() {return 1.;} |
---|
33 | |
---|
34 | // Select flavour. No colour or anticolour. |
---|
35 | virtual void setIdColAcol() {setId( -11, 11, 999999); |
---|
36 | setColAcol( 0, 0, 0, 0, 0, 0);} |
---|
37 | |
---|
38 | // Info on the subprocess. |
---|
39 | virtual string name() const {return "GenericResonance";} |
---|
40 | virtual int code() const {return 9001;} |
---|
41 | virtual string inFlux() const {return "ffbarSame";} |
---|
42 | |
---|
43 | }; |
---|
44 | |
---|
45 | //========================================================================== |
---|
46 | |
---|
47 | int main() { |
---|
48 | |
---|
49 | // Pythia generator. |
---|
50 | Pythia pythia; |
---|
51 | |
---|
52 | // A class to generate the fictitious resonance initial state. |
---|
53 | SigmaProcess* sigma1GenRes = new Sigma1GenRes(); |
---|
54 | |
---|
55 | // Hand pointer to Pythia. |
---|
56 | pythia.setSigmaPtr( sigma1GenRes); |
---|
57 | |
---|
58 | // Read in the rest of the settings and data from a separate file. |
---|
59 | pythia.readFile("main07.cmnd"); |
---|
60 | |
---|
61 | // Initialization. |
---|
62 | pythia.init(); |
---|
63 | |
---|
64 | // Extract settings to be used in the main program. |
---|
65 | int nEvent = pythia.mode("Main:numberOfEvents"); |
---|
66 | int nAbort = pythia.mode("Main:timesAllowErrors"); |
---|
67 | |
---|
68 | // Histogram particle spectra. |
---|
69 | Hist eGamma("energy spectrum of photons", 100, 0., 250.); |
---|
70 | Hist eE( "energy spectrum of e+ and e-", 100, 0., 250.); |
---|
71 | Hist eP( "energy spectrum of p and pbar", 100, 0., 250.); |
---|
72 | Hist eNu( "energy spectrum of neutrinos", 100, 0., 250.); |
---|
73 | Hist eRest( "energy spectrum of rest particles", 100, 0., 250.); |
---|
74 | |
---|
75 | // Begin event loop. |
---|
76 | int iAbort = 0; |
---|
77 | for (int iEvent = 0; iEvent < nEvent; ++iEvent) { |
---|
78 | |
---|
79 | // Generate events. Quit if many failures. |
---|
80 | if (!pythia.next()) { |
---|
81 | if (++iAbort < nAbort) continue; |
---|
82 | cout << " Event generation aborted prematurely, owing to error!\n"; |
---|
83 | break; |
---|
84 | } |
---|
85 | |
---|
86 | // Loop over all particles and analyze the final-state ones. |
---|
87 | for (int i = 0; i < pythia.event.size(); ++i) |
---|
88 | if (pythia.event[i].isFinal()) { |
---|
89 | int idAbs = pythia.event[i].idAbs(); |
---|
90 | double eI = pythia.event[i].e(); |
---|
91 | if (idAbs == 22) eGamma.fill(eI); |
---|
92 | else if (idAbs == 11) eE.fill(eI); |
---|
93 | else if (idAbs == 2212) eP.fill(eI); |
---|
94 | else if (idAbs == 12 || idAbs == 14 || idAbs == 16) eNu.fill(eI); |
---|
95 | else { |
---|
96 | eRest.fill(eI); |
---|
97 | cout << " Error: stable id = " << pythia.event[i].id() << endl; |
---|
98 | } |
---|
99 | } |
---|
100 | |
---|
101 | // End of event loop. |
---|
102 | } |
---|
103 | |
---|
104 | // Final statistics and histograms. |
---|
105 | pythia.stat(); |
---|
106 | cout << eGamma << eE << eP << eNu << eRest; |
---|
107 | |
---|
108 | // Done. |
---|
109 | return 0; |
---|
110 | } |
---|