-
Notifications
You must be signed in to change notification settings - Fork 6
/
Copy pathbtilde.f
378 lines (373 loc) · 13 KB
/
btilde.f
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
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
function btilde(xx,www0,ifirst,imode,
1 retval,retval0)
c retval is the function return value
c retval0 is an 'avatar' function the has similar value, but is much
c easier to compute (i.e. the Born term in this case)
c imode = 0 compute only the Born contribution (retval = retval0 in this case)
c imode = 1 compute everything
c return value: 0: success; 1: retval0 was not computed
c (for functions that do not support an avatar function)
c here we always return 0.
implicit none
include 'nlegborn.h'
include 'pwhg_flst.h'
include 'pwhg_kn.h'
include 'pwhg_rad.h'
include 'pwhg_flg.h'
include 'pwhg_math.h'
c independent variables for real graph: number of final state
c legs times 3, take away 4 for 4-momentum conservation, add 2
c for x_1 and x_2, and take away an overall azimuth
real * 8 xx(ndiminteg),www0,retval,retval0
real * 8 xrad(3)
real * 8 xborn(ndiminteg-3)
integer btilde,ifirst,imode,iret
real * 8 resborn(maxprocborn),resvirt(maxprocborn),
# resreal(maxprocborn),rescoll(maxprocborn)
real * 8 results(maxprocborn)
real * 8 tmp,suppfact,www,wwwtot
integer j
save results,resborn,resvirt,wwwtot,suppfact
real * 8 seconds
real *8 totborn,totvirt,ptotborn
logical pwhg_isfinite
external pwhg_isfinite
real * 8 powheginput
logical ini,btildebornon,btildevirton,btildecollon,btilderealon
data ini/.true./
save ini,btildebornon,btildevirton,btildecollon,btilderealon
if(ini) then
btildebornon = .not.(powheginput("#btildeborn").eq.0)
btildevirton = .not.(powheginput("#btildevirt").eq.0)
btildecollon = .not.(powheginput("#btildecoll").eq.0)
btilderealon = .not.(powheginput("#btildereal").eq.0)
ini = .false.
endif
btilde=0
www=www0*hc2
do j=1,ndiminteg-3
xborn(j)=xx(j)
enddo
do j=1,3
xrad(j)=xx(ndiminteg-3 + j)
enddo
if(ifirst.eq.0) then
wwwtot=www
c sets born momenta in kin. common block
call gen_born_phsp(xborn)
call born_suppression(suppfact)
c set scales
call setscalesbtilde
call allborn
c sets xscaled, y, phi in kinematics common block
call btildeborn(resborn)
if(.not.btildebornon) resborn = 0
if (.not.flg_bornonly.and..not.imode.eq.0) then
call reset_timer
if(btildevirton .and. kn_jacborn /= 0) then
call btildevirt(resvirt)
else
resvirt = 0
endif
call get_timer(seconds)
call addtocnt('virt time (sec)',seconds)
call reset_timer
if(btildecollon .and. kn_jacborn /= 0) then
call btildecoll(xrad,rescoll,www)
else
rescoll = 0
endif
if(btilderealon .and. kn_jacborn /= 0) then
call btildereal(xrad,resreal,www)
else
resreal = 0
endif
call get_timer(seconds)
call addtocnt('real time (sec)',seconds)
endif
c accumulate values
retval=0
do j=1,flst_nborn
c jacobians are already included in rescoll and resreal
tmp=resborn(j)
if (.not.flg_bornonly.and..not.imode.eq.0) then
tmp = tmp
# + resvirt(j)
# + rescoll(j)
# + resreal(j)
endif
c initial value in results
results(j)=tmp*www*suppfact
retval=retval+tmp*www*suppfact
enddo
elseif(ifirst.eq.1) then
c subsequent calls:
c In case of folding the call to btildeborn and btildevirt can be
c avoided, since results are the same.
c If the NLO calculation is performed also (flg_nlotest is set)
c we need to accumulate all weight within a single folding sequence
c in order to later output the correct Born and Virtual contribution
c to the NLO analysis routine.
wwwtot=wwwtot+www
if (.not.flg_bornonly.and..not.imode.eq.0) then
c btildecoll and btildereal take care themselves to invoke the NLO
c analysis if required.
call reset_timer
c in case btlscalereal is set, we need to reset the scales
c to the underlying Born value for the computation of the
c collinear remnants.
call setscalesbtilde
if(btildecollon .and. kn_jacborn /= 0) then
call btildecoll(xrad,rescoll,www)
else
rescoll = 0
endif
if(btilderealon .and. kn_jacborn /= 0) then
call btildereal(xrad,resreal,www)
else
resreal = 0
endif
call get_timer(seconds)
call addtocnt('real time (sec)',seconds)
endif
retval=0
do j=1,flst_nborn
tmp=resborn(j)
if (.not.flg_bornonly.and..not.imode.eq.0) then
tmp = tmp
# + resvirt(j)
# + rescoll(j)
# + resreal(j)
endif
c accumulate values in results
results(j)=results(j)+tmp*www*suppfact
retval=retval+tmp*www*suppfact
enddo
elseif(ifirst.eq.2) then
totborn=0d0
c compute Born, to return in retval0
ptotborn=0
do j=1,flst_nborn
totborn=totborn+resborn(j)
ptotborn=ptotborn+abs(resborn(j))
enddo
totborn=totborn*wwwtot
ptotborn=ptotborn*wwwtot
if (.not.pwhg_isfinite(totborn)) then
totborn = 0d0
ptotborn = 0d0
resborn = 0d0
endif
if(flg_nlotest) then
c output Born
if(.not.imode.eq.0) then
call analysis_extrainfo('born',flst_nborn,resborn,wwwtot)
call analysis_driver(totborn,0)
endif
if(.not.flg_bornonly.and..not.imode.eq.0) then
c output virtual
totvirt=0d0
do j=1,flst_nborn
totvirt=totvirt+resvirt(j)
enddo
totvirt=totvirt*wwwtot
if (.not.pwhg_isfinite(totvirt)) then
totvirt = 0d0
resvirt = 0d0
endif
call analysis_extrainfo('virt',flst_nborn,resvirt,wwwtot)
call analysis_driver(totvirt,0)
endif
call pwhgaccumup
endif
c Make the born part of the result available;
c (to test, if bornonly is set, should equal the output btilde when ifirst=2)
retval0=ptotborn*suppfact
btilde=0
c closing call to end a sequence of correlated events in the
c analysis routines.
c closing call: accumulate values with correct signs
retval=0
call adduptotals(results,flst_nborn)
do j=1,flst_nborn
c this is only useful if withnegweights on (i.e. =1 in powheg.input,
c logical true here). However, better set a default (Les Houches
c interface will simply output this sign for the event.)
rad_btilde_sign(j)=1
if(flg_withnegweights) then
if(results(j).lt.0) then
results(j)=-results(j)
rad_btilde_sign(j)=-1
endif
else
if(results(j).lt.0) then
results(j)=0
endif
endif
retval=retval+results(j)
c Transfer all flavour components of btilde to the array
c in common block; will be used to decide the underlying
c flavour of the event
rad_btilde_arr(j)=results(j)
enddo
else
write(*,*) 'wrong value of ifirst in btilde => ',ifirst
call exit(-1)
endif
end
subroutine adduptotals(results,n)
implicit none
include 'nlegborn.h'
integer n
real * 8 results(n)
real * 8 tot,totabs,totpos,totneg,etot,etotabs,etotpos,etotneg
real * 8 totj(maxprocborn),totabsj(maxprocborn),
1 totposj(maxprocborn),totnegj(maxprocborn),
2 etotj(maxprocborn),etotabsj(maxprocborn),
3 etotposj(maxprocborn),etotnegj(maxprocborn)
integer nentries
common/cadduptotals/tot,totabs,totpos,totneg,etot,etotabs,
1 etotpos,etotneg,totj,totabsj,totposj,totnegj,
1 etotj,etotabsj,etotposj,etotnegj,nentries
real * 8 dtot,dtotabs,dtotpos,dtotneg
integer j
nentries=nentries+1
dtot=0
dtotabs=0
dtotpos=0
dtotneg=0
do j=1,n
dtot=dtot+results(j)
dtotabs=dtotabs+abs(results(j))
if(results(j).gt.0) then
dtotpos=dtotpos+results(j)
else
dtotneg=dtotneg-results(j)
endif
enddo
tot=tot+dtot
totabs=totabs+dtotabs
totpos=totpos+dtotpos
totneg=totneg+dtotneg
etot=etot+dtot**2
etotabs=etotabs+dtotabs**2
etotpos=etotpos+dtotpos**2
etotneg=etotneg+dtotneg**2
c j contributions
do j=1,n
dtot=results(j)
dtotabs=abs(results(j))
if(results(j).gt.0) then
dtotpos=results(j)
else
dtotneg=-results(j)
endif
totj(j)=totj(j)+dtot
totabsj(j)=totabsj(j)+dtotabs
totposj(j)=totposj(j)+dtotpos
totnegj(j)=totnegj(j)+dtotneg
etotj(j)=etotj(j)+dtot**2
etotabsj(j)=etotabsj(j)+dtotabs**2
etotposj(j)=etotposj(j)+dtotpos**2
etotnegj(j)=etotnegj(j)+dtotneg**2
enddo
end
subroutine resettotals
implicit none
include 'nlegborn.h'
include 'pwhg_flst.h'
real * 8 tot,totabs,totpos,totneg,etot,etotabs,etotpos,etotneg
real * 8 totj(maxprocborn),totabsj(maxprocborn),
1 totposj(maxprocborn),totnegj(maxprocborn),
2 etotj(maxprocborn),etotabsj(maxprocborn),
3 etotposj(maxprocborn),etotnegj(maxprocborn)
integer nentries
common/cadduptotals/tot,totabs,totpos,totneg,etot,etotabs,
1 etotpos,etotneg,totj,totabsj,totposj,totnegj,
1 etotj,etotabsj,etotposj,etotnegj,nentries
integer j
nentries=0
tot=0
etot=0
totabs=0
etotabs=0
totpos=0
etotpos=0
totneg=0
etotneg=0
do j=1,flst_nborn
totj(j)=0
etotj(j)=0
totabsj(j)=0
etotabsj(j)=0
totposj(j)=0
etotposj(j)=0
totnegj(j)=0
etotnegj(j)=0
enddo
end
subroutine finaltotals
implicit none
include 'nlegborn.h'
include 'pwhg_flst.h'
include 'pwhg_rad.h'
real * 8 tot,totabs,totpos,totneg,etot,etotabs,etotpos,etotneg
real * 8 totj(maxprocborn),totabsj(maxprocborn),
1 totposj(maxprocborn),totnegj(maxprocborn),
2 etotj(maxprocborn),etotabsj(maxprocborn),
3 etotposj(maxprocborn),etotnegj(maxprocborn)
integer nentries
common/cadduptotals/tot,totabs,totpos,totneg,etot,etotabs,
1 etotpos,etotneg,totj,totabsj,totposj,totnegj,
1 etotj,etotabsj,etotposj,etotnegj,nentries
integer n,j,k
character * 80 format
real * 8 tmp_totbtlj,tmp_etotbtlj,tmp_totabsbtlj,
1 tmp_etotabsbtlj,tmp_totposbtlj,tmp_etotposbtlj,
2 tmp_totnegbtlj,tmp_etotnegbtlj
real * 8 tmp
integer iun
character * 20 pwgprefix
integer lprefix
common/cpwgprefix/pwgprefix,lprefix
real * 8 powheginput
external powheginput
n=nentries
rad_totbtl=tot/n
rad_etotbtl=sqrt((etot/n-(tot/n)**2)/n)
rad_totabsbtl=totabs/n
rad_etotabsbtl=sqrt((etotabs/n-(totabs/n)**2)/n)
rad_totposbtl=totpos/n
rad_etotposbtl=sqrt((etotpos/n-(totpos/n)**2)/n)
rad_totnegbtl=totneg/n
rad_etotnegbtl=sqrt((etotneg/n-(totneg/n)**2)/n)
write(*,*) 'tot:',rad_totbtl,'+-',rad_etotbtl
write(*,*) 'abs:',rad_totabsbtl,'+-',rad_etotabsbtl
write(*,*) 'pos:',rad_totposbtl,'+-',rad_etotposbtl
write(*,*) 'neg:',rad_totnegbtl,'+-',rad_etotnegbtl
c
if(powheginput('#ubsigmadetails').eq.1) then
call newunit(iun)
open(iun,file=pwgprefix(1:lprefix)//'ubsigma.dat')
format='( (i8,1x),4(a,1x,d10.4,a,d7.1))'
write(format(2:4),'(i3)') nlegborn
tmp=0
do j=1,flst_nborn
tmp_totbtlj=totj(j)/n
tmp_etotbtlj=sqrt((etotj(j)/n-(totj(j)/n)**2)/n)
tmp_totabsbtlj=totabsj(j)/n
tmp_etotabsbtlj=sqrt((etotabsj(j)/n-(totabsj(j)/n)**2)/n)
tmp_totposbtlj=totposj(j)/n
tmp_etotposbtlj=sqrt((etotposj(j)/n-(totposj(j)/n)**2)/n)
tmp_totnegbtlj=totnegj(j)/n
tmp_etotnegbtlj=sqrt((etotnegj(j)/n-(totnegj(j)/n)**2)/n)
write(iun,format) (flst_born(k,j),k=1,nlegborn),
1 'tot:',tmp_totbtlj,' +- ',tmp_etotbtlj,
2 '; abs:',tmp_totabsbtlj,' +- ',tmp_etotabsbtlj,
3 '; pos:',tmp_totposbtlj,' +- ',tmp_etotposbtlj,
4 '; neg:',tmp_totnegbtlj,' +- ',tmp_etotnegbtlj
tmp=tmp+tmp_totbtlj
enddo
write(iun,*) tmp
endif
end