]>
Commit | Line | Data |
---|---|---|
9419eeef | 1 | // main29.cc is a part of the PYTHIA event generator. |
2 | // Copyright (C) 2010 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 | ||
14 | #include "Pythia.h" | |
15 | ||
16 | using namespace Pythia8; | |
17 | ||
18 | //========================================================================== | |
19 | ||
20 | // A derived class for (e+ e- ->) GenericResonance -> various final states. | |
21 | ||
22 | class Sigma1GenRes : public Sigma1Process { | |
23 | ||
24 | public: | |
25 | ||
26 | // Constructor. | |
27 | Sigma1GenRes() {} | |
28 | ||
29 | // Evaluate sigmaHat(sHat): dummy unit cross section. | |
30 | virtual double sigmaHat() {return 1.;} | |
31 | ||
32 | // Select flavour. No colour or anticolour. | |
33 | virtual void setIdColAcol() {setId( -11, 11, 999999); | |
34 | setColAcol( 0, 0, 0, 0, 0, 0);} | |
35 | ||
36 | // Info on the subprocess. | |
37 | virtual string name() const {return "GenericResonance";} | |
38 | virtual int code() const {return 9001;} | |
39 | virtual string inFlux() const {return "ffbarSame";} | |
40 | ||
41 | }; | |
42 | ||
43 | //========================================================================== | |
44 | ||
45 | int main() { | |
46 | ||
47 | // Pythia generator. | |
48 | Pythia pythia; | |
49 | ||
50 | // A class to generate the fictitious resonance initial state. | |
51 | SigmaProcess* sigma1GenRes = new Sigma1GenRes(); | |
52 | ||
53 | // Hand pointer to Pythia. | |
54 | pythia.setSigmaPtr( sigma1GenRes); | |
55 | ||
56 | // Read in the rest of the settings and data from a separate file. | |
57 | pythia.readFile("main29.cmnd"); | |
58 | ||
59 | // Initialization. | |
60 | pythia.init(); | |
61 | ||
62 | // Extract settings to be used in the main program. | |
63 | int nEvent = pythia.mode("Main:numberOfEvents"); | |
64 | int nList = pythia.mode("Main:numberToList"); | |
65 | int nShow = pythia.mode("Main:timesToShow"); | |
66 | int nAbort = pythia.mode("Main:timesAllowErrors"); | |
67 | bool showCS = pythia.flag("Main:showChangedSettings"); | |
68 | bool showCPD = pythia.flag("Main:showChangedParticleData"); | |
69 | ||
70 | // List changes. | |
71 | if (showCS) pythia.settings.listChanged(); | |
72 | if (showCPD) pythia.particleData.listChanged(); | |
73 | ||
74 | // Histogram particle spectra. | |
75 | Hist eGamma("energy spectrum of photons", 100, 0., 250.); | |
76 | Hist eE( "energy spectrum of e+ and e-", 100, 0., 250.); | |
77 | Hist eP( "energy spectrum of p and pbar", 100, 0., 250.); | |
78 | Hist eNu( "energy spectrum of neutrinos", 100, 0., 250.); | |
79 | Hist eRest( "energy spectrum of rest particles", 100, 0., 250.); | |
80 | ||
81 | // Begin event loop. | |
82 | int nPace = max(1, nEvent / max(1, nShow) ); | |
83 | int iAbort = 0; | |
84 | for (int iEvent = 0; iEvent < nEvent; ++iEvent) { | |
85 | if (nShow > 0 && iEvent%nPace == 0) | |
86 | cout << " Now begin event " << iEvent << endl; | |
87 | ||
88 | // Generate events. Quit if many failures. | |
89 | if (!pythia.next()) { | |
90 | if (++iAbort < nAbort) continue; | |
91 | cout << " Event generation aborted prematurely, owing to error!\n"; | |
92 | break; | |
93 | } | |
94 | ||
95 | // List first few events. | |
96 | if (iEvent < nList) { | |
97 | pythia.info.list(); | |
98 | pythia.event.list(); | |
99 | } | |
100 | ||
101 | // Loop over all particles and analyze the final-state ones. | |
102 | for (int i = 0; i < pythia.event.size(); ++i) | |
103 | if (pythia.event[i].isFinal()) { | |
104 | int idAbs = pythia.event[i].idAbs(); | |
105 | double eI = pythia.event[i].e(); | |
106 | if (idAbs == 22) eGamma.fill(eI); | |
107 | else if (idAbs == 11) eE.fill(eI); | |
108 | else if (idAbs == 2212) eP.fill(eI); | |
109 | else if (idAbs == 12 || idAbs == 14 || idAbs == 16) eNu.fill(eI); | |
110 | else { | |
111 | eRest.fill(eI); | |
112 | cout << " Error: stable id = " << pythia.event[i].id() << endl; | |
113 | } | |
114 | } | |
115 | ||
116 | // End of event loop. | |
117 | } | |
118 | ||
119 | // Final statistics and histograms. | |
120 | pythia.statistics(); | |
121 | cout << eGamma << eE << eP << eNu << eRest; | |
122 | ||
123 | // Done. | |
124 | return 0; | |
125 | } |