anchor.diff:  #endif
   
     if(firstrun) then
  -
  -    if(iturn == 1) then
  -      sum_ax(ie) = zero
  -      sum_ay(ie) = zero
  -    end if
  -
  -    do j = 1, napx
  -      xj  = (xv1(j)-torbx(ie)) /c1e3
  -      xpj = (yv1(j)-torbxp(ie))/c1e3
  -      yj  = (xv2(j)-torby(ie)) /c1e3
  -      ypj = (yv2(j)-torbyp(ie))/c1e3
  -      pj  = ejv(j)/c1e3
  -
  -      if(tbetax(ie) > zero) then
  -        gammax = (one + talphax(ie)**2)/tbetax(ie)
  -        gammay = (one + talphay(ie)**2)/tbetay(ie)
  -      else
  -        gammax = (one + talphax(ie-1)**2)/tbetax(ie-1)
  -        gammay = (one + talphay(ie-1)**2)/tbetay(ie-1)
  -      endif
  -
  -      if(part_abs_pos(j) == 0 .and. part_abs_turn(j) == 0) then
  -        if(tbetax(ie) > 0.) then
  -          nspx = sqrt(abs( gammax*(xj)**2 + two*talphax(ie)*xj*xpj +   tbetax(ie)*xpj**2 )/c_emitx0_collgap)
  -          nspy = sqrt(abs( gammay*(yj)**2 + two*talphay(ie)*yj*ypj +   tbetay(ie)*ypj**2 )/c_emity0_collgap)
  -        else
  -          nspx = sqrt(abs( gammax*(xj)**2 + two*talphax(ie-1)*xj*xpj + tbetax(ie-1)*xpj**2 )/c_emitx0_collgap)
  -          nspy = sqrt(abs( gammay*(yj)**2 + two*talphay(ie-1)*yj*ypj + tbetay(ie-1)*ypj**2 )/c_emity0_collgap)
  -        end if
  -
  -        sum_ax(ie)   = sum_ax(ie) + nspx
  -        sqsum_ax(ie) = sqsum_ax(ie) + nspx**2
  -        sum_ay(ie)   = sum_ay(ie) + nspy
  -        sqsum_ay(ie) = sqsum_ay(ie) + nspy**2
  -        nampl(ie)    = nampl(ie) + 1
  -      else
  -        nspx = zero
  -        nspy = zero
  -      end if
  -    end do
  +    call coll_computeStats
     end if
   
     ! Note: Not in firstrun
     if(dowritetracks) then
  +    gammax = (one + talphax(ie)**2)/tbetax(ie)
  +    gammay = (one + talphay(ie)**2)/tbetay(ie)
  +
       do j=1,napx
         xj  = (xv1(j)-torbx(ie)) /c1e3
         xpj = (yv1(j)-torbxp(ie))/c1e3
  
positive.code: 
    integer j,fUnit,ieff,ieffdpop
    real(kind=fPrec) gammax,gammay,xdisp,dnormx,dnormy,driftx,drifty,xnorm,xpnorm,xangle,ynorm,ypnorm,&
      yangle,c_dpop,dpopmin,dpopmax
  
  #ifdef HDF5
    ! For tracks2
  
negative.code: end subroutine collimate_end_turn
  
  ! ================================================================================================ !
  !  Find the smallest gap, and also write sigmasettings.out
  !  Updated: 2019-10-10
  ! ================================================================================================ !
  subroutine coll_getMinGapID(minGapID)
  