anchor.diff:    integer(kind=int16) :: nnuc0,nnuc1
   #endif
   
  -!-----------------------------------------------------------------------
  -!GRD NEW COLLIMATION PARAMETERS
  -!-----------------------------------------------------------------------
  -!++  Get the aperture from the beta functions and emittance
  -!++  A simple estimate of beta beating can be included that
  -!++  has twice the betatron phase advance
  +  ! Get the aperture from the beta functions and emittance
  +  ! A simple estimate of beta beating can be included that has twice the betatron phase advance
     if(.not. cdb_doNSig) then
       nsig = cdb_cNSig(icoll)
     end if
  
positive.code:   use geant4
  #endif
  
    implicit none
  
    real(kind=fPrec), intent(in) :: stracki
  
    integer j, iSlice, nSlices
    logical onesided, linside(napx), isAbs, isHit
    real(kind=fPrec) jawLength,jawAperture,jawOffset,jawTilt(2),x_Dump,xpDump,y_Dump,ypDump,s_Dump,   &
      xmax,ymax,calc_aperture,ldrift,c_nex2,c_ney2,Nap1pos,Nap2pos,Nap1neg,Nap2neg,tiltOffsPos1,      &
      tiltOffsPos2,tiltOffsNeg1,tiltOffsNeg2,beamsize1,beamsize2,betax1,betax2,betay1,betay2,alphax1, &
      alphax2,alphay1,alphay2,minAmpl,zpj,xmax_pencil,ymax_pencil,xmax_nom,ymax_nom,nom_aperture,     &
      scale_bx,scale_by,c_tilt(2),c_offset,c_aperture,c_rotation,cry_bendangle,cry_tilt,cry_tilt0
  
  #ifdef G4COLLIMATION
  
