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
|
//////////////////////////////////////////////////////////////////////////
// Matt.Dobbs@Cern.CH, October 2002
// example of generating events with Herwig using HepMC/HerwigWrapper.h
// Events are read into the HepMC event record from the FORTRAN HEPEVT
// common block using the IO_HERWIG strategy.
//////////////////////////////////////////////////////////////////////////
/// To Compile: go to the HepMC directory and type:
/// gmake examples/example_MyHerwig.exe
///
/// In this example the precision and number of entries for the HEPEVT
/// fortran common block are explicitly defined to correspond to those
/// used in the Herwig version of the HEPEVT common block.
/// If you get funny output from HEPEVT in your own code, probably you have
/// set these values incorrectly!
///
#include <iostream>
#include "HepMC/HerwigWrapper.h"
#include "HepMC/IO_HERWIG.h"
#include "HepMC/IO_GenEvent.h"
#include "HepMC/GenEvent.h"
#include "HepMC/HEPEVT_Wrapper.h"
extern "C" void hwaend_() {}
int main() {
//
//........................................HEPEVT
// Herwig 6.4 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);
//
//.......................................INITIALIZATIONS
hwproc.PBEAM1 = 7000.; // energy of beam1
hwproc.PBEAM2 = 7000.; // energy of beam2
// 1610 = gg->H--> WW, 1706 = qq-->ttbar, 2510 = ttH -> ttWW
hwproc.IPROC = 1706; // qq -> ttbar production
hwproc.MAXEV = 100; // number of events
// tell it what the beam particles are:
for ( unsigned int i = 0; i < 8; ++i ) {
hwbmch.PART1[i] = (i < 1) ? 'P' : ' ';
hwbmch.PART2[i] = (i < 1) ? 'P' : ' ';
}
hwigin(); // INITIALISE OTHER COMMON BLOCKS
hwevnt.MAXPR = 1; // number of events to print
hwuinc(); // compute parameter-dependent constants
hweini(); // initialise elementary process
//........................................HepMC INITIALIZATIONS
//
// Instantiate an IO strategy for reading from HEPEVT.
HepMC::IO_HERWIG hepevtio;
// Instantiate an IO strategy to write the data to file
HepMC::IO_GenEvent ascii_io("example_MyHerwig.dat",std::ios::out);
//
//........................................EVENT LOOP
for ( int i = 1; i <= hwproc.MAXEV; i++ ) {
if ( i%50==1 ) std::cout << "Processing Event Number "
<< i << std::endl;
// initialise event
hwuine();
// generate hard subprocess
hwepro();
// generate parton cascades
hwbgen();
// do heavy object decays
hwdhob();
// do cluster formation
hwcfor();
// do cluster decays
hwcdec();
// do unstable particle decays
hwdhad();
// do heavy flavour hadron decays
hwdhvy();
// add soft underlying event if needed
hwmevt();
// finish event
hwufne();
HepMC::GenEvent* evt = hepevtio.read_next_event();
// define the units (Herwig uses GeV and mm)
evt->use_units(HepMC::Units::GEV, HepMC::Units::MM);
// set cross section information
evt->set_cross_section( HepMC::getHerwigCrossSection(i) );
// add some information to the event
evt->set_event_number(i);
evt->set_signal_process_id(20);
if (i<=hwevnt.MAXPR) {
std::cout << "\n\n This is the FIXED version of HEPEVT as "
<< "coded in IO_HERWIG " << std::endl;
HepMC::HEPEVT_Wrapper::print_hepevt();
evt->print();
}
// write the event to the ascii file
ascii_io << evt;
// we also need to delete the created event from memory
delete evt;
}
//........................................TERMINATION
hwefin();
return 0;
}
|
