diff --git a/main.f90 b/main.f90
index f07cca56ebb25e17eca1f7f3486c8492e6a8ce90..d73b2e28c780e2829b2fdcf16a43a991526654ab 100755
--- a/main.f90
+++ b/main.f90
@@ -88,6 +88,8 @@ integer :: NumArgs,NArg
read(arg(10:20),*) VegasNc0
elseif( arg(1:9).eq."VegasNc1=" ) then
read(arg(10:20),*) VegasNc1
+ elseif( arg(1:10).eq."Partition=" ) then
+ read(arg(11:12),*) Partition
endif
enddo
@@ -180,6 +182,9 @@ double precision :: epsrel, epsabs, flatness
call qlinit
MuRen = 100d0*GeV
+
+
+
if( NPoint.eq.3 ) then
print *, "Three-point integral:"
shat = M_H**2
@@ -300,59 +305,16 @@ endif
-! k0min = -300d0*GeV
-! k0max = +300d0*GeV
-! k0steps = 100 + 1
+k0min = +172d0*GeV + 0*GeV
+k0max = +172d0*GeV
+k0steps = 0 + 1
-! iinit=0
-! do i = 1,k0steps
-!
-! kEfix = k0min + (k0max-k0min)*(i-1)/dble(k0steps-1)
-! print *, i, "k0 FIX=",kEfix*100d0
-!
-! ncall = VegasNc0
-! itmx = VegasIt0
-! iinit=0
-! #if UseMPIVegas==1
-! call vegas_mpi(yrange(1:2*ndim),ndim,VegasIntegrand1,iinit,ncall,itmx,nprn,NUMFUNCTIONS,PDIM,WORKERS,MC_Result,MC_Error,MC_Chi2)
-! #else
-! call nvegas(yrange(1:2*ndim),ndim,VegasIntegrand1,iinit,ncall,itmx,nprn,NUMFUNCTIONS,PDIM,WORKERS,MC_Result,MC_Error,MC_Chi2)
-! #endif
-!
-!
-!
-! ncall = VegasNc1
-! itmx = VegasIt1
-! iinit=1
-! #if UseMPIVegas==1
-! call vegas_mpi(yrange(1:2*ndim),ndim,VegasIntegrand1,iinit,ncall,itmx,nprn,NUMFUNCTIONS,PDIM,WORKERS,MC_Result,MC_Error,MC_Chi2)
-! #else
-! call nvegas(yrange(1:2*ndim),ndim,VegasIntegrand1,iinit,ncall,itmx,nprn,NUMFUNCTIONS,PDIM,WORKERS,MC_Result,MC_Error,MC_Chi2)
-! #endif
-!
-!
-! write(30,*) kEFix*100d0,MC_Result/dble(k0steps),MC_Error/dble(k0steps)
-!
-! enddo
+iinit=0
+do i = 1,k0steps
-
-
-
-
-
-
-
-
-! k0steps = 400 + 1
-! k0min = 0.0000000001d0
-! k0max = 0.9999999999d0
-!
-! do i = 1,k0steps
-! x1fix = k0min + (k0max-k0min)*(i-1)/dble(k0steps-1)
+ kEfix = k0min !+ (k0max-k0min)*(i-1)/dble(k0steps-1)
+ print *, i, "k0 FIX=",kEfix*100d0
- k0steps=1
-! x1fix = 0.212313d0
-
ncall = VegasNc0
itmx = VegasIt0
iinit=0
@@ -362,6 +324,8 @@ endif
call nvegas(yrange(1:2*ndim),ndim,VegasIntegrand1,iinit,ncall,itmx,nprn,NUMFUNCTIONS,PDIM,WORKERS,MC_Result,MC_Error,MC_Chi2)
#endif
+
+
ncall = VegasNc1
itmx = VegasIt1
iinit=1
@@ -370,12 +334,53 @@ endif
#else
call nvegas(yrange(1:2*ndim),ndim,VegasIntegrand1,iinit,ncall,itmx,nprn,NUMFUNCTIONS,PDIM,WORKERS,MC_Result,MC_Error,MC_Chi2)
#endif
-
-
- write(32,"(4E24.12)") x1fix,kEfix*100d0,MC_Result/dble(k0steps),MC_Error/dble(k0steps)
+
+write(30,*) kEFix*100d0,MC_Result/dble(k0steps),MC_Error/dble(k0steps)
+
+enddo
return
-! enddo
+
+
+
+
+
+
+
+
+! ! k0steps = 400 + 1
+! ! k0min = 0.0000000001d0
+! ! k0max = 0.9999999999d0
+! !
+! ! do i = 1,k0steps
+! ! x1fix = k0min + (k0max-k0min)*(i-1)/dble(k0steps-1)
+!
+! k0steps=1
+! ! x1fix = 0.212313d0
+!
+! ncall = VegasNc0
+! itmx = VegasIt0
+! iinit=0
+! #if UseMPIVegas==1
+! call vegas_mpi(yrange(1:2*ndim),ndim,VegasIntegrand1,iinit,ncall,itmx,nprn,NUMFUNCTIONS,PDIM,WORKERS,MC_Result,MC_Error,MC_Chi2)
+! #else
+! call nvegas(yrange(1:2*ndim),ndim,VegasIntegrand1,iinit,ncall,itmx,nprn,NUMFUNCTIONS,PDIM,WORKERS,MC_Result,MC_Error,MC_Chi2)
+! #endif
+!
+! ncall = VegasNc1
+! itmx = VegasIt1
+! iinit=1
+! #if UseMPIVegas==1
+! call vegas_mpi(yrange(1:2*ndim),ndim,VegasIntegrand1,iinit,ncall,itmx,nprn,NUMFUNCTIONS,PDIM,WORKERS,MC_Result,MC_Error,MC_Chi2)
+! #else
+! call nvegas(yrange(1:2*ndim),ndim,VegasIntegrand1,iinit,ncall,itmx,nprn,NUMFUNCTIONS,PDIM,WORKERS,MC_Result,MC_Error,MC_Chi2)
+! #endif
+!
+!
+! write(32,"(4E24.12)") x1fix,kEfix*100d0,MC_Result/dble(k0steps),MC_Error/dble(k0steps)
+!
+! return
+! ! enddo
diff --git a/mod_Integrand.f90 b/mod_Integrand.f90
index 31562646a8c38641a2a151ea6c3b6edf8eec85da..92c37cdc0f09d78dd8f317a367f0a412d3946a59 100755
--- a/mod_Integrand.f90
+++ b/mod_Integrand.f90
@@ -86,31 +86,37 @@ real(8) :: test,ki(1:4,1:3)
endif
call Generate1LoopMom(yRnd(1:4),kLoopJac)
+
+
call Init1LoopEvent(kExt)
! if( LoopMomInInterior ) return ! only integrating Exterior
! if( .not. LoopMomInInterior ) return ! only integrating Interior
+
call GetContourJacobian1(ConJac)
call Eval1LoopIntegrand(LoopIntegrand)
+
Res(1) = kLoopJac * dreal(LoopIntegrand * ConJac)
-
+
+! print *, "kLoopJac",kLoopJac
+! print *, "ConJac",ConJac
+! print *, "res(1)",res(1)
+! pause
+
Res(1) = Res(1) / dreal(FinalRes)
-! write(18,*) kLoop1(1),res(1)
! ki(1:4,1) = kLoop1(1:4) - qi(1:4,1)
! ki(1:4,2) = kLoop1(1:4) - qi(1:4,2)
! ki(1:4,3) = kLoop1(1:4)!- qi(1:4,3)
-! write(19,"(10F16.8)") kloop1(4),Res(1),cdabs(ConJac),&
-! (ki(1:4,1).dot.kappa),1d0/((ki(1:4,1).dot.ki(1:4,1))-msqi(1)),&
-! (ki(1:4,2).dot.kappa),1d0/((ki(1:4,2).dot.ki(1:4,2))-msqi(2)),&
-! (ki(1:4,3).dot.kappa),1d0/((ki(1:4,3).dot.ki(1:4,3))-msqi(3))
-!
-
-
+! write(19,"(10F16.8)") kloop1(1),kloop1(4),& ! 1 2
+! Res(1),& ! 3
+! 0d0
+
+
VegasIntegrand1=0
RETURN
END FUNCTION
diff --git a/mod_LoopSpace.f90 b/mod_LoopSpace.f90
index 617ee192fa54dd2ff02491b40bdc18650e8e93b4..c2c76e8fcb3bd51c2d8cef1015a2582f2cbac817 100755
--- a/mod_LoopSpace.f90
+++ b/mod_LoopSpace.f90
@@ -14,7 +14,7 @@ public :: Eval1LoopIntegrand,Eval2LoopIntegrand,Eval1LoopDenominator_ggH
real(8),public :: kappa(1:4),kLoop1(1:4),kLoop2(1:4)
complex(8), public :: kLoop1Full(1:4)
-integer,public,parameter :: NPoint=4
+integer,public,parameter :: NPoint=3
integer,public,parameter :: NPointMAX=6
real(8),private ::mu1sq,mu2sq,E_soft,gamma1,gamma2
@@ -59,9 +59,9 @@ real(8) :: cutout
! urnd(3) = 0.0000001d0; print *, "fixing theta to zero"
! urnd(3) = 0.9999999d0; print *, "fixing theta to one"
-! uRnd(1) = 0.99999d0 ! going to UV
-! uRnd(2) = 0.5d0*( 1d0-1d0/dsqrt(2d0) ) ! this sets zeta to pi/4.
-! uRnd(3) = 1d0
+! uRnd(1) = 0.999999999d0 ! going to UV
+! uRnd(2) = 0.231d0!0.5d0*( 1d0-1d0/dsqrt(2d0) ) ! this sets zeta to pi/4.
+! uRnd(3) = 0.78d0
! uRnd(4) = 0d0
@@ -75,51 +75,58 @@ real(8) :: cutout
! endif
+
+! !Weinzierl's spherical parameterization
+! kE = dsqrt( mu_Loop1_sq*dtan( 0.5d0*Pi*uRnd(1) ) )! 0..inf
+! kE_sq=kE**2
+! zeta = dacos(1d0-2d0*uRnd(2))! 0..Pi
+! theta = dacos(1d0-2d0*uRnd(3))! 0..Pi
+! phi = 2d0*Pi*uRnd(4) ! 0..2Pi
+! call sincos(phi,sinP,cosP)
+! call sincos(theta,sinT,cosT)
+! call sincos(zeta,sinZ,cosZ)
+!
+! kLoop1(1) = kE * cosZ
+! kLoop1(2) = kE * sinZ * sinT * sinP
+! kLoop1(3) = kE * sinZ * sinT * cosP
+! kLoop1(4) = kE * sinZ * cosT
+!
+! Jacobian = 2d0*Pi**2*kE_sq/mu_Loop1_sq * ( kE_sq**2+mu_Loop1_sq**2 )*sinZ
+! !note: for zeta=1/4,3/4Pi we have sinZ^2=cosZ^2 such that kLoop1.kLoop1 = kE^2 ( sinZ^2 - cosZ^2 ) = 0, without constraint on kE
+
+! uRnd(1) = x1fix
+
+
+ !My cylindrical parameterization
+ k0 = mu_Loop1*dtan( Pi*(uRnd(1)-0.5d0) ) ! -inf..inf
- !Weinzierl's spherical parameterization
- kE = dsqrt( mu_Loop1_sq*dtan( 0.5d0*Pi*uRnd(1) ) )! 0..inf
- kE_sq=kE**2
- zeta = dacos(1d0-2d0*uRnd(2))! 0..Pi
+! ! kEfix = k0
+! k0 = kEfix
+! ! uRnd(3) = 0d0
+! uRnd(4) = 0.32d0
+
+ kE = mu_Loop1*dsqrt( dtan(0.5d0*Pi*uRnd(2)) ) ! 0..inf
theta = dacos(1d0-2d0*uRnd(3))! 0..Pi
phi = 2d0*Pi*uRnd(4) ! 0..2Pi
call sincos(phi,sinP,cosP)
call sincos(theta,sinT,cosT)
- call sincos(zeta,sinZ,cosZ)
- kLoop1(1) = kE * cosZ
- kLoop1(2) = kE * sinZ * sinT * sinP
- kLoop1(3) = kE * sinZ * sinT * cosP
- kLoop1(4) = kE * sinZ * cosT
-
- Jacobian = 2d0*Pi**2*kE_sq/mu_Loop1_sq * ( kE_sq**2+mu_Loop1_sq**2 )*sinZ
- !note: for zeta=1/4,3/4Pi we have sinZ^2=cosZ^2 such that kLoop1.kLoop1 = kE^2 ( sinZ^2 - cosZ^2 ) = 0, without constraint on kE
-
+ kLoop1(1) = k0
+ kLoop1(2) = kE * sinT * sinP
+ kLoop1(3) = kE * sinT * cosP
+ kLoop1(4) = kE * cosT
+
+ Jacobian = Pi**3 * kE*mu_Loop1**3 * ((k0/mu_Loop1)**2+1d0) * ( (kE/mu_Loop1)**4+1d0 )
+
+
+
+
+
+
+
-
-! uRnd(1) = x1fix
-
-
-
-! !My cylindrical parameterization
-! k0 = mu_Loop1*dtan( Pi*(uRnd(1)-0.5d0) ) ! -inf..inf
-!
-! ! kEfix = k0
-! ! k0 = kEfix
-!
-! kE = mu_Loop1*dsqrt( dtan(0.5d0*Pi*uRnd(2)) ) ! 0..inf
-! theta = dacos(1d0-2d0*uRnd(3))! 0..Pi
-! phi = 2d0*Pi*uRnd(4) ! 0..2Pi
-! call sincos(phi,sinP,cosP)
-! call sincos(theta,sinT,cosT)
-!
-! kLoop1(1) = k0
-! kLoop1(2) = kE * sinT * sinP
-! kLoop1(3) = kE * sinT * cosP
-! kLoop1(4) = kE * cosT
-!
-! Jacobian = Pi**3 * kE*mu_Loop1**3 * ((k0/mu_Loop1)**2+1d0) * ( (kE/mu_Loop1)**4+1d0 )
@@ -278,8 +285,6 @@ real(8) :: Q(1:4),EHat
kappa(:) = Get_kappa(kLoop1)
kLoop1Full(:) = kLoop1(:) + cI*kappa(:)
-
-
RETURN
END SUBROUTINE
@@ -537,9 +542,11 @@ SUBROUTINE Eval1LoopIntegrand(LoopIntegrand)
use ModMisc
use ModParameters
implicit none
-complex(8) :: LoopIntegrand,D(1:NPointMAX),rho(1:NPointMAX),PNorm
+complex(8) :: LoopIntegrand,D(1:NPointMAX),Dex(1:NPointMAX),rho(0:NPointMAX),PNorm
complex(8),parameter :: LoopIntNorm=1d0/cI/Pi**2
-real(8), parameter :: alph=1.2d0
+real(8), parameter :: alph=1d0
+
+ rho(0) = 1d0
if( NPoint.eq.3 ) then
D(1) = Get_MInv2(kLoop1Full-qi(1:4,1))-msqi(1)
@@ -547,14 +554,19 @@ real(8), parameter :: alph=1.2d0
D(3) = Get_MInv2(kLoop1Full-qi(1:4,3))-msqi(3)
LoopIntegrand = LoopIntNorm/(D(1)*D(2)*D(3))
-! ! partitioning factors such that rho(1)+rho(2)+rho(3)=1
-! ! PNorm =! 0 and rho(i) cancels two propagators
-! PNorm = D(1)**alph*D(2)**alph + D(1)**alph*D(3)**alph + D(2)**alph*D(3)**alph
-! rho(1) = D(2)**alph * D(3)**alph/PNorm
-! rho(2) = D(1)**alph * D(3)**alph/PNorm
-! rho(3) = D(1)**alph * D(2)**alph/PNorm
-! LoopIntegrand = LoopIntegrand * rho(3)
-
+ ! partitioning factors such that rho(1)+rho(2)+rho(3)=1
+ ! PNorm =! 0 and rho(i) cancels two propagators
+! D(1) = Get_MInv2(kLoop1-qi(1:4,1))-msqi(1)
+! D(2) = Get_MInv2(kLoop1-qi(1:4,2))-msqi(2)
+! D(3) = Get_MInv2(kLoop1-qi(1:4,3))-msqi(3)
+
+ Dex(:) = D(:)**alph
+ PNorm = Dex(1)*Dex(2) + Dex(1)*Dex(3) + Dex(2)*Dex(3)
+ rho(1) = Dex(2) * Dex(3)/PNorm
+ rho(2) = Dex(1) * Dex(3)/PNorm
+ rho(3) = Dex(1) * Dex(2)/PNorm
+ LoopIntegrand = LoopIntegrand * rho(Partition)
+
elseif( NPoint.eq.4 ) then
D(1) = Get_MInv2(kLoop1Full-qi(1:4,1))-msqi(1)
D(2) = Get_MInv2(kLoop1Full-qi(1:4,2))-msqi(2)
@@ -564,11 +576,11 @@ real(8), parameter :: alph=1.2d0
! ! partitioning factors such that rho(1)+rho(2)+rho(3)+rho(4)=1
! ! PNorm =! 0 and rho(i) cancels two propagators
-! PNorm = D(1)**alph*D(2)**alph*D(3)**alph + D(1)**alph*D(2)**alph*D(4)**alph + D(1)**alph*D(3)**alph*D(4)**alph + D(2)**alph*D(3)**alph*D(4)**alph
-! rho(1) = D(2)**alph * D(3)**alph * D(4)**alph/PNorm
-! rho(2) = D(1)**alph * D(3)**alph * D(4)**alph/PNorm
-! rho(3) = D(1)**alph * D(2)**alph * D(4)**alph/PNorm
-! rho(4) = D(1)**alph * D(2)**alph * D(3)**alph/PNorm
+! PNorm = Dex(1)*Dex(2)*Dex(3) + Dex(1)*Dex(2)*Dex(4) + Dex(1)*Dex(3)*Dex(4) + Dex(2)*Dex(3)*Dex(4)
+! rho(1) = Dex(2) * Dex(3) * Dex(4)/PNorm
+! rho(2) = Dex(1) * Dex(3) * Dex(4)/PNorm
+! rho(3) = Dex(1) * Dex(2) * Dex(4)/PNorm
+! rho(4) = Dex(1) * Dex(2) * Dex(3)/PNorm
! LoopIntegrand = LoopIntegrand * rho(1)
elseif( NPoint.eq.5 ) then
@@ -592,6 +604,22 @@ real(8), parameter :: alph=1.2d0
D(6) = Get_MInv2(kLoop1Full-qi(1:4,6))-msqi(6)
LoopIntegrand = LoopIntNorm/(D(1)*D(2)*D(3)*D(4)*D(5)*D(6))
+
+ Dex(:) = D(:)**alph
+ PNorm = Dex(1)*Dex(2)*Dex(3)*Dex(4)*Dex(5) &
+ + Dex(1)*Dex(2)*Dex(3)*Dex(4)*Dex(6) &
+ + Dex(1)*Dex(2)*Dex(3)*Dex(5)*Dex(6) &
+ + Dex(1)*Dex(2)*Dex(4)*Dex(5)*Dex(6) &
+ + Dex(1)*Dex(3)*Dex(4)*Dex(5)*Dex(6) &
+ + Dex(2)*Dex(3)*Dex(4)*Dex(5)*Dex(6)
+ rho(1) = Dex(2)*Dex(3)*Dex(4)*Dex(5)*Dex(6)/PNorm
+ rho(2) = Dex(1)*Dex(3)*Dex(4)*Dex(5)*Dex(6)/PNorm
+ rho(3) = Dex(1)*Dex(2)*Dex(4)*Dex(5)*Dex(6)/PNorm
+ rho(4) = Dex(1)*Dex(2)*Dex(3)*Dex(5)*Dex(6)/PNorm
+ rho(5) = Dex(1)*Dex(2)*Dex(3)*Dex(4)*Dex(6)/PNorm
+ rho(6) = Dex(1)*Dex(2)*Dex(3)*Dex(4)*Dex(5)/PNorm
+ LoopIntegrand = LoopIntegrand * rho(Partition)
+
! print *, "2",Get_MInv2(qi(1:4,1)), Get_MInv2(qi(1:4,2)-qi(1:4,1)), Get_MInv2(qi(1:4,3)-qi(1:4,2)), Get_MInv2(qi(1:4,4)-qi(1:4,3)), Get_MInv2(qi(1:4,5)-qi(1:4,4)), Get_MInv2(qi(1:4,5)),Get_MInv2(qi(1:4,2)),Get_MInv2(qi(1:4,3)-qi(1:4,1)),Get_MInv2(qi(1:4,4)-qi(1:4,2)),Get_MInv2(qi(1:4,5)-qi(1:4,3)),Get_MInv2(qi(1:4,4)),Get_MInv2(qi(1:4,1)-qi(1:4,5)), Get_MInv2(qi(1:4,3)),Get_MInv2(qi(1:4,4)-qi(1:4,1)),Get_MInv2(qi(1:4,5)-qi(1:4,2))
! pause
@@ -599,8 +627,6 @@ real(8), parameter :: alph=1.2d0
endif
-
-
RETURN
END SUBROUTINE
@@ -851,7 +877,7 @@ implicit none
real(8) :: Get_kappa(1:4),kTilde(:)
if( NPoint.eq.3 ) then
- Get_kappa = Get_kappa_Tri(kTilde)
+ Get_kappa = Get_kappa_Tri(kTilde)
elseif( NPoint.eq.4 ) then
Get_kappa = Get_kappa_Quad(kTilde)
elseif( NPoint.eq.5 ) then
@@ -871,6 +897,7 @@ END FUNCTION
FUNCTION Get_kappa_Tri(kTilde)
use ModMisc
+use ModParameters
implicit none
real(8) :: Get_kappa_Tri(1:4),kTilde(:),ki(1:4,1:3)
real(8) :: c_pl,c_mi
@@ -894,31 +921,57 @@ complex(8) :: D(1:3)
P_pl = get_Zpl(qi(:,1),+qi(:,1))
P_mi = get_Zmi(qi(:,2),-qi(:,2))
-
kappa_Ext(:)=0d0
kappa_Int(:)=0d0
-! damping = 1d0
-! kE_sq = VectorTimes(kTilde,kTilde)
-! damping = dexp(-((ki(1:4,1).dot.ki(1:4,1))-msqi(1))**2/1d0) + dexp(-((ki(1:4,2).dot.ki(1:4,2))-msqi(2))**2/1d0) + dexp(-((ki(1:4,3).dot.ki(1:4,3))-msqi(3))**2/1d0)
-! damping = damping + kE_sq/(kE_sq+1d7 * 0.01d0)
-
! diagnostics
LoopMomInInterior = .false.
if( InsideFWDLightcone(P_pl,kTilde,0d0) .and. InsideBWDLightcone(P_mi,kTilde,0d0) ) LoopMomInInterior = .not. LoopMomInInterior
-
+
+! !---------------------- TEST -----------------------------
+!
+! Get_kappa_Tri(1) = +2d0*ki(1,Partition)
+! Get_kappa_Tri(2) = -2d0*ki(2,Partition)
+! Get_kappa_Tri(3) = -2d0*ki(3,Partition)
+! Get_kappa_Tri(4) = -2d0*ki(4,Partition)
+!
+! D(1) = Get_MInv2(ki(1:4,1))-msqi(1)
+! D(2) = Get_MInv2(ki(1:4,2))-msqi(2)
+! D(3) = Get_MInv2(ki(1:4,3))-msqi(3)
+!
+! damp_width=50d0*kTilde(1)+10*GeV
+! ! damping = dexp(-dreal(D(Partition))**2/damp_width**2)
+! ! Get_kappa_Tri= Get_kappa_Tri * damping
+!
+! ! Get_kappa_Tri(1) = dsqrt(Get_kappa_Tri(2)**2+Get_kappa_Tri(3)**2+Get_kappa_Tri(4)**2)
+!
+!
+!
+!
+!
+! Get_kappa_Tri(1) = +2d0*ki(1,1)*dexp(-dreal(D(1))**2/damp_width**2) +2d0*ki(1,2)*dexp(-dreal(D(2))**2/damp_width**2) !+2d0*ki(1,3)*dexp(-dreal(D(3))**2/damp_width**2)
+! Get_kappa_Tri(2) = -2d0*ki(2,1)*dexp(-dreal(D(1))**2/damp_width**2) -2d0*ki(2,2)*dexp(-dreal(D(2))**2/damp_width**2) !-2d0*ki(2,3)*dexp(-dreal(D(3))**2/damp_width**2)
+! Get_kappa_Tri(3) = -2d0*ki(3,1)*dexp(-dreal(D(1))**2/damp_width**2) -2d0*ki(3,2)*dexp(-dreal(D(2))**2/damp_width**2) !-2d0*ki(3,3)*dexp(-dreal(D(3))**2/damp_width**2)
+! Get_kappa_Tri(4) = -2d0*ki(4,1)*dexp(-dreal(D(1))**2/damp_width**2) -2d0*ki(4,2)*dexp(-dreal(D(2))**2/damp_width**2) !-2d0*ki(4,3)*dexp(-dreal(D(3))**2/damp_width**2)
+!
+!
+! RETURN
+! !------------------ END TEST -----------------------------
+
+
! compute kappa_Ext
call get_hdel_pm(ki(:,1),msqi(1),mu1sq,hdel1_pl,hdel1_mi)
call get_hdel_pm(ki(:,2),msqi(2),mu1sq,hdel2_pl,hdel2_mi)
call get_hdel_pm(ki(:,3),msqi(3),mu1sq,hdel3_pl,hdel3_mi)
- c_pl = hdel1_mi * hdel2_mi * hdel3_mi
- c_mi = hdel1_pl * hdel2_pl * hdel3_pl
+ c_pl = hdel1_mi * hdel2_mi * hdel3_mi
+ c_mi = hdel1_pl * hdel2_pl * hdel3_pl
+
k_pl(1:4) = kTilde(1:4)-P_pl(1:4)
k_mi(1:4) = kTilde(1:4)-P_mi(1:4)
@@ -927,9 +980,7 @@ complex(8) :: D(1:3)
kappa_Ext(2:4) =-c_pl * k_pl(2:4) - c_mi * k_mi(2:4)
-
-
-! !compute kappa_Int
+! ! !compute kappa_Int
vij(:,ij12) = 0.5d0*( qi(:,1)+qi(:,2) - (dsqrt(msqi(1))-dsqrt(msqi(2)))/(Get_MInv(qi(1:4,1)-qi(1:4,2))+1d-14) *(qi(:,1)-qi(:,2)) )
vij(:,ij13) = 0.5d0*( qi(:,1)+qi(:,3) - (dsqrt(msqi(1))-dsqrt(msqi(3)))/(Get_MInv(qi(1:4,1)-qi(1:4,3))+1d-14) *(qi(:,1)-qi(:,3)) )
vij(:,ij23) = 0.5d0*( qi(:,2)+qi(:,3) - (dsqrt(msqi(2))-dsqrt(msqi(3)))/(Get_MInv(qi(1:4,2)-qi(1:4,3))+1d-14) *(qi(:,2)-qi(:,3)) )
@@ -949,25 +1000,22 @@ complex(8) :: D(1:3)
call get_kappa_soft(gfac,kappa_soft)
kappa_soft(:) = 0d0! removing kappa_soft for now...
- kappa_Int(:) = - c_i(1)*ki(:,1) - c_ij(ij12)*kij(:,ij12) - c_ij(ij13)*kij(:,ij13) &
- - c_i(2)*ki(:,2) - c_ij(ij23)*kij(:,ij23) &
+ kappa_Int(:) = - c_i(1)*ki(:,1) &
+ - c_i(2)*ki(:,2) &
- c_i(3)*ki(:,3) &
+ - c_ij(ij12)*kij(:,ij12) &
+ - c_ij(ij13)*kij(:,ij13) &
+ - c_ij(ij23)*kij(:,ij23) &
+ kappa_soft(:)
-! kappa_Ext = 0d0 ; call PrintOnce("removing kappaExt")
-! kappa_Int = 0d0 ; call PrintOnce("removing kappaInt")
-
+
! compute scaling function
kappa_0(:) = kappa_Int(:)+kappa_Ext(:)
kappa0_sq = Get_MInv2(kappa_0)
- ! removing kappa^2 terms by making it lightlike
-! kappa_0(1) = dsqrt( dabs( kappa_0(2)**2+kappa_0(3)**2+kappa_0(4)**2 ) )
-! call printonce("removing kappa^2 terms by making it lightlike")
-
lambdai_sq(:) = 1d0
do j=1,3
@@ -1085,7 +1133,7 @@ complex(8) :: D(1:3)
- damping = 1d0
+! damping = 1d0
! damp_width = 0.5d0
@@ -1103,7 +1151,7 @@ complex(8) :: D(1:3)
- Get_kappa_Tri(1:4) = Get_kappa_Tri(1:4) !* damping
+! Get_kappa_Tri(1:4) = Get_kappa_Tri(1:4) * damping
diff --git a/mod_Misc.f90 b/mod_Misc.f90
index 87b04346a5d90f728a5cf65c9422e0e588b59237..bdbd24506b2520d7799b94a4959ce4455125d92a 100755
--- a/mod_Misc.f90
+++ b/mod_Misc.f90
@@ -251,7 +251,7 @@ real(8) :: RefMom(1:4),TestMom(1:4)
real(8) :: Qsq,RefScale
real(8),optional :: Msq
! *******************************************
-real(8),parameter :: RelSpaceing = 1d-3 !*
+real(8),parameter :: RelSpaceing = 1d-1 !*
! *******************************************
diff --git a/mod_Parameters.f90 b/mod_Parameters.f90
index 9a0eebc71ae8e48d47235cd396ead01d2346c9b6..09255c4127c1080a1bb14b07282336c690d2ecc8 100755
--- a/mod_Parameters.f90
+++ b/mod_Parameters.f90
@@ -11,7 +11,7 @@ integer(8), public, save :: EvalCounter=0
integer(8), public, save :: SkipCounter=0
logical, public :: Seed_random,Init
-integer(8) :: ItMx,NCall,NPrn,It
+integer(8) :: ItMx,NCall,NPrn,It,Partition=0
integer(8), parameter :: MXDIM=25! has to match pvegas_mpi.h
#if compiler==1
integer, parameter :: NumSeeds=2
@@ -43,8 +43,8 @@ real(8), public,save :: dummy(1:10) = 0d0
real(8), public, parameter :: GF = (1.16639d-5)/GeV**2
-real(8), public, parameter :: m_Top = 172d0*GeV !50d0*GeV
-real(8), public, parameter :: m_H = 125d0*GeV ! 200d0*GeV
+real(8), public, parameter :: m_Top = 172d0*GeV !50d0*GeV !
+real(8), public, parameter :: m_H = 125d0*GeV ! 200d0*GeV !
real(8), public, parameter :: m_Z = 91.1876*GeV
real(8), public, parameter :: m_W = 80.399d0*GeV