|
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
|
#include <algorithm>
#include <iostream>
#include <iterator>
#include <sstream>
#include <string>
#include <memory>
#include <cassert>
#include "GeneratorInterface/Pythia8Interface/plugins/LHAupLesHouches.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
using namespace Pythia8;
bool LHAupLesHouches::setInit() {
if (!runInfo)
return false;
const lhef::HEPRUP &heprup = *runInfo->getHEPRUP();
setBeamA(heprup.IDBMUP.first, heprup.EBMUP.first, heprup.PDFGUP.first, heprup.PDFSUP.first);
setBeamB(heprup.IDBMUP.second, heprup.EBMUP.second, heprup.PDFGUP.second, heprup.PDFSUP.second);
setStrategy(heprup.IDWTUP);
for (int i = 0; i < heprup.NPRUP; i++)
addProcess(heprup.LPRUP[i], heprup.XSECUP[i], heprup.XERRUP[i], heprup.XMAXUP[i]);
//hadronisation->onInit().emit();
//runInfo.reset();
//fill SLHA header information if available
std::vector<std::string> slha = runInfo->findHeader("slha");
if (!slha.empty()) {
std::string slhaheader;
for (std::vector<std::string>::const_iterator iter = slha.begin(); iter != slha.end(); ++iter) {
slhaheader.append(*iter);
}
infoPtr->setHeader("slha", slhaheader);
}
//will be used to work around missing initialization inside pythia8
if (!fEvAttributes) {
fEvAttributes = new std::map<std::string, std::string>;
} else {
fEvAttributes->clear();
}
return true;
}
bool LHAupLesHouches::setEvent(int inProcId) {
if (!event)
return false;
if (event->getReadAttempts() > 0)
return false; // record already tried
const lhef::HEPEUP &hepeup = *event->getHEPEUP();
if (!hepeup.NUP)
return false;
setProcess(hepeup.IDPRUP, hepeup.XWGTUP, hepeup.SCALUP, hepeup.AQEDUP, hepeup.AQCDUP);
const std::vector<float> &scales = event->scales();
unsigned int iscale = 0;
for (int i = 0; i < hepeup.NUP; i++) {
//retrieve scale corresponding to each particle
double scalein = -1.;
//handle clustering scales if present,
//applies to outgoing partons only
if (setScalesFromLHEF_ && !scales.empty() && hepeup.ISTUP[i] == 1) {
if (iscale >= scales.size()) {
edm::LogError("InvalidLHEInput") << "Pythia8 requires"
<< "cluster scales for all outgoing partons or for none" << std::endl;
}
scalein = scales[iscale];
++iscale;
}
addParticle(hepeup.IDUP[i],
hepeup.ISTUP[i],
hepeup.MOTHUP[i].first,
hepeup.MOTHUP[i].second,
hepeup.ICOLUP[i].first,
hepeup.ICOLUP[i].second,
hepeup.PUP[i][0],
hepeup.PUP[i][1],
hepeup.PUP[i][2],
hepeup.PUP[i][3],
hepeup.PUP[i][4],
hepeup.VTIMUP[i],
hepeup.SPINUP[i],
scalein);
}
if (!infoPtr->eventAttributes) {
fEvAttributes->clear();
infoPtr->eventAttributes = fEvAttributes;
} else {
infoPtr->eventAttributes = fEvAttributes; // make sure still there
infoPtr->eventAttributes->clear();
}
//fill parton multiplicities if available
int npLO = event->npLO();
int npNLO = event->npNLO();
//default value of -99 indicates tags were not present in the original LHE file
//don't pass to pythia in this case to emulate pythia internal lhe reader behaviour
if (npLO != -99) {
char buffer[100];
snprintf(buffer, 100, "%i", npLO);
(*infoPtr->eventAttributes)["npLO"] = buffer;
}
if (npNLO != -99) {
char buffer[100];
snprintf(buffer, 100, "%i", npNLO);
(*infoPtr->eventAttributes)["npNLO"] = buffer;
}
//add #rwgt info from comments
const std::vector<std::string> &comments = event->getComments();
for (unsigned i = 0; i < comments.size(); i++) {
if (comments[i].rfind("#rwgt", 0) == 0)
(*infoPtr->eventAttributes)["#rwgt"] = comments[i];
}
const lhef::LHEEvent::PDF *pdf = event->getPDF();
if (pdf) {
this->setPdf(pdf->id.first,
pdf->id.second,
pdf->x.first,
pdf->x.second,
pdf->scalePDF,
pdf->xPDF.first,
pdf->xPDF.second,
true);
} else {
this->setPdf(hepeup.IDUP[0],
hepeup.IDUP[1],
hepeup.PUP[0][3] / runInfo->getHEPRUP()->EBMUP.first,
hepeup.PUP[1][3] / runInfo->getHEPRUP()->EBMUP.second,
0.,
0.,
0.,
false);
}
event->attempted();
return true;
}
|