diff --git a/mod_LoopSpace.f90 b/mod_LoopSpace.f90
index 7d1764669c0e52118bad246146b377be9acaa468..c8db663b805296f57ffef0eb58a9f73a7a6c62f6 100755
--- a/mod_LoopSpace.f90
+++ b/mod_LoopSpace.f90
@@ -35,9 +35,10 @@ real(8) :: sinZ,cosZ,sinT,cosT,sinP,cosP
 
 
 !   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"
+
 
     kE    = dsqrt( mu_Loop1_sq*dtan( 0.5d0*Pi*uRnd(1) ) )
     kE_sq=kE**2
@@ -52,7 +53,7 @@ real(8) :: sinZ,cosZ,sinT,cosT,sinP,cosP
     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 
 
@@ -561,7 +562,7 @@ integer :: j
     kappa_Ext(1)   = c_pl * k_pl(1)   + c_mi * k_mi(1)
     kappa_Ext(2:4) =-c_pl * k_pl(2:4) - c_mi * k_mi(2:4)
     
-    
+   
 !     print *, "kTilde",kTilde
 !     print *, "absvec kTilde",dsqrt(kTilde(2)**2+kTilde(3)**2+kTilde(4)**2),1000d0*0.01d0
 !     print *, "kappa_Ext",kappa_Ext
@@ -569,8 +570,8 @@ integer :: j
 !     print *, "ki(1:4,2).dot.kappa_Ext",ki(1:4,2).dot.kappa_Ext(1:4),(ki(1:4,2).dot.ki(1:4,2))-msqi(2)
     
         
-    kappa_Int(:)=0d0 ! removing kappa_int for now...
-!   compute kappa_Int   
+!     kappa_Int(:)=0d0 ! removing kappa_int for now...
+!   !compute kappa_Int   
 !     vij(:,ij12) = 0.5d0*( qi(:,1)+qi(:,2) - (dsqrt(msqi(1))-dsqrt(msqi(2)))*(qi(:,1)-qi(:,2))/Get_MInv(qi(1:4,1)-qi(1:4,2)) )
 !     vij(:,ij13) = 0.5d0*( qi(:,1)+qi(:,3) - (dsqrt(msqi(1))-dsqrt(msqi(3)))*(qi(:,1)-qi(:,3))/Get_MInv(qi(1:4,1)-qi(1:4,3)) )
 !     vij(:,ij23) = 0.5d0*( qi(:,2)+qi(:,3) - (dsqrt(msqi(2))-dsqrt(msqi(3)))*(qi(:,2)-qi(:,3))/Get_MInv(qi(1:4,2)-qi(1:4,3)) )
@@ -632,6 +633,46 @@ integer :: j
 !  pause
  
 
+ 
+ 
+!  if( ((kTilde).dot.(Get_kappa) ).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(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)
+!     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
 END FUNCTION
 
@@ -640,10 +681,10 @@ END FUNCTION
 
 
 
-FUNCTION IsInInterior(kTilde)
+FUNCTION IsInInterior(k)
 use ModMisc
 implicit none
-real(8) :: kTilde(1:4),P_pl(1:4),P_mi(1:4)
+real(8) :: k(1:4),P_pl(1:4),P_mi(1:4)
 logical :: IsInInterior
 
     IsInInterior=.false.
@@ -651,7 +692,7 @@ logical :: IsInInterior
     P_pl = get_Zpl(qi(:,1),+qi(:,1))
     P_mi = get_Zmi(qi(:,2),-qi(:,2))   
 
-    if( InsideFWDLightcone(P_pl,kTilde) .and. InsideBWDLightcone(P_mi,kTilde) ) IsInInterior=.true.
+    if( InsideFWDLightcone(P_pl,k) .and. InsideBWDLightcone(P_mi,k) ) IsInInterior=.true.
 
 RETURN
 END FUNCTION 
diff --git a/mod_Misc.f90 b/mod_Misc.f90
index 40701eeb7f9e620f33b9c7a985a6a374875fdcbc..5a38f3b5bfadc5ef43a08040e6e3a1818a7f2d33 100755
--- a/mod_Misc.f90
+++ b/mod_Misc.f90
@@ -235,7 +235,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-2    !*
 ! *******************************************
     
     
@@ -1419,6 +1419,19 @@ END SUBROUTINE
 
 
 
+
+SUBROUTINE printMom(p)
+implicit none
+real(8) :: p(1:4)
+
+  write(*,"(A,1F22.16,A,1F22.16,A,1F22.16,A,1F22.16)") "E=",p(1),"    pX=",p(2),"    pY=",p(3),"    pZ=",p(4)
+
+
+
+END SUBROUTINE
+
+
+
 SUBROUTINE SwitchEnergyComponent(p)
 implicit none
 real(8) :: p(:,:),px,py,pz,E