← Back to branch summary

~maddevelopers/mg5amcnlo/3.0.2

« back to all changes in this revision

Viewing changes to Template/NLO/SubProcesses/driver_mintFO.f

  • Committer: Rikkert Frederix
  • Date: 2019-08-15 14:34:57 UTC
  • mfrom: (956.1.11 mint_module)
  • Revision ID: frederix@physik.uzh.ch-20190815143457-e0kkm72494mh2c6r
Major rewrite of the mint-integrator package. Converted it to a module
for easier extension with additional ideas.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/SubProcesses/driver_mintFO.f
3
3
c     This is the driver for the whole calculation
4
4
c**************************************************************************
5
5
      use extra_weights
 
6
      use mint_module
6
7
      implicit none
7
8
C
8
9
C     CONSTANTS
11
12
      parameter       (ZERO = 0d0)
12
13
      include 'nexternal.inc'
13
14
      include 'genps.inc'
14
 
      INTEGER    ITMAX,   NCALL
15
 
 
16
 
      common/citmax/itmax,ncall
17
15
C
18
16
C     LOCAL
19
17
C
20
 
      integer i,j,k,l,l1,l2
 
18
      integer i,j,k,l,l1,l2,kchan
21
19
      character*130 buf
22
20
c
23
21
c     Global
27
25
      include 'coupl.inc'
28
26
      
29
27
c     Vegas stuff
30
 
      integer         ndim
31
 
      common/tosigint/ndim
 
28
      integer         nndim
 
29
      common/tosigint/nndim
32
30
 
33
31
      real*8 sigint
34
32
      external sigint
43
41
      double precision xratmax
44
42
      common/ccheckcnt/i_momcmp_count,xratmax
45
43
 
 
44
      character*4      abrv
 
45
      common /to_abrv/ abrv
46
46
      integer n_mp, n_disc
47
 
c For MINT:
48
 
      include "mint.inc"
49
 
      integer nhits_in_grids(maxchannels)
50
 
      real* 8 xgrid(0:nintervals,ndimmax,maxchannels),ymax(nintervals
51
 
     $     ,ndimmax,maxchannels),ymax_virt(0:maxchannels),ans(nintegrals
52
 
     $     ,0:maxchannels),unc(nintegrals,0:maxchannels),chi2(nintegrals
53
 
     $     ,0:maxchannels),x(ndimmax),itmax_fl
54
47
      integer ixi_i,iphi_i,iy_ij,vn
55
 
      integer ifold(ndimmax) 
56
 
      common /cifold/ifold
57
 
      integer ifold_energy,ifold_phi,ifold_yij
58
 
      common /cifoldnumbers/ifold_energy,ifold_phi,ifold_yij
59
48
      logical putonshell
60
 
      integer imode,dummy
61
49
      logical unwgt
62
50
      double precision evtsgn
63
51
      common /c_unwgt/evtsgn,unwgt
72
60
 
73
61
      double precision virtual_over_born
74
62
      common/c_vob/virtual_over_born
75
 
      double precision average_virtual(0:n_ave_virt,maxchannels)
76
 
     $     ,virtual_fraction(maxchannels)
77
 
      common/c_avg_virt/average_virtual,virtual_fraction
78
63
      include 'orders.inc'
79
 
      integer              n_ord_virt
80
 
      common /c_n_ord_virt/n_ord_virt
81
64
 
82
65
c timing statistics
83
66
      include "timing_variables.inc"
95
78
      integer ntot_granny,derntot,ncase(0:6)
96
79
      common /c_granny_counters/ ntot_granny,ncase,derntot,deravg,derstd
97
80
     &     ,dermax,xi_i_fks_ev_der_max,y_ij_fks_ev_der_max
98
 
      logical              fixed_order,nlo_ps
99
 
      common /c_fnlo_nlops/fixed_order,nlo_ps
100
81
 
 
82
      logical useitmax
 
83
      common/cuseitmax/useitmax
101
84
 
102
85
C-----
103
86
C  BEGIN CODE
104
 
C-----  
 
87
C-----
 
88
      useitmax=.false. ! to be overwritten in open_output_files.f if need be
105
89
c
106
90
c     Setup the timing variable
107
91
c
112
96
c     Read general MadFKS parameters
113
97
c
114
98
      call FKSParamReader(paramFileName,.TRUE.,.FALSE.)
 
