|
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
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
|
#include <iostream>
#include <cstring>
#include <vector>
#include <memory>
#include <string>
#include <HepMC/GenEvent.h>
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h"
#include "GeneratorInterface/LHEInterface/interface/LHEEvent.h"
#include "GeneratorInterface/PartonShowerVeto/interface/JetMatchingAlpgen.h"
namespace gen {
extern "C" {
// Put here all the functions for interfacing Fortran code.
// extern void alinit_();
extern void alsetp_();
// extern void alevnt_();
extern void alveto_(int* ipveto);
extern void dbpart_();
extern void pyupre_();
void alshcd_(char csho[3]);
void alshen_();
} // extern "C"
// Constructor
JetMatchingAlpgen::JetMatchingAlpgen(const edm::ParameterSet& params)
: JetMatching(params), applyMatching(params.getParameter<bool>("applyMatching")), runInitialized(false) {
ahopts_.etclus = params.getParameter<double>("etMin");
ahopts_.rclus = params.getParameter<double>("drMin");
ahopts_.iexc = params.getParameter<bool>("exclusive");
}
// Destructor
JetMatchingAlpgen::~JetMatchingAlpgen() {}
std::set<std::string> JetMatchingAlpgen::capabilities() const {
std::set<std::string> result;
result.insert("psFinalState");
result.insert("hepevt");
result.insert("pythia6");
return result;
}
// Implements the Alpgen MLM method - use alpg_match.F
void JetMatchingAlpgen::init(const lhef::LHERunInfo* runInfo) {
// Read Alpgen run card stored in the LHERunInfo object.
std::vector<std::string> headerLines = runInfo->findHeader("AlpgenUnwParFile");
if (headerLines.empty())
throw cms::Exception("Generator|PartonShowerVeto") << "In order to use Alpgen jet matching, "
"the input file has to contain the corresponding "
"Alpgen headers."
<< std::endl;
// Parse the header using its bultin function.
header.parse(headerLines.begin(), headerLines.end());
// I don't want to print this right now.
// std::cout << "Alpgen header" << std::endl;
// std::cout << "========================" << std::endl;
// std::cout << "\tihrd = " << header.ihrd << std::endl;
// std::cout << "\tmc = " << header.masses[AlpgenHeader::mc]
// << ", mb = " << header.masses[AlpgenHeader::mb]
// << ", mt = " << header.masses[AlpgenHeader::mt]
// << ", mw = " << header.masses[AlpgenHeader::mw]
// << ", mz = " << header.masses[AlpgenHeader::mz]
// << ", mh = " << header.masses[AlpgenHeader::mh]
// << std::endl;
// for(std::map<AlpgenHeader::Parameter, double>::const_iterator iter =
// header.params.begin(); iter != header.params.end(); ++iter)
// std::cout << "\t" << AlpgenHeader::parameterName(iter->first)
// << " = " << iter->second << std::endl;
// std::cout << "\txsec = " << header.xsec
// << " +-" << header.xsecErr << std::endl;
// std::cout << "========================" << std::endl;
// Here we pass a few header variables to common block and
// call Alpgen init routine to do the rest.
// The variables passed first are the ones directly that
// need to be set up "manually": IHRD and the masses.
// (ebeam is set just to mimic the original code)
// Afterwards, we pass the full spectrum of Alpgen
// parameters directly into the AHPARS structure, to be
// treated by AHSPAR which is called inside alsetp_().
std::copy(header.masses, header.masses + AlpgenHeader::MASS_MAX, ahppara_.masses);
ahppara_.ihrd = header.ihrd;
ahppara_.ebeam = header.ebeam;
for (std::map<AlpgenHeader::Parameter, double>::const_iterator iter = header.params.begin();
iter != header.params.end();
++iter) {
if (iter->first <= 0 || iter->first >= (int)AHPARS::nparam - 1)
continue;
ahpars_.parval[(int)iter->first - 1] = iter->second;
}
// Run the rest of the setup.
alsetp_();
// When we reach this point, the run is fully initialized.
runInitialized = true;
}
void JetMatchingAlpgen::beforeHadronisation(const lhef::LHEEvent* event) {
// We can't continue if the run has not been initialized.
if (!runInitialized)
throw cms::Exception("Generator|PartonShowerVeto") << "Run not initialized in JetMatchingAlpgen" << std::endl;
// We are called just after LHEInterface has filled in
// the Fortran common block (and Pythia6 called UPEVNT).
// Possibly not interesting for us.
// (except perhaps for debugging?)
// pyupre_();
// dbpart_();
eventInitialized = true;
}
/*
int JetMatchingAlpgen::match(const HepMC::GenEvent *partonLevel,
const HepMC::GenEvent *finalState,
bool showeredFinalState)
*/
int JetMatchingAlpgen::match(const lhef::LHEEvent* partonLevel, const std::vector<fastjet::PseudoJet>* jetInput) {
/*
if (!showeredFinalState)
throw cms::Exception("Generator|PartonShowerVeto")
<< "Alpgen matching expected parton shower "
"final state." << std::endl;
*/
if (!runInitialized)
throw cms::Exception("Generator|PartonShowerVeto") << "Run not initialized in JetMatchingAlpgen" << std::endl;
if (!eventInitialized)
throw cms::Exception("Generator|PartonShowerVeto") << "Event not initialized in JetMatchingAlpgen" << std::endl;
// If matching not required (e.g., icckw = 0), don't run the
// FORTRAN veto code.
if (!applyMatching)
return 0;
// Call the Fortran veto code.
int veto = 0;
alveto_(&veto);
eventInitialized = false;
// If event was vetoed, the variable veto will contain the number 1.
// In this case, we must return 1 - that will be used as the return value from UPVETO.
// If event was accepted, the variable veto will contain the number 0.
// In this case, we must return 0 - that will be used as the return value from UPVETO.
return veto ? 1 : 0;
}
void alshcd_(char csho[3]) {
std::memcpy(csho, "PYT", 3); // or "HER"
}
void alshen_() {}
} // end namespace gen
|