Matrix element integration over 2 final state neutrinos with MG5 standalone

Dear MG5 experts,

I generated a standalone for WW production, with the Ws decaying leptonically into first and second generations, and I'm trying to integrate over the 2 final state neutrinos momentum components.

I'm using scipy.nquad and the integration is taking forever, returning the warning "Integral is probably divergent or slowly convergent".

Do you have any tips to speed up the computations? You can find my piece of code below.

Thank you, hope this was the right place where to ask this question, even if not involving MG5 problems.

Matteo

Here's the part of code doing the integration:

         proc_id = 0 # e+ processes
         pdgf = [-11,12,13,-14]
         if e_pdg < 0: # e- processes
             proc_id = 1
             pdgf = [11,-12,-13,14]

         pdgin_list = [[1,-1],[2,-2],[3,-3],[4,-4],[5,-5]]

         def integrand (pvey,pvex,pvmz,pvez):

             # Neutrinos momenta
             ve_mom = [(pvex**2+pvey**2+pvez**2)**.5,pvex,pvey,pvez]
             pvmx = -(e_mom[1]+mu_mom[1]+pvex)
             pvmy = -(e_mom[2]+mu_mom[2]+pvey)
             vm_mom = [(pvmx**2+pvmy**2+pvmz**2)**.5,pvmx,pvmy,pvmz]

             # Initial state
             pfz = e_mom[3]+mu_mom[3]+ve_mom[3]+vm_mom[3]
             Ef = e_mom[0]+mu_mom[0]+ve_mom[0]+vm_mom[0]
             in1 = [(Ef+pfz)/2,0.,0.,(Ef+pfz)/2]
             in2 = [(Ef-pfz)/2,0.,0.,(pfz-Ef)/2]

             # Matrix element computation
             tmp = 0.

             for pin in [[in1,in2],[in2,in1]]:
                 p_all = [pin[0],pin[1],e_mom,ve_mom,mu_mom,vm_mom]
                 pinv_all = myf.invert_momenta (p_all)

                 for pdgin in pdgin_list:
                     pdgs = pdgin+pdgf
                     pdfw = pdf.xfxQ (pdgin[0],min(p_all[0][0]/ELHC,1),FactScale) * pdf.xfxQ (pdgin[1],min(p_all[1][0]/ELHC,1),FactScale)
                     me = allmatrix2py.smatrixhel (pdgs,proc_id,pinv_all,.09,0.,-1) # Sum over all helicities
                     tmp += pdfw*me

             return tmp

         # Integration limits
         pzmax = ((1.3e4-e_mom[0]-mu_mom[0])**2-(e_mom[1]+mu_mom[1])**2-(e_mom[2]+mu_mom[2])**2)**.5

         def limz (pvez):
             return [-pzmax-pvez,pzmax-pvez]

         def limx (pvmz,pvez):
             pxmax = ((1.3e4-e_mom[0]-mu_mom[0])**2-(e_mom[2]+mu_mom[2])**2-(pvez+pvmz)**2)**.5
             return [-pxmax,pxmax]

         def limy (pvex,pvmz,pvez):
             pymax = ((1.3e4-e_mom[0]-mu_mom[0])**2-(e_mom[1]+mu_mom[1])**2-(pvez+pvmz)**2)**.5
             return [-pymax,pymax]

         try:
             meint = integrate.nquad (integrand,[limy,limx,limz,[-pzmax,pzmax]])
         except:
             print (myrank,pzmax)

         delta_sm = weight*np.sign(meint[0])
         sigma_meas += delta_sm