5
real * 8 dy,dpt,dr,dptzoom,dm,dmzoom
4
real * 8 dy,dpt,dptzoom,dm,dmzoom,ss,dshat
16
call bookupeqbins('ycm', dy,-5d0, 5d0)
18
call bookupeqbins('y_top', dy,-5d0, 5d0)
19
call bookupeqbins('eta_top',dy,-5d0, 5d0)
20
call bookupeqbins('pt_top', dpt,0d0,500d0)
21
call bookupeqbins('pt_top_zoom', dptzoom,0d0,100d0)
23
call bookupeqbins('y_jet', dy,-5d0, 5d0)
24
call bookupeqbins('eta_jet',dy,-5d0, 5d0)
25
call bookupeqbins('pt_jet', dpt,0d0,500d0)
26
call bookupeqbins('m_jet', dm,0d0,500d0)
27
call bookupeqbins('pt_jet_zoom', dpt,0d0,50d0)
28
call bookupeqbins('m_jet_zoom', dm,0d0,50d0)
30
C call bookupeqbins('y_w', dy,-5d0, 5d0)
31
call bookupeqbins('eta_w', dy,-5d0, 5d0)
32
call bookupeqbins('pt_w', dpt,0d0,500d0)
33
call bookupeqbins('m_w', dm,0d0,500d0)
34
call bookupeqbins('pt_w_zoom', dptzoom,0d0,100d0)
35
call bookupeqbins('m_w_zoom', dmzoom,0d0,100d0)
37
call bookupeqbins('dy', dy,-10d0, 10d0)
38
call bookupeqbins('deta', dy,-10d0, 10d0)
39
call bookupeqbins('dphi', pi/20d0,0d0, pi)
40
call bookupeqbins('dr', dr, 0d0, 10d0)
13
ss = 14000d0 ! to be changed ...
16
c Book histograms as usual but with the reweight plugin
17
call reweight_book("total", 1d0, 0d0, 1d0)
18
call reweight_book('ycm', dy,-5d0, 5d0)
19
call reweight_book('sshat', dshat,0d0, ss)
20
call reweight_book('y_top', dy,-5d0, 5d0)
21
call reweight_book('eta_top',dy,-5d0, 5d0)
22
call reweight_book('pt_top', dpt,0d0,500d0)
23
call reweight_book('pt_top_zoom', dptzoom,0d0,200d0)
24
call reweight_book('eta_w', dy,-5d0, 5d0)
25
call reweight_book('pt_w', dpt,0d0,500d0)
26
call reweight_book('m_w', dm,0d0,500d0)
27
call reweight_book('pt_w_zoom', dptzoom,0d0,100d0)
28
call reweight_book('m_w_zoom', dmzoom,0d0,100d0)
44
32
subroutine analysis(dsig0)
49
37
include 'nlegborn.h'
50
38
include 'pwhg_flst.h'
53
41
include 'pwhg_flg.h'
54
42
include 'LesHouches.h'
55
43
include 'pwhg_weights.h'
68
56
external powheginput
69
57
C -----------------------------------
70
58
real * 8 p(0:3,nlegreal), q2merge
71
integer nfinal,jmerge, kmerge, mergedfl, flav(nlegreal)
59
integer nfinal,jmerge, kmerge, mergedfl, flav(nlegreal)
73
61
parameter (onem=1000000)
74
real *8 p_top(0:3),pt_top,y_top,eta_top,m_top
75
real *8 p_jet(0:3),pt_jet,y_jet,eta_jet,m_jet
76
real *8 p_bot(0:3),p_w(0:3),pt_w,y_w,eta_w,m_w
77
real *8 ycm, dy,deta,dphi,dr
62
real *8 tMom(0:3),pt_top,y_top,eta_top,m_top,q2w
63
real *8 bmom(0:3), Wmom(0:3), jMom(0:3)
64
real *8 p_jet(0:3),pt_jet,y_jet,eta_jet,m_jet
65
real *8 p_bot(0:3),p_w(0:3),pt_w,y_w,eta_w,m_w
66
real *8 ycm, dy,deta,dphi,dr,sshat
80
69
if(weights_num.eq.0) then
97
86
if(sum(abs(dsig)).eq.0) return
102
if (nhep > nlegreal) then
91
if (nhep > nlegreal) then
103
92
stop 'pwhg_analysis-phib: nhep> nlegreal'
105
C make a copy of phep with energy in the standard position
106
C get flavour and assign 0 to gluon
107
flav(1:nhep) = idhep(1:nhep)
95
C make a copy of phep with energy in the standard position
96
C get flavour and assign 0 to gluon
97
flav(1:nhep) = idhep(1:nhep)
109
99
if(abs(idhep(ihep)).eq.6.and.isthep(ihep).eq.1) then
112
102
p(0,ihep) = phep(4,ihep)
113
103
p(1:3,ihep) = phep(1:3,ihep)
114
if (flav(ihep) == 21) flav(ihep) = 0 ! set flavour of gluon to 0
104
if (flav(ihep) == 21) flav(ihep) = 0 ! set flavour of gluon to 0
117
C do clustering until two partons are left in the final
119
do while (nfinal > 2)
120
call findNearestNeighbours(p,flav,jmerge,kmerge,mergedfl,q2merge)
121
if(q2merge.lt.1d10) then
123
p(:,jmerge)=p(:,jmerge)+p(:,kmerge)
125
p(3,jmerge)=p(3,jmerge)-p(3,kmerge)
126
p(0,jmerge)=abs(p(3,jmerge))
129
flav(jmerge)=mergedfl
132
write(*,*) 'flav', flav
133
stop 'no merging took place'
137
! sanity check that we end up in t-channel single-top Born configuration
138
if (.not. any(abs(flav(1:2)) ==5 )) stop 'no b in initial state'
139
if (.not. abs(flav(3)) ==6) stop 'no top in position 3'
108
call ClusterEvent(p,flav, q2merge,tMom,bMom,WMom,jMom)
110
! sanity check that we end up in t-channel single-top Born configuration
111
if (.not. abs(flav(3)) ==6) stop 'no top in position 3'
142
if (abs(flav(i)) > 0 .and. abs(flav(i)) .ne. onem) then
115
if (abs(flav(i)) > 0 .and. abs(flav(i)) .ne. onem) then
145
C find top, bottom and light jet momentum
146
if (abs(flav(i)) == 6) p_top = p(:,i)
147
if (abs(flav(i)) == 5) p_bot = p(:,i)
148
if (i.gt.4 .and. flav(i).ne. onem) p_jet = p(:,i)
151
C sanity check that we end up with 4 quarks...
119
C sanity check that we end up with 4 quarks...
120
if (nq /= 4 .and. nq /= 6) then
153
121
write(*,*) 'flav', flav
155
stop ' number of quarks different from 4'
123
stop 'number of quarks different from 4 or 6'
158
C now compute some Born variables
160
call getyetaptmass(p_top,y_top,eta_top,pt_top,m_top)
161
call getyetaptmass(p_jet,y_jet,eta_jet,pt_jet,m_jet)
162
call getyetaptmass(p_w,y_w,eta_w,pt_w,m_w)
163
call getdydetadphidr(p_top,p_jet,dy,deta,dphi,dr)
165
C write(*,*) 'p_w', p_w
169
ycm = 0.5*log(p(0,1)/p(0,2))
171
call filld('ycm', ycm, dsig)
173
call filld('y_top', y_top, dsig)
174
call filld('eta_top', eta_top, dsig)
175
call filld('pt_top', pt_top, dsig)
176
call filld('pt_top_zoom', pt_top, dsig)
178
call filld('y_jet', y_jet, dsig)
179
call filld('eta_jet', eta_jet, dsig)
180
call filld('pt_jet', pt_jet, dsig)
181
call filld('m_jet', m_jet, dsig)
182
call filld('pt_jet_zoom', pt_jet, dsig)
183
call filld('m_jet_zoom', m_jet, dsig)
185
C call filld('y_w', y_w, dsig)
186
call filld('eta_w', eta_w, dsig)
187
call filld('pt_w', pt_w, dsig)
188
call filld('m_w', m_w, dsig)
189
call filld('pt_w_zoom', pt_w, dsig)
190
call filld('m_w_zoom', m_w, dsig)
192
call filld('dy', dy, dsig)
193
call filld('deta', deta, dsig)
194
call filld('dphi', dphi, dsig)
195
call filld('dr', dr, dsig)
126
C now compute some Born variables
127
call getyetaptmass(tMom,y_top,eta_top,pt_top,m_top)
128
call getyetaptmass(jMom,y_jet,eta_jet,pt_jet,m_jet)
129
call getyetaptmass(WMom,y_w,eta_w,pt_w,m_w)
132
ycm = 0.5*log(p(0,1)/p(0,2))
133
sshat = sqrt(p(0,1)*p(0,2)-p(1,1)*p(1,2)-p(2,1)*p(2,2)-p(3,1)*p(3,2))
135
C Fill histograms as usual
136
call reweight_fill('total', 0.5d0, dsig)
137
call reweight_fill('ycm', ycm, dsig)
138
call reweight_fill('sshat', sshat, dsig)
139
call reweight_fill('y_top', y_top, dsig)
140
call reweight_fill('eta_top', eta_top, dsig)
141
call reweight_fill('pt_top', pt_top, dsig)
142
call reweight_fill('pt_top_zoom', pt_top, dsig)
143
call reweight_fill('eta_w', eta_w, dsig)
144
call reweight_fill('pt_w', pt_w, dsig)
145
call reweight_fill('m_w', m_w, dsig)
146
call reweight_fill('pt_w_zoom', pt_w, dsig)
147
call reweight_fill('m_w_zoom', m_w, dsig)
149
C Collect the log term and dump to file
150
Q2W=abs(WMom(0)**2-WMom(1)**2-WMom(2)**2-WMom(3)**2)
151
call reweight_store(dsig, log(q2merge/Q2W))