diff --git a/QCDLoop-1.9/ff/ffinit_mine.f b/QCDLoop-1.9/ff/ffinit_mine.f
old mode 100644
new mode 100755
index 38c1d759f92388301c7ab70a4e185023e613c0c9..058930560ef56fa4be51b3c8309875e024b3ecb6
--- a/QCDLoop-1.9/ff/ffinit_mine.f
+++ b/QCDLoop-1.9/ff/ffinit_mine.f
@@ -621,7 +621,7 @@
+ 'ffwarn: warning: illegal value for ierr'
xlosti(i) = 0
1 continue
- call ffopen(ifile,'ffwarn.dat',ier)
+ call ffopen(ifile,'./ff/ffwarn.dat',ier)
if ( ier.ne.0 ) goto 100
rewind(ifile)
read(ifile,'(a)')warn1
@@ -783,8 +783,8 @@ C path = '/user/gj/lib/
30 continue
* second try - the system directory
C path = '/usr/local/ff/'
- path = '/afs/cern.ch/user/m/maschulz/projects
- +/TOPAZ/TTBZ/QCDLoop-1.9/ff/'
+ path='/users/pep/mschulze/projects/NumTwoLoop
+ +/Fortran/QCDLoop-1.9/'
fullname = path(1:index(path,' ')-1)//name
open(ifile,file=fullname,status='OLD',err=40)
return
diff --git a/QCDLoop-1.9/ff/makefile b/QCDLoop-1.9/ff/makefile
old mode 100644
new mode 100755
index c69731ebfe4b23db5a0a83849839c27616a0f853..4ff6f5a45417273a986aec36bb8b0efe3110a9ed
--- a/QCDLoop-1.9/ff/makefile
+++ b/QCDLoop-1.9/ff/makefile
@@ -1,4 +1,5 @@
FC = ifort -O2
+# FC = gfortran -O2
FFLAGS =
LFLAGS = $(FFLAGS)
diff --git a/QCDLoop-1.9/ql/makefile b/QCDLoop-1.9/ql/makefile
old mode 100644
new mode 100755
index ee04dff777ccd42b64f58b0fec1c1012336142cb..70fb6d976a068ba29e5a9c6b6d932d7adb4702b3
--- a/QCDLoop-1.9/ql/makefile
+++ b/QCDLoop-1.9/ql/makefile
@@ -2,6 +2,7 @@ FFLAGS = -O2
HERE = $(PWD)
LIBDIR = $(HERE)/../ff/
FC = ifort
+# FC = gfortran
LIBRARY = libqcdloop.a
diff --git a/main.f90 b/main.f90
index 83e6a78a6f2d3b71f6955946c8b4b08ac1e93d18..0a55d13ee15931be8e41f7ffd84e090c7de7a752 100755
--- a/main.f90
+++ b/main.f90
@@ -10,18 +10,19 @@
PROGRAM NumTwoLoop
use ModParameters
use ModKinematics
+#if compiler==1
use ifport
+#endif
implicit none
real(8) :: VG_Result,VG_Error,chi2
-!DEC$ IF(_UseMPIVegas .EQ.1)
+#if UseMPIVegas==1
include 'mpif.h'
integer ::ierror
call MPI_INIT(ierror)
call MPI_COMM_RANK(MPI_COMM_WORLD,MPI_Rank,ierror)
-!DEC$ ELSE
+#else
MPI_Rank=0
-!DEC$ ENDIF
-
+#endif
call GetCommandlineArgs()
call SetFileNames()
@@ -43,10 +44,9 @@ integer ::ierror
endif
call CloseFiles()
-!DEC$ IF(_UseMPIVegas .EQ.1)
+#if UseMPIVegas==1
call MPI_FINALIZE(ierror)
-!DEC$ ENDIF
-
+#endif
END PROGRAM
@@ -71,8 +71,11 @@ integer :: NumArgs,NArg
VegasIt0=-1
VegasIt1=-1
-
+#if compiler==1
NumArgs = NArgs()-1
+#elif compiler==2
+ NumArgs = COMMAND_ARGUMENT_COUNT()
+#endif
do NArg=1,NumArgs
call GetArg(NArg,arg)
if( arg(1:9).eq."VegasIt0=" ) then
@@ -100,7 +103,6 @@ END SUBROUTINE
SUBROUTINE SetFileNames()
use ModParameters
-use ifport
implicit none
@@ -135,9 +137,9 @@ use ModParameters
use modLoopSpace
implicit none
real(8) :: VG_Result,VG_Error,VG_Chi2
-!DEC$ IF(_UseMPIVegas .EQ.1)
+#if UseMPIVegas==1
include 'mpif.h'
-!DEC$ ENDIF
+#endif
integer i,init,NDim
real(8) :: yrange(1:2*MXDIM)
complex(8) :: qlI3,qlI4,Int
@@ -151,10 +153,9 @@ real(8) :: shat,that,uhat
yrange(i+ndim)=1d0
enddo
-!DEC$ IF(_UseMPIVegas .EQ.1)
+#if UseMPIVegas==1
call ClearRedHisto()
-!DEC$ ENDIF
-
+#endif
call qlinit
MuRen = 100d0*GeV
@@ -165,7 +166,6 @@ if( NPoint.eq.3 ) then
Int = qlI3(shat,0d0,0d0,m_top**2,m_top**2,m_top**2,MuRen**2,-1); print *, "QCD Lib: C0 1/ep^1",Int
Int = qlI3(shat,0d0,0d0,m_top**2,m_top**2,m_top**2,MuRen**2,+0); print *, "QCD Lib: C0 1/ep^0",Int
elseif( NPoint.eq.4 ) then
-
print *, "Four-point integral:"
shat = +429.839225955798d0
that = -121.141609980866d0
@@ -178,22 +178,21 @@ endif
ncall = VegasNc0
itmx = VegasIt0
init=0
-!DEC$ IF(_UseMPIVegas .EQ.1)
+#if UseMPIVegas==1
call vegas_mpi(yrange(1:2*ndim),ndim,VegasIntegrand1,init,ncall,itmx,nprn,NUMFUNCTIONS,PDIM,WORKERS,VG_Result,VG_Error,VG_Chi2)
-!DEC$ ELSE
+#else
call vegas(yrange(1:2*ndim),ndim,VegasIntegrand1,init,ncall,itmx,nprn,NUMFUNCTIONS,PDIM,WORKERS,VG_Result,VG_Error,VG_Chi2)
-!DEC$ ENDIF
-
+#endif
ncall = VegasNc1
itmx = VegasIt1
init=1
-!DEC$ IF(_UseMPIVegas .EQ.1)
+#if UseMPIVegas==1
call vegas_mpi(yrange(1:2*ndim),ndim,VegasIntegrand1,init,ncall,itmx,nprn,NUMFUNCTIONS,PDIM,WORKERS,VG_Result,VG_Error,VG_Chi2)
-!DEC$ ELSE
+#else
call vegas(yrange(1:2*ndim),ndim,VegasIntegrand1,init,ncall,itmx,nprn,NUMFUNCTIONS,PDIM,WORKERS,VG_Result,VG_Error,VG_Chi2)
-!DEC$ ENDIF
+#endif
@@ -355,10 +354,11 @@ integer, dimension(:), allocatable :: gen_seed
integer :: n,i,sclock,SeedSize
+
if( seed_random ) then
call random_seed()
call random_seed(size=SeedSize)
- if(SeedSize.ne.2) call Error("ifort SeedSize has changed from 2 to ",SeedSize)
+ if(SeedSize.ne.2) call Error("SeedSize has changed from 2 to ",SeedSize)
Seeds(0) = SeedSize
call random_seed(get=Seeds(1:SeedSize))
else
diff --git a/makefile b/makefile
index 43a1779c34fb3f4ec71f45f5415a35317efee0e0..7776a26957b25f5b6e12d1ed9b0a11945a5c2cb0 100755
--- a/makefile
+++ b/makefile
@@ -20,21 +20,27 @@ useMPI = No
ifeq ($(useMPI),Yes)
Exec = ./NumInt_MPI
- F95compiler = mpif90-ifort -lpthread -D_UseMPIVegas=1
+ F95compiler = mpif90-ifort -lpthread -DUseMPIVegas=1
ccomp = mpicc -lpthread -lm
else
Exec = ./NumInt
- F95compiler = ifort -D_UseMPIVegas=0
+ F95compiler = ifort -fpp -DUseMPIVegas=0 -Dcompiler=1
ccomp = icc -O2
+
+# F95compiler = gfortran -cpp -DUseMPIVegas=0 -Dcompiler=2
+# ccomp = gcc -O3
endif
ifeq ($(Opt),Yes)
- IfortOpts = -O2 -fpp -I$(Here)/colors -I$(VegasDir) -module $(ModuleDir)
+ IfortOpts = -O2 -I$(VegasDir) -module $(ModuleDir)
+
+# IfortOpts = -O3 -I$(VegasDir) -I$(ModuleDir) -J$(ModuleDir) -ffree-line-length-none
+
else
- IfortOpts = -O0 -fpp -implicitnone -check bounds -check pointer -warn interfaces -ftrapuv -debug extended -g -traceback -fpe0 -check uninit -I$(Here)/colors -I$(VegasDir) -module $(ModuleDir)
+ IfortOpts = -O0 -implicitnone -check bounds -check pointer -warn interfaces -ftrapuv -debug extended -g -traceback -fpe0 -check uninit -I$(VegasDir) -module $(ModuleDir)
endif
fcomp = $(F95compiler) $(IfortOpts) @$(ConfigFile)
diff --git a/mod_Integrand.f90 b/mod_Integrand.f90
old mode 100644
new mode 100755
index d5816242ef4881542534a054e8269de2a9b1d9af..2ab754259df11e522872ca1f3ee93f038cd96536
--- a/mod_Integrand.f90
+++ b/mod_Integrand.f90
@@ -20,16 +20,18 @@ real(8) :: kLoopJac
complex(8) :: LoopIntegrand,ConJac
-! EHat=m_H
-! call EvalPhasespace_2to1(EHat,kExt,PSWgt)
-
- EHat = 1.9d0*( M_H+M_Z )
- call EvalPhasespace_2to2(EHat,(/M_H,M_Z/),(/0.212d0,0.632d0/),kExt,PSWgt)
-
- call GenerateLoopMom(yRnd(1:4),kLoopJac)
- call InitLoopEvent(kExt)
- call GetContourJacobian(ConJac)
- call EvalLoopIntegrand(LoopIntegrand)
+ if( NPoint.eq.3 ) then
+ EHat=m_H
+ call EvalPhasespace_2to1(EHat,kExt,PSWgt)
+ elseif( NPoint.eq.4 ) then
+ EHat = 1.9d0*( M_H+M_Z )
+ call EvalPhasespace_2to2(EHat,(/M_H,M_Z/),(/0.212d0,0.632d0/),kExt,PSWgt)
+ endif
+
+ call Generate1LoopMom(yRnd(1:4),kLoopJac)
+ call Init1LoopEvent(kExt)
+ call GetContourJacobian1(ConJac)
+ call Eval1LoopIntegrand(LoopIntegrand)
Res(1) = kLoopJac * dreal(LoopIntegrand * ConJac)
@@ -43,6 +45,36 @@ END FUNCTION
+FUNCTION VegasIntegrand2(yRnd,VgsWgt,Res)
+use ModParameters
+use ModMisc
+use ModKinematics
+use ModLoopSpace
+implicit none
+integer :: VegasIntegrand2
+real(8) :: yRnd(*),Res(*),VgsWgt
+real(8) :: Ehat,kExt(1:4,1:4),PSWgt
+real(8) :: kLoopJac
+complex(8) :: LoopIntegrand,ConJac
+
+
+ EHat=m_H
+
+ call Generate2LoopMom(yRnd(1:8),kLoopJac)
+ call Init2LoopEvent(kExt)
+ call GetContourJacobian2(ConJac)
+ call Eval2LoopIntegrand(LoopIntegrand)
+
+
+ Res(1) = kLoopJac * dreal(LoopIntegrand * ConJac)
+
+
+ VegasIntegrand2=0
+RETURN
+END FUNCTION
+
+
+
END MODULE ModIntegrand
diff --git a/mod_LoopSpace.f90 b/mod_LoopSpace.f90
index 15c0bc5ba01caff3a38e3e0c64f65188cf98cac8..3315d81915f8db2a585294a165f8faea57b52203 100755
--- a/mod_LoopSpace.f90
+++ b/mod_LoopSpace.f90
@@ -6,12 +6,16 @@ implicit none
private
+public :: Generate1LoopMom,Generate2LoopMom
+public :: Init1LoopEvent,Init2LoopEvent
+public :: GetContourJacobian1,GetContourJacobian2
+public :: Eval1LoopIntegrand,Eval2LoopIntegrand
-public :: InitLoopEvent,GetContourJacobian,GenerateLoopMom,EvalLoopIntegrand
-real(8),public :: kappa(1:4),kLoop1(1:4)
+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
+
real(8),private :: msqi(1:NPoint),mu1sq,mu2sq,E_soft,gamma1,gamma2
real(8),private :: qi(1:4,1:NPoint)
real(8),private :: kbar(1:4),ImMuUVsq
@@ -22,7 +26,7 @@ real(8),private :: kbar(1:4),ImMuUVsq
-SUBROUTINE GenerateLoopMom(uRnd,Jacobian)
+SUBROUTINE Generate1LoopMom(uRnd,Jacobian)
use ModParameters
use modmisc
implicit none
@@ -165,8 +169,60 @@ END SUBROUTINE
+SUBROUTINE Generate2LoopMom(uRnd,Jacobian)
+use ModParameters
+use modmisc
+implicit none
+real(8) :: uRnd(:),Jacobian
+real(8) :: kE,kE_sq,zeta,theta,phi
+real(8),parameter :: mu_Loop1_sq = (100d0*GeV)**2
+real(8),parameter :: mu_Loop2_sq = (100d0*GeV)**2
+real(8) :: sinZ,cosZ,sinT,cosT,sinP,cosP
+
+
+
+ kE = dsqrt( mu_Loop1_sq*dtan( 0.5d0*Pi*uRnd(1) ) )
+ kE_sq=kE**2
+ zeta = dacos(1d0-2d0*uRnd(2))
+ theta = dacos(1d0-2d0*uRnd(3))
+ phi = 2d0*Pi*uRnd(4)
+ 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
+
+
-SUBROUTINE InitLoopEvent(kExt)
+ kE = dsqrt( mu_Loop2_sq*dtan( 0.5d0*Pi*uRnd(5) ) )
+ kE_sq=kE**2
+ zeta = dacos(1d0-2d0*uRnd(6))
+ theta = dacos(1d0-2d0*uRnd(7))
+ phi = 2d0*Pi*uRnd(8)
+ call sincos(phi,sinP,cosP)
+ call sincos(theta,sinT,cosT)
+ call sincos(zeta,sinZ,cosZ)
+
+ kLoop2(1) = kE * cosZ
+ kLoop2(2) = kE * sinZ * sinT * sinP
+ kLoop2(3) = kE * sinZ * sinT * cosP
+ kLoop2(4) = kE * sinZ * cosT
+
+ Jacobian = Jacobian * 2d0*Pi**2*kE_sq/mu_Loop2_sq * ( kE_sq**2+mu_Loop2_sq**2 )*sinZ
+
+
+RETURN
+END SUBROUTINE
+
+
+
+SUBROUTINE Init1LoopEvent(kExt)
use ModMisc
use ModParameters
implicit none
@@ -224,8 +280,46 @@ END SUBROUTINE
+SUBROUTINE Init2LoopEvent(kExt)
+use ModMisc
+use ModParameters
+implicit none
+real(8) :: kExt(:,:)
+real(8) :: Q(1:4),EHat
+
+
+
+ EHat = 2d0*kExt(1,1)
+ mu1sq = (0.035d0*EHat)**2
+ mu2sq = (0.7d0*EHat)**2
+ E_soft = 0.03d0*EHat
+ gamma1 = 0.7d0
+ gamma2 = 0.008d0
+
+ Q(1:4) = (/1000d0,1000d0,1000d0,1000d0/)*GeV
+ ImMuUVsq = -2500d0*GeV
+ kbar(:) = kLoop1(:) - Q(:) ! defining kbar real i.e. with k_tilde and not k
+
+
+ qi(1:4,1) = -kExt(1:4,1)
+ qi(1:4,2) = -kExt(1:4,1)+kExt(1:4,3)
+ qi(1:4,3) = 0d0
+
+ msqi(1) = m_Top**2
+ msqi(2) = m_Top**2
+
+
+ kappa(:) = Get_kappa(kLoop1)
+ kLoop1Full(:) = kLoop1(:) + cI*kappa(:)
+
+
+RETURN
+END SUBROUTINE
+
+
+
-SUBROUTINE GetContourJacobian(Jacobian)
+SUBROUTINE GetContourJacobian1(Jacobian)
!
! J_ij = delta_ij + cI * d(kappa^i)/d(kLoop^j)
!
@@ -256,6 +350,35 @@ END SUBROUTINE
+
+SUBROUTINE GetContourJacobian2(Jacobian)
+!
+! J_ij = delta_ij + cI * d(kappa^i)/d(kLoop^j)
+!
+use modMisc
+use ModParameters
+implicit none
+complex(8) :: Jacobian !,JacobiMatrixC(1:8,1:8)
+real(8) :: JacobiMatrix(1:8,1:8)
+
+
+
+! call GetJacobiMatrix_kappa(kLoop1,kLoop2,JacobiMatrix)
+! Jacobian = GetDet8x8OnePlusID(JacobiMatrix)
+
+
+! write(*,"(4d14.3)") JacobiMatrix(1,1:4)
+! write(*,"(4d14.3)") JacobiMatrix(2,1:4)
+! write(*,"(4d14.3)") JacobiMatrix(3,1:4)
+! write(*,"(4d14.3)") JacobiMatrix(4,1:4)
+! print *, Jacobian
+! pause
+
+RETURN
+END SUBROUTINE
+
+
+
SUBROUTINE GetJacobiMatrix_kappa(kLoop,der) ! this is d(kappa)/d(kLoop) at kLoop
@@ -350,7 +473,7 @@ END SUBROUTINE
-SUBROUTINE EvalLoopIntegrand(LoopIntegrand)
+SUBROUTINE Eval1LoopIntegrand(LoopIntegrand)
use ModMisc
use ModParameters
implicit none
@@ -376,6 +499,31 @@ END SUBROUTINE
+
+
+SUBROUTINE Eval2LoopIntegrand(LoopIntegrand)
+use ModMisc
+use ModParameters
+implicit none
+complex(8) :: LoopIntegrand,D(1:5)
+complex(8),parameter :: LoopIntNorm=1d0/(4d0*Pi)**4
+
+
+! D(1) = Get_MInv2(kLoop1Full-qi(1:4,1))-msqi(1)
+! D(2) = Get_MInv2(kLoop1Full-qi(1:4,2))-msqi(2)
+! D(3) = Get_MInv2(kLoop1Full-qi(1:4,3))-msqi(3)
+! D(4) = Get_MInv2(kLoop1Full-qi(1:4,4))-msqi(4)
+! D(5) = Get_MInv2(kLoop1Full-qi(1:4,5))-msqi(5)
+! LoopIntegrand = LoopIntNorm/(D(1)*D(2)*D(3)*D(4)*D(5))
+
+
+RETURN
+END SUBROUTINE
+
+
+
+
+
SUBROUTINE get_hdel_pm(k,msq,M1sq,hdel_pl,hdel_mi)
implicit none
real(8),intent(in) :: k(:),msq,M1sq
diff --git a/mod_Misc.f90 b/mod_Misc.f90
index 4a6d741781917b8d923575f8f26b7cc1c1460184..63c30a3edc81cbc1e25e46873a23b22ae2b061a0 100755
--- a/mod_Misc.f90
+++ b/mod_Misc.f90
@@ -166,6 +166,22 @@ END FUNCTION
+FUNCTION GetDet8x8OnePlusID(d) ! determinant of 1 + i*d(i,j)
+implicit none
+real(8) :: d(:,:)
+complex(8) :: GetDet8x8OnePlusID
+complex(8),parameter :: cI=(0d0,1d0)
+
+ GetDet8x8OnePlusID = 0d0
+
+
+RETURN
+END FUNCTION
+
+
+
+
+
FUNCTION InsideFWDLightcone(RefMom,TestMom,Msq)
implicit none
@@ -1382,7 +1398,7 @@ complex(16) :: mult, tsum
work(1:N,1:N+1) = matrix(1:N,1:N+1)
do i = 1,N
- if (cqabs(work(i,i)) <= 1.0e-6) then
+ if (abs(work(i,i)) <= 1.0e-6) then
! print *, i,N
! print *, "go_Gauss, Zero pivot element: ",cdabs(work(i,i))
! call Error("go_Gauss, Zero pivot element")
diff --git a/mod_Parameters.f90 b/mod_Parameters.f90
index a3a8e31be09d1a2a8cc2c689690284a7aa24ffff..6fa09af38114bf2a2561e2e13f96f09fbcdfa700 100755
--- a/mod_Parameters.f90
+++ b/mod_Parameters.f90
@@ -13,7 +13,12 @@ integer(8), public, save :: SkipCounter=0
logical, public :: Seed_random,WarmUp
integer(8) :: ItMx,NCall,NPrn,It
integer(8), parameter :: MXDIM=25! has to match pvegas_mpi.h
-integer, public :: TheSeeds(0:2) = (/2,700470849,476/)! only used if seed_random=.false., the first entry is the total number of seeds
+#if compiler==1
+integer(8), parameter :: NumSeeds=2
+#else
+integer(8), parameter :: NumSeeds=12
+#endif
+integer, public :: TheSeeds(0:12) = (/NumSeeds,700470849,476,123123,756,21434,6754654,12313,54353,12313,65464,12313,65464/)! only used if seed_random=.false., the first entry is the total number of seeds
real(8), public :: IntegrandMax=-1d30
! PVegas (MPI) part