Merge lp:~maddm/maddm/maddm_resonances into lp:maddm

[Olivier Mattelaer]

fix a couple of issues including a compatibility issue with 2.6.6

[lp:~maddm/maddm/maddm_resonances]

35. By olivier-mattelaer on 2019-09-04

fix an issue with compatibility with madgraph

[Chiara Arina (carina)]

Salut Olivier,

je ne sais pas si je dois faire quelque chose mais j’ai regarde le merging que tu proposes et c’est ok pour moi. Dit moi si je dois confirmer quelque part sur le site web..
Chiara

[lp:~maddm/maddm/maddm_resonances]

36. By olivier-mattelaer on 2019-09-05

revert chiara QED->QCD change but allow to test [noborn=QED] if no diagram are generated for [noborn=QCD]

=== modified file 'Templates/matrix_elements/makefile'
--- Templates/matrix_elements/makefile	2016-05-14 13:38:36 +0000
+++ Templates/matrix_elements/makefile	2019-09-05 08:03:10 +0000
@@ -1,9 +1,10 @@
-fort= gfortran
+include ../Source/make_opts
+fort= $(FC)
 
 libname = matrix_elements
 incdir = ../include
 
-flags = -O -fno-align-commons -ffixed-line-length-132 -I$(incdir)
+flags = $(FFLAGS) -I$(incdir)
 
 objs = smatrix.o $(patsubst %.f,%.o,$(wildcard *.f))
 

=== modified file 'Templates/src/directional_detection.f'
--- Templates/src/directional_detection.f	2017-07-27 15:49:35 +0000
+++ Templates/src/directional_detection.f	2019-09-05 08:03:10 +0000
@@ -1831,7 +1831,7 @@
       integer flag
       double precision x, eff 
 
-
+      efficiency = 0d0
       if(flag .eq. 0) then
          efficiency = 1.d0 ! for 100% detection efficiency
 

=== modified file 'Templates/src/init.f'
--- Templates/src/init.f	2017-01-23 13:00:58 +0000
+++ Templates/src/init.f	2019-09-05 08:03:10 +0000
@@ -107,9 +107,13 @@
 c Y_eq = 45/(4*pi^4) g_i/g_*S x^2 K_2(x), with x=m_i/T
 c K_2(x) is expanded for large values of x ( > 10)
       if (mass.gt.5.d0*temp) then
+        if (mass/temp.gt.706)then
+           get_Y_eq = 0d0
+        else
         get_Y_eq = 45.d0/(4.d0*pi**4)*(dof/rel_dofs)*(mass/temp)**(1.5d0)
      .        * dsqrt(pi/2.d0)*dexp(-mass/temp) * (1.d0 + 15.d0/(8.d0*(mass/temp))
      .        + 105.d0/(128.d0*(mass/temp)**2) - 315.d0/(1024.d0*(mass/temp)**3))
+      endif
       else if (mass.gt.temp/5.d0) then
         get_Y_eq = 45.d0/(4.d0*pi**4)*(dof/rel_dofs)*(mass/temp)**2
      .        * bessk(2,mass/temp)
@@ -194,6 +198,11 @@
       double precision Wij_hold, Wij_total, beta_step
       logical flag1, flag2, flag3
 
+      integer i, ii
+      double precision next_beta(nres), tmp_beta
+      integer res_step(nres)
+      integer nres_points_local, max_width_factor
+
 c output to the screen to show progress
       if (print_out) then
         write(*,*) 'Calculating all the Wijs'
@@ -202,9 +211,28 @@
         flag3 = .true.
       endif
 
+      nres_points_local = 100
+      max_width_factor = 6 ! this is e^N-1 times the width for the further point
+
+      do ii=1, nres
+         res_step(ii) = - nres_points_local
+         tmp_beta = -1d0
+         if (beta_res(ii).lt.0d0) then
+            res_step(ii) = nres_points_local +1
+            cycle
+         endif
+         do while (tmp_beta.lt.0d0)
+            tmp_beta = beta_res(ii) - beta_res_width(ii)*(DEXP(1d0*res_step(ii)*max_width_factor/nres_points_local)-1)
+            res_step(ii) = res_step(ii) + 1
+         enddo
+c         write(*,*) 'resonances',ii,beta_res(ii), beta_res_width(ii), res_step(ii)
+      enddo
+
+
 c setting up step sizes
