|
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
|
#include <memory>
#include "GeneratorInterface/Core/interface/GeneratorFilter.h"
#include "GeneratorInterface/ExternalDecays/interface/ExternalDecayDriver.h"
#include "GeneratorInterface/Pythia8Interface/interface/Py8GunBase.h"
namespace gen {
class Py8MassGun : public Py8GunBase {
public:
Py8MassGun(edm::ParameterSet const&);
~Py8MassGun() override {}
bool generatePartonsAndHadronize() override;
const char* classname() const override;
private:
// PtGun particle(s) characteristics
double fMinEta;
double fMaxEta;
double fMinP;
double fMaxP;
double fMinPt;
double fMaxPt;
double fMinM;
double fMaxM;
int fMomMode;
};
// implementation
//
Py8MassGun::Py8MassGun(edm::ParameterSet const& ps) : Py8GunBase(ps) {
// ParameterSet defpset ;
edm::ParameterSet pgun_params = ps.getParameter<edm::ParameterSet>("PGunParameters"); // , defpset ) ;
fMinEta = pgun_params.getParameter<double>("MinEta"); // ,-2.2);
fMaxEta = pgun_params.getParameter<double>("MaxEta"); // , 2.2);
fMinP = pgun_params.getParameter<double>("MinP"); // , 0.);
fMaxP = pgun_params.getParameter<double>("MaxP"); // , 0.);
fMinPt = pgun_params.getParameter<double>("MinPt"); // , 0.);
fMaxPt = pgun_params.getParameter<double>("MaxPt"); // , 0.);
fMinM = pgun_params.getParameter<double>("MinM"); // , 0.);
fMaxM = pgun_params.getParameter<double>("MaxM"); // , 0.);
fMomMode = pgun_params.getParameter<int>("MomMode"); // , 1);
}
bool Py8MassGun::generatePartonsAndHadronize() {
fMasterGen->event.reset();
size_t pSize = fPartIDs.size();
if (pSize > 2)
return false;
int NTotParticles = fPartIDs.size();
// energy below is dummy, it is not used
(fMasterGen->event).append(990, -11, 0, 0, 2, 1 + NTotParticles, 0, 0, 0., 0., 0., 15000., 15000.);
// Pick a flat mass range
double phi, eta, the, ee, pp;
double m0 = (fMaxM - fMinM) * randomEngine().flat() + fMinM;
// Global eta
eta = (fMaxEta - fMinEta) * randomEngine().flat() + fMinEta;
if (pSize == 2) {
// Masses.
double m1 = fMasterGen->particleData.m0(fPartIDs[0]);
double m2 = fMasterGen->particleData.m0(fPartIDs[1]);
// Energies and absolute momentum in the rest frame.
if (m1 + m2 > m0)
return false;
double e1 = 0.5 * (m0 * m0 + m1 * m1 - m2 * m2) / m0;
double e2 = 0.5 * (m0 * m0 + m2 * m2 - m1 * m1) / m0;
double pAbs = 0.5 * sqrt((m0 - m1 - m2) * (m0 + m1 + m2) * (m0 + m1 - m2) * (m0 - m1 + m2)) / m0;
// Isotropic angles in rest frame give three-momentum.
double cosTheta = 2. * randomEngine().flat() - 1.;
double sinTheta = sqrt(1. - cosTheta * cosTheta);
phi = 2. * M_PI * randomEngine().flat();
double pX = pAbs * sinTheta * cos(phi);
double pY = pAbs * sinTheta * sin(phi);
double pZ = pAbs * cosTheta;
(fMasterGen->event).append(fPartIDs[0], 1, 1, 0, 0, 0, 0, 0, pX, pY, pZ, e1, m1);
(fMasterGen->event).append(fPartIDs[1], 1, 1, 0, 0, 0, 0, 0, -pX, -pY, -pZ, e2, m2);
} else {
(fMasterGen->event).append(fPartIDs[0], 1, 1, 0, 0, 0, 0, 0, 0.0, 0.0, 0.0, m0, m0);
}
//now the boost (from input params)
if (fMomMode == 0) {
pp = (fMaxP - fMinP) * randomEngine().flat() + fMinP;
} else {
double pT = (fMaxPt - fMinPt) * randomEngine().flat() + fMinPt;
pp = pT * cosh(eta);
}
ee = sqrt(m0 * m0 + pp * pp);
//the boost direction (from input params)
//
phi = (fMaxPhi - fMinPhi) * randomEngine().flat() + fMinPhi;
the = 2. * atan(exp(-eta));
double betaX = pp / ee * std::sin(the) * std::cos(phi);
double betaY = pp / ee * std::sin(the) * std::sin(phi);
double betaZ = pp / ee * std::cos(the);
// boost all particles
//
(fMasterGen->event).bst(betaX, betaY, betaZ);
if (!fMasterGen->next())
return false;
event() = std::make_unique<HepMC::GenEvent>();
return toHepMC.fill_next_event(fMasterGen->event, event().get());
}
const char* Py8MassGun::classname() const { return "Py8MassGun"; }
typedef edm::GeneratorFilter<gen::Py8MassGun, gen::ExternalDecayDriver> Pythia8MassGun;
} // namespace gen
using gen::Pythia8MassGun;
DEFINE_FWK_MODULE(Pythia8MassGun);
|