negative.code: 
  end subroutine coll_endTurn
  
  
  ! ================================================================================================ !
  !  R. Bruce
  !
  !  Create distribution where the normalized distance between jaw and beam is the smallest.
  !  This is where particles will first impact:
  !  Without imperfections, it is:
  !   - at the face of the collimator for the case of beta'<0 (POSITIVE alpha - beam converging) and
  !   - at the exit of the collimator for the case of beta'>0 (NEGATIVE alpha beam diverging)
  !  With imperfections:
  !   - include errors on gap, tilt and offset. We have to calculate the normalized distance to each
  !     corner separately!
  !
  ! Calculate optical parameters at start and end of collimator (half a collimator length upstream
  ! and downstream of present s-position). Assuming a purely vertical or horizontal halo - need to
  ! add more conditions for other cases!
  !
  ! Using standard twiss transfer matrix for a drift at start of collimator
  ! ================================================================================================ !
  subroutine coll_matchedHalo(c_tilt,c_offset,c_aperture)
  
    use crcoall
    use coll_dist
    use mod_common
    use mod_common_main
    use mod_common_track
    use mathlib_bouncer
  
    real(kind=fPrec), intent(in) :: c_tilt(2)
    real(kind=fPrec), intent(in) :: c_offset
    real(kind=fPrec), intent(in) :: c_aperture
  
    integer j
    real(kind=fPrec) Nap1pos,Nap2pos,Nap1neg,Nap2neg,tiltOffsPos1,tiltOffsPos2,tiltOffsNeg1,     &
      tiltOffsNeg2,beamsize1,beamsize2,minAmpl,ldrift,c_nex2,c_ney2,betax1,betax2,betay1,betay2, &
      alphax1,alphax2,alphay1,alphay2
  
    ! Assign the drift length over which the optics functions are propagated
    ldrift = -c_length/two
    betax1 = tbetax(ie) - two*ldrift*talphax(ie) + (ldrift**2 * (one+talphax(ie)**2))/tbetax(ie)
    betay1 = tbetay(ie) - two*ldrift*talphay(ie) + (ldrift**2 * (one+talphay(ie)**2))/tbetay(ie)
  
    alphax1 = talphax(ie) - (ldrift*(1+talphax(ie)**2))/tbetax(ie)
    alphay1 = talphay(ie) - (ldrift*(1+talphay(ie)**2))/tbetay(ie)
  
    ! At end of collimator:
    ldrift = c_length/two
    betax2 = tbetax(ie) - two*ldrift*talphax(ie) + (ldrift**2 * (one+talphax(ie)**2))/tbetax(ie)
    betay2 = tbetay(ie) - two*ldrift*talphay(ie) + (ldrift**2 * (one+talphay(ie)**2))/tbetay(ie)
  
    alphax2 = talphax(ie) - (ldrift*(1+talphax(ie)**2))/tbetax(ie)
    alphay2 = talphay(ie) - (ldrift*(1+talphay(ie)**2))/tbetay(ie)
  
    ! Calculate beam size at start and end of collimator. account for collimation plane
    if(cdist_ampX > zero .and. cdist_ampY == zero) then ! Horizontal halo
      beamsize1 = sqrt(betax1 * c_emitx0_collgap)
      beamsize2 = sqrt(betax2 * c_emitx0_collgap)
    else if(cdist_ampX == zero .and. cdist_ampY > zero) then ! Vertical halo
      beamsize1 = sqrt(betay1 * c_emity0_collgap)
      beamsize2 = sqrt(betay2 * c_emity0_collgap)
    else
      write(lerr,"(a)") "COLL> ERROR Attempting to use a halo not purely in the horizontal "//&
        "or vertical plane with pencil_dist=3 - abort."
      call prror
    end if
  
    ! Calculate offset from tilt of positive and negative jaws, at start and end
    ! Remember: tilt angle is defined such that one corner stays at nominal position, the other corner is more open
  
    ! Jaw in positive x (or y):
    if(c_tilt(1) >= zero) then
      tiltOffsPos1 = zero
      tiltOffsPos2 = abs(sin_mb(c_tilt(1))) * c_length
    else
      tiltOffsPos1 = abs(sin_mb(c_tilt(1))) * c_length
      tiltOffsPos2 = zero
    end if
  
    ! Jaw in negative x (or y):
    if(c_tilt(2) >= zero) then
      tiltOffsNeg1 = abs(sin_mb(c_tilt(2))) * c_length
      tiltOffsNeg2 = zero
    else
      tiltOffsNeg1 = zero
      tiltOffsNeg2 = abs(sin_mb(c_tilt(2))) * c_length
    end if
  
    ! Calculate half distance from jaws to beam center (in units of beam sigma) at the beginning of the collimator,
    ! Positive and neg jaws.
    Nap1pos = ((c_aperture/two + c_offset) + tiltOffsPos1)/beamsize1
    Nap2pos = ((c_aperture/two + c_offset) + tiltOffsPos2)/beamsize2
    Nap1neg = ((c_aperture/two - c_offset) + tiltOffsNeg1)/beamsize1
    Nap2neg = ((c_aperture/two - c_offset) + tiltOffsNeg2)/beamsize2
  
    ! Minimum normalized distance from jaw to beam center - this is the n_sigma at which the halo should be generated
    minAmpl = min(Nap1pos,Nap2pos,Nap1neg,Nap2neg)
  
    ! Assign amplitudes in x and y for the halo generation function
    if(cdist_ampX > zero .and. cdist_ampY == zero) then ! Horizontal halo
      c_nex2 = minAmpl
    else if(cdist_ampX == zero .and. cdist_ampY > zero) then ! Vertical halo
      c_ney2 = minAmpl
    end if ! Other cases taken care of above - in these cases, program has already stopped
  
    ! Assign optics parameters to use for the generation of the starting halo - at start or end of collimator
    if(minAmpl == Nap1pos .or. minAmpl == Nap1neg) then ! min normalized distance occurs at start of collimator
      c_betax  = betax1
      c_betay  = betay1
      c_alphax = alphax1
      c_alphay = alphay1
      ldrift   = -c_length/two
    else ! Min normalized distance occurs at end of collimator
      c_betax  = betax2
      c_betay  = betay2
      c_alphax = alphax2
      c_alphay = alphay2
      ldrift   = c_length/two
    end if
  
    ! create new pencil beam distribution with spread at start or end of collimator at the minAmpl
    ! note: if imperfections are active, equal amounts of particles are still generated on the two jaws.
    ! but it might be then that only one jaw is hit on the first turn, thus only by half of the particles
    ! the particle generated on the other side will then hit the same jaw several turns later, possibly smearing the impact parameter
    ! This could possibly be improved in the future.
    call cdist_makeDist_coll(c_alphax,c_alphay,c_betax,c_betay,c_nex2,c_ney2)
  
    do j=1,napx
      xv1(j)  = c1e3*xv1(j) + torbx(ie)
      yv1(j)  = c1e3*yv1(j) + torbxp(ie)
      xv2(j)  = c1e3*xv2(j) + torby(ie)
      yv2(j)  = c1e3*yv2(j) + torbyp(ie)
  
      ! as main routine will track particles back half a collimator length (to start of jaw),
      ! track them now forward (if generated at face) or backward (if generated at end)
      ! 1/2 collimator length to center of collimator (ldrift pos or neg)
      xv1(j)  = xv1(j) - ldrift*yv1(j)
      xv2(j)  = xv2(j) - ldrift*yv2(j)
    end do
  
  end subroutine coll_matchedHalo
  
  ! ================================================================================================ !
  !  Find the smallest gap, and also write sigmasettings.out
  !  Updated: 2019-10-10
  