-      beta_step = 0.01d0
-
+      beta_step = 1e-3
+      beta_step_min =  1e-6
+      beta_step_max = 1e-2
 c starting the calculation of the Wij's
       nWij = 1
       betas(nWij) = 0.001d0
@@ -213,11 +241,32 @@
 c continue to do beta steps until we are near beta = 1
       do while (betas(nWij).le.0.99d0-beta_step)
         nWij = nWij+1
-        betas(nWij) = betas(nWij-1) + beta_step
+        tmp_beta = betas(nWij-1) + beta_step
+        do ii=1, nres
+           if (res_step(ii).le.0)then
+              next_beta(ii) = beta_res(ii) - beta_res_width(ii)*(DEXP(1d0*res_step(ii)*max_width_factor/nres_points_local)-1)
+           elseif(res_step(ii).le.nres_points_local)then
+               next_beta(ii) = beta_res(ii) + beta_res_width(ii)*(DEXP(1d0*res_step(ii)*max_width_factor/nres_points_local)-1)
+           else
+              next_beta(ii) = 1d0
+            endif
+            tmp_beta = min(tmp_beta, next_beta(ii),1d0)
+         enddo
+        betas(nWij) = tmp_beta
         Wij_total = get_Wij_ann(betas(nWij))
 
+
 c if the new Wij is 100% larger than the previous Wij, then reduce the stepsize
         if (((dabs(Wij_total-Wij_hold)/Wij_hold).gt.1.d0).and.(beta_step.gt.beta_step_min)) then
+c           write(*,*) 'reduce', betas(nWij), beta_step, ((dabs(Wij_total-Wij_hold)/Wij_hold)), Wij_total, Wij_hold
+
+           if (((dabs(Wij_total-Wij_hold)/Wij_hold).gt.2.d0))then
+c              write(*,*) 'rewind',  betas(nWij-1)
+              nWij = nWij-1
+              beta_step = beta_step/2.d0
+             cycle
+           endif
+
           beta_step = beta_step/2.d0
           if (beta_step .lt. beta_step_min) then
             beta_step = beta_step_min
@@ -226,12 +275,24 @@
 
 c if the new Wij is within 10% of the previous Wij, then increase the stepsize
         if (((dabs(Wij_total-Wij_hold)/Wij_hold).lt.0.1d0).and.(beta_step.lt.beta_step_max)) then
-          beta_step = beta_step*2.d0
+          beta_step = beta_step*2d0
           if (beta_step .gt. beta_step_max) then
             beta_step = beta_step_max
           endif
+c          write(*,*) 'increase', betas(nWij), beta_step
         endif
 
+         do ii=1, nres
+            if (next_beta(ii).eq.tmp_beta) then
+               res_step(ii) = res_step(ii)+1
+               beta_step = max(beta_step,1e-3)
+c               if (beta_step.eq.1e-3) then
+c                  write(*,*) 'reset'
+c               endif
+            endif
+         enddo
+
+
 c Hold the Wij_value so that if the Wij integral is too small we have a stored value ot use
         Wij_hold = Wij_total
 
@@ -256,7 +317,7 @@
       betas(nWij) = 0.999d0
       Wij_total = get_Wij_ann(betas(nWij))
 
-      if (print_out) write(*,*) 'Done!'
+      if (print_out)  write(*,*) 'Done!', nWij
       Wij_ann_calc = .true.
 
       return
@@ -284,7 +345,6 @@
 
 c set beta_pass to pass beta to Wij_integrand
       beta_pass = beta
-
 c loop over all the dm particles
       do x1=1, ndmparticles
         do x2=x1, ndmparticles
@@ -311,6 +371,54 @@
       return
       end
 
