1 | // main91.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 | // This is a simple test program illustrating how to link in Pythia 6.4 |
---|
7 | // for the generation of hard processes. That part is now considered |
---|
8 | // obsolete, but for some debug work this example has been collected |
---|
9 | // from various code pieces that were separately available up until |
---|
10 | // Pythia 8.125. In addition to modifying the below code to fit your |
---|
11 | // needs, you also have to modify examples/Makefile to link properly |
---|
12 | // for main91, including the Pythia 6.4xx library to be used (where |
---|
13 | // xx is the current subversion number). |
---|
14 | |
---|
15 | // All hard PYTHIA 6.4 processes should be available for full generation |
---|
16 | // in PYTHIA 8, at least to the extent that they are defined for p p, |
---|
17 | // pbar p or e+ e-collisions. Soft processes, i.e. elastic and diffractive |
---|
18 | // scattering, as well as minimum-bias events, require a different |
---|
19 | // kinematics machinery, and can only be generated with the internal |
---|
20 | // PYTHIA 8 processes. |
---|
21 | |
---|
22 | // PYTHIA 6.4 does its own rejection of events internally, according to |
---|
23 | // the strategy option 3. However, the current runtime interface does not |
---|
24 | // take cross-section information from PYTHIA 6.4. This means that both |
---|
25 | // the initial maxima and the final cross sections printed by the PYTHIA 8 |
---|
26 | // routines are irrelevant in this case. Instead you have to study the |
---|
27 | // standard PYTHIA 6.4 initialization printout and call on pystat(...) |
---|
28 | // at the end of the run to obtain this information. It also means that |
---|
29 | // you cannot mix with internally generated PYTHIA 8 events. |
---|
30 | |
---|
31 | //========================================================================== |
---|
32 | |
---|
33 | #include "Pythia.h" |
---|
34 | #include "LHAFortran.h" |
---|
35 | |
---|
36 | using namespace Pythia8; |
---|
37 | |
---|
38 | //========================================================================== |
---|
39 | |
---|
40 | // Declare the Fortran subroutines that may be used. |
---|
41 | // This code section is generic. |
---|
42 | |
---|
43 | #ifdef _WIN32 |
---|
44 | #define pygive_ PYGIVE |
---|
45 | #define pyinit_ PYINIT |
---|
46 | #define pyupin_ PYUPIN |
---|
47 | #define pyupev_ PYUPEV |
---|
48 | #define pylist_ PYLIST |
---|
49 | #define pystat_ PYSTAT |
---|
50 | #endif |
---|
51 | |
---|
52 | extern "C" { |
---|
53 | #ifdef _WIN32 |
---|
54 | extern void pyinit_(const char*, int, const char*, int, const char*, |
---|
55 | int, double&); |
---|
56 | #else |
---|
57 | extern void pyinit_(const char*, const char*, const char*, double&, |
---|
58 | int, int, int); |
---|
59 | #endif |
---|
60 | } |
---|
61 | |
---|
62 | extern "C" { |
---|
63 | extern void pygive_(const char*, int); |
---|
64 | extern void pyupin_(); |
---|
65 | extern void pyupev_(); |
---|
66 | extern void pylist_(int&); |
---|
67 | extern void pystat_(int&); |
---|
68 | } |
---|
69 | |
---|
70 | //========================================================================== |
---|
71 | |
---|
72 | // Interfaces to the above routines, to make the C++ calls similar to Fortran. |
---|
73 | // This code section is generic. |
---|
74 | |
---|
75 | class Pythia6Interface { |
---|
76 | |
---|
77 | public: |
---|
78 | |
---|
79 | // Give in a command to change a setting. |
---|
80 | static void pygive(const string cmnd) { |
---|
81 | const char* cstring = cmnd.c_str(); int len = cmnd.length(); |
---|
82 | pygive_(cstring, len); |
---|
83 | } |
---|
84 | |
---|
85 | // Initialize the generation for the given beam confiuration. |
---|
86 | static void pyinit(const string frame, const string beam, |
---|
87 | const string target, double wIn) { |
---|
88 | const char* cframe = frame.c_str(); int lenframe = frame.length(); |
---|
89 | const char* cbeam = beam.c_str(); int lenbeam = beam.length(); |
---|
90 | const char* ctarget = target.c_str(); int lentarget = target.length(); |
---|
91 | #ifdef _WIN32 |
---|
92 | pyinit_(cframe, lenframe, cbeam, lenbeam, ctarget, lentarget, wIn); |
---|
93 | #else |
---|
94 | pyinit_(cframe, cbeam, ctarget, wIn, lenframe, lenbeam, lentarget); |
---|
95 | #endif |
---|
96 | } |
---|
97 | |
---|
98 | // Fill the initialization information in the HEPRUP commonblock. |
---|
99 | static void pyupin() {pyupin_();} |
---|
100 | |
---|
101 | // Generate the next hard process and |
---|
102 | // fill the event information in the HEPEUP commonblock |
---|
103 | static void pyupev() {pyupev_();} |
---|
104 | |
---|
105 | // List the event at the process level. |
---|
106 | static void pylist(int mode) {pylist_(mode);} |
---|
107 | |
---|
108 | // Print statistics on the event generation process. |
---|
109 | static void pystat(int mode) {pystat_(mode);} |
---|
110 | |
---|
111 | }; |
---|
112 | |
---|
113 | //========================================================================== |
---|
114 | |
---|
115 | // Implement initialization fillHepRup method for Pythia6 example. |
---|
116 | // This code section is specific to the kind of precesses you generate. |
---|
117 | |
---|
118 | // Of all parameters that could be set with pygive, only those that |
---|
119 | // influence the generation of the hard processes have any impact, |
---|
120 | // since this is the only part of the Fortran code that is used. |
---|
121 | |
---|
122 | bool LHAupFortran::fillHepRup() { |
---|
123 | |
---|
124 | // Set process to generate. |
---|
125 | // Example 1: QCD production; must set pTmin. |
---|
126 | //Pythia6Interface::pygive("msel = 1"); |
---|
127 | //Pythia6Interface::pygive("ckin(3) = 20."); |
---|
128 | // Example 2: t-tbar production. |
---|
129 | //Pythia6Interface::pygive("msel = 6"); |
---|
130 | // Example 3: Z0 production; must set mMin. |
---|
131 | Pythia6Interface::pygive("msel = 11"); |
---|
132 | Pythia6Interface::pygive("ckin(1) = 50."); |
---|
133 | |
---|
134 | // Speed up initialization: multiparton interactions only in C++ code. |
---|
135 | Pythia6Interface::pygive("mstp(81)=0"); |
---|
136 | |
---|
137 | // Initialize for 14 TeV pp collider. |
---|
138 | Pythia6Interface::pyinit("cms","p","p",14000.); |
---|
139 | |
---|
140 | // Fill initialization information in HEPRUP. |
---|
141 | Pythia6Interface::pyupin(); |
---|
142 | |
---|
143 | // Done. |
---|
144 | return true; |
---|
145 | |
---|
146 | } |
---|
147 | |
---|
148 | //========================================================================== |
---|
149 | |
---|
150 | // Implement event generation fillHepEup method for Pythia 6 example. |
---|
151 | // This code section is generic. |
---|
152 | |
---|
153 | bool LHAupFortran::fillHepEup() { |
---|
154 | |
---|
155 | // Generate and fill the next Pythia6 event in HEPEUP. |
---|
156 | Pythia6Interface::pyupev(); |
---|
157 | |
---|
158 | // Done. |
---|
159 | return true; |
---|
160 | |
---|
161 | } |
---|
162 | |
---|
163 | //========================================================================== |
---|
164 | |
---|
165 | // The main program. |
---|
166 | // This code section is specific to the physics study you want to do. |
---|
167 | |
---|
168 | int main() { |
---|
169 | |
---|
170 | // Generator. Shorthand for the event and for settings. |
---|
171 | Pythia pythia; |
---|
172 | Event& event = pythia.event; |
---|
173 | Settings& settings = pythia.settings; |
---|
174 | |
---|
175 | // Set Pythia8 generation aspects. Examples only. |
---|
176 | pythia.readString("BeamRemnants:primordialKThard = 2."); |
---|
177 | pythia.readString("MultipartonInteractions:bProfile = 3"); |
---|
178 | pythia.readString("Next:numberShowInfo = 0"); |
---|
179 | pythia.readString("Next:numberShowProcess = 0"); |
---|
180 | pythia.readString("Next:numberShowEvent = 0"); |
---|
181 | |
---|
182 | // Initialize to access Pythia6 generator by Les Houches interface. |
---|
183 | pythia.readString("Beams:frameType = 5"); |
---|
184 | LHAupFortran pythia6; |
---|
185 | pythia.setLHAupPtr( &pythia6); |
---|
186 | pythia.init(); |
---|
187 | |
---|
188 | // Set some generation values. |
---|
189 | int nEvent = 100; |
---|
190 | int nList = 1; |
---|
191 | int nAbort = 10; |
---|
192 | |
---|
193 | // List changed settings data. |
---|
194 | settings.listChanged(); |
---|
195 | |
---|
196 | // Histograms. |
---|
197 | double eCM = 14000.; |
---|
198 | double epTol = 1e-7 * eCM; |
---|
199 | Hist epCons("deviation from energy-momentum conservation",100,0.,epTol); |
---|
200 | Hist nFinal("final particle multiplicity",100,-0.5,1599.5); |
---|
201 | Hist nChg("final charged multiplicity",100,-0.5,799.5); |
---|
202 | |
---|
203 | // Begin event loop. |
---|
204 | int iAbort = 0; |
---|
205 | for (int iEvent = 0; iEvent < nEvent; ++iEvent) { |
---|
206 | |
---|
207 | // Generate events. Quit if too many failures. |
---|
208 | if (!pythia.next()) { |
---|
209 | if (++iAbort < nAbort) continue; |
---|
210 | cout << " Event generation aborted prematurely, owing to error!\n"; |
---|
211 | break; |
---|
212 | } |
---|
213 | |
---|
214 | // List first few events, both hard process and complete events. |
---|
215 | if (iEvent < nList) { |
---|
216 | pythia.info.list(); |
---|
217 | // This call to Pythia6 is superfluous, but shows it can be done. |
---|
218 | Pythia6Interface::pylist(1); |
---|
219 | pythia.process.list(); |
---|
220 | event.list(); |
---|
221 | } |
---|
222 | |
---|
223 | // Reset quantities to be summed over event. |
---|
224 | int nfin = 0; |
---|
225 | int nch = 0; |
---|
226 | Vec4 pSum = - (event[1].p() + event[2].p()); |
---|
227 | |
---|
228 | // Loop over final particles in the event. |
---|
229 | for (int i = 0; i < event.size(); ++i) |
---|
230 | if (event[i].isFinal()) { |
---|
231 | ++nfin; |
---|
232 | if (event[i].isCharged()) ++nch; |
---|
233 | pSum += event[i].p(); |
---|
234 | } |
---|
235 | |
---|
236 | // Fill summed quantities. |
---|
237 | double epDev = abs(pSum.e()) + abs(pSum.px()) + abs(pSum.py()) |
---|
238 | + abs(pSum.pz()); |
---|
239 | epCons.fill(epDev); |
---|
240 | nFinal.fill(nfin); |
---|
241 | nChg.fill(nch); |
---|
242 | |
---|
243 | // End of event loop. |
---|
244 | } |
---|
245 | |
---|
246 | // Final statistics. Must do call to Pythia6 explicitly. |
---|
247 | pythia.stat(); |
---|
248 | Pythia6Interface::pystat(1); |
---|
249 | |
---|
250 | // Histogram output. |
---|
251 | cout << epCons << nFinal<< nChg; |
---|
252 | |
---|
253 | // Done. |
---|
254 | return 0; |
---|
255 | } |
---|