« back to all changes in this revision
Viewing changes to pwhg_analysis-phiB.f
- browse files at revision 116
- compare with another revision
- download diff
- download tarball
- view history from revision 116
- Committer: Rikkert Frederix
- Date: 2017-06-19 18:34:45 UTC
- mfrom: (112.1.26 PWHG_STJ)
- Revision ID: frederix@physik.uzh.ch-20170619183445-5832xjj3q30gdr9f
merge with latest PWHG_STJ
- files added:
- auxiliary-progress.f
- reweight
- reweight/obscorrelations
- files modified:
- Born_phsp.f
added
removed
pwhg_analysis-phiB.f
1
1
subroutine init_hist
2
2
implicit none
3
include 'pwhg_math.h'
4
integer j,k,i
5
real * 8 dy,dpt,dr,dptzoom,dm,dmzoom
3
include 'pwhg_math.h'
4
real * 8 dy,dpt,dptzoom,dm,dmzoom,ss,dshat
6
5
7
6
call inihists
8
7
9
8
dy=0.5d0
10
dr =0.5d0
11
9
dpt=10d0
12
10
dptzoom=2d0
13
11
dm=10d0
14
12
dmzoom=1d0
15
16
call bookupeqbins('ycm', dy,-5d0, 5d0)
17
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)
22
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)
29
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)
36
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 ...
14
dshat = 14000d0/100d0
15
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)
41
29
42
30
end
43
31
44
32
subroutine analysis(dsig0)
45
use auxiliary
33
use auxiliary
46
34
implicit none
47
35
real * 8 dsig0
48
36
include 'hepevt.h'
49
37
include 'nlegborn.h'
50
38
include 'pwhg_flst.h'
51
include 'pwhg_math.h'
52
include 'pwhg_rad.h'
39
include 'pwhg_math.h'
40
include 'pwhg_rad.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)
72
integer onem,nq
59
integer nfinal,jmerge, kmerge, mergedfl, flav(nlegreal)
60
integer onem,nq
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
78
67
79
68
if(ini) then
80
69
if(weights_num.eq.0) then
96
85
97
86
if(sum(abs(dsig)).eq.0) return
98
87
99
C initialization
100
p = 0d0
101
flav = onem
102
if (nhep > nlegreal) then
88
C initialization
89
p = 0d0
90
flav = onem
91
if (nhep > nlegreal) then
103
92
stop 'pwhg_analysis-phib: nhep> nlegreal'
104
93
endif
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)
94
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)
108
98
do ihep=1,nhep
109
99
if(abs(idhep(ihep)).eq.6.and.isthep(ihep).eq.1) then
110
p_top=phep(1:4,ihep)
100
tMom=phep(1:4,ihep)
111
101
endif
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
115
105
enddo
116
117
C do clustering until two partons are left in the final
118
nfinal = nhep-2
119
do while (nfinal > 2)
120
call findNearestNeighbours(p,flav,jmerge,kmerge,mergedfl,q2merge)
121
if(q2merge.lt.1d10) then
122
if(jmerge.gt.2) then
123
p(:,jmerge)=p(:,jmerge)+p(:,kmerge)
124
else
125
p(3,jmerge)=p(3,jmerge)-p(3,kmerge)
126
p(0,jmerge)=abs(p(3,jmerge))
127
endif
128
flav(kmerge)=onem
129
flav(jmerge)=mergedfl
130
nfinal = nfinal-1
131
else
132
write(*,*) 'flav', flav
133
stop 'no merging took place'
134
endif
135
end do
136
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'
140
nq = 0
106
107
! clustering event
108
call ClusterEvent(p,flav, q2merge,tMom,bMom,WMom,jMom)
109
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'
112
113
nq = 0
141
114
do i=1,nhep
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
143
116
nq = nq+1
144
117
endif
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)
149
118
enddo
150
151
C sanity check that we end up with 4 quarks...
152
if (nq /= 4 ) then
119
C sanity check that we end up with 4 quarks...
120
if (nq /= 4 .and. nq /= 6) then
153
121
write(*,*) 'flav', flav
154
write(*,*) 'nq', nq
155
stop ' number of quarks different from 4'
122
write(*,*) 'nq', nq
123
stop 'number of quarks different from 4 or 6'
156
124
endif
157
125
158
C now compute some Born variables
159
p_w = p_top-p_bot
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)
164
165
C write(*,*) 'p_w', p_w
166
C write(*,*)
167
168
C CM rapidity
169
ycm = 0.5*log(p(0,1)/p(0,2))
170
171
call filld('ycm', ycm, dsig)
172
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)
177
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)
184
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)
191
192
call filld('dy', dy, dsig)
193
call filld('deta', deta, dsig)
194
call filld('dphi', dphi, dsig)
195
call filld('dr', dr, dsig)
196
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)
130
131
C CM rapidity
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))
134
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)
148
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))
152
197
153
end
198
199
200
201
Loggerhead is a web-based interface for Breezy
Version: 2.0.1