+c-------------------------------------------------------------------------c
+      function get_Wij_ann_nogrid(beta)
+c-------------------------------------------------------------------------c
+c                                                                         c
+c  Given the input velocity beta, this function integrates the matrix     c
+c  elements generated by MadGraph over the solid angle for all            c
+c  combinations of \chi_i, \chi_j to all final states.                    c
+c                                                                         c
+c-------------------------------------------------------------------------c
+      implicit none
+      include 'maddm.inc'
+
+c input parameters
+      double precision beta
+      double precision Wij_ann_nogrid(maxdms,maxdms)
+c external function used in this subroutine
+      external Wij_ann_integrand
+
+c set beta_pass to pass beta to Wij_integrand
+      beta_pass = beta
+
+      get_Wij_ann_nogrid =  0d0
+c loop over all the dm particles
+      do x1=1, ndmparticles
+        do x2=x1, ndmparticles
+          dmi = x1
+          dmj = x2
+
+c integrate the 2 -> 2 process over x = (cos(theta)+1)/2.d0
+          call romberg(Wij_ann_integrand, 0.d0, 1.d0, Wij_ann_nogrid(x1,x2), eps_wij, iter_wij)
+
+c the Wijs for (DMi, DMj) are the same for (DMj, DMi)
+          if (x1.ne.x2) Wij_ann_nogrid(x2,x1) = Wij_ann_nogrid(x1,x2)
+        enddo
+      enddo
+
+c the return of this function will be the sum of all the Wij's 
+c so the next beta step can be determined
+      get_Wij_ann = 0.d0
+      do x1=1, ndmparticles
+        do x2=x1, ndmparticles
+          get_Wij_ann_nogrid = get_Wij_ann_nogrid + Wij_ann_nogrid(x1,x2)
+        enddo
+      enddo
+
+      return
+      end
+
 
 
 c-------------------------------------------------------------------------c
@@ -483,4 +591,4 @@
       read(*,*)
 
       return
-      end
\ No newline at end of file
+      end

=== modified file 'Templates/src/integrate.f'
--- Templates/src/integrate.f	2017-08-10 20:21:07 +0000
+++ Templates/src/integrate.f	2019-09-05 08:03:10 +0000
@@ -89,6 +89,8 @@
       double precision left, right, whole, end_pt, start_pt, width
 
       integer ii, jj, kk, grid_pos
+      integer nres_points_local,max_width_factor
+      
 
 c    Set up the integration grid
 c    Initialize to 10*b, so that we can identify the relevant points
@@ -97,7 +99,7 @@
           grid=10.d0*b
 
 c         write(*,*) 'width: ', beta_res_width(1)
-
+          call find_resonances()
 c    Make sure that the array is large enough to store the grid
          if (grid_npts.lt.ngrid_init) then
              write(*,*) 'Error: grid array not large enough!'
@@ -116,14 +118,18 @@
                 grid_pos = grid_pos+1
             endif
        enddo
+       
+       nres_points_local = 20
+       max_width_factor = 2 ! allow to go as far as e^N-1 times the width for the extra point
 
 c then add more points around the resonance
        do ii=1, nres
           if (beta_res(ii).ge.a.and.beta_res(ii).lt.b.and.beta_res(ii).ge.0.d0) then
-              do jj=1, nres_points
+              do jj=1, nres_points_local
 c                pts_to_add_adaptive = ceiling(real(pts_to_add_adaptive / 2))
 c                do kk = 1, pts_to_add_adaptive
-                    additional_pt = beta_res(ii) + (exp(real(beta_res_width(ii))/10.d0*exp(real(jj)))-1.d0)
+                    additional_pt = beta_res(ii) + beta_res_width(ii)*(DEXP(1d0*jj**max_width_factor/nres_points_local)-1)
+c                    additional_pt = beta_res(ii) + (exp(real(beta_res_width(ii))/10.d0*exp(real(jj)))-1.d0)
                     if (additional_pt.lt.b) then
                         grid(grid_pos) = additional_pt
                         grid_pos=grid_pos+1
@@ -133,7 +139,7 @@
                         endif
                     endif
 
-                    additional_pt = width - (exp(real(beta_res_width(ii))/10.d0*exp(real(jj)))-1.d0)
+                    additional_pt = beta_res(ii) - beta_res_width(ii)*(DEXP(1d0*jj**max_width_factor/nres_points_local)-1)
                     if (additional_pt.gt.a) then
                         grid(grid_pos) = additional_pt
                         grid_pos=grid_pos+1
@@ -196,7 +202,7 @@
       external func
       integer gr_npts, ii
       double precision grid_simp(gr_npts)
-
+      double precision tmp1, tmp2, tmp3
       simpson = 0.d0
       ii = 1
 
@@ -215,8 +221,16 @@
 c            write(101,*) start_pt, func(start_pt)
 c           write(*,*) 'start, end', start_pt, end_pt
 c           write(*,*) 'func: ', func(0.d0), func(0.00001d0)
