1
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2
c MINT Integrator Package
3
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4
c Original version by Paolo Nason (for POWHEG (BOX))
5
c Modified by Rikkert Frederix (for aMC@NLO)
6
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
7
c subroutine mint(fun,ndim,ncalls0,nitmax,imode,
8
c ndim=number of dimensions
9
c ncalls0=# of calls per iteration
10
c nitmax =# of iterations
11
c fun(xx,www,ifirst): returns the function to be integrated multiplied by www;
12
c xx(1:ndim) are the variables of integration
13
c ifirst=0: normal behaviour
17
c same as imode=0 as far as this routine is concerned, except for the
18
c fact that a grid is read at the beginning (rather than initialized).
19
c The return value of imode will be zero.
22
c When called with imode=0 the routine integrates the absolute value of
23
c the function and sets up a grid xgrid(0:50,ndim) such that in each
24
c ndim-1 dimensional slice (i.e. xgrid(m-1,n)<xx(n)<xgrid(m,n)) the
25
c contribution of the integral is the same the array xgrid is setup at
26
c this stage; ans and err are the integral and its error
28
c imode=1 (in fact #0)
29
c When called with imode=1, the routine performs the integral of the
30
c function fun using the grid xgrid. If some number in the array ifold,
31
c (say, ifold(n)) is different from 1, it must be a divisor of 50, and
32
c the 50 intervals xgrid(0:50,n) are grouped into ifold(n) groups, each
33
c group containing 50/ifold(n) nearby intervals. For example, if
34
c ifold(1)=5, the 50 intervals for the first dimension are divided in 5
35
c groups of 10. The integral is then performed by folding on top of each
36
c other these 5 groups. Suppose, for example, that we choose a random
37
c point in xx(1) = xgrid(2,1)+x*(xgrid(3,1)-xgrid(2,1)), in the group of
38
c the first 5 interval. we sum the contribution of this point to the
39
c contributions of points
40
c xgrid(2+m*10,1)+x*(xgrid(3+m*10,1)-xgrid(2+m*10,1)), with m=1,...,4.
41
c In the sequence of calls to the function fun, the call for the first
42
c point is performed with ifirst=0, and that for all subsequent points
43
c with ifirst=1, so that the function can avoid to compute quantities
44
c that only depend upon dimensions that have ifold=1, and do not change
45
c in each group of folded call. The values returned by fun in a sequence
46
c of folded calls with ifirst=0 and ifirst=1 are not used. The function
47
c itself must accumulate the values, and must return them when called
50
c Added the posibility to keep track of more than one integral:
52
c nintegrals=1 : the function that is used to update the grids. This is
53
c the ABS cross section. If imode.eq.1, this does not contain the
54
c virtual corrections because for them a separate maximum is kept using (5).
55
c nintegrals=2 : the actual cross section. This includes virtual corrections.
56
c nintegrals=3 : the cross section from the M_Virt/M_Born ratio alone:
57
c this defines the average virtual that is added to each phase-space
59
c nintegrals=4 : the cross section of the actual virtual minus the
60
c average virtual. This is used to determine the fraction of phase-space
61
c points for which we include the virtual.
62
c nintegrals=5 : abs of 3
63
c nintegrals=6 : born*alpha_S/2Pi
65
subroutine mint(fun,ndim,ncalls0,nitmax,imode,xgrid,ymax,ymax_virt
66
$ ,ans,unc,chi2,nhits_in_grids)
67
c imode= 0: integrate and adapt the grid
68
c imode= 1: frozen grid, compute the integral and the upper bounds
69
c imode=-1: same as imode=0, but use previously generated grids
70
c others: same as 1 (for now)
73
include "FKSParams.inc"
74
integer i,j,ncalls0,ndim,nitmax,imode
75
real * 8 fun,xgrid(0:nintervals,ndimmax,maxchannels)
76
$ ,xint(maxchannels),ymax(nintervals,ndimmax,maxchannels)
77
$ ,ans(nintegrals,0:maxchannels),unc(nintegrals,0:maxchannels)
78
$ ,ans3(nintegrals,3),unc3(nintegrals,3),ans_l3(nintegrals)
79
$ ,unc_l3(nintegrals),chi2_l3(nintegrals),vol_chan
80
real * 8 xint_virt(maxchannels),ymax_virt(0:maxchannels)
81
$ ,ans_chan(0:maxchannels)
82
real * 8 x(ndimmax),vol
83
real * 8 xacc(0:nintervals,ndimmax,maxchannels)
84
integer icell(ndimmax),ncell(ndimmax),ncell_virt
85
integer ifold(ndimmax),kfold(ndimmax)
87
integer nhits(nintervals,ndimmax,maxchannels)
88
$ ,nhits_in_grids(maxchannels)
89
real * 8 rand(ndimmax),HwU_values(2)
90
real * 8 dx(ndimmax),f(nintegrals),vtot(nintegrals,0:maxchannels)
91
$ ,etot(nintegrals,0:maxchannels),prod,f1(nintegrals)
92
$ ,chi2(nintegrals,0:maxchannels),efrac(nintegrals),dummy
93
integer kdim,kint,kpoint,nit,ncalls,iret,nintcurr,nintcurr_virt
94
$ ,ifirst,nit_included,kpoint_iter,non_zero_point(nintegrals)
95
$ ,ntotcalls(nintegrals),nint_used,nint_used_virt,min_it,np
98
logical even,double_events,bad_iteration,regridded(maxchannels)
100
double precision average_virtual(maxchannels)
101
$ ,virtual_fraction(maxchannels)
102
common/c_avg_virt/average_virtual,virtual_fraction
103
character*13 title(nintegrals)
105
common /c_new_point/ new_point
106
data title(1)/'ABS integral '/
107
data title(2)/'Integral '/
108
data title(3)/'Virtual '/
109
data title(4)/'Virtual ratio'/
110
data title(5)/'ABS virtual '/
111
data title(6)/'Born*ao2pi '/
114
common /for_applgrid/ iappl
115
logical fixed_order,nlo_ps
116
common /c_fnlo_nlops/fixed_order,nlo_ps
118
double precision cross_section
119
common /for_FixedOrder_lhe/ cross_section,npoints
122
c if ncalls0 is greater than 0, use the default running, i.e. do not
123
c double the events after each iteration as well as use a fixed number
124
c of intervals in the grids.
125
if (ncalls0.gt.0) then
126
double_events=.false.
128
nint_used_virt=nintervals_virt
130
c if ncalls0.le.0, reset it and double the events per iteration
131
ncalls0=80*ndim*(nchans/3+1)
133
if (imode.eq.1 .or. imode.eq.-1) then
135
nint_used_virt=nintervals_virt
138
nint_used_virt=min_inter
141
bad_iteration=.false.
143
ncalls=0 ! # PS points (updated below)
145
c Grids read from file
154
ans_chan(kchan)=ans(1,kchan)
155
ans_chan(0)=ans_chan(0)+ans(1,kchan)
157
elseif(imode.eq.0) then
165
xgrid(kint,kdim,kchan)=dble(kint)/nint_used
167
nhits(kint,kdim,kchan)=0
170
regridded(kchan)=.true.
173
do kchan=1,maxchannels
174
nhits_in_grids(kchan)=0
176
call init_ave_virt(nint_used_virt,ndim)
180
if (double_events) then
181
c when double events, start with the very first channel
182
c only. For the first iteration, we compute each channel
186
ncalls0=ncalls0/nchans
188
elseif(imode.eq.1) then
189
c Initialize upper bounding envelope
191
xint(kchan)=ans(1,kchan)
192
xint_virt(kchan)=ans(5,kchan)
197
nintcurr=nint_used/ifold(kdim)
198
nintcurr_virt=nint_used_virt/ifold(kdim)
199
if(nintcurr*ifold(kdim).ne.nint_used .or.
200
& nintcurr_virt*ifold(kdim).ne.nint_used_virt) then
201
write(*,*) 'mint: the values in the ifold array'/
202
$ /'shoud be divisors of',nint_used,'and'
208
ymax(kint,kdim,kchan)=xint(kchan)**(1d0/ndim)
213
ymax_virt(kchan)=xint_virt(kchan)
217
ans_chan(kchan)=ans(1,kchan)
218
ans_chan(0)=ans_chan(0)+ans(1,kchan)
221
cross_section=ans_chan(0)
237
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
238
c Main loop over the iterations
243
np=np+nhits(kint,1,kchan)
245
write (*,250) 'channel',kchan,':',iconfigs(kchan)
246
$ ,regridded(kchan),np,nhits_in_grids(kchan),ans_chan(kchan)
247
$ ,ans(2,kchan),virtual_fraction(kchan)
250
if(nit.ge.nitmax) then
251
c We did enough iterations, update arguments and return
253
if (nit_included.ge.2) then
254
chi2(1,kchan)=chi2(1,kchan)/dble(nit_included-1)
259
write (*,*) '-------'
260
ncalls0=ncalls*kpoint_iter ! return number of points used
261
if (double_events) then
266
cross_section=ans(2,0)
268
if (regridded(kchan)) then
271
np=np+nhits(kint,1,kchan)
273
nhits_in_grids(kchan)=np ! set equal to number of points
274
! used for last update
280
write (*,*) '------- iteration',nit
281
if (even .and. ncalls.ne.ncalls0) then
282
c Uses more evenly distributed random numbers. This overwrites the
284
call initialize_even_random_numbers(ncalls0,ndim,ncalls)
285
write (*,*) 'Update # PS points (even): ',ncalls0,' --> '
287
elseif (ncalls0.ne.ncalls) then
289
write (*,*) 'Update # PS points: ',ncalls0,' --> ',ncalls
292
c Reset the accumulated results for grid updating
295
if (regridded(kchan).or.reset) then ! only reset if grids were
296
! updated (or forced reset)
297
if (regridded(kchan) .and. .not. reset) then
300
np=np+nhits(kint,1,kchan)
302
nhits_in_grids(kchan)=np ! set equal to number of
303
! points used for last update
304
elseif (regridded(kchan) .and. reset) then
305
nhits_in_grids(kchan)=0
309
xacc(kint,kdim,kchan)=0
311
nhits(kint,kdim,kchan)=0
329
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
330
c Loop over PS points
331
2 kpoint_iter=kpoint_iter+1
334
c first determine which integration channel we want to pick this point from
335
call get_channel(ans_chan,vol_chan)
336
c find random x, and its random cell
339
c if(even), we should compute the ncell and the rand from the ran3()
341
rand(kdim)=ran3(even)
342
ncell(kdim)= min(int(rand(kdim)*nint_used)+1,
344
rand(kdim)=rand(kdim)*nint_used-(ncell(kdim)-1)
346
ncell(kdim)=min(int(nint_used/ifold(kdim)*ran3(even))+1,
348
rand(kdim)=ran3(even)
356
vol=1d0/vol_chan * wgt_mult
357
c c convert 'flat x' ('rand') to 'vegas x' ('x') and include jacobian ('vol')
359
nintcurr=nint_used/ifold(kdim)
360
icell(kdim)=ncell(kdim)+(kfold(kdim)-1)*nintcurr
361
dx(kdim)=xgrid(icell(kdim),kdim,ichan)-xgrid(icell(kdim)-1
363
vol=vol*dx(kdim)*nintcurr
364
x(kdim)=xgrid(icell(kdim)-1,kdim,ichan)+rand(kdim)*dx(kdim)
365
if(imode.eq.0) nhits(icell(kdim),kdim,ichan)=
366
$ nhits(icell(kdim),kdim,ichan)+1
368
call get_ave_virt(x,nint_used_virt,ndim,average_virtual)
369
c contribution to integral
371
dummy=fun(x,vol,ifirst,f1)
376
c this accumulated value will not be used
377
dummy=fun(x,vol,ifirst,f1)
382
call nextlexi(ndim,ifold,kfold,iret)
384
c closing call: accumulated value with correct sign
385
dummy=fun(x,vol,2,f1)
392
c accumulate the function in xacc(icell(kdim),kdim) to adjust the grid later
394
xacc(icell(kdim),kdim,ichan)=
395
$ xacc(icell(kdim),kdim,ichan) + f(1)
397
c Set the Born contribution (to compute the average_virtual) to zero if
398
c the virtual was not computed for this phase-space point. Compensate by
399
c including the virtual_fraction.
400
if (f(3).ne.0d0) then
401
f(6)=f(6)/virtual_fraction(ichan)
402
call fill_ave_virt(x,nint_used_virt,ndim,f(3),f(6))
407
c update the upper bounding envelope total rate
410
prod=prod*ymax(ncell(kdim),kdim,ichan)
413
if (prod.gt.1d0) then
414
c Weight for this PS point is larger than current upper bound. Increase
415
c the bound so that it is equal to the current max weight. If the new
416
c point is more than twice as large as current upper bound, increase
417
c bound by factor 2 only to prevent a single unstable points to
418
c completely screw up the efficiency
420
prod=prod**(1d0/dble(ndim))
422
ymax(ncell(kdim),kdim,ichan)=
423
$ ymax(ncell(kdim),kdim,ichan)*prod
426
c Update the upper bounding envelope virtual. Do not include the
427
c enhancement due to the virtual_fraction. (And again limit by factor 2
429
if (f(5)*virtual_fraction(ichan).gt.ymax_virt(ichan))
430
$ ymax_virt(ichan)=min(f(5)*virtual_fraction(ichan)
431
$ ,ymax_virt(ichan)*2d0)
432
c for consistent printing in the log files (in particular when doing LO
433
c runs), set also f(6) to zero when imode.eq.1 and the virtuals are not
435
if (f(3).eq.0) f(6)=0d0
438
if (f(i).ne.0d0) non_zero_point(i)=non_zero_point(i)+1
440
c Add the PS point to the result of this iteration
442
vtot(i,ichan)=vtot(i,ichan)+f(i)
443
etot(i,ichan)=etot(i,ichan)+f(i)**2
445
if (f(1).ne.0d0) call HwU_add_points
448
c Number of phase-space points used
449
ntotcalls(i)=ncalls*kpoint_iter
450
c Special for the computation of the 'computed virtual'
451
if (i.eq.4 .and. non_zero_point(i).ne.0 )
452
& ntotcalls(i) = non_zero_point(i)
454
if (ntotcalls(1).gt.max_points .and. non_zero_point(1).lt.25 .and.
455
& double_events) then
456
write (*,*) 'ERROR: INTEGRAL APPEARS TO BE ZERO.'
457
write (*,*) 'TRIED',ntotcalls(1),'PS POINTS AND ONLY '
458
& ,non_zero_point(1),' GAVE A NON-ZERO INTEGRAND.'
461
c Goto beginning of loop over PS points until enough points have found
463
if (non_zero_point(1).lt.int(0.99*ncalls)
464
& .and. double_events) goto 2
465
c This is the loop over all the channels for the very first iteration.
466
if (imode.eq.0 .and. nit.eq.1 .and. double_events) then
468
if (ans_chan(kchan).eq.1d0) then
470
vtot(i,kchan)=vtot(i,kchan)/dble(ntotcalls(i))
471
etot(i,kchan)=etot(i,kchan)/dble(ntotcalls(i))
472
etot(i,kchan)=sqrt(abs(etot(i,kchan)-vtot(i,kchan)**2)
473
$ /dble(ntotcalls(i)))
475
if (kchan.eq.nchans) exit
477
ans_chan(kchan+1)=1d0
486
c set the total result for the first iteration as the sum over all the channels
489
vtot(i,0)=vtot(i,0)+vtot(i,kchan)
490
etot(i,0)=etot(i,0)+etot(i,kchan)**2
492
etot(i,0)=sqrt(etot(i,0))
494
ncalls0=ncalls0*nchans
496
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
497
c Iteration done. Update the accumulated results and print them to the
500
if (imode.eq.0 .and. nit.eq.1 .and. double_events) then
505
vtot(i,0)=vtot(i,0)+vtot(i,kchan)
506
etot(i,0)=etot(i,0)+etot(i,kchan)
508
vtot(i,kchan)=vtot(i,kchan)/dble(ntotcalls(i))
509
etot(i,kchan)=etot(i,kchan)/dble(ntotcalls(i))
510
c the abs is to avoid tiny negative values
511
etot(i,kchan)=sqrt(abs(etot(i,kchan)-vtot(i,kchan)**2)
512
$ /dble(ntotcalls(i)))
515
if (vtot(i,0).ne.0d0) then
516
efrac(i)=abs(etot(i,0)/vtot(i,0))
522
write(*,'(a,1x,e10.4,1x,a,1x,e10.4,1x,a,1x,f7.3,1x,a)')
523
$ title(i)//' =',vtot(i,0),' +/- ',etot(i,0),' (',efrac(i)
526
C If there was a large fluctation in this iteration, be careful with
527
C including it in the accumalated results and plots.
528
if (efrac(1).gt.0.3d0 .and. iappl.eq.0 .and. nit.gt.3) then
529
c Do not include the results in the plots
530
if (fixed_order) call accum(.false.)
531
if (fixed_order) call HwU_accum_iter(.false.,ntotcalls(1)
534
if (efrac(1).gt.0.3d0 .and. nit.gt.3 .and. iappl.eq.0) then
535
c Do not include the results in the updating of the grids.
536
write (*,*) 'Large fluctuation ( >30 % ).'
537
& //'Not including iteration in results.'
538
c empty the accumulated results in the MC over integers
539
call empty_MC_integer
540
c empty the accumalated results for the MINT grids
542
c emptying accum. results is done above when the iteration starts
544
elseif (imode.eq.1) then
545
c Cannot really skip the increase of the upper bounding envelope. So,
546
c simply continue here. Note that no matter how large the integrand for
547
c the PS point, the upper bounding envelope is at most increased by a
548
c factor 2, so this should be fine.
551
c double the number of points for the next iteration
552
if (double_events) ncalls0=ncalls0*2
553
if (bad_iteration .and. imode.eq.0 .and. double_events) then
554
c 2nd bad iteration is a row. Reset grids
555
write (*,*)'2nd bad iteration in a row. '/
556
& /'Resetting grids and starting from scratch...'
557
if (double_events) then
558
if (imode.eq.0) nint_used=min_inter ! reset number of intervals
559
ncalls0=ncalls0/8 ! Start with larger number
563
c Reset the MINT grids
568
xgrid(kint,kdim,kchan)=dble(kint)/nint_used
570
regridded(kchan)=.true.
573
do kchan=1,maxchannels
574
nhits_in_grids(kchan)=0
576
call init_ave_virt(nint_used_virt,ndim)
580
if (double_events) then
583
ncalls0=ncalls0/nchans
585
elseif (imode.eq.1) then
587
nintcurr=nint_used/ifold(kdim)
590
ymax(kint,kdim,kchan)=xint(kchan)**(1d0/ndim)
593
nintcurr_virt=nint_used_virt/ifold(kdim)
596
ymax_virt(kchan)=xint_virt(kchan)
599
call reset_MC_grid ! reset the grid for the integers
600
if (fixed_order) call initplot ! Also reset all the plots
615
bad_iteration=.false.
621
bad_iteration=.false.
623
HwU_values(1)=etot(1,0)
624
HwU_values(2)=unc(1,0)
628
ans(i,kchan)=vtot(i,kchan)
629
unc(i,kchan)=etot(i,kchan)
631
ans_chan(kchan)=ans(1,kchan)
633
write (*,'(a,1x,e10.4)') 'Chi^2 per d.o.f.',0d0
635
c prevent annoying division by zero for nearly zero
637
do kchan=nchans,0,-1 ! go backwards so that kchan=0 goes last
638
! (this makes sure central value is
639
! correctly updated).
641
if(etot(i,0).eq.0.and.unc(i,0).eq.0) then
642
if(ans(i,0).eq.vtot(i,0) .and. i.eq.1) then
643
c double the number of points for the next iteration
644
if (double_events) ncalls0=ncalls0*2
647
unc(i,kchan)=abs(vtot(i,kchan)-ans(i,kchan))
648
etot(i,kchan)=abs(vtot(i,kchan)-ans(i,kchan))
650
elseif(etot(i,0).eq.0) then
651
etot(i,kchan)=unc(i,kchan)
652
elseif(unc(i,0).eq.0) then ! 1st iteration; set to a large value
653
unc(i,kchan)=etot(i,kchan)*1d99
655
if (i.ne.1 .and. (etot(i,0).eq.0 .or. unc(i,0).eq.0)) then
660
ans(i,kchan)=(ans(i,kchan)/unc(i,0)+vtot(i,kchan)/etot(i
661
$ ,0))/(1/unc(i,0)+1/etot(i,0))
662
unc(i,kchan)=1/sqrt(1/unc(i,kchan)**2+1/etot(i,kchan)**2)
663
chi2(i,kchan)=chi2(i,kchan)+(vtot(i,kchan)-ans(i,kchan))
664
$ **2/etot(i,kchan)**2
667
ans_chan(kchan)=ans(1,kchan)
669
write (*,'(a,1x,e10.4)') 'Chi^2=',(vtot(1,0)-ans(1,0))**2
672
nit_included=nit_included+1
674
if (ans(i,0).ne.0d0) then
675
efrac(i)=abs(unc(i,0)/ans(i,0))
679
write(*,'(a,1x,e10.4,1x,a,1x,e10.4,1x,a,1x,f7.3,1x,a)')
680
$ 'accumulated results '//title(i)//' =',ans(i,0),' +/- '
681
$ ,unc(i,0) ,' (',efrac(i)*100d0,'%)'
683
cross_section=ans(1,0)
684
if (nit_included.le.1) then
685
write (*,'(a,1x,e10.4)') 'accumulated result Chi^2 per DoF ='
688
write (*,'(a,1x,e10.4)') 'accumulated result Chi^2 per DoF =',
689
& chi2(1,0)/dble(nit_included-1)
692
c Update the fraction of the events for which we include the virtual corrections
695
virtual_fraction(kchan)=max(min(virtual_fraction(kchan)
696
$ *max(min(2d0*etot(3,kchan)/etot(1,kchan),2d0),0.25d0)
697
$ ,1d0),Min_virt_fraction)
700
c Update the results of the last tree iterations
703
ans3(i,j)=ans3(i,j+1)
704
unc3(i,j)=unc3(i,j+1)
709
c Compute the results of the last three iterations
710
if (nit_included.ge.4) then
713
unc_l3(i)=ans3(i,1)*1d99
716
if (i.ne.1 .and. (unc_l3(i).eq.0d0 .or. unc3(i
722
ans_l3(i)=(ans_l3(i)/unc_l3(i)+ans3(i,j)/unc3(i,j))/
723
& (1/unc_l3(i)+1/unc3(i,j))
724
unc_l3(i)=1/sqrt(1/unc_l3(i)**2+1/unc3(i,j)**2)
725
chi2_l3(i)=chi2_l3(i)+
726
& (ans3(i,j)-ans_l3(i))**2/unc3(i,j)**2
729
chi2_l3(i)=chi2_l3(i)/2d0 ! three iterations, so 2 degrees of freedom
732
if (ans_l3(i).ne.0d0) then
733
efrac(i)=abs(unc_l3(i)/ans_l3(i))
737
write(*,'(a,1x,e10.4,1x,a,1x,e10.4,1x,a,1x,f7.3,1x,a)')
738
$ 'accumulated results last 3 iterations '//title(i)/
739
$ /' =' ,ans_l3(i),' +/- ',unc_l3(i) ,' (',efrac(i)*100d0
742
write(*,'(a,1x,e10.4)')
743
$ 'accumulated result last 3 iterrations Chi^2 per DoF ='
747
c Iteration is finished; now rearrange the grid
750
call regrid(xacc(0,kdim,kchan),xgrid(0,kdim,kchan)
751
$ ,nhits(1,kdim,kchan),nint_used,nhits_in_grids(kchan)
755
call regrid_ave_virt(nint_used_virt,ndim)
756
c Regrid the MC over integers (used for the MC over FKS dirs)
757
call regrid_MC_integer
759
c Quit if the desired accuracy has been reached
760
if (nit_included.ge.min_it .and. accuracy.gt.0d0) then
761
if (unc(1,0)/ans(1,0)*max(1d0,chi2(1,0)/dble(nit_included-1))
763
write (*,*) 'Found desired accuracy'
765
c Improve the stats in the plots
766
if (fixed_order) call accum(.true.)
767
if (fixed_order) call HwU_accum_iter(.true.,ntotcalls(1)
770
elseif(unc_l3(1)/ans_l3(1)*max(1d0,chi2_l3(1)).lt.accuracy)
773
& 'Found desired accuracy in last 3 iterations'
778
chi2(i,0)=chi2_l3(i)*dble(nit_included-1)
780
c Improve the stats in the plots
781
if (fixed_order) call accum(.true.)
782
if (fixed_order) call HwU_accum_iter(.true.,ntotcalls(1)
787
c Double the number of intervals in the grids if not yet reach the maximum
788
if (2*nint_used.le.nintervals .and. double_events) then
791
call double_grid(xgrid(0,kdim,kchan),nhits(1,kdim,kchan)
795
nint_used=2*nint_used
797
if (2*nint_used_virt.le.nintervals_virt .and. double_events) then
798
call double_ave_virt(nint_used_virt,ndim)
799
nint_used_virt=2*nint_used_virt
802
c double the number of points for the next iteration
803
if (double_events) ncalls0=ncalls0*2
804
c Also improve stats in plots
805
if (fixed_order) call accum(.true.)
806
if (fixed_order) call HwU_accum_iter(.true.,ntotcalls(1)
811
250 format(a7,i5,1x,a1,1x,i5,1x,l,1x,i8,1x,i8,2x,e10.4,2x,e10.4,2x
815
subroutine double_grid(xgrid,nhits,ninter)
818
integer ninter,nhits(nintervals)
819
real * 8 xgrid(0:nintervals)
823
xgrid(i*2-1)=(xgrid(i)+xgrid(i-1))/2d0
824
nhits(i*2)=nhits(i)/2
825
nhits(i*2-1)=nhits(i)-nhits(i*2)
831
subroutine regrid(xacc,xgrid,nhits,ninter,nhits_in_grid,regridded)
835
integer ninter,nhits(nintervals),np,nhits_in_grid
836
real * 8 xacc(0:nintervals),xgrid(0:nintervals)
837
real * 8 xn(nintervals),r,tiny,xl,xu,nl,nu,sum
838
parameter ( tiny=1d-8 )
840
c compute total number of points and update grids if large
846
if (np.lt.nhits_in_grid) return
849
c Use the same smoothing as in VEGAS uses for the grids, i.e. use the
850
c average of the central and the two neighbouring grid points: (Only do
851
c this if we are already at the maximum intervals, because the doubling
852
c of the grids also includes a smoothing).
853
if (ninter.eq.nintervals) then
859
nhits(1)=nint((nl+nu)/2d0)
864
xacc(kint)=(xacc(kint)+xu)/3d0
868
nhits(kint)=nint((nhits(kint)+nu)/3d0)
870
xacc(ninter)=(xu+xl)/2d0
871
nhits(ninter)=nint((nu+nl)/2d0)
876
if (nhits(kint).ne.0) then
877
xacc(kint)=abs(xacc(kint))/nhits(kint)
883
c xacc (xerr) already contains a factor equal to the interval size
884
c Thus the integral of rho is performed by summing up
885
if(nhits(kint).ne.0) then
886
c take logarithm to help convergence (taken from LO dsample.f)
887
if (xacc(kint).ne.sum) then
888
xacc(kint)=((xacc(kint)/sum-1d0)/
889
& log(xacc(kint)/sum))**1.5
893
xacc(kint)= xacc(kint-1) + abs(xacc(kint))
895
xacc(kint)=xacc(kint-1)
899
xacc(kint)=xacc(kint)/xacc(ninter)
901
c Check that we have a reasonable result and update the accumulated
904
if (xacc(kint).lt.(xacc(kint-1)+tiny)) then
905
xacc(kint)=xacc(kint-1)+tiny
908
c it could happen that the change above yielded xacc() values greater
909
c than 1: one more update needed
912
if (xacc(ninter-kint).gt.(xacc(ninter-kint+1)-tiny)) then
913
xacc(ninter-kint)=1d0-dble(kint)*tiny
920
r=dble(kint)/dble(ninter)
922
if(r.lt.xacc(jint)) then
923
xn(kint)=xgrid(jint-1)+(r-xacc(jint-1))
924
# /(xacc(jint)-xacc(jint-1))*(xgrid(jint)-xgrid(jint-1))
928
if(jint.ne.ninter+1.and.kint.ne.ninter) then
929
write(*,*) ' error',jint,ninter
940
subroutine nextlexi(ndim,iii,kkk,iret)
941
c kkk: array of integers 1 <= kkk(j) <= iii(j), j=1,ndim
942
c at each call iii is increased lexicographycally.
943
c for example, starting from ndim=3, kkk=(1,1,1), iii=(2,3,2)
944
c subsequent calls to nextlexi return
945
c kkk(1) kkk(2) kkk(3) iret
960
integer ndim,iret,kkk(ndim),iii(ndim)
964
if(kkk(k).lt.iii(k)) then
980
subroutine gen(fun,ndim,xgrid,ymax,ymax_virt,imode,x,vn)
981
c imode=0 to initialize
982
c imode=1 to generate
983
c imode=3 store generation efficiency in x(1)
987
real * 8 fun,xgrid(0:nintervals,ndimmax,maxchannels)
988
$ ,ymax(nintervals,ndimmax,maxchannels)
989
$ ,ymax_virt(0:maxchannels),x(ndimmax)
990
real * 8 dx(ndimmax),xx(ndimmax),vol_chan,dummy
991
integer icell(ndimmax),ncell(ndimmax),ncell_virt
992
integer ifold(ndimmax),kfold(ndimmax)
994
real * 8 r,f(nintegrals),f1(nintegrals),ubound,vol,ran3
995
$ ,xmmm(nintervals,ndimmax,maxchannels)
996
real * 8 rand(ndimmax)
998
integer icalls,mcalls,kdim,kint,nintcurr,nintcurr_virt,iret,ifirst
999
$ ,i,vn,icalls_virt,mcalls_virt,icalls_nz,icalls_virt_nz
1000
double precision average_virtual(maxchannels)
1001
$ ,virtual_fraction(maxchannels)
1002
common/c_avg_virt/average_virtual,virtual_fraction
1003
save icalls,mcalls,icalls_virt,mcalls_virt,xmmm,icalls_nz
1006
common /c_new_point/ new_point
1010
nintcurr=nintervals/ifold(kdim)
1011
xmmm(1,kdim,kchan)=ymax(1,kdim,kchan)
1013
xmmm(kint,kdim,kchan)=xmmm(kint-1,kdim,kchan)
1014
$ +ymax(kint,kdim,kchan)
1017
xmmm(kint,kdim,kchan)=xmmm(kint,kdim,kchan)
1018
$ /xmmm(nintcurr,kdim,kchan)
1033
elseif(imode.eq.3) then
1034
if(icalls.gt.0) then
1035
x(1)=dble(mcalls)/dble(icalls)
1036
x(2)=xx(2)/dble(icalls)
1037
x(3)=xx(3)/dble(icalls)
1043
if(icalls_virt.gt.0) then
1044
x(4)=dble(mcalls_virt)/dble(icalls_virt)
1045
x(5)=xx(5)/dble(icalls_virt)
1046
x(6)=xx(6)/dble(icalls_virt)
1052
call increasecnt(' ',imode)
1056
mcalls_virt=mcalls_virt+1
1057
elseif(vn.eq.2 .or. vn.eq.3) then
1060
write (*,*) 'vn not correct in mint-integrator2.f',vn,imode
1064
call get_channel(ymax_virt,vol_chan)
1067
icalls_virt=icalls_virt+1
1068
elseif(vn.eq.2 .or. vn.eq.3) then
1074
ncell(kdim)=min(int(ran3(.false.)*nintcurr)+1,nintcurr)
1075
rand(kdim)=ran3(.false.)
1077
elseif(vn.eq.2 .or. vn.eq.3) then
1079
nintcurr=nintervals/ifold(kdim)
1082
if(r.lt.xmmm(kint,kdim,ichan)) then
1087
rand(kdim)=ran3(.false.)
1090
if (vn.eq.2 .or. vn.eq.3) then
1093
ubound=ubound*ymax(ncell(kdim),kdim,ichan)
1106
nintcurr=nintervals/ifold(kdim)
1107
nintcurr_virt=nintervals_virt/ifold(kdim)
1108
icell(kdim)=ncell(kdim)+(kfold(kdim)-1)*nintcurr
1109
dx(kdim)=xgrid(icell(kdim),kdim,ichan)-xgrid(icell(kdim)-1,kdim
1111
vol=vol*dx(kdim)*nintervals/ifold(kdim)
1112
x(kdim)=xgrid(icell(kdim)-1,kdim,ichan)+rand(kdim)*dx(kdim)
1114
call get_ave_virt(x,nintcurr_virt,ndim,average_virtual)
1116
ubound=ymax_virt(ichan)
1118
dummy=fun(x,vol,ifirst,f1)
1123
call nextlexi(ndim,ifold,kfold,iret)
1124
if(iret.eq.0) goto 5
1125
c get final value (x and vol not used in this call)
1126
dummy=fun(x,vol,2,f1)
1130
if (vn.eq.2 .or. vn.eq.3) then
1137
call increasecnt('another call to the function',imode)
1138
if (f(1).eq.0d0) then
1139
call increasecnt('failed generation cuts',imode)
1142
icalls_virt_nz=icalls_virt_nz+1
1143
elseif(vn.eq.2 .or.vn.eq.3) then
1144
icalls_nz=icalls_nz+1
1148
write(*,*) 'gen: non positive function'
1151
if(f(1).gt.ubound) then
1153
call increasecnt('ubound fail novi',imode)
1154
elseif (vn.eq.1) then
1155
call increasecnt('ubound fail virt',imode)
1156
elseif (vn.eq.3) then
1157
call increasecnt('ubound fail born',imode)
1160
ubound=ubound*ran3(.false.)
1161
if(ubound.gt.f(1)) then
1163
& ('vetoed calls in inclusive cross section',imode)
1167
call increasecnt('events gen novi',imode)
1168
elseif (vn.eq.1) then
1169
call increasecnt('events gen virt',imode)
1170
elseif (vn.eq.3) then
1171
call increasecnt('events gen born',imode)
1176
c Dummy subroutine (normally used with vegas when resuming plots)
1181
subroutine increasecnt(argument,imode)
1182
c Be careful, argument should be at least 15 characters
1183
c long for this subroutine to work properly
1185
character*(*) argument
1186
character*15 list(100)
1187
integer ilist(0:100),i,j,imode
1189
data firsttime/.true./
1203
do while (i.le.ilist(0))
1204
if(i.eq.ilist(0)) then
1205
list(i)=argument(1:15)
1210
if (argument(1:15).eq.list(i)) then
1216
write (*,*) 'error #1 in increasecnt'
1218
write (*,*) list(j),ilist(j)
1226
write (*,*) list(i),ilist(i)
1231
double precision function ran3(even)
1233
double precision ran2,get_ran
1244
subroutine initialize_even_random_numbers(ncalls0,ndim,ncalls)
1245
c Recompute the number of calls. Uses the algorithm from VEGAS
1247
integer ncalls0,ndim,ncalls,i
1248
integer dim,ng,npg,k
1250
common /even_ran/dim,ng,npg,k,firsttime
1251
c Make sure that hypercubes are newly initialized
1253
c Number of dimension of the integral
1255
c Number of elements in which we can split one dimension
1256
ng=(ncalls0/2.)**(1./ndim)
1257
c Total number of hypercubes
1259
c Number of PS points in each hypercube (at least 2)
1260
npg=max(ncalls0/k,2)
1261
c Number of PS points for this iteration
1267
double precision function get_ran()
1269
double precision ran2,dng
1271
integer dim,ng,npg,k
1273
common /even_ran/dim,ng,npg,k,firsttime
1275
parameter (maxdim=100)
1276
integer iii(maxdim),kkk(maxdim),i,iret
1278
save current_dim,dng,kkk,iii
1288
current_dim=mod(current_dim,dim)+1
1289
c This is the random number in the hypercube 'k' for current_dim
1290
get_ran=dng*(ran2()+dble(kkk(current_dim)-1))
1291
c Got random numbers for all dimensions, update kkk() for the next call
1292
if (current_dim.eq.dim) then
1293
call nextlexi(dim,iii,kkk,iret)
1295
call nextlexi(dim,iii,kkk,iret)
1302
subroutine init_ave_virt(ninter,ndim)
1305
integer kdim,ndim,ninter,i
1309
nvirt(i,kdim,kchan)=0
1310
ave_virt(i,kdim,kchan)=0d0
1311
nvirt_acc(i,kdim,kchan)=0
1312
ave_virt_acc(i,kdim,kchan)=0d0
1313
ave_born_acc(i,kdim,kchan)=0d0
1320
subroutine get_ave_virt(x,ninter,ndim,average_virtual)
1323
integer kdim,ndim,ninter,ncell
1324
double precision x(ndimmax),average_virtual(maxchannels)
1325
average_virtual(ichan)=0d0
1327
ncell=min(int(x(kdim)*ninter)+1,ninter)
1328
average_virtual(ichan)=average_virtual(ichan)+
1329
& ave_virt(ncell,kdim,ichan)
1331
average_virtual(ichan)=average_virtual(ichan)/ndim
1335
subroutine fill_ave_virt(x,ninter,ndim,virtual,born)
1338
integer kdim,ndim,ninter,ncell
1339
double precision x(ndimmax),virtual,born
1341
ncell=min(int(x(kdim)*ninter)+1,ninter)
1342
nvirt_acc(ncell,kdim,ichan)=nvirt_acc(ncell,kdim,ichan)+1
1343
ave_virt_acc(ncell,kdim,ichan)=ave_virt_acc(ncell,kdim,ichan)
1345
ave_born_acc(ncell,kdim,ichan)=ave_born_acc(ncell,kdim,ichan)
1351
subroutine regrid_ave_virt(ninter,ndim)
1354
integer ninter,ndim,kdim,i
1355
c need to solve for k_new = (virt+k_old*born)/born = virt/born + k_old
1359
if (ave_born_acc(i,kdim,kchan).eq.0d0) cycle
1360
if (ave_virt(i,kdim,kchan).eq.0d0) then ! i.e. first iteration
1361
ave_virt(i,kdim,kchan)= ave_virt_acc(i,kdim,kchan)
1362
$ /ave_born_acc(i,kdim,kchan)+ave_virt(i,kdim,kchan)
1363
else ! give some importance to the iterations already
1364
! done: simply base it on the number of points already
1365
! computed in that bin. Give half the weight to old iterations
1366
! k_new = [N_new*(V/B+k_old) + (N_old/2)*k_old] / [N_new+(N_old/2)]
1367
ave_virt(i,kdim,kchan)=(nvirt_acc(i,kdim,kchan)
1368
$ *ave_virt_acc(i,kdim,kchan)/ave_born_acc(i,kdim
1369
$ ,kchan)+(nvirt_acc(i,kdim,kchan)+nint(nvirt(i,kdim
1370
$ ,kchan)/2d0))*ave_virt(i,kdim,kchan))
1371
$ /dble(nvirt_acc(i,kdim,kchan)+nint(nvirt(i,kdim
1376
c reset the acc values
1379
nvirt(i,kdim,kchan)=nvirt(i,kdim,kchan)+nvirt_acc(i,kdim
1381
nvirt_acc(i,kdim,kchan)=0
1382
ave_born_acc(i,kdim,kchan)=0d0
1383
ave_virt_acc(i,kdim,kchan)=0d0
1391
subroutine double_ave_virt(ninter,ndim)
1394
integer kdim,ndim,i,ninter
1398
ave_virt(i*2,kdim,kchan)=ave_virt(i,kdim,kchan)
1399
if (nvirt(i,kdim,kchan).ne.0) then
1400
nvirt(i*2,kdim,kchan)=max(nvirt(i,kdim,kchan)/2,1)
1402
nvirt(i*2,kdim,kchan)=0
1405
ave_virt(i*2-1,kdim,kchan)=(ave_virt(i,kdim,kchan)
1406
$ +ave_virt(i-1,kdim,kchan))/2d0
1407
if (nvirt(i,kdim,kchan)+nvirt(i-1,kdim,kchan).ne.0) then
1408
nvirt(i*2-1,kdim,kchan)=
1409
& max((nvirt(i,kdim,kchan)+nvirt(i-1,kdim,kchan))/4,1)
1411
nvirt(i*2-1,kdim,kchan)=0
1414
if (nvirt(1,kdim,kchan).ne.0) then
1415
nvirt(1,kdim,kchan)=max(nvirt(1,kdim,kchan)/2,1)
1417
nvirt(1,kdim,kchan)=0
1428
subroutine get_channel(ans,vol_chan)
1429
c Picks one random 'ichan' among the 'nchans' integration channels and
1430
c fills the channels common block in mint.inc.
1433
double precision ans(0:maxchannels),vol_chan,ran3,target,sum
1435
if (nchans.eq.1) then
1437
iconfig=iconfigs(ichan)
1439
elseif (nchans.gt.1) then
1440
if (ans(0).le.0d0) then
1441
c pick one at random (flat)
1442
ichan=int(ran3(.false.)*nchans)+1
1443
iconfig=iconfigs(ichan)
1444
vol_chan=1d0/dble(nchans)
1446
c pick one at random (weighted by cross section)
1451
if (abs(sum-ans(0))/(sum+ans(0)).gt.1d-8) then
1452
write (*,*) 'ERROR: sum should be equal to ans',
1456
target=ans(0)*ran3(.false.)
1459
do while (sum.lt.target)
1463
if (ichan.eq.0 .or. ichan.gt.nchans) then
1464
write (*,*) 'ERROR: ichan cannot be zero or'/
1465
$ /' larger than nchans',ichan,nchans
1468
iconfig=iconfigs(ichan)
1469
vol_chan=ans(ichan)/ans(0)
added
removed
Template/NLO/SubProcesses/mint-integrator2.f
