From c60c25747f1c411e4c239c873503563346020230 Mon Sep 17 00:00:00 2001
From: Markus Schulze <markus.schulze@physik.hu-berlin.de>
Date: Fri, 26 Feb 2021 16:14:48 +0100
Subject: [PATCH] some cleanup
---
main.f90 | 107 --------------
makefile | 0
mod_Amplitudes.f90 | 0
mod_Integrand.f90 | 18 ---
mod_Kinematics.f90 | 0
mod_LoopSpace.f90 | 341 +--------------------------------------------
mod_Misc.f90 | 0
mod_Parameters.f90 | 0
8 files changed, 3 insertions(+), 463 deletions(-)
mode change 100755 => 100644 main.f90
mode change 100755 => 100644 makefile
mode change 100755 => 100644 mod_Amplitudes.f90
mode change 100755 => 100644 mod_Integrand.f90
mode change 100755 => 100644 mod_Kinematics.f90
mode change 100755 => 100644 mod_LoopSpace.f90
mode change 100755 => 100644 mod_Misc.f90
mode change 100755 => 100644 mod_Parameters.f90
diff --git a/main.f90 b/main.f90
old mode 100755
new mode 100644
index d73b2e2..8c57218
--- a/main.f90
+++ b/main.f90
@@ -243,12 +243,6 @@ elseif( NPoint.eq.6 ) then
print *, "Six-point integral:"
call EvalPhasespace_2to4(1.9d0*( 2*M_H+2*M_Z ),(/M_H,M_H,M_Z,M_Z/),(/0.212d0,0.632d0,0.275d0,0.820d0,0.123d0,0.452d0,0.821d0,0.231d0/),kExt,PSWgt)
-! shat = Get_MInv2(-kExt(1:4,1)-kExt(1:4,2))
-! uhat = Get_MInv2(-kExt(1:4,1)+kExt(1:4,3))
-! vhat = Get_MInv2( kExt(1:4,3)+kExt(1:4,4))
-! what = Get_MInv2( kExt(1:4,4)+kExt(1:4,5))
-! zhat = Get_MInv2(-kExt(1:4,1)+kExt(1:4,3)+kExt(1:4,4))
-
#if useCollier==1
Nmax = 6
@@ -348,62 +342,6 @@ return
-! ! 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
-
-
-
-
-
-
-
-! ncall = VegasNc0*VegasIt0+VegasNc1*VegasIt1*10
-! iinit=0
-! call DCUHRE(ndim,NUMFUNCTIONS,yrange(1:4),yrange(5:8),ncall/10,ncall,VegasIntegrand1_DCUHRE,0d0,1d-6,0,NW,iinit,MC_result,MC_error,NEVAL,IFAIL,WORK)
-!
-!
-! print *, "cuhre ",i, MC_result,MC_error,NEVAL,IFAIL
-
-
-
-
-
-
-
-
-
-
@@ -473,51 +411,6 @@ use modMisc
p(i,3) = MomVec(3,i)
end do
-! ! use momentum set from Racoon
-! p(1,0) = -0.9500000000000000D+02 *000000d0
-! p(1,1) = 0.0000000000000000D+00 *000000d0
-! p(1,2) = 0.0000000000000000D+00 *000000d0
-! p(1,3) = 0.9500000000000000D+02 *000000d0
-!
-! p(2,0) = -0.9500000000000000D+02 *000000d0
-! p(2,1) = 0.0000000000000000D+00 *000000d0
-! p(2,2) = 0.0000000000000000D+00 *000000d0
-! p(2,3) = -0.9500000000000000D+02 *000000d0
-!
-! p(3,0) = 0.2046111001757171d+02 *000000d0
-! p(3,1) = 0.1057734233089455D+02 *000000d0
-! p(3,2) = -0.2324961261504543D+01 *000000d0
-! p(3,3) = 0.1736005205921753D+02 *000000d0
-!
-! p(4,0) = 0.3558305163378869D+01 *000000d0
-! p(4,1) = 0.1436222934374051d+01 *000000d0
-! p(4,2) = -0.2174258125294355D+01 *000000d0
-! p(4,3) = -0.2423097382091398D+01 *000000d0
-!
-! p(5,0) = 0.8154540918019539D+02 *000000d0
-! p(5,1) = -0.5230395944682889D+02 *000000d0
-! p(5,2) = 0.3083642435466509D+02 *000000d0
-! p(5,3) = 0.5443403822581044D+02 *000000d0
-!
-! p(6,0) = 0.8443517563885433D+02 *000000d0
-! p(6,1) = 0.4029039418156027D+02 *000000d0
-! p(6,2) = -0.2633720496786619D+02 *000000d0
-! p(6,3) = -0.6937099290293661d+02 *000000d0
-
-
-! do i=N+1,6
-! p(N,0) = p(N,0) + p(i,0)
-! p(N,1) = p(N,1) + p(i,1)
-! p(N,2) = p(N,2) + p(i,2)
-! p(N,3) = p(N,3) + p(i,3)
-! end do
-
-! if (present(MomVec)) then
-! MomVec(0:3,1) = p(1,0:3)
-! do i=2,N-1
-! MomVec(0:3,i) = MomVec(0:3,i-1)+p(i,0:3)
-! end do
-! end if
! define invariants
do k=0,N-1
diff --git a/makefile b/makefile
old mode 100755
new mode 100644
diff --git a/mod_Amplitudes.f90 b/mod_Amplitudes.f90
old mode 100755
new mode 100644
diff --git a/mod_Integrand.f90 b/mod_Integrand.f90
old mode 100755
new mode 100644
index 92c37cd..71d56d8
--- a/mod_Integrand.f90
+++ b/mod_Integrand.f90
@@ -90,32 +90,14 @@ real(8) :: test,ki(1:4,1:3)
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)
-
-! 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(1),kloop1(4),& ! 1 2
-! Res(1),& ! 3
-! 0d0
-
VegasIntegrand1=0
RETURN
diff --git a/mod_Kinematics.f90 b/mod_Kinematics.f90
old mode 100755
new mode 100644
diff --git a/mod_LoopSpace.f90 b/mod_LoopSpace.f90
old mode 100755
new mode 100644
index c2c76e8..d319e32
--- a/mod_LoopSpace.f90
+++ b/mod_LoopSpace.f90
@@ -53,28 +53,6 @@ real(8) :: sinZ,cosZ,sinT,cosT,sinP,cosP
real(8) :: cutout
-
-! urnd(1:4)=(/0.9999999999d0, 0.146446609406726d0, 0.37581005831787d0, 0.422403883500291d0 /)
-! urnd(2) = 0.5d0; print *, "fixing zeta to zero"
-! urnd(3) = 0.0000001d0; print *, "fixing theta to zero"
-! urnd(3) = 0.9999999d0; print *, "fixing theta to one"
-
-! 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
-
-
-! call PrintOnce("Fixing yRnd")
-! uRnd(2) = 0.5d0 ! cosZ=0d0 kLoop(1) =0
-! ! uRnd(3) = 1d-16 ! cosT = -1d0
-! if( uRnd(3).lt.0.5d0 ) then
-! uRnd(3)=1d-15
-! else
-! uRnd(3)=1d0
-! endif
-
-
! !Weinzierl's spherical parameterization
! kE = dsqrt( mu_Loop1_sq*dtan( 0.5d0*Pi*uRnd(1) ) )! 0..inf
@@ -100,13 +78,7 @@ real(8) :: cutout
!My cylindrical parameterization
- k0 = mu_Loop1*dtan( Pi*(uRnd(1)-0.5d0) ) ! -inf..inf
-
-! ! kEfix = k0
-! k0 = kEfix
-! ! uRnd(3) = 0d0
-! uRnd(4) = 0.32d0
-
+ k0 = mu_Loop1*dtan( Pi*(uRnd(1)-0.5d0) ) ! -inf..inf
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
@@ -120,42 +92,6 @@ real(8) :: cutout
Jacobian = Pi**3 * kE*mu_Loop1**3 * ((k0/mu_Loop1)**2+1d0) * ( (kE/mu_Loop1)**4+1d0 )
-
-
-
-
-
-
-
-
-
-
-
-! kloop1(1) = 200d0*GeV; call printonce("fixing kE")
-
-! cutout = m_Top*1d-2 ! cutting out a point in loop space
-! if( dabs(kLoop1(1)-(+000d0*GeV)).lt.cutout*GeV .and. &
-! dabs(kLoop1(2)-(+000d0*GeV)).lt.cutout*GeV .and. &
-! dabs(kLoop1(3)-(+000d0*GeV)).lt.cutout*GeV .and. &
-! dabs(kLoop1(4)-(-100d0*GeV)).lt.cutout*GeV ) Jacobian = 0d0
-
-! dummy(1) = dummy(1) + 0.5d0*GeV
-! print *, "Warning: fixing kloop"
-! kLoop1(1) = 0d0
-! kLoop1(2) = 0d0
-! kLoop1(3) = 0d0
-! kLoop1(4) = -0.133975d0
-! kLoop1(4) = -100d0*GeV - 86.603d0*GeV
-! kLoop1(4) = -250d0*GeV + dummy(1)
-
-
-! print *, "ke",dsqrt(kE_sq)
-! print *, "k.k",(kLoop1.dot.kLoop1)
-! pause
-
-
-
-
RETURN
END SUBROUTINE
@@ -349,13 +285,6 @@ real(8) :: JacobiMatrix(1:4,1:4)
Jacobian = GetDet4x4OnePlusID(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
@@ -377,14 +306,6 @@ 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
@@ -417,25 +338,6 @@ logical :: okay=.false.
der(1:4,j) = 0.5d0/hh*dkap(1:4)
-
-! canc(1:4) = dabs(dkap(1:4)/dmax1(fmi(1:4),fpl(1:4)))
-! if( dmax1( canc(1),canc(2),canc(3),canc(4) ) .gt.1d-2 ) then
-! print *, "warning"
-! print *, canc(1:4)
-! endif
-
-
-! do i=1,4
-! if( dabs(der(i,j)).gt.100d0 ) then
-! print *, "xx",der(i,j)
-! fpl(1:4) = Get_kappa(xpl)
-! fmi(1:4) = Get_kappa(xmi)
-! print *, dkap(i),fpl(i),fmi(i)
-! print *, kLoop(j),hh
-! pause
-! endif
-! enddo
-
enddo
@@ -498,43 +400,6 @@ END SUBROUTINE
-
-
-! SUBROUTINE GetJacobiMatrix(kLoop,res)
-! use ModParameters
-! implicit none
-! complex(8) :: res(:,:)
-! real(8) :: fpl(1:4),fmi(1:4)
-! real(8) :: kLoop(:),xpl(1:4),xmi(1:4),dkap(1:4)
-! integer(kind=1) :: j
-! real(8),parameter :: h=1d-10
-!
-!
-!
-! do j=1,4
-!
-! xpl(1:4) = kLoop(:)
-! xmi(1:4) = kLoop(:)
-! xpl(j) = xpl(j)+h
-! xmi(j) = xpl(j)-h
-!
-! fpl(1:4) = Get_kappa(xpl)
-! fmi(1:4) = Get_kappa(xmi)
-! dkap(1:4) = fpl(1:4)-fmi(1:4)
-! res(1:4,j) = 0.5d0/h*dkap(1:4) * cI
-!
-! res(j,j) = 1d0 + res(j,j)! adding the d(kLoop)/d(kTilde) contribution
-!
-! if( dabs(1d0-fpl(1)/fmi(1)).lt.1d-14 ) print *, "Problem in GetJacobiMatrix 1,",j
-! if( dabs(1d0-fpl(2)/fmi(2)).lt.1d-14 ) print *, "Problem in GetJacobiMatrix 2,",j
-! if( dabs(1d0-fpl(3)/fmi(3)).lt.1d-14 ) print *, "Problem in GetJacobiMatrix 3,",j
-! if( dabs(1d0-fpl(4)/fmi(4)).lt.1d-14 ) print *, "Problem in GetJacobiMatrix 4,",j
-!
-! enddo
-!
-! END SUBROUTINE
-
-
@@ -590,10 +455,7 @@ real(8), parameter :: alph=1d0
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))
-
-
-! 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,4)),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,3)),Get_MInv2(qi(1:4,1)-qi(1:4,4))
-! pause
+
elseif( NPoint.eq.6 ) then
D(1) = Get_MInv2(kLoop1Full-qi(1:4,1))-msqi(1)
@@ -619,10 +481,7 @@ real(8), parameter :: alph=1d0
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
+
endif
@@ -932,38 +791,6 @@ complex(8) :: D(1:3)
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)
@@ -1048,169 +875,7 @@ complex(8) :: D(1:3)
-! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !
-!
-!
-! j=1
-! kappa_0(1) = SmoothSignFunc(+ki(1,j),1d0/dabs(kTilde(1)) )
-! kappa_0(2) = SmoothSignFunc(-ki(2,j),1d0/dabs(kTilde(1)) )/dsqrt(3d0)
-! kappa_0(3) = SmoothSignFunc(-ki(3,j),1d0/dabs(kTilde(1)) )/dsqrt(3d0)
-! kappa_0(4) = SmoothSignFunc(-ki(4,j),1d0/dabs(kTilde(1)) )/dsqrt(3d0)
-!
-! kappa_0(2:4) = -2d0*ki(2:4,j)
-! kappa_0(1) = SmoothSignFunc(+ki(1,j),1d0/dabs(kTilde(1)))*dsqrt(kappa_0(2)**2+kappa_0(3)**2+kappa_0(4)**2 )
-!
-!
-! k_pl(1:4) = ki(1:4,j)
-! k_pl(1) = k_pl(1) +30*GeV*(1d0+kTilde(1)/50d0*GeV)
-! D(j) = Get_MInv2(k_pl)-msqi(j)
-! damping = SmoothStepFunc( +dreal(D(j)) )
-! damping2= SmoothStepFunc( -dreal(D(j)) )
-!
-! k_pl(1:4) = ki(1:4,j)
-! k_pl(1) = k_pl(1) -30*GeV*(1d0+kTilde(1)/50d0*GeV)
-! D(j) = Get_MInv2(k_pl)-msqi(j)
-! damping = damping * SmoothStepFunc(-dreal(D(j)) )
-! damping2= damping2* SmoothStepFunc(+dreal(D(j)) )
-!
-! kappa_0(:) = kappa_0(:) * (damping+damping2)
-! Get_kappa_Tri(1:4) = kappa_0(:)
-!
-!
-!
-! j=2
-! kappa_0(1) = SmoothSignFunc(+ki(1,j),1d0/dabs(kTilde(1)) )
-! kappa_0(2) = SmoothSignFunc(-ki(2,j),1d0/dabs(kTilde(1)) )/dsqrt(3d0)
-! kappa_0(3) = SmoothSignFunc(-ki(3,j),1d0/dabs(kTilde(1)) )/dsqrt(3d0)
-! kappa_0(4) = SmoothSignFunc(-ki(4,j),1d0/dabs(kTilde(1)) )/dsqrt(3d0)
-!
-! kappa_0(2:4) = -2d0*ki(2:4,j)
-! kappa_0(1) = SmoothSignFunc(+ki(1,j),1d0/dabs(kTilde(1)))*dsqrt(kappa_0(2)**2+kappa_0(3)**2+kappa_0(4)**2 )
-!
-! k_pl(1:4) = ki(1:4,j)
-! k_pl(1) = k_pl(1) +30*GeV*(1d0+kTilde(1)/50d0*GeV)
-! D(j) = Get_MInv2(k_pl)-msqi(j)
-! damping = SmoothStepFunc( +dreal(D(j)) )
-! damping2= SmoothStepFunc( -dreal(D(j)) )
-!
-! k_pl(1:4) = ki(1:4,j)
-! k_pl(1) = k_pl(1) -30*GeV*(1d0+kTilde(1)/50d0*GeV)
-! D(j) = Get_MInv2(k_pl)-msqi(j)
-! damping = damping * SmoothStepFunc(-dreal(D(j)) )
-! damping2= damping2* SmoothStepFunc(+dreal(D(j)) )
-!
-! kappa_0(:) = kappa_0(:) * (damping+damping2)
-! Get_kappa_Tri(1:4) = Get_kappa_Tri(1:4) + kappa_0(:)
-!
-!
-!
-! j=3
-! kappa_0(1) = SmoothSignFunc(+ki(1,j),1d0/dabs(kTilde(1)) )
-! kappa_0(2) = SmoothSignFunc(-ki(2,j),1d0/dabs(kTilde(1)) )/dsqrt(3d0)
-! kappa_0(3) = SmoothSignFunc(-ki(3,j),1d0/dabs(kTilde(1)) )/dsqrt(3d0)
-! kappa_0(4) = SmoothSignFunc(-ki(4,j),1d0/dabs(kTilde(1)) )/dsqrt(3d0)
-!
-! kappa_0(2:4) = -2d0*ki(2:4,j)
-! kappa_0(1) = SmoothSignFunc(+ki(1,j),1d0/dabs(kTilde(1)))*dsqrt(kappa_0(2)**2+kappa_0(3)**2+kappa_0(4)**2 )
-!
-! k_pl(1:4) = ki(1:4,j)
-! k_pl(1) = k_pl(1) +30*GeV*(1d0+kTilde(1)/50d0*GeV)
-! D(j) = Get_MInv2(k_pl)-msqi(j)
-! damping = SmoothStepFunc( +dreal(D(j)) )
-! damping2= SmoothStepFunc( -dreal(D(j)) )
-!
-! k_pl(1:4) = ki(1:4,j)
-! k_pl(1) = k_pl(1) -30*GeV*(1d0+kTilde(1)/50d0*GeV)
-! D(j) = Get_MInv2(k_pl)-msqi(j)
-! damping = damping * SmoothStepFunc(-dreal(D(j)) )
-! damping2= damping2* SmoothStepFunc(+dreal(D(j)) )
-!
-! kappa_0(:) = kappa_0(:) * (damping+damping2)
-! Get_kappa_Tri(1:4) = Get_kappa_Tri(1:4) + kappa_0(:)
-!
-! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !
-
-
-
-
-! damping = 1d0
-
-
-! damp_width = 0.5d0
-! kE_sq = VectorTimes(kTilde,kTilde)
-! D(1) = Get_MInv2(kTilde-qi(1:4,1))-msqi(1)
-! D(2) = Get_MInv2(kTilde-qi(1:4,2))-msqi(2)
-! D(3) = Get_MInv2(kTilde-qi(1:4,3))-msqi(3)
-!
-! damping = dmax1( dexp(-dreal(D(1))**2/damp_width**2)/dexp(0d0), &
-! dexp(-dreal(D(2))**2/damp_width**2)/dexp(0d0), &
-! dexp(-dreal(D(3))**2/damp_width**2)/dexp(0d0), &
-! smooth_theta( (kE_sq)**0.5-(1d4*GeV)**0.5 ) )
-
-! damping = dexp(-dreal(D(2))**2/damp_width**2)/dexp(0d0)
-
-
-
-! Get_kappa_Tri(1:4) = Get_kappa_Tri(1:4) * damping
-
-
-
-
-
-
-
-
-
-
-
-! if( ((kTilde).dot.(Get_kappa_Tri) ).lt.0d0 ) then
-!
-! print *, ""
-! print *, "wrong deformation",((kTilde).dot.(Get_kappa_Tri) )
-! print *, "Interior?",IsInInterior(kTilde)
-!
-!
-! endif
-
-
-
-! if( ((kTilde).dot.(Get_kappa_Tri) ).lt.0d0 .and. CloseToLightcone(qi(1:4,1),kTilde,msqi(1)) ) then ! if wrong deformation and close to lightcone of q1
-!
-! kLoop1Full(1:4) = kLoop1(1:4) + (0d0,1d0)*Get_kappa_Tri(1:4)
-! print *, "-----"
-! print *, "close to lightcone q1"
-! print *, "D1tilde",Get_MInv2(kTilde-qi(1:4,1))-msqi(1)
-! print *, ""
-! print *, "inside lightcone q1",InsideLightcone(qi(1:4,1),kTilde,msqi(1))
-! print *, "kappa.ktilde",(kTilde).dot.(Get_kappa_Tri)
-! print *, ""
-! print *, "-----"
-! pause
-!
-! endif
-
-
-
-! if( IsInInterior(kTilde) ) then
-!
-! print *, "in interior"
-! call printMom( kTilde(1:4) )
-! pause
-! else
-! print *, "in exterior"
-!
-! call printMom( kTilde(1:4) )
-!
-! endif
-
-
-! if( kTilde(4).lt.-5d0 ) then
-! print *, "IsInInterior",IsInInterior(kTilde)
-! call printMom( kTilde(1:4) )
-! pause
-! endif
-
RETURN
diff --git a/mod_Misc.f90 b/mod_Misc.f90
old mode 100755
new mode 100644
diff --git a/mod_Parameters.f90 b/mod_Parameters.f90
old mode 100755
new mode 100644
--
GitLab