-           simpson = simpson+ (end_pt - start_pt)/6.d0*(func(start_pt) + 4.d0*func(0.5d0*
+
+           tmp1 = (func(start_pt) + 4.d0*func(0.5d0*
      .                                  (start_pt+end_pt)) + func(end_pt) )
+           if (tmp1.lt.1e-20*simpson)then
+              exit
+           else if (tmp1.gt.1d-300) then
+           simpson = simpson+ (end_pt - start_pt)/6.d0*tmp1
+           else 
+              exit
+           endif
 c            simpson = simpson+ (end_pt - start_pt)/8.d0*(max(0.d0, func(start_pt))
 c     .             + 3.d0*max(0.d0,func(0.33d0*(2.d0*start_pt+end_pt)))
 c     .             + 3.d0*max(0.d0,func(0.33d0*(start_pt+2.d0*end_pt)))  + max(0.d0,func(end_pt)))
@@ -471,8 +485,10 @@
           endif
           return
         endif
-c        if(dabs(hnext) .lt. hmin) pause
-c     &     'stepsize smaller than minimum in odeint'
+        if(dabs(hnext) .lt. hmin) then
+           write(*,*) 'too small step'
+           stop 1
+        endif 
         h = hnext
         enddo
 c      pause 'too many steps in odeint'
@@ -682,7 +698,7 @@
 
 c Set up for stratification
         if (mds.ne.0) then
-          ng = (dble(ncall)/2.d0+0.25d0)**(1.d0/dble(ndim))
+          ng = int((dble(ncall)/2.d0+0.25d0)**(1.d0/dble(ndim)))
           mds = 1
           if ((2*ng-ndmx).ge.0) then
             mds = -1

=== modified file 'Templates/src/interpolate.f'
--- Templates/src/interpolate.f	2015-11-26 03:18:05 +0000
+++ Templates/src/interpolate.f	2019-09-05 08:03:10 +0000
@@ -21,6 +21,7 @@
 c parameters used to find the index
       integer index, lowindex, highindex, mid
 
+      get_Wij_value = 0d0
 c If the xvalue is outside the given array then set the index to the boundary
       if (beta .lt. betas(1)) then
         index = 1
@@ -42,6 +43,8 @@
 
 c interpolate the array with a line
       if (group.eq.1) then
+c         write(*,*) 'bypass interpolation'
+c        get_Wij_value = get_Wij_ann_nogrid(beta)
         get_Wij_value = (Wij_ann(i,j,index+1)-Wij_ann(i,j,index))/(betas(index+1)-betas(index))
      .        *(beta-betas(index))+Wij_ann(i,j,index)
       else if (group.eq.2) then
@@ -80,7 +83,7 @@
 
 c parameters used to find the index
       integer x_index, lowindex, highindex, mid
-
+      get_taacs_value = 0d0
 c If the xvalue is outside the given array then set the index to the boundary
       if (x .lt. x_taacs(1)) then
         x_index = 1
@@ -392,4 +395,4 @@
       endif
 
       return
-      end
\ No newline at end of file
+      end

=== modified file 'Templates/src/makefile'
--- Templates/src/makefile	2017-08-17 10:04:18 +0000
+++ Templates/src/makefile	2019-09-05 08:03:10 +0000
@@ -1,10 +1,11 @@
+include ../Source/make_opts
 prog = maddm
-fort= gfortran
+fort= $(FC)
 
 libdir = ../lib/
 incdir = ../include
 
-flags = -O -fno-align-commons -ffixed-line-length-132 -L$(libdir) -I$(incdir) -fbounds-check
+flags = $(FFLAGS) -L$(libdir) -I$(incdir)
 
 libs = -lmatrix_elements -ldhelas -lmodel
 

=== modified file 'Templates/src/relic_canon.f'
--- Templates/src/relic_canon.f	2017-08-10 20:21:07 +0000
+++ Templates/src/relic_canon.f	2019-09-05 08:03:10 +0000
@@ -81,7 +81,7 @@
             endif
 
 			x_1 = x_1 - dx_step
-            if (x_1.le.1d0) then
+           if (x_1.le.1d0) then
 			    write(*,*) 'Freezeout occurs too early! Relic density too big to comprehend. Seeting Omegah^2 = -1.0 \n'
 			    relic_canon=-1.0
 			    return
@@ -255,10 +255,16 @@
       enddo
 
 c integrand for the thermally averaged annihilation cross section
+      if (2.d0*x_pass*(1.d0-1.d0/(1.d0-beta**2)**(0.5d0)).lt.-695)then
+         taacs_integrand_canon =mdm(1)**2*beta/(1.d0-beta**2)**2 * Wij_sum / (2.d0*T*K2_sum**2)
+     .    * 0d0
+     .    * (1.d0 + 3.d0*T/(8.d0*dsqrt(s)) - 15.d0*T**2/(128.d0*s) + 105.d0*T**3/(1024.d0*s**1.5d0))
+     .    / (1.d0 + 15.d0*T/(8.d0*mdm(1)) + 105.d0*T**2/(128.d0*mdm(1)**2) - 315.d0*T**3/(1024.d0*mdm(1)**3))**2
+      else
       taacs_integrand_canon = mdm(1)**2*beta/(1.d0-beta**2)**2 * Wij_sum / (2.d0*T*K2_sum**2)
      .    * dsqrt(2.d0/(pi*dsqrt(s)*T))*mdm(1)*dexp(2.d0*x_pass*(1.d0-1.d0/(1.d0-beta**2)**(0.5d0)))
      .    * (1.d0 + 3.d0*T/(8.d0*dsqrt(s)) - 15.d0*T**2/(128.d0*s) + 105.d0*T**3/(1024.d0*s**1.5d0))
      .    / (1.d0 + 15.d0*T/(8.d0*mdm(1)) + 105.d0*T**2/(128.d0*mdm(1)**2) - 315.d0*T**3/(1024.d0*mdm(1)**3))**2
-
+      endif
       return
       end

=== modified file 'Templates/src/relic_density.f'
Binary files Templates/src/relic_density.f	2016-11-04 11:35:41 +0000 and Templates/src/relic_density.f	2019-09-05 08:03:10 +0000 differ
=== modified file 'Templates/src/tests.f'
--- Templates/src/tests.f	2017-04-11 13:16:20 +0000
+++ Templates/src/tests.f	2019-09-05 08:03:10 +0000
@@ -637,6 +637,7 @@
         write(*,*) 'Invalid value of p_or_E'
         write(*,*) 'p_or_E = 1 for initial momentum'
         write(*,*) 'p_or_E = 0 for initial energy'
+      cross_check_process = 0d0
         return
       endif
 
@@ -768,6 +769,10 @@
           msq = smatrix_ann(p_ext,i_pass,j_pass,k_pass)
         else if (group_pass.eq.2) then
           msq = smatrix_dm2dm(p_ext,i_pass,j_pass,k_pass)
+        else
+           msq = 0d0
+           write(*,*) 'undefined msq'
+           stop 1
 C         else if (group_pass.eq.3) then
 C           msq = smatrix_scattering(p_ext,i_pass,j_pass,k_pass)
         endif

=== modified file 'maddm_interface.py'
--- maddm_interface.py	2019-01-10 10:53:52 +0000
+++ maddm_interface.py	2019-09-05 08:03:10 +0000
@@ -139,8 +139,8 @@
         super(MadDM_interface, self).preloop(*args, **opts)
         self.prompt = 'MadDM>'
     
-    def change_principal_cmd(self, name):
-        out = super(MadDM_interface, self).change_principal_cmd(name)
+    def change_principal_cmd(self, name, *args, **opts):
+        out = super(MadDM_interface, self).change_principal_cmd(name,*args, **opts)
         self.prompt = 'MadDM>'
         return out
     
@@ -1087,8 +1087,11 @@
                 except (self.InvalidCmd, diagram_generation.NoDiagramException), error:
                     if allow_loop_induce:
                         proc = '%s %s > %s %s [noborn=QCD] @ID ' % (name, antiname, ' '.join(argument), coupling)
-                        self.do_add('process %s' % proc)
-  
+                        try:                    
+                            self.do_add('process %s' % proc)
+                        except (self.InvalidCmd, diagram_generation.NoDiagramException), error:
+                            proc = '%s %s > %s %s [noborn=QED] @ID ' % (name, antiname, ' '.join(argument), coupling)
+                            self.do_add('process %s' % proc)
   
     def install_indirect(self):
         """Code to install the reduction library if needed"""

=== modified file 'maddm_run_interface.py'
--- maddm_run_interface.py	2019-01-10 10:53:52 +0000
+++ maddm_run_interface.py	2019-09-05 08:03:10 +0000
@@ -74,7 +74,7 @@
 
     def __init__(self):
 
-        self._allowed_final_states = ['qqx', 'ccx', 'gg', 'bbx', 'ttx', 'e+e-', 'mu+mu-', 'ta+ta-', 'w+w-', 'zz', 'hh', 
+        self._allowed_final_states = ['qqx', 'ccx', 'gg', 'bbx', 'ttx', 'emep', 'mummup', 'tamtap', 'wpwm', 'zz', 'hh', 
                                       'hess2013','hess2016', 'aaER16','aaIR90','aaNFWcR3','aaNFWR41']
 
         self._oh2_planck = 0.1198
@@ -88,10 +88,10 @@
                                'gg'    :pjoin(MDMDIR, 'ExpData', 'MadDM_Fermi_Limit_gg.dat'),
                                'bbx'   :pjoin(MDMDIR, 'ExpData', 'MadDM_Fermi_Limit_bb.dat'),
                                'ttx'   :pjoin(MDMDIR, 'ExpData', 'MadDM_Fermi_Limit_tt.dat'),
-                               'e+e-'  :pjoin(MDMDIR, 'ExpData', 'MadDM_Fermi_Limit_ee.dat'),
-                               'mu+mu-':pjoin(MDMDIR, 'ExpData', 'MadDM_Fermi_Limit_mumu.dat'),
-                               'ta+ta-':pjoin(MDMDIR, 'ExpData', 'MadDM_Fermi_Limit_tautau.dat'),
-                               'w+w-'  :pjoin(MDMDIR, 'ExpData', 'MadDM_Fermi_Limit_WW.dat'),
+                               'emep'  :pjoin(MDMDIR, 'ExpData', 'MadDM_Fermi_Limit_ee.dat'),
+                               'mummup':pjoin(MDMDIR, 'ExpData', 'MadDM_Fermi_Limit_mumu.dat'),
+                               'tamtap':pjoin(MDMDIR, 'ExpData', 'MadDM_Fermi_Limit_tautau.dat'),
+                               'wpwm'  :pjoin(MDMDIR, 'ExpData', 'MadDM_Fermi_Limit_WW.dat'),
                                'zz'    :pjoin(MDMDIR, 'ExpData', 'MadDM_Fermi_Limit_ZZ.dat'),
                                'hh'    :pjoin(MDMDIR, 'ExpData', 'MadDM_Fermi_Limit_hh.dat'),
 
@@ -103,8 +103,8 @@
                                'aaNFWcR3':pjoin(MDMDIR,'ExpData', 'Fermi_lines_2015_NFWcontracted_R3.dat'),
                                'aaNFWR41':pjoin(MDMDIR,'ExpData', 'Fermi_lines_2015_NFW_R41.dat')}
 
-        self._id_limit_vel = {'qqx':2.0E-5,   'ccx':2.0E-5, 'gg':2.0E-5, 'bbx':2.0E-5,'ttx':2.0E-5,'e+e-':2.0E-5,'mu+mu-':2.0E-5,'ta+ta-':2.0E-5,
-                              'w+w-':2.0E-5,  'zz':2.0E-5,  'hh':2.0E-5,
+        self._id_limit_vel = {'qqx':2.0E-5,   'ccx':2.0E-5, 'gg':2.0E-5, 'bbx':2.0E-5,'ttx':2.0E-5,'emep':2.0E-5,'mummup':2.0E-5,'tamtap':2.0E-5,
+                              'wpwm':2.0E-5,  'zz':2.0E-5,  'hh':2.0E-5,
                               'aaER16':1.0E-3,'aaIR90':1.0E-3,'aaNFWcR3':1.0E-3,'aaNFWR41':1.0E-3 ,
                               'hess2013': 1.0E-3 , 'hess2016': 1.0E-3 } 
 
@@ -243,7 +243,7 @@
         self.channels      = ['ee', 'mumu', 'tautau', 'qq', 'cc', 'bb', 'tt', 'ZZ', 'WW', 'hh', 'gammagamma', 'gg']
 
         self.map_allowed_final_state_PPPC = {'qqx':'qq', 'ccx':'cc', 'gg':'gg', 'bbx':'bb', 'ttx':'tt',
-                                             'e+e-':'ee', 'mu+mu-':'mumu', 'ta+ta-':'tautau', 'w+w-':'WW', 'zz':'ZZ', 'hh':'hh' }
+                                             'emep':'ee', 'mummup':'mumu', 'tamtap':'tautau', 'wpwm':'WW', 'zz':'ZZ', 'hh':'hh' }
 
     def check_mass(self,mdm):
         if (mdm < 5.0 or mdm > 100000):
@@ -265,13 +265,13 @@
             return
 
         if not corr:
-            dic =  np.load(PPPCDIR+'/PPPC_Tables_noEW.npy').item()
+            dic =  np.load(PPPCDIR+'/PPPC_Tables_noEW.npy', allow_pickle = True).item()
             logger.info('PPPC4DMID Spectra at source loaded')
             Spectra.PPPCdata = dic
             Spectra.PPPC_type = ('source', corr)
             return dic
         elif corr == 'ew':
-            dic =  np.load(PPPCDIR+'/PPPC_Tables_noEW.npy').item()
+            dic =  np.load(PPPCDIR+'/PPPC_Tables_noEW.npy', allow_pickle= True).item()
             logger.info('PPPC4DMID Spectra at source (with EW corrections) loaded')
             Spectra.PPPCdata = dic
             Spectra.PPPC_type = ('source', corr)
@@ -287,7 +287,7 @@
             return
         
         else:
-            dic = np.load(PPPCDIR+'/PPPC_Tables_epEarth_'+prof+'.npy').item()
+            dic = np.load(PPPCDIR+'/PPPC_Tables_epEarth_'+prof+'.npy' , allow_pickle= True).item()
             Spectra.PPPCdata = dic
             Spectra.PPPC_type = ('earth', prof)
             
@@ -934,10 +934,10 @@
  
             # *** Indirect detection
             if self.mode['indirect']:
-                halo_vel = self.maddm_card['vave_indirect']
-                if halo_vel > (3*10**(-6)) and halo_vel < ( 1.4*10**(-4) ):  # cannot calculate Fermi limits
+                    halo_vel = self.maddm_card['vave_indirect']
+                    #if halo_vel > (3*10**(-6)) and halo_vel < ( 1.4*10**(-4) ):  # cannot calculate Fermi limits
                  
-                #if not self._two2twoLO:
+                    #if not self._two2twoLO:
                     detailled_keys = [k for k in self.last_results if k.startswith('taacsID#') ]
                     if len(detailled_keys)>1:
                         for key in detailled_keys:
@@ -945,7 +945,8 @@
                             clean_key = clean_key_list[0]+"_"+clean_key_list[1]
 
                             order +=[clean_key]
-                            order +=['lim_'+clean_key]
+                            if halo_vel > (3*10**(-6)) and halo_vel < ( 1.4*10**(-4) ):
+                                    order +=['lim_'+clean_key]
 
                     order.append('taacsID')
                     order.append('tot_SM_xsec')
@@ -1746,7 +1747,7 @@
 
           detailled_keys = [k for k in self.last_results.keys() if k.startswith('taacsID#')]
           logger.info('<sigma v> method: %s ' % self.maddm_card['sigmav_method'] )
-          logger.info('DM particle halo velocity: %s/c ' % halo_vel) # the velocity should be the same                                                           
+          logger.info('DM particle halo velocity: %s/c ' % halo_vel) # the velocity should be the same
 
           skip = []
           if halo_vel > (3*10**(-6)) and halo_vel < ( 1.4*10**(-4) ): # range of validity of Fermi limits
@@ -1770,7 +1771,7 @@
                         skip.append(f_original)
                         continue 
                     if finalstate in self.limits._allowed_final_states :
-                        s_ul   = self.limits.ID_max(mdm, finalstate)    
+                        s_ul   = self.limits.ID_max(mdm, finalstate)
                         self.last_results['Fermi_lim_'+key] = self.limits.ID_max(mdm, finalstate)
                     else:  s_ul   = '-1'
 
@@ -1910,7 +1911,7 @@
             out_dir = dir_point # all the various output must be saved here ==  pjoin(events, run_name )     
    
             if 'off' in save_switch:
-                if 'inclusive' not in self.maddm_card['sigmav_method'] : # here, need to remove everyhting i.e. the whole run_xx folder
+                if 'inclusive' not in self.maddm_card['sigmav_method'] : # here, need to remove everything i.e. the whole run_xx folder
                     shutil.rmtree(out_dir)
                                       
             elif 'all' in save_switch :
@@ -1921,7 +1922,7 @@
                         shutil.move(f_path , pjoin(self.dir_path , 'output', run_name, F) )
 
                  
-                if 'inclusive' in self.maddm_card['sigmav_method']: # saving spectra onyl in inclusive mode since madevent/resh have their own pythia8 spectra             
+                if 'inclusive' in self.maddm_card['sigmav_method']: # saving spectra only in inclusive mode since madevent/resh have their own pythia8 spectra             
                     self.save_spec_flux(out_dir = pjoin(self.dir_path , 'output', run_name), \
                                             spec_source = spec_source, flux_source = flux_source, flux_earth = flux_earth)
                  
@@ -2540,9 +2541,17 @@
         pythia8_card_path = pjoin(opts['mother_interface'].dir_path, 'Cards', 'pythia8_card.dat')
 
         cards = [param_card_path, 'maddm', 'multinest', pythia8_card_path]
+        
+        tmp_allow_arg = list(self.allow_arg)
+        opts = dict(opts)
+        opts['allow_arg'] = []
+        
         common_run.AskforEditCard.__init__(self, question, cards,
-                                            *args, **opts)
+                                            *args, **opts )
+
+        self.allow_arg += tmp_allow_arg
         self.question = self.create_question()
+
                
     def write_switch(self, path=None):
         """store value of the switch for the default at next run"""
@@ -2822,7 +2831,7 @@
 
             self.maddm.do_help(args[start])
             
-        ### CHECK if a help_xxx exist (only for those wihtout a do_xxx)
+        ### CHECK if a help exist (only for those wihtout a do_xxx)
         if len(args) == 1:
             if not hasattr(self, 'do_%s' % args[0]) and hasattr(self, 'help_%s' % args[0]):
                 getattr(self, 'help_%s' %args[0])()
@@ -3185,7 +3194,7 @@
         self.add_param('dx_step', 1.0, hidden=True)
         
         self.add_param('ngrid_init', 50, hidden=True, comment="Initial number of points in the grid to integrate over (of dm velocity)")
-        self.add_param('nres_points', 50, hidden=True, comment="Number of points to add for one width around each resonance peak. (the code will add points which exponentially increase in distance from the pole)")
+        self.add_param('nres_points', 100, hidden=True, comment="Number of points to add for one width around each resonance peak. (the code will add points which exponentially increase in distance from the pole)")
         self.add_param('eps_wij', 0.01, hidden=True, comment="precision of the romberg integration for Wij (irrelevant if simpson's rule used)")
         self.add_param('iter_wij', 2, hidden=True, comment="minimum number of iterations in the romberg integration algorithm for both Wij (irrelevant if simpson's rule used)")
         

=== modified file 'python_templates/maddm.inc'
--- python_templates/maddm.inc	2017-08-10 20:21:07 +0000
+++ python_templates/maddm.inc	2019-09-05 08:03:10 +0000
@@ -3,6 +3,7 @@
       double precision simpson, simpson_adaptive
 c init.f
       double precision get_Y_eq, get_gstar, get_Wij_ann, Wij_ann_integrand
+      double precision get_Wij_ann_nogrid
 c interpolate.f
       integer getindex
       double precision get_Wij_value, get_taacs_value, lineint, bessk, bessk1
@@ -80,13 +81,17 @@
       common/dmparms_logical/ solved_ODE, found_xf, dm_sm_scattering, print_out, dmsm_warning
 
 c relic density integrals
+      integer fitsize
+      parameter (fitsize=2500)
       integer nWij
       common/num_Wijs/ nWij
-      double precision betas(1000), beta_step_min, beta_step_max
+      double precision betas(fitsize), beta_step_min, beta_step_max
       common/beta_steps/ betas, beta_step_min,beta_step_max
-      double precision Wij_ann(maxdms,maxdms,1000), Wij_dm2dm(maxdms,maxdms,max_dm2dm,1000)
-      double precision Wij_scattering(maxdms,maxdms,max_dm2dm,1000), x_taacs(1000), taacs_ann(maxdms,maxdms,1000)
-      double precision taacs_dm2dm(maxdms,maxdms,max_dm2dm,1000), taacs_scattering(maxdms,maxdms,max_dm2dm,1000)
+      double precision Wij_ann(maxdms,maxdms,fitsize), Wij_dm2dm(maxdms,maxdms,max_dm2dm,fitsize)
+      double precision Wij_scattering(maxdms,maxdms,max_dm2dm,fitsize), x_taacs(fitsize)
+      double precision taacs_ann(maxdms,maxdms,fitsize)
+      double precision taacs_dm2dm(maxdms,maxdms,max_dm2dm,fitsize)
+      double precision taacs_scattering(maxdms,maxdms,max_dm2dm,fitsize)
       common/dmintegrals/ Wij_ann, Wij_dm2dm, Wij_scattering, x_taacs, taacs_ann, taacs_dm2dm, taacs_scattering
 
 c number of annihilation particles and associated diagrams