1 | ////////////////////////////////////////////////////////////////////////// |
---|
2 | // Matt.Dobbs@Cern.CH, December 1999 |
---|
3 | // November 2000, updated to use Pythia 6.1 |
---|
4 | // |
---|
5 | ////////////////////////////////////////////////////////////////////////// |
---|
6 | /// example of generating events with Pythia using HepMC/PythiaWrapper.h |
---|
7 | /// Events are read into the HepMC event record from the FORTRAN HEPEVT |
---|
8 | /// common block using the IO_HEPEVT strategy |
---|
9 | /// |
---|
10 | /// To Compile: go to the HepMC directory and type: |
---|
11 | /// gmake examples/example_MyPythia.exe |
---|
12 | /// |
---|
13 | /// In this example the precision and number of entries for the HEPEVT |
---|
14 | /// fortran common block are explicitly defined to correspond to those |
---|
15 | /// used in the Pythia version of the HEPEVT common block. |
---|
16 | /// |
---|
17 | /// If you get funny output from HEPEVT in your own code, probably you have |
---|
18 | /// set these values incorrectly! |
---|
19 | /// |
---|
20 | ////////////////////////////////////////////////////////////////////////// |
---|
21 | /// |
---|
22 | /// pythia_out(): |
---|
23 | /// Events are read into the HepMC event record from the FORTRAN HEPEVT |
---|
24 | /// common block using the IO_HEPEVT strategy and then output to file in |
---|
25 | /// ascii format using the IO_GenEvent strategy. |
---|
26 | /// |
---|
27 | /// pythia_particle_out(): |
---|
28 | /// Events are read into the HepMC event record from the FORTRAN HEPEVT |
---|
29 | /// common block using the IO_HEPEVT strategy and then output to file in |
---|
30 | /// ascii format using the IO_AsciiParticles strategy. |
---|
31 | /// This is identical to pythia_out() except for the choice of output format. |
---|
32 | /// |
---|
33 | /// event_selection(): |
---|
34 | /// Events are read into the HepMC event record from the FORTRAN HEPEVT |
---|
35 | /// common block using the IO_HEPEVT strategy and then a very simple event |
---|
36 | /// selection is performed. |
---|
37 | /// |
---|
38 | /// pythia_in(): |
---|
39 | /// Read the file created by pythia_out(). |
---|
40 | /// |
---|
41 | /// pythia_in_out(): |
---|
42 | /// generate events with Pythia, write a file, and read the resulting output |
---|
43 | /// Notice that we use scope to explicitly close the ouput files. |
---|
44 | /// The two output files should be identical. |
---|
45 | /// |
---|
46 | |
---|
47 | |
---|
48 | #include <iostream> |
---|
49 | #include "HepMC/PythiaWrapper.h" |
---|
50 | #include "HepMC/IO_HEPEVT.h" |
---|
51 | #include "HepMC/IO_GenEvent.h" |
---|
52 | #include "HepMC/IO_AsciiParticles.h" |
---|
53 | #include "HepMC/GenEvent.h" |
---|
54 | #include "PythiaHelper.h" |
---|
55 | |
---|
56 | //! example class |
---|
57 | |
---|
58 | /// \class IsGoodEventMyPythia |
---|
59 | /// event selection predicate. returns true if the event contains |
---|
60 | /// a photon with pT > 25 GeV |
---|
61 | class IsGoodEventMyPythia { |
---|
62 | public: |
---|
63 | /// returns true if event is "good" |
---|
64 | bool operator()( const HepMC::GenEvent* evt ) { |
---|
65 | for ( HepMC::GenEvent::particle_const_iterator p |
---|
66 | = evt->particles_begin(); p != evt->particles_end(); ++p ){ |
---|
67 | if ( (*p)->pdg_id() == 22 && (*p)->momentum().perp() > 25. ) { |
---|
68 | //std::cout << "Event " << evt->event_number() |
---|
69 | // << " is a good event." << std::endl; |
---|
70 | //(*p)->print(); |
---|
71 | return 1; |
---|
72 | } |
---|
73 | } |
---|
74 | return 0; |
---|
75 | } |
---|
76 | }; |
---|
77 | |
---|
78 | |
---|
79 | void pythia_out(); |
---|
80 | void pythia_in(); |
---|
81 | void pythia_in_out(); |
---|
82 | void event_selection(); |
---|
83 | void pythia_particle_out(); |
---|
84 | |
---|
85 | int main() { |
---|
86 | // example to generate events and write output |
---|
87 | pythia_out(); |
---|
88 | // example to generate events and perform simple event selection |
---|
89 | event_selection(); |
---|
90 | // example to read the file written by pythia_out |
---|
91 | pythia_in(); |
---|
92 | // example to generate events, write them, and read them back |
---|
93 | pythia_in_out(); |
---|
94 | |
---|
95 | return 0; |
---|
96 | } |
---|
97 | |
---|
98 | |
---|
99 | void pythia_out() |
---|
100 | { |
---|
101 | std::cout << std::endl; |
---|
102 | std::cout << "Begin pythia_out()" << std::endl; |
---|
103 | //........................................HEPEVT |
---|
104 | // Pythia 6.1 uses HEPEVT with 4000 entries and 8-byte floating point |
---|
105 | // numbers. We need to explicitly pass this information to the |
---|
106 | // HEPEVT_Wrapper. |
---|
107 | // |
---|
108 | HepMC::HEPEVT_Wrapper::set_max_number_entries(4000); |
---|
109 | HepMC::HEPEVT_Wrapper::set_sizeof_real(8); |
---|
110 | // |
---|
111 | //........................................PYTHIA INITIALIZATIONS |
---|
112 | initPythia(); |
---|
113 | |
---|
114 | //........................................HepMC INITIALIZATIONS |
---|
115 | // |
---|
116 | // Instantiate an IO strategy for reading from HEPEVT. |
---|
117 | HepMC::IO_HEPEVT hepevtio; |
---|
118 | // |
---|
119 | { // begin scope of ascii_io |
---|
120 | // Instantiate an IO strategy to write the data to file |
---|
121 | HepMC::IO_GenEvent ascii_io("example_MyPythia.dat",std::ios::out); |
---|
122 | // |
---|
123 | //........................................EVENT LOOP |
---|
124 | for ( int i = 1; i <= 100; i++ ) { |
---|
125 | if ( i%50==1 ) std::cout << "Processing Event Number " |
---|
126 | << i << std::endl; |
---|
127 | call_pyevnt(); // generate one event with Pythia |
---|
128 | // pythia pyhepc routine converts common PYJETS in common HEPEVT |
---|
129 | call_pyhepc( 1 ); |
---|
130 | HepMC::GenEvent* evt = hepevtio.read_next_event(); |
---|
131 | // define the units (Pythia uses GeV and mm) |
---|
132 | evt->use_units(HepMC::Units::GEV, HepMC::Units::MM); |
---|
133 | // add some information to the event |
---|
134 | evt->set_event_number(i); |
---|
135 | evt->set_signal_process_id(20); |
---|
136 | // set number of multi parton interactions |
---|
137 | evt->set_mpi( pypars.msti[31-1] ); |
---|
138 | // set cross section information |
---|
139 | evt->set_cross_section( HepMC::getPythiaCrossSection() ); |
---|
140 | // write the event out to the ascii files |
---|
141 | ascii_io << evt; |
---|
142 | // we also need to delete the created event from memory |
---|
143 | delete evt; |
---|
144 | } |
---|
145 | //........................................TERMINATION |
---|
146 | // write out some information from Pythia to the screen |
---|
147 | call_pystat( 1 ); |
---|
148 | } // end scope of ascii_io |
---|
149 | } |
---|
150 | |
---|
151 | |
---|
152 | void event_selection() |
---|
153 | { |
---|
154 | std::cout << std::endl; |
---|
155 | std::cout << "Begin event_selection()" << std::endl; |
---|
156 | //........................................HEPEVT |
---|
157 | // Pythia 6.1 uses HEPEVT with 4000 entries and 8-byte floating point |
---|
158 | // numbers. We need to explicitly pass this information to the |
---|
159 | // HEPEVT_Wrapper. |
---|
160 | // |
---|
161 | HepMC::HEPEVT_Wrapper::set_max_number_entries(4000); |
---|
162 | HepMC::HEPEVT_Wrapper::set_sizeof_real(8); |
---|
163 | // |
---|
164 | //........................................PYTHIA INITIALIZATIONS |
---|
165 | initPythia(); |
---|
166 | // |
---|
167 | //........................................HepMC INITIALIZATIONS |
---|
168 | // Instantiate an IO strategy for reading from HEPEVT. |
---|
169 | HepMC::IO_HEPEVT hepevtio; |
---|
170 | // declare an instance of the event selection predicate |
---|
171 | IsGoodEventMyPythia is_good_event; |
---|
172 | //........................................EVENT LOOP |
---|
173 | int icount=0; |
---|
174 | int num_good_events=0; |
---|
175 | for ( int i = 1; i <= 100; i++ ) { |
---|
176 | icount++; |
---|
177 | if ( i%50==1 ) std::cout << "Processing Event Number " |
---|
178 | << i << std::endl; |
---|
179 | call_pyevnt(); // generate one event with Pythia |
---|
180 | // pythia pyhepc routine convert common PYJETS in common HEPEVT |
---|
181 | call_pyhepc( 1 ); |
---|
182 | HepMC::GenEvent* evt = hepevtio.read_next_event(); |
---|
183 | // define the units (Pythia uses GeV and mm) |
---|
184 | evt->use_units(HepMC::Units::GEV, HepMC::Units::MM); |
---|
185 | // set number of multi parton interactions |
---|
186 | evt->set_mpi( pypars.msti[31-1] ); |
---|
187 | // set cross section information |
---|
188 | evt->set_cross_section( HepMC::getPythiaCrossSection() ); |
---|
189 | // do event selection |
---|
190 | if ( is_good_event(evt) ) { |
---|
191 | std::cout << "Good Event Number " << i << std::endl; |
---|
192 | ++num_good_events; |
---|
193 | } |
---|
194 | // we also need to delete the created event from memory |
---|
195 | delete evt; |
---|
196 | } |
---|
197 | //........................................TERMINATION |
---|
198 | // write out some information from Pythia to the screen |
---|
199 | call_pystat( 1 ); |
---|
200 | //........................................PRINT RESULTS |
---|
201 | std::cout << num_good_events << " out of " << icount |
---|
202 | << " processed events passed the cuts. Finished." << std::endl; |
---|
203 | } |
---|
204 | |
---|
205 | void pythia_in() |
---|
206 | { |
---|
207 | std::cout << std::endl; |
---|
208 | std::cout << "Begin pythia_in()" << std::endl; |
---|
209 | std::cout << "reading example_MyPythia.dat" << std::endl; |
---|
210 | //........................................define an input scope |
---|
211 | { |
---|
212 | // open input stream |
---|
213 | std::ifstream istr( "example_MyPythia.dat" ); |
---|
214 | if( !istr ) { |
---|
215 | std::cerr << "example_ReadMyPythia: cannot open example_MyPythia.dat" << std::endl; |
---|
216 | exit(-1); |
---|
217 | } |
---|
218 | HepMC::IO_GenEvent ascii_in(istr); |
---|
219 | // open output stream (alternate method) |
---|
220 | HepMC::IO_GenEvent ascii_out("example_MyPythia2.dat",std::ios::out); |
---|
221 | // now read the file |
---|
222 | int icount=0; |
---|
223 | HepMC::GenEvent* evt = ascii_in.read_next_event(); |
---|
224 | while ( evt ) { |
---|
225 | icount++; |
---|
226 | if ( icount%50==1 ) std::cout << "Processing Event Number " << icount |
---|
227 | << " its # " << evt->event_number() |
---|
228 | << std::endl; |
---|
229 | // write the event out to the ascii file |
---|
230 | ascii_out << evt; |
---|
231 | delete evt; |
---|
232 | ascii_in >> evt; |
---|
233 | } |
---|
234 | //........................................PRINT RESULT |
---|
235 | std::cout << icount << " events found. Finished." << std::endl; |
---|
236 | } // ascii_out and istr destructors are called here |
---|
237 | } |
---|
238 | |
---|
239 | void pythia_in_out() |
---|
240 | { |
---|
241 | std::cout << std::endl; |
---|
242 | std::cout << "Begin pythia_in_out()" << std::endl; |
---|
243 | //........................................HEPEVT |
---|
244 | // Pythia 6.3 uses HEPEVT with 4000 entries and 8-byte floating point |
---|
245 | // numbers. We need to explicitly pass this information to the |
---|
246 | // HEPEVT_Wrapper. |
---|
247 | // |
---|
248 | HepMC::HEPEVT_Wrapper::set_max_number_entries(4000); |
---|
249 | HepMC::HEPEVT_Wrapper::set_sizeof_real(8); |
---|
250 | // |
---|
251 | //........................................PYTHIA INITIALIZATIONS |
---|
252 | initPythia(); |
---|
253 | |
---|
254 | //........................................HepMC INITIALIZATIONS |
---|
255 | // |
---|
256 | // Instantiate an IO strategy for reading from HEPEVT. |
---|
257 | HepMC::IO_HEPEVT hepevtio; |
---|
258 | // |
---|
259 | //........................................define the output scope |
---|
260 | { |
---|
261 | // Instantial an IO strategy to write the data to file |
---|
262 | HepMC::IO_GenEvent ascii_io("example_MyPythiaRead.dat",std::ios::out); |
---|
263 | // |
---|
264 | //........................................EVENT LOOP |
---|
265 | for ( int i = 1; i <= 100; i++ ) { |
---|
266 | if ( i%50==1 ) std::cout << "Processing Event Number " |
---|
267 | << i << std::endl; |
---|
268 | call_pyevnt(); // generate one event with Pythia |
---|
269 | // pythia pyhepc routine converts common PYJETS in common HEPEVT |
---|
270 | call_pyhepc( 1 ); |
---|
271 | HepMC::GenEvent* evt = hepevtio.read_next_event(); |
---|
272 | // define the units (Pythia uses GeV and mm) |
---|
273 | evt->use_units(HepMC::Units::GEV, HepMC::Units::MM); |
---|
274 | // set cross section information |
---|
275 | evt->set_cross_section( HepMC::getPythiaCrossSection() ); |
---|
276 | // add some information to the event |
---|
277 | evt->set_event_number(i); |
---|
278 | evt->set_signal_process_id(20); |
---|
279 | // write the event out to the ascii file |
---|
280 | ascii_io << evt; |
---|
281 | // we also need to delete the created event from memory |
---|
282 | delete evt; |
---|
283 | } |
---|
284 | //........................................TERMINATION |
---|
285 | // write out some information from Pythia to the screen |
---|
286 | call_pystat( 1 ); |
---|
287 | } // ascii_io destructor is called here |
---|
288 | // |
---|
289 | //........................................define an input scope |
---|
290 | { |
---|
291 | // now read the file we wrote |
---|
292 | HepMC::IO_GenEvent ascii_in("example_MyPythiaRead.dat",std::ios::in); |
---|
293 | HepMC::IO_GenEvent ascii_io2("example_MyPythiaRead2.dat",std::ios::out); |
---|
294 | int icount=0; |
---|
295 | HepMC::GenEvent* evt = ascii_in.read_next_event(); |
---|
296 | while ( evt ) { |
---|
297 | icount++; |
---|
298 | if ( icount%50==1 ) std::cout << "Processing Event Number " << icount |
---|
299 | << " its # " << evt->event_number() |
---|
300 | << std::endl; |
---|
301 | // write the event out to the ascii file |
---|
302 | ascii_io2 << evt; |
---|
303 | delete evt; |
---|
304 | ascii_in >> evt; |
---|
305 | } |
---|
306 | //........................................PRINT RESULT |
---|
307 | std::cout << icount << " events found. Finished." << std::endl; |
---|
308 | } // ascii_io2 and ascii_in destructors are called here |
---|
309 | } |
---|
310 | |
---|
311 | void pythia_particle_out() |
---|
312 | { |
---|
313 | std::cout << std::endl; |
---|
314 | std::cout << "Begin pythia_particle_out()" << std::endl; |
---|
315 | //........................................HEPEVT |
---|
316 | // Pythia 6.1 uses HEPEVT with 4000 entries and 8-byte floating point |
---|
317 | // numbers. We need to explicitly pass this information to the |
---|
318 | // HEPEVT_Wrapper. |
---|
319 | // |
---|
320 | HepMC::HEPEVT_Wrapper::set_max_number_entries(4000); |
---|
321 | HepMC::HEPEVT_Wrapper::set_sizeof_real(8); |
---|
322 | // |
---|
323 | //........................................PYTHIA INITIALIZATIONS |
---|
324 | initPythia(); |
---|
325 | |
---|
326 | //........................................HepMC INITIALIZATIONS |
---|
327 | // |
---|
328 | // Instantiate an IO strategy for reading from HEPEVT. |
---|
329 | HepMC::IO_HEPEVT hepevtio; |
---|
330 | // |
---|
331 | { // begin scope of ascii_io |
---|
332 | // Instantiate an IO strategy to write the data to file |
---|
333 | HepMC::IO_AsciiParticles ascii_io("example_PythiaParticle.dat",std::ios::out); |
---|
334 | // |
---|
335 | //........................................EVENT LOOP |
---|
336 | for ( int i = 1; i <= 100; i++ ) { |
---|
337 | if ( i%50==1 ) std::cout << "Processing Event Number " |
---|
338 | << i << std::endl; |
---|
339 | call_pyevnt(); // generate one event with Pythia |
---|
340 | // pythia pyhepc routine converts common PYJETS in common HEPEVT |
---|
341 | call_pyhepc( 1 ); |
---|
342 | HepMC::GenEvent* evt = hepevtio.read_next_event(); |
---|
343 | // define the units (Pythia uses GeV and mm) |
---|
344 | evt->use_units(HepMC::Units::GEV, HepMC::Units::MM); |
---|
345 | // set cross section information |
---|
346 | evt->set_cross_section( HepMC::getPythiaCrossSection() ); |
---|
347 | // add some information to the event |
---|
348 | evt->set_event_number(i); |
---|
349 | evt->set_signal_process_id(20); |
---|
350 | // write the event out to the ascii file |
---|
351 | ascii_io << evt; |
---|
352 | // we also need to delete the created event from memory |
---|
353 | delete evt; |
---|
354 | } |
---|
355 | //........................................TERMINATION |
---|
356 | // write out some information from Pythia to the screen |
---|
357 | call_pystat( 1 ); |
---|
358 | } // end scope of ascii_io |
---|
359 | } |
---|
360 | |
---|