-
Notifications
You must be signed in to change notification settings - Fork 6
/
Copy pathmaxrat.f
706 lines (679 loc) · 22.5 KB
/
maxrat.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
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
subroutine do_maxrat(mcalls,icalls,imode,iret)
c imode = -1: load or generate
c imode = 0: force generating
c imode = 1: force load
implicit none
integer imode,iret
include 'pwhg_math.h'
include 'nlegborn.h'
include 'pwhg_flst.h'
include 'pwhg_flg.h'
include 'pwhg_kn.h'
include 'pwhg_rad.h'
include 'pwhg_rnd.h'
integer mcalls,icalls
character * 20 pwgprefix
integer lprefix
common/cpwgprefix/pwgprefix,lprefix
character * 4 chseed
integer nubound,ios,iun,nynormsold,ncsinormsold,
1 icsi,iy,j,k,jfile,nfiles,fl1,fl2,flemitter
real * 8 csiynorms(rad_ncsiynormsmx,
1 rad_ncsiynormsmx,nlegborn-1,maxprocborn)
logical lpresent,loaded,manyfiles,is_coloured,is_charged
real * 8 powheginput,random
external powheginput,random,is_coloured,is_charged
call zerohistnorms
loaded=.false.
nubound=powheginput('nubound')
rad_nynorms=powheginput('#iymax')
if(rad_nynorms.lt.0) rad_nynorms=1
rad_ncsinorms=powheginput('#icsimax')
if(rad_ncsinorms.lt.0) rad_ncsinorms=1
if(rad_nynorms.gt.rad_ncsiynormsmx.or.
# rad_ncsinorms.gt.rad_ncsiynormsmx) then
write(*,*)
# ' error in poweg.input: rad_nynorms or rad_ncsinorms too large'
stop
endif
rad_normfact=powheginput('#xupbound')
if(rad_normfact.lt.0) rad_normfact=2
if(nubound.eq.0) then
write(*,*) ' ubound set =0, cannot proceed'
call write_counters
call exit(1)
endif
call newunit(iun)
c force generation
if(imode.eq.0) goto 111
c force load
if(imode.eq.1) then
call newunit(iun)
open(unit=iun,file=pwgprefix(1:lprefix)//'ubound.dat',
1 status='old',iostat=ios)
if(ios.eq.0) then
nfiles=1
manyfiles=.false.
else
nfiles=9999
manyfiles=.true.
endif
goto 100
endif
c imode = -1, load or generate
if(powheginput('use-old-ubound').eq.1) then
call newunit(iun)
open(unit=iun,file=pwgprefix(1:lprefix)//'ubound.dat',
1 status='old',iostat=ios)
if(ios.eq.0) then
nfiles=1
manyfiles=.false.
else
if(rnd_cwhichseed.ne.'none') then
c Are we required to generate a grid file or to load a bunch of them?v-
c See if the grid file exists already
inquire(file=pwgprefix(1:lprefix)//'ubound-'//
1 rnd_cwhichseed//'.dat',exist=lpresent)
if(.not.lpresent) then
c we must generate it: go to grid computation
goto 111
endif
endif
c try sequence of number files
nfiles=9999
manyfiles=.true.
endif
else
c new grid required: go to grid computation
goto 111
endif
100 continue
do jfile=1,nfiles
if(nfiles.ne.1) then
c sequence of number files
write(chseed,'(i4)') jfile
do j=1,4
if(chseed(j:j).eq.' ') chseed(j:j)='0'
enddo
inquire(file=pwgprefix(1:lprefix)//'ubound-'//
1 chseed//'.dat',exist=lpresent)
c This one is not present: skip it and go to the next
if(.not.lpresent) goto 110
open(unit=iun,file=pwgprefix(1:lprefix)//'ubound-'//
1 chseed//'.dat', status='old',iostat=ios)
c error in loading; assume that grids have to be redone
if(ios.ne.0) goto 111
write(*,*) ' opened ',pwgprefix(1:lprefix)//'ubound-'//
1 chseed//'.dat'
endif
c rad_ncsinorms, rad_nynorms will be overridden by nynormsold,ncsinormsold
c if all goes well
read(iun,fmt=*,iostat=ios,end=111) nynormsold,ncsinormsold
do j=1,flst_nborn
read(iun,fmt=*,iostat=ios,end=111)
do iy=1,nynormsold
read(iun,fmt=*,iostat=ios,end=111)
do icsi=1,ncsinormsold
read(iun,fmt=*,iostat=ios,end=111)
read(iun,fmt=*,iostat=ios,end=111)
1 (csiynorms(icsi,iy,k,j),k=1,rad_nkinreg)
enddo
enddo
enddo
close(iun)
do iy=1,nynormsold
do icsi=1,ncsinormsold
do j=1,flst_nborn
do k=1,rad_nkinreg
if(jfile.eq.1) then
rad_csiynorms(icsi,iy,k,j)=
1 csiynorms(icsi,iy,k,j)
else
rad_csiynorms(icsi,iy,k,j)=max(
1 rad_csiynorms(icsi,iy,k,j),
2 csiynorms(icsi,iy,k,j))
endif
enddo
enddo
enddo
enddo
c all went well: override rad_ncsinorms, rad_nynorms
rad_ncsinorms=ncsinormsold
rad_nynorms=nynormsold
loaded=.true.
110 continue
enddo
if(loaded) then
write(*,*)
1 ' normalization of upper bounding function for radiation',
2 ' successfully loaded'
iret = 0
goto 998
elseif(imode.eq.1) then
iret = -1
return
endif
c Compute and store normalization grid
111 continue
write(*,*)
write(*,*)
# ' POWHEG: computing normalization of upper'//
# ' bounding function for radiation'
do j=1,flst_nborn
do k=1,rad_nkinreg
do iy=1,rad_nynorms
do icsi=1,rad_ncsinorms
rad_csiynorms(icsi,iy,k,j)=0
enddo
enddo
rad_norms(k,j)=0
enddo
enddo
do j=1,nubound
call gen_born(mcalls,icalls)
call gen_uborn_idx
kn_csitilde=random()
kn_y=2*random()-1
kn_azi=2*pi*random()
do rad_kinreg=1,rad_nkinreg
if(rad_kinreg.eq.1.and.rad_kinreg_on(1)) then
c initial state radiation
fl1=flst_born(1,rad_ubornidx)
fl2=flst_born(2,rad_ubornidx)
if((.not.is_coloured(fl1).and..not.is_coloured(fl2))
1 .and.(is_charged(fl1).or.is_charged(fl2))) then
flg_em_rad = .true.
else
flg_em_rad = .false.
endif
c better set kn_emitter<=2 to avoid troubles ...
kn_emitter = 0
call gen_real_phsp_isr_rad0
call inc_norms
elseif(rad_kinreg_on(rad_kinreg)) then
c final state radiation
kn_emitter=flst_lightpart+rad_kinreg-2
flemitter=flst_born(kn_emitter,rad_ubornidx)
if(.not.is_coloured(flemitter).and.is_charged(flemitter))
1 then
flg_em_rad = .true.
else
flg_em_rad = .false.
endif
call gen_real_phsp_fsr_rad0
call inc_norms
endif
enddo
enddo
call evenmaxrat
call newunit(iun)
if(rnd_cwhichseed.eq.'none') then
open(unit=iun,file=pwgprefix(1:lprefix)//'ubound.dat',
1 status='unknown')
else
open(unit=iun,file=pwgprefix(1:lprefix)//'ubound-'
1 //rnd_cwhichseed//'.dat',status='unknown')
endif
write(iun,*) rad_nynorms,rad_ncsinorms, ' n-y and n-csi grid'
do j=1,flst_nborn
write(iun,fmt=*,iostat=ios) j,' Born index'
do iy=1,rad_nynorms
write(iun,fmt=*,iostat=ios) iy, ' i-y'
do icsi=1,rad_ncsinorms
write(iun,fmt=*,iostat=ios) icsi, ' i-csi'
write(iun,fmt=*,iostat=ios)
#(rad_csiynorms(icsi,iy,k,j),k=1,rad_nkinreg)
enddo
enddo
enddo
close(iun)
write(*,*) ' Normalization of upper bounding function'//
1 ' for radiation computed and stored'
call printhistnorms
998 do j=1,flst_nborn
do k=1,rad_nkinreg
do iy=1,rad_nynorms
do icsi=1,rad_ncsinorms
rad_csiynorms(icsi,iy,k,j)=
# rad_csiynorms(icsi,iy,k,j)*rad_normfact
rad_norms(k,j)=max(rad_norms(k,j),
# rad_csiynorms(icsi,iy,k,j))
enddo
enddo
enddo
enddo
c this subroutine identifies the upper bounds for identical
c contributions
c call maxratident
end
subroutine evenmaxrat
implicit none
include 'nlegborn.h'
include 'pwhg_flst.h'
include 'pwhg_flg.h'
include 'pwhg_rad.h'
integer equivto(maxprocborn)
common/cequivtoborn/equivto
integer j,k,iy,icsi,ir
if(.not.flg_evenmaxrat) return
do j=1,flst_nborn
if(equivto(j).eq.-1) then
c first store in j maximum of norms in equivalente amplitudes
do k=1,flst_nborn
if(equivto(k).eq.j) then
do iy=1,rad_nynorms
do icsi=1,rad_ncsinorms
do ir=1,rad_nkinreg
rad_csiynorms(icsi,iy,ir,j)=
1 max(rad_csiynorms(icsi,iy,ir,j),
2 rad_csiynorms(icsi,iy,ir,k))
enddo
enddo
enddo
endif
enddo
c Now set all equivalent amplitudes equal to j
do k=1,flst_nborn
if(equivto(k).eq.j) then
do iy=1,rad_nynorms
do icsi=1,rad_ncsinorms
do ir=1,rad_nkinreg
rad_csiynorms(icsi,iy,ir,k)=
1 rad_csiynorms(icsi,iy,ir,j)
enddo
enddo
enddo
endif
enddo
endif
enddo
end
subroutine inc_norms
implicit none
include 'pwhg_math.h'
include 'nlegborn.h'
include 'pwhg_flst.h'
include 'pwhg_kn.h'
include 'pwhg_flg.h'
include 'pwhg_rad.h'
real * 8 ptsq,born,sig,xnorm
real * 8 histnorms(nlegborn-1,maxprocborn,100)
real * 8 efficiencies(nlegborn-1,maxprocborn)
common/chistnorms/histnorms,efficiencies
integer icsi,iy,ind,k
logical refuse_pdf
real * 8 pwhg_pt2,pwhg_upperb_rad
external pwhg_pt2,pwhg_upperb_rad
logical debug
parameter (debug=.false.)
ccccccccccccccccccccccccccccccccc
logical verbose
parameter (verbose=.false.)
include 'pwhg_st.h'
include 'pwhg_pdf.h'
include 'pwhg_br.h'
real * 8 pdf1(-pdf_nparton:pdf_nparton),
1 pdf2(-pdf_nparton:pdf_nparton)
real*8 ptsqmin
real *8 kt2minqed
common/showerqed/kt2minqed
cccccccccccccccccccccccccccccccc
ptsq=pwhg_pt2()
c.....mod
if(flg_em_rad) then
ptsqmin = kt2minqed
else
ptsqmin = rad_ptsqmin
endif
c.....mod
if(ptsq.gt.ptsqmin) then
call set_rad_scales(ptsq)
if(.not.refuse_pdf()) then
call sigborn_rad(born)
call sigreal_rad(sig)
sig=sig*kn_jacreal
if(born.le.0) then
if(verbose) then
write(*,*) ' *****************************'
write(*,*) ' inc_norms: Warning, born zero'
write(*,*) ' RAD REG: ', rad_kinreg
write(*,*) ' inc_norms: x1,x2', kn_xb1,kn_xb2
write(*,*) ' muf= ',sqrt(st_mufact2)
write(*,*) ' Flavour: ', (flst_born(ind,rad_ubornidx),ind
$ =1,5)
write(*,*) ' Born w/o pdf: ',br_born(rad_ubornidx)
$ ,' jac: ',kn_jacborn
call pdfcall(1,kn_xb1,pdf1)
call pdfcall(2,kn_xb2,pdf2)
write(*,*) " PDF's ",pdf1(flst_born(1
$ ,rad_ubornidx)),pdf2(flst_born(2,rad_ubornidx))
if ((abs(flst_born(1
$ ,rad_ubornidx)).ge.4).or.(abs(flst_born(2
$ ,rad_ubornidx)).ge.4)) then
write(*,*) ' inc_norms: hvq in initial state!'
else
write(*,*) ' inc_norms: more serious problem!'
endif
write(*,*) ' *****************************'
endif
return
endif
xnorm=sig/born/pwhg_upperb_rad()
if (debug) then
if(xnorm.gt.0.1) then
write(*,*) ' **********************'
write(*,*)
$ ' event with abnormally high ratio to bound:'
$ ,xnorm
call printevent
do k=1,10000
write(*,*) ' enter csitilde,y'
read(*,*) kn_csitilde, kn_y
call gen_real_phsp_isr_rad0
ptsq=pwhg_pt2()
if(ptsq.gt.ptsqmin) then
call set_rad_scales(ptsq)
call sigreal_rad(sig)
sig=sig*kn_jacreal
xnorm=sig/born/pwhg_upperb_rad()
write(*,*) 'xnorm:',
$ xnorm
call printevent
endif
enddo
endif
if(xnorm.gt.3e-5.and.xnorm.lt.3.1e-5) then
write(*,*) ' **********************'
write(*,*) ' typical event:',
$ xnorm
call printevent
endif
endif
else
return
endif
else
return
endif
c fill the norm array
iy=abs(kn_y)*rad_nynorms+1
icsi=log(1/(1-kn_csi))/log(kn_sbeams/max(1d0,kn_minmass**2))
1 *rad_ncsinorms+1
c
if(iy.lt.1) iy=1
if(iy.gt.rad_nynorms) iy=rad_nynorms
if(icsi.lt.1) icsi=1
if(icsi.gt.rad_ncsinorms) icsi=rad_ncsinorms
c
rad_csiynorms(icsi,iy,rad_kinreg,rad_ubornidx)=
# max(xnorm,rad_csiynorms(icsi,iy,rad_kinreg,rad_ubornidx))
rad_norms(rad_kinreg,rad_ubornidx)=
# max(xnorm,rad_norms(rad_kinreg,rad_ubornidx))
if(xnorm.gt.0) then
ind=log(xnorm)+50
else
ind=-100
endif
ind=max(ind,1)
ind=min(ind,100)
efficiencies(rad_kinreg,rad_ubornidx)=
# efficiencies(rad_kinreg,rad_ubornidx)+xnorm
histnorms(rad_kinreg,rad_ubornidx,ind)=
# histnorms(rad_kinreg,rad_ubornidx,ind)+1
end
subroutine printevent
implicit none
include 'pwhg_math.h'
include 'nlegborn.h'
include 'pwhg_flst.h'
include 'pwhg_kn.h'
include 'pwhg_rad.h'
real * 8 pwhg_pt2
external pwhg_pt2
integer j,k,alr,mu
write(*,*) 'kinematic region (1=isr,>2=fsr)',rad_kinreg
write(*,*) ' born flavour:'
write(*,*)
# (flst_born(j,rad_ubornidx),j=1,nlegborn)
write(*,*) ' real flavour:'
do k=1,flst_born2alr(0,rad_ubornidx)
alr=flst_born2alr(k,rad_ubornidx)
write(*,*) (flst_alr(j,alr),j=1,nlegreal)
enddo
write(*,*) ' pt', sqrt(pwhg_pt2())
write(*,*) ' x1,x2 = ', kn_x1, kn_x2
write(*,*) ' underlying born x1,x2', kn_xb1, kn_xb2
write(*,*) ' y, csi:', kn_y,kn_csi
c write(*,*) ' born pt', sqrt(kn_born_pt2)
c write(*,*) ' cosine decay angle', kn_cthdec
write(*,*) ' born cm kinematics'
do k=1,nlegborn
write(*,100) (kn_cmpborn(mu,k),mu=0,3)
enddo
100 format('[',3(d14.8,','),d14.8,'],')
end
subroutine printhistnorms
implicit none
include 'nlegborn.h'
include 'pwhg_flst.h'
include 'pwhg_rad.h'
real * 8 ptsq,born,sig,xnorm
real * 8 histnorms(nlegborn-1,maxprocborn,100)
real * 8 efficiencies(nlegborn-1,maxprocborn)
common/chistnorms/histnorms,efficiencies
integer iun,ireg,iborn,j
real * 8 x1,x2,xm,dx
real * 8 powheginput
external powheginput
call newunit(iun)
open(iun,status='unknown',file='pwghistnorms.top')
do ireg=1,rad_nkinreg
do iborn=1,flst_nborn
write(iun,*)'set scale x log'
write(iun,*)'set order x dx y'
write(iun,*)'set limits x 1e-14 1e14'
write(iun,*)'title top "region=',ireg,', born=',iborn,'"'
do j=1,100
x1=exp(j-50d0)
x2=exp(j+1-50d0)
xm=(x1+x2)/2
dx=x2-xm
write(iun,*) xm,dx,histnorms(ireg,iborn,j)
enddo
write(iun,*)' hist'
if (rad_norms(ireg,iborn).ne.0d0) then
write(iun,*) 'title bottom " efficiency',
1 efficiencies(ireg,iborn)/powheginput('nubound')/
2 rad_norms(ireg,iborn),'"'
endif
write(iun,*)' newplot'
enddo
enddo
c do ireg=1,rad_nkinreg
c do iborn=1,flst_nborn
c write(iun,*)'( efficiencies: region=',ireg,', born=',iborn
c write(iun,*) '(', efficiencies(ireg,iborn)/
c 1 powheginput('nubound')/rad_norms(ireg,iborn)
c enddo
c enddo
close(iun)
end
subroutine zerohistnorms
implicit none
include 'nlegborn.h'
include 'pwhg_flst.h'
include 'pwhg_rad.h'
real * 8 ptsq,born,sig,xnorm
real * 8 histnorms(nlegborn-1,maxprocborn,100)
real * 8 efficiencies(nlegborn-1,maxprocborn)
common/chistnorms/histnorms,efficiencies
integer ireg,iborn,j
do ireg=1,rad_nkinreg
do iborn=1,flst_nborn
efficiencies(ireg,iborn)=0
do j=1,100
histnorms(ireg,iborn,j)=0
enddo
enddo
enddo
end
function refuse_pdf()
c Some pdf have (inconsistently) vanishing quark density at x values
c when the gluon density is not zero; cannot use these points to compute
c maxrat
implicit none
logical refuse_pdf
include 'nlegborn.h'
include 'pwhg_flst.h'
include 'pwhg_kn.h'
include 'pwhg_rad.h'
include 'pwhg_st.h'
include 'pwhg_pdf.h'
real * 8 pdf(-pdf_nparton:pdf_nparton),sn,suppfc,x
integer k,j,fl
suppfc(x)=(1-x)*x
c in FKS, for final state radiation, pdf always cancel in real/born
if(rad_kinreg.ne.1) then
refuse_pdf=.false.
return
endif
if(st_mufact2.lt.2.or.kn_x1.gt.0.9.or.kn_x2.gt.0.9) then
refuse_pdf=.true.
return
endif
do k=1,2
if(k.eq.1) then
fl=flst_born(1,rad_ubornidx)
x=kn_xb1
call pdfcall(1,x,pdf)
else
fl=flst_born(2,rad_ubornidx)
x=kn_xb2
call pdfcall(2,x,pdf)
endif
sn=0
do j=1,5
sn=sn+pdf(j)+pdf(-j)
enddo
if(fl.eq.0) then
if(sn*suppfc(x).gt.30*pdf(0)) then
refuse_pdf=.true.
return
endif
else
if(pdf(0)*suppfc(x).gt.30*pdf(fl)) then
refuse_pdf=.true.
return
endif
endif
enddo
refuse_pdf=.false.
end
subroutine uboundfct(fct,csi,y)
implicit none
include 'pwhg_math.h'
include 'nlegborn.h'
include 'pwhg_flst.h'
include 'pwhg_kn.h'
include 'pwhg_rad.h'
real * 8 fct,csi,y
real * 8 unorm
integer iy,icsi
iy=abs(y)*rad_nynorms+1
icsi=log(1/(1-csi))/log(kn_sbeams/kn_minmass**2)*rad_ncsinorms+1
c
if(iy.lt.1) iy=1
if(iy.gt.rad_nynorms) iy=rad_nynorms
if(icsi.lt.1) icsi=1
if(icsi.gt.rad_ncsinorms) icsi=rad_ncsinorms
c
unorm=rad_norms(rad_kinreg,rad_ubornidx)
if(unorm.ne.0) then
fct=rad_csiynorms(icsi,iy,rad_kinreg,rad_ubornidx)/unorm
else
fct=1
endif
end
subroutine gen_born(mcalls,icalls)
implicit none
include 'nlegborn.h'
include 'pwhg_flst.h'
include 'pwhg_flg.h'
integer mcalls,icalls
integer ndimbornint
parameter (ndimbornint=ndiminteg-3)
real * 8 xgrid(0:50,ndimbornint),ymax(50,ndimbornint),
1 ymaxrat(50,ndimbornint),xacc(0:50,ndimbornint),
1 xmmm(0:50,ndimbornint),xint,sigbtl,errbtl
integer ifold(ndimbornint),nhits(0:50,ndimbornint),ncall1,ncall2,
1 itmx1,itmx2,imode,j,k,iun
character * 20 pwgprefix
integer lprefix
common/cpwgprefix/pwgprefix,lprefix
logical savenlotest,save_storemintupb,save_fastbtlbound,ini,
1 save_monitorubound
data ini/.true./
real * 8 xx(ndimbornint)
real * 8 btilde,powheginput
external btilde,powheginput
save xgrid,ymax,xmmm,ini,ifold
if(flg_bornonly) then
call gen_btilde(mcalls,icalls)
return
endif
flg_bornonly=.true.
save_storemintupb = flg_storemintupb
save_fastbtlbound = flg_fastbtlbound
save_monitorubound = flg_monitorubound
flg_storemintupb = .false.
flg_fastbtlbound = .false.
flg_monitorubound = .false.
c
savenlotest=flg_nlotest
flg_nlotest=.false.
if(ini) then
write(*,'(/,a,/,a)')
1 ' Initializing grids for the generation of the underlying',
1 ' Born configurations according to the Born cross section'
c set up the grids
do j=1,ndimbornint
do k=0,50
xgrid(k,j)=0
xmmm(k,j)=0
enddo
do k=1,50
ymax(k,j)=0
enddo
ifold(j)=1
enddo
call initxgrid(xgrid,ndimbornint)
call newunit(iun)
open(unit=iun,file=pwgprefix(1:lprefix)//'borngrid.top',
1 status='unknown')
imode=0
ncall1=powheginput('ncall1')
itmx1=powheginput('itmx1')
c This is not really needed; the accumulated totals
c will not be used;
call resettotals
call mint(btilde,ndimbornint,ncall1,itmx1,ifold,imode,iun,
1 xgrid,xint,xacc,nhits,ymax,ymaxrat,sigbtl,errbtl)
close(iun)
imode=1
ncall2=powheginput('ncall2')
itmx2=powheginput('itmx2')
call mint(btilde,ndimbornint,ncall2,itmx2,ifold,imode,iun,
1 xgrid,xint,xacc,nhits,ymax,ymaxrat,sigbtl,errbtl)
call gen(btilde,ndimbornint,xgrid,ymax,ymaxrat,xmmm,ifold,0,
1 mcalls,icalls,xx)
ini=.false.
endif
C - N.B. if -fno-automatic isn't used the ifold array may need reset to 1 here.
call gen(btilde,ndimbornint,xgrid,ymax,ymaxrat,xmmm,ifold,1,
# mcalls,icalls,xx)
flg_nlotest=savenlotest
flg_bornonly=.false.
flg_storemintupb = save_storemintupb
flg_fastbtlbound = save_fastbtlbound
flg_monitorubound = save_monitorubound
end