|
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
|
// Hook for setting shower scale in top and W resonances
// for Powheg ttb_NLO_dec and b_bbar_4l processes
// C++ port of algorithm by Jezo et. al. (arXiv:1607.04538, Appendix B.2)
#include "Pythia8/Pythia.h"
using namespace Pythia8;
#include "GeneratorInterface/Pythia8Interface/plugins/PowhegResHook.h"
double PowhegResHook::scaleResonance(const int iRes, const Event& event) {
calcScales_ = settingsPtr->flag("POWHEGres:calcScales");
double scale = 0.;
int nDau = event[iRes].daughterList().size();
if (!calcScales_ or nDau == 0) {
// No resonance found, set scale to high value
// Pythia will shower any MC generated resonance unrestricted
scale = 1e30;
}
else if (nDau < 3) {
// No radiating resonance found
scale = 0.8;
}
else if (abs(event[iRes].id()) == 6) {
// Find top daughters
int idw = -1, idb = -1, idg = -1;
for (int i = 0; i < nDau; i++) {
int iDau = event[iRes].daughterList()[i];
if (abs(event[iDau].id()) == 24)
idw = iDau;
if (abs(event[iDau].id()) == 5)
idb = iDau;
if (abs(event[iDau].id()) == 21)
idg = iDau;
}
// Get daughter 4-vectors in resonance frame
Vec4 pw(event[idw].p());
pw.bstback(event[iRes].p());
Vec4 pb(event[idb].p());
pb.bstback(event[iRes].p());
Vec4 pg(event[idg].p());
pg.bstback(event[iRes].p());
// Calculate scale
scale = sqrt(2 * pg * pb * pg.e() / pb.e());
}
else if (abs(event[iRes].id()) == 24) {
// Find W daughters
int idq = -1, ida = -1, idg = -1;
for (int i = 0; i < nDau; i++) {
int iDau = event[iRes].daughterList()[i];
if (event[iDau].id() == 21)
idg = iDau;
else if (event[iDau].id() > 0)
idq = iDau;
else if (event[iDau].id() < 0)
ida = iDau;
}
// Get daughter 4-vectors in resonance frame
Vec4 pq(event[idq].p());
pq.bstback(event[iRes].p());
Vec4 pa(event[ida].p());
pa.bstback(event[iRes].p());
Vec4 pg(event[idg].p());
pg.bstback(event[iRes].p());
// Calculate scale
Vec4 pw = pq + pa + pg;
double q2 = pw * pw;
double csi = 2 * pg.e() / sqrt(q2);
double yq = 1 - pg * pq / (pg.e() * pq.e());
double ya = 1 - pg * pa / (pg.e() * pa.e());
scale = sqrt(min(1 - yq, 1 - ya) * pow2(csi) * q2 / 2);
}
return scale;
}
|