|
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
|
#include <memory>
#include "GeneratorInterface/Pythia8Interface/interface/Py8GunBase.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Concurrency/interface/SharedResourceNames.h"
// EvtGen plugin
//
#include "Pythia8Plugins/EvtGen.h"
using namespace Pythia8;
const std::vector<std::string> gen::Py8GunBase::p8SharedResources = {edm::SharedResourceNames::kPythia8};
namespace gen {
Py8GunBase::Py8GunBase(edm::ParameterSet const& ps) : Py8InterfaceBase(ps) {
// PGun specs
//
edm::ParameterSet pgun_params = ps.getParameter<edm::ParameterSet>("PGunParameters");
// although there's the method ParameterSet::empty(),
// it looks like it's NOT even necessary to check if it is,
// before trying to extract parameters - if it is empty,
// the default values seem to be taken
//
fPartIDs = pgun_params.getParameter<std::vector<int> >("ParticleID");
fMinPhi = pgun_params.getParameter<double>("MinPhi"); // ,-3.14159265358979323846);
fMaxPhi = pgun_params.getParameter<double>("MaxPhi"); // , 3.14159265358979323846);
}
// specific to Py8GunHad !!!
//
bool Py8GunBase::initializeForInternalPartons() {
// NO MATTER what was this setting below, override it before init
// - this is essencial for the PGun mode
// Key requirement: switch off ProcessLevel, and thereby also PartonLevel.
fMasterGen->readString("ProcessLevel:all = off");
fMasterGen->readString("ProcessLevel::resonanceDecays=on");
fMasterGen->init();
// init decayer
fDecayer->readString("ProcessLevel:all = off"); // Same trick!
fDecayer->readString("ProcessLevel::resonanceDecays=on");
fDecayer->init();
if (useEvtGen) {
edm::LogInfo("Pythia8Interface") << "Creating and initializing pythia8 EvtGen plugin";
evtgenDecays = std::make_shared<EvtGenDecays>(fMasterGen.get(), evtgenDecFile, evtgenPdlFile);
for (unsigned int i = 0; i < evtgenUserFiles.size(); i++)
evtgenDecays->readDecayFile(evtgenUserFiles.at(i));
}
return true;
}
bool Py8GunBase::residualDecay() {
Event* pythiaEvent = &(fMasterGen->event);
int NPartsBeforeDecays = pythiaEvent->size() - 1; // do NOT count the very 1st "system" particle
// in Pythia8::Event record; it does NOT even
// get translated by the HepMCInterface to the
// HepMC::GenEvent record !!!
int NPartsAfterDecays = event().get()->particles_size();
if (NPartsAfterDecays == NPartsBeforeDecays)
return true;
bool result = true;
for (int ipart = NPartsAfterDecays; ipart > NPartsBeforeDecays; ipart--) {
HepMC::GenParticle* part = event().get()->barcode_to_particle(ipart);
if (part->status() == 1 && (fDecayer->particleData).canDecay(part->pdg_id())) {
fDecayer->event.reset();
Particle py8part(part->pdg_id(),
93,
0,
0,
0,
0,
0,
0,
part->momentum().x(),
part->momentum().y(),
part->momentum().z(),
part->momentum().t(),
part->generated_mass());
HepMC::GenVertex* ProdVtx = part->production_vertex();
py8part.vProd(
ProdVtx->position().x(), ProdVtx->position().y(), ProdVtx->position().z(), ProdVtx->position().t());
py8part.tau((fDecayer->particleData).tau0(part->pdg_id()));
fDecayer->event.append(py8part);
int nentries = fDecayer->event.size();
if (!fDecayer->event[nentries - 1].mayDecay())
continue;
fDecayer->next();
int nentries1 = fDecayer->event.size();
if (nentries1 <= nentries)
continue; //same number of particles, no decays...
part->set_status(2);
result = toHepMC.fill_next_event(*(fDecayer.get()), event().get(), -1, true, part);
}
}
return result;
}
void Py8GunBase::finalizeEvent() {
//******** Verbosity ********
if (maxEventsToPrint > 0 && (pythiaPylistVerbosity || pythiaHepMCVerbosity || pythiaHepMCVerbosityParticles)) {
maxEventsToPrint--;
if (pythiaPylistVerbosity) {
fMasterGen->info.list();
fMasterGen->event.list();
}
if (pythiaHepMCVerbosity) {
std::cout << "Event process = " << fMasterGen->info.code() << "\n"
<< "----------------------" << std::endl;
event()->print();
}
if (pythiaHepMCVerbosityParticles) {
std::cout << "Event process = " << fMasterGen->info.code() << "\n"
<< "----------------------" << std::endl;
ascii_io->write_event(event().get());
}
}
return;
}
void Py8GunBase::statistics() {
fMasterGen->stat();
double xsec = fMasterGen->info.sigmaGen(); // cross section in mb
xsec *= 1.0e9; // translate to pb (CMS/Gen "convention" as of May 2009)
runInfo().setInternalXSec(xsec);
return;
}
void Py8GunBase::evtGenDecay() {
if (evtgenDecays.get())
evtgenDecays->decay();
}
} // namespace gen
|