|
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
|
#include "Pythia8/Pythia.h"
#include "TF1.h"
class PtHatReweightUserHook : public Pythia8::UserHooks {
public:
PtHatReweightUserHook(double _pt = 15, double _power = 4.5) : pt(_pt), power(_power) {}
~PtHatReweightUserHook() override {}
bool canBiasSelection() override { return true; }
double biasSelectionBy(const Pythia8::SigmaProcess* sigmaProcessPtr,
const Pythia8::PhaseSpace* phaseSpacePtr,
bool inEvent) override {
//the variable selBias of the base class should be used;
if ((sigmaProcessPtr->nFinal() == 2)) {
selBias = pow(phaseSpacePtr->pTHat() / pt, power);
return selBias;
}
selBias = 1.;
return selBias;
}
private:
double pt, power;
};
class PtHatEmpReweightUserHook : public Pythia8::UserHooks {
public:
PtHatEmpReweightUserHook(const std::string& tuneName = "") {
if (tuneName == "CP5" || tuneName == "CP5Run3")
p = {7377.94700788, 8.38168461349, -4.70983112392, -0.0310148108446, -0.028798537937, 925.335472326};
//Default reweighting - works good for tune CUEPT8M1
else
p = {5.3571961909810e+13,
1.0907678218282e+01,
-2.5898069229451e+00,
-5.1575514014931e-01,
5.5951279807561e-02,
3.5e+02};
const double ecms = (tuneName == "CP5Run3" ? 13600. : 13000.);
sigma = [this, ecms](double x) -> double {
return (p[0] * pow(x, p[2] + p[3] * log(0.01 * x) + p[4] * pow(log(0.01 * x), 2)) *
pow(1 - 2 * x / (ecms + p[5]), p[1])) *
x;
};
}
~PtHatEmpReweightUserHook() override {}
bool canBiasSelection() override { return true; }
double biasSelectionBy(const Pythia8::SigmaProcess* sigmaProcessPtr,
const Pythia8::PhaseSpace* phaseSpacePtr,
bool inEvent) override {
//the variable selBias of the base class should be used;
if ((sigmaProcessPtr->nFinal() == 2)) {
selBias = 1.0 / sigma(phaseSpacePtr->pTHat());
return selBias;
}
selBias = 1.;
return selBias;
}
private:
std::vector<double> p;
std::function<double(double)> sigma;
};
class RapReweightUserHook : public Pythia8::UserHooks {
public:
RapReweightUserHook(const std::string& _yLabsigma_func,
double _yLab_power,
const std::string& _yCMsigma_func,
double _yCM_power,
double _pTHatMin,
double _pTHatMax)
: yLabsigma_func(_yLabsigma_func),
yCMsigma_func(_yCMsigma_func),
yLab_power(_yLab_power),
yCM_power(_yCM_power),
pTHatMin(_pTHatMin),
pTHatMax(_pTHatMax) {
// empirical parametrizations defined in configuration file
yLabsigma = TF1("yLabsigma", yLabsigma_func.c_str(), pTHatMin, pTHatMax);
yCMsigma = TF1("yCMsigma", yLabsigma_func.c_str(), pTHatMin, pTHatMax);
}
~RapReweightUserHook() override {}
bool canBiasSelection() override { return true; }
double biasSelectionBy(const Pythia8::SigmaProcess* sigmaProcessPtr,
const Pythia8::PhaseSpace* phaseSpacePtr,
bool inEvent) override {
//the variable selBias of the base class should be used;
if ((sigmaProcessPtr->nFinal() == 2)) {
double x1 = phaseSpacePtr->x1();
double x2 = phaseSpacePtr->x2();
double yLab = 0.5 * log(x1 / x2);
double yCM = 0.5 * log(phaseSpacePtr->tHat() / phaseSpacePtr->uHat());
double pTHat = phaseSpacePtr->pTHat();
double sigmaLab = yLabsigma.Eval(pTHat);
double sigmaCM = yCMsigma.Eval(pTHat);
// empirical reweighting function
selBias = exp(pow(fabs(yLab), yLab_power) / (2 * sigmaLab * sigmaLab) +
pow(fabs(yCM), yCM_power) / (2 * sigmaCM * sigmaCM));
return selBias;
}
selBias = 1.;
return selBias;
}
private:
std::string yLabsigma_func, yCMsigma_func;
double yLab_power, yCM_power, pTHatMin, pTHatMax;
TF1 yLabsigma, yCMsigma;
};
class PtHatRapReweightUserHook : public Pythia8::UserHooks {
public:
PtHatRapReweightUserHook(const std::string& _yLabsigma_func,
double _yLab_power,
const std::string& _yCMsigma_func,
double _yCM_power,
double _pTHatMin,
double _pTHatMax,
double _pt = 15,
double _power = 4.5)
: yLabsigma_func(_yLabsigma_func),
yCMsigma_func(_yCMsigma_func),
yLab_power(_yLab_power),
yCM_power(_yCM_power),
pTHatMin(_pTHatMin),
pTHatMax(_pTHatMax),
pt(_pt),
power(_power) {
// empirical parametrizations defined in configuration file
yLabsigma = TF1("yLabsigma", yLabsigma_func.c_str(), pTHatMin, pTHatMax);
yCMsigma = TF1("yCMsigma", yLabsigma_func.c_str(), pTHatMin, pTHatMax);
}
~PtHatRapReweightUserHook() override {}
bool canBiasSelection() override { return true; }
double biasSelectionBy(const Pythia8::SigmaProcess* sigmaProcessPtr,
const Pythia8::PhaseSpace* phaseSpacePtr,
bool inEvent) override {
//the variable selBias of the base class should be used;
if ((sigmaProcessPtr->nFinal() == 2)) {
double x1 = phaseSpacePtr->x1();
double x2 = phaseSpacePtr->x2();
double yLab = 0.5 * log(x1 / x2);
double yCM = 0.5 * log(phaseSpacePtr->tHat() / phaseSpacePtr->uHat());
double pTHat = phaseSpacePtr->pTHat();
double sigmaLab = yLabsigma.Eval(pTHat);
double sigmaCM = yCMsigma.Eval(pTHat);
// empirical reweighting function
selBias = pow(pTHat / pt, power) * exp(pow(fabs(yLab), yLab_power) / (2 * sigmaLab * sigmaLab) +
pow(fabs(yCM), yCM_power) / (2 * sigmaCM * sigmaCM));
return selBias;
}
selBias = 1.;
return selBias;
}
private:
std::string yLabsigma_func, yCMsigma_func;
double yLab_power, yCM_power, pTHatMin, pTHatMax, pt, power;
TF1 yLabsigma, yCMsigma;
};
|