1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
|
//////////////////////////////////////////////////////////////////////////
// example_PythiaStreamIO.cc
//
// garren@fnal.gov, May 2009
//
//////////////////////////////////////////////////////////////////////////
/// example of generating events with Pythia using HepMC/PythiaWrapper.h
/// Events are read into the HepMC event record from the FORTRAN HEPEVT
/// common block using the IO_HEPEVT strategy
///
/// To Compile: go to the HepMC example directory and type:
/// make example_PythiaStreamIO.exe
///
/// This example uses streaming I/O
/// writePythiaStreamIO() sets the cross section in GenRun
/// readPythiaStreamIO() reads the file written by writePythiaStreamIO()
///
//////////////////////////////////////////////////////////////////////////
#include <fstream>
#include <iostream>
#include "HepMC/PythiaWrapper.h"
#include "HepMC/IO_HEPEVT.h"
#include "HepMC/GenEvent.h"
#include "PythiaHelper.h"
void writePythiaStreamIO();
void readPythiaStreamIO();
int main() {
writePythiaStreamIO();
readPythiaStreamIO();
return 0;
}
void writePythiaStreamIO() {
// example to generate events and write output
std::cout << std::endl;
std::cout << "Begin pythia_out()" << std::endl;
//........................................HEPEVT
// Pythia 6.1 uses HEPEVT with 4000 entries and 8-byte floating point
// numbers. We need to explicitly pass this information to the
// HEPEVT_Wrapper.
//
HepMC::HEPEVT_Wrapper::set_max_number_entries(4000);
HepMC::HEPEVT_Wrapper::set_sizeof_real(8);
//
//........................................PYTHIA INITIALIZATIONS
initPythia();
//........................................HepMC INITIALIZATIONS
//
// Instantiate an IO strategy for reading from HEPEVT.
HepMC::IO_HEPEVT hepevtio;
//
{ // begin scope of ascii_io
// declare an output stream
const char outfile[] = "example_PythiaStreamIO_write.dat";
std::ofstream ascii_io( outfile );
if( !ascii_io ) {
std::cerr << "cannot open " << outfile << std::endl;
exit(-1);
}
// use the default IO_GenEvent precision
ascii_io.precision(16);
// write the line that defines the beginning of a GenEvent block
HepMC::write_HepMC_IO_block_begin( ascii_io );
//
//........................................EVENT LOOP
for ( int i = 1; i <= 100; i++ ) {
if ( i%50==1 ) std::cout << "Processing Event Number "
<< i << std::endl;
call_pyevnt(); // generate one event with Pythia
// pythia pyhepc routine converts common PYJETS in common HEPEVT
call_pyhepc( 1 );
HepMC::GenEvent* evt = hepevtio.read_next_event();
// define the units (Pythia uses GeV and mm)
evt->use_units(HepMC::Units::GEV, HepMC::Units::MM);
// add some information to the event
evt->set_event_number(i);
evt->set_signal_process_id(20);
// set number of multi parton interactions
evt->set_mpi( pypars.msti[31-1] );
// set cross section information
evt->set_cross_section( HepMC::getPythiaCrossSection() );
// write the event out to the ascii files
ascii_io << (*evt);;
// we also need to delete the created event from memory
delete evt;
}
// write the line that defines the end of a GenEvent block
HepMC::write_HepMC_IO_block_end( ascii_io );
//........................................TERMINATION
// write out some information from Pythia to the screen
call_pystat( 1 );
} // end scope of ascii_io
}
void readPythiaStreamIO() {
// example to read events written by writePythiaStreamIO
// and write them back out
std::cout << std::endl;
// input units are GeV and mm
const char infile[] = "example_PythiaStreamIO_write.dat";
std::ifstream is( infile );
if( !is ) {
std::cerr << "cannot open " << infile << std::endl;
exit(-1);
}
//
{ // begin scope of ascii_io
// declare an output stream
const char outfile[] = "example_PythiaStreamIO_read.dat";
std::ofstream ascii_io( outfile );
if( !ascii_io ) {
std::cerr << "cannot open " << outfile << std::endl;
exit(-1);
}
ascii_io.precision(16);
HepMC::write_HepMC_IO_block_begin( ascii_io );
//
//........................................EVENT LOOP
HepMC::GenEvent evt;
int i = 0;
while ( is ) {
evt.read( is );
// make sure we have a valid event
if( evt.is_valid() ) {
++i;
if ( i%50==1 ) std::cout << "Processing Event Number "
<< i << std::endl;
if ( i%25==2 ) {
// write the cross section if it exists
if( evt.cross_section() ) {
std::cout << "cross section at event " << i << " is "
<< evt.cross_section()->cross_section()
<< std::endl;
}
}
// write the event out to the ascii files
evt.write( ascii_io );
}
}
//........................................TERMINATION
HepMC::write_HepMC_IO_block_end( ascii_io );
} // end scope of ascii_io
}
|