99
      min_virt_fraction_mint=min_virt_fraction
115
100
      do kchan=1,maxchannels
116
101
         do i=0,n_ave_virt
117
102
            average_virtual(i,kchan)=0d0
153
138
c     Get user input
154
139
c
155
140
      write(*,*) "getting user params"
156
 
      call get_user_params(ncall,itmax,imode)
 
141
      call get_user_params(ncalls0,itmax,imode)
157
142
      if(imode.eq.0)then
158
143
        flat_grid=.true.
159
144
      else
162
147
      ndim = 3*(nexternal-nincoming)-4
163
148
      if (abs(lpp(1)) .ge. 1) ndim=ndim+1
164
149
      if (abs(lpp(2)) .ge. 1) ndim=ndim+1
 
150
      nndim=ndim
165
151
c Don't proceed if muF1#muF2 (we need to work out the relevant formulae
166
152
c at the NLO)
167
153
      if( ( fixed_fac_scale .and.
172
158
        write(*,*)'NLO computations require muF1=muF2'
173
159
        stop
174
160
      endif
175
 
      write(*,*) "about to integrate ", ndim,ncall,itmax
 
161
      write(*,*) "about to integrate ", ndim,ncalls0,itmax
176
162
c APPLgrid
177
163
      if (imode.eq.0) iappl=0 ! overwrite when starting completely fresh
178
164
      if(iappl.ne.0) then
184
170
         call find_iproc_map
185
171
         write(6,*) "   ... done."
186
172
      endif
 
173
      if (abrv(1:4).eq.'virt') then
 
174
         only_virt=.true.
 
175
      else
 
176
         only_virt=.false.
 
177
      endif
 
178
c     Prepare the MINT folding
 
179
      do j=1,ndimmax
 
180
         if (j.le.ndim) then
 
181
            ifold(j)=1
 
182
         else
 
183
            ifold(j)=0
 
184
         endif
 
185
      enddo
 
186
      ifold_energy=ndim-2
 
187
      ifold_yij=ndim-1
 
188
      ifold_phi=ndim
 
189
c      
187
190
      i_momcmp_count=0
188
191
      xratmax=0.d0
189
192
      unwgt=.false.
194
197
            doreweight=.false.
195
198
            do_rwgt_scale=.false.
196
199
            do_rwgt_pdf=.false.
197
 
            do kchan=1,nchans
198
 
               do i=1,ndimmax
199
 
                  do j=0,nintervals
200
 
                     xgrid(j,i,kchan)=0.d0
201
 
                  enddo
202
 
               enddo
203
 
            enddo
204
200
         else
205
201
            doreweight=do_rwgt_scale.or.do_rwgt_pdf
206
 
c to restore grids:
207
 
            open (unit=12, file='mint_grids',status='old')
208
 
            ans(1,0)=0d0
209
 
            unc(1,0)=0d0
210
 
            do kchan=1,nchans
211
 
               do j=0,nintervals
212
 
                  read (12,*) (xgrid(j,i,kchan),i=1,ndim)
213
 
               enddo
214
 
               do j=1,nintervals_virt
215
 
                  do k=0,n_ord_virt
216
 
                     read (12,*) (ave_virt(j,i,k,kchan),i=1,ndim)
217
 
                  enddo
218
 
               enddo
219
 
               read(12,*) ans(1,kchan),unc(1,kchan),dummy,dummy
220
 
     $              ,nhits_in_grids(kchan)
221
 
               read(12,*) virtual_fraction(kchan),average_virtual(0
222
 
     $              ,kchan)
223
 
               ans(1,0)=ans(1,0)+ans(1,kchan)
224
 
               unc(1,0)=unc(1,0)+unc(1,kchan)**2
225
 
            enddo
226
 
            unc(1,0)=sqrt(unc(1,0))
227
 
            close (12)
228
 
            write (*,*) "Update iterations and points to",itmax,ncall
229
202
         endif
230
203
c
231
204
         write (*,*) 'imode is ',imode
236
209
               virtual_fraction(kchan)=1d0
237
210
            enddo
238
211
         endif
239
 
C check for zero cross-section
240
 
C if restoring grids corresponding to sigma=0, just terminate the run
241
 
         if (imode.ne.0.and.ans(1,0).eq.0d0.and.unc(1,0).eq.0d0) then
242
 
            call initplot()
243
 
            call close_run_zero_res(ncall, itmax, ndim)
244
 
            stop
245
 
         endif
246
 
         call mint(sigint,ndim,ncall,itmax,imode,xgrid,ymax
247
 
     $        ,ymax_virt,ans,unc,chi2,nhits_in_grids)
 
212
         call mint(sigint)
248
213
         call topout
249
214
         call deallocate_weight_lines
250
 
         write(*,*)'Final result [ABS]:',ans(1,0),' +/-',unc(1,0)
251
 
         write(*,*)'Final result:',ans(2,0),' +/-',unc(2,0)
252
 
         write(*,*)'chi**2 per D.o.F.:',chi2(1,0)
253
 
         open(unit=58,file='results.dat',status='unknown')
254
 
         do kchan=0,nchans
255
 
            write(58,*) ans(1,kchan),unc(2,kchan),0d0,0,0,0,0,0d0,0d0
256
 
     $           ,ans(2,kchan)
257
 
         enddo
258
 
         close(58)
259
215
c
260
 
c to save grids:
261
 
         open (unit=12, file='mint_grids',status='unknown')
262
 
         do kchan=1,nchans
263
 
            do j=0,nintervals
264
 
               write (12,*) (xgrid(j,i,kchan),i=1,ndim)
265
 
            enddo
266
 
            do j=1,nintervals_virt
267
 
               do k=0,n_ord_virt
268
 
                  write (12,*) (ave_virt(j,i,k,kchan),i=1,ndim)
269
 
               enddo
270
 
            enddo
271
 
            write (12,*) ans(1,kchan),unc(1,kchan),ncall,itmax
272
 
     $           ,nhits_in_grids(kchan)
273
 
            write (12,*) virtual_fraction(kchan),average_virtual(0
274
 
     $           ,kchan)
275
 
         enddo
276
 
         close (12)
277
216
      else
278
217
         write (*,*) 'Unknown imode',imode
279
218
         stop
351
290
      open (unit=12, file='res.dat',status='unknown')
352
291
      do kchan=0,nchans
353
292
         write (12,*)ans(1,kchan),unc(1,kchan),ans(2,kchan),unc(2,kchan)
354
 
     $        ,itmax,ncall,tTot
 
293
     $        ,itmax,ncalls0,tTot
355
294
      enddo
356
295
      close(12)
357
296
 
392
331
      double precision function sigint(xx,vegas_wgt,ifl,f)
393
332
      use weight_lines
394
333
      use extra_weights
 
334
      use mint_module
395
335
      implicit none
396
336
      include 'nexternal.inc'
397
 
      include 'mint.inc'
398
337
      include 'nFKSconfigs.inc'
399
338
      include 'run.inc'
400
339
      include 'orders.inc'
410
349
      integer             ini_fin_fks(maxchannels)
411
350
      common/fks_channels/ini_fin_fks
412
351
      data sum /.false./
413
 
      integer         ndim
414
 
      common/tosigint/ndim
 
352
      integer         nndim
 
353
      common/tosigint/nndim
415
354
      logical       nbody
416
355
      common/cnbody/nbody
417
356
      double precision p1_cnt(0:3,nexternal,-2:2),wgt_cnt(-2:2)
419
358
      common/counterevnts/p1_cnt,wgt_cnt,pswgt_cnt,jac_cnt
420
359
      double precision p_born(0:3,nexternal-1)
421
360
      common /pborn/   p_born
422
 
      double precision           virt_wgt_mint(0:n_ave_virt),
423
 
     &                           born_wgt_mint(0:n_ave_virt)
424
 
      common /virt_born_wgt_mint/virt_wgt_mint,born_wgt_mint
425
361
      double precision virtual_over_born
426
362
      common/c_vob/virtual_over_born
427
363
      logical                calculatedBorn
434
370
      common /c_wgt_ME_tree/ wgt_ME_born,wgt_ME_real
435
371
      integer ini_fin_fks_map(0:2,0:fks_configs)
436
372
      save ini_fin_fks_map
437
 
      logical new_point
438
 
      common /c_new_point/ new_point
439
373
      if (new_point .and. ifl.ne.2) then
440
374
         pass_cuts_check=.false.
441
375
      endif
484
418
      else
485
419
         jac=0.5d0
486
420
      endif
487
 
      call generate_momenta(ndim,iconfig,jac,x,p)
 
421
      call generate_momenta(nndim,iconfig,jac,x,p)
488
422
      if (p_born(0,1).lt.0d0) goto 12
489
423
      call compute_prefactors_nbody(vegas_wgt)
490
424
      call set_cms_stuff(izero)
526
460
         wgt_me_real=0d0
527
461
         jac=MC_int_wgt
528
462
         call update_fks_dir(iFKS)
529
 
         call generate_momenta(ndim,iconfig,jac,x,p)
 
463
         call generate_momenta(nndim,iconfig,jac,x,p)
530
464
         if (p_born(0,1).lt.0d0) cycle
531
465
         call compute_prefactors_n1body(vegas_wgt,jac)
532
466
         call set_cms_stuff(izero)
711
645
      end
712
646
 
713
647
      subroutine update_vegas_x(xx,x)
 
648
      use mint_module
714
649
      implicit none
715
 
      include 'mint.inc'
716
650
      integer i
717
651
      double precision xx(ndimmax),x(99),ran2
718
652
      external ran2
719
 
      integer ndim
720
 
      common/tosigint/ndim
721
 
      character*4 abrv
 
653
      integer         nndim
 
654
      common/tosigint/nndim
 
655
      character*4      abrv
722
656
      common /to_abrv/ abrv
723
657
      do i=1,99
724
658
         if (abrv.eq.'born'.or.abrv(1:2).eq.'vi') then
725
 
            if(i.le.ndim-3)then
 
659
            if(i.le.nndim-3)then
726
660
               x(i)=xx(i)
727
 
            elseif(i.le.ndim) then
 
661
            elseif(i.le.nndim) then
728
662
               x(i)=ran2()      ! Choose them flat when not including real-emision
729
663
            else
730
664
               x(i)=0.d0
731
665
            endif
732
666
         else
733
 
            if(i.le.ndim)then
 
667
            if(i.le.nndim)then
734
668
               x(i)=xx(i)
735
669
            else
736
670
               x(i)=0.d0
741
675
      end
742
676
 
743
677
c
744
 
      subroutine get_user_params(ncall,itmax,irestart)
 
678
      subroutine get_user_params(ncall,nitmax,irestart)
745
679
c**********************************************************************
746
680
c     Routine to get user specified parameters for run
747
681
c**********************************************************************
 
682
      use mint_module
748
683
      implicit none
749
684
c
750
685
c     Constants
754
689
      include 'nFKSconfigs.inc'
755
690
      include 'fks_info.inc'
756
691
      include 'run.inc'
757
 
      include 'mint.inc'
758
692
      include 'orders.inc'
759
693
c
760
694
c     Arguments
761
695
c
762
 
      integer ncall,itmax
 
696
      integer ncall,nitmax
763
697
c
764
698
c     Local
765
699
c
766
 
      integer i, j
 
700
      integer i, j, kchan
767
701
      double precision dconfig(maxchannels)
768
702
c
769
703
c     Global
797
731
      character * 70 idstring
798
732
      logical savegrid
799
733
 
800
 
      logical usexinteg,mint
801
 
      common/cusexinteg/usexinteg,mint
802
734
      logical unwgt
803
735
      double precision evtsgn
804
736
      common /c_unwgt/evtsgn,unwgt
813
745
c-----
814
746
c  Begin Code
815
747
c-----
816
 
      mint=.true.
817
748
      unwgt=.false.
818
749
      open (unit=83,file='input_app.txt',status='old')
819
750
      done=.false.
825
756
            read(buffer,*) ncall
826
757
            write (*,*) 'Number of phase-space points per iteration:',ncall
827
758
         elseif(buffer(1:11).eq.'NITERATIONS') then
828
 
            read(buffer(14:),*) itmax
829
 
            write (*,*) 'Maximum number of iterations is:',itmax
 
759
            read(buffer(14:),*) nitmax
 
760
            write (*,*) 'Maximum number of iterations is:',nitmax
830
761
         elseif(buffer(1:8).eq.'ACCURACY') then
831
762
            read(buffer(11:),*) accuracy
832
763
            write (*,*) 'Desired accuracy is:',accuracy