anchor.diff:    sqsum    = zero
   
   #ifdef G4COLLIMATION
  -  do j = 1, napx
  -      if (stracki.eq.0.) then
  -        if(iexact .eqv. .false.) then
  -!          write(lout,*) 'iexact = 0', rcxp(j), rcyp(j)
  -          rcx(j)  = rcx(j) - half*c_length*rcxp(j)
  -          rcy(j)  = rcy(j) - half*c_length*rcyp(j)
  -        else
  -!          write(lout,*) 'iexact = 1', rcxp(j), rcyp(j)
  -          zpj=sqrt(one-rcxp(j)**2-rcyp(j)**2)
  +  do j=1,napx
  +    if(stracki == zero) then
  +      if(iexact .eqv. .false.) then
  +        rcx(j) = rcx(j) - half*c_length*rcxp(j)
  +        rcy(j) = rcy(j) - half*c_length*rcyp(j)
  +      else
  +        zpj    = sqrt(one-rcxp(j)**2-rcyp(j)**2)
  +        rcx(j) = rcx(j) - half*c_length*(rcxp(j)/zpj)
  +        rcy(j) = rcy(j) - half*c_length*(rcyp(j)/zpj)
  +      end if
  +    end if
  +
  +    ! Now copy data back to original verctor
  +    xv1(j) = rcx(j)  * c1e3 + torbx(ie)
  +    yv1(j) = rcxp(j) * c1e3 + torbxp(ie)
  +    xv2(j) = rcy(j)  * c1e3 + torby(ie)
  +    yv2(j) = rcyp(j) * c1e3 + torbyp(ie)
  +    ejv(j) = rcp(j)  * c1e3
  +
  +    ! Update mtc and other arrays.
  +    ejfv    (j) = sqrt(ejv(j)**2-nucm(j)**2)
  +    rvv     (j) = (ejv(j)*e0f)/(e0*ejfv(j))
  +    dpsv    (j) = (ejfv(j)*(nucm0/nucm(j))-e0f)/e0f
  +    oidpsv  (j) = one/(one+dpsv(j))
  +    dpsv1   (j) = (dpsv(j)*c1e3)*oidpsv(j)
  +    mtc     (j) = (nqq(j)*nucm0)/(qq0*nucm(j))
  +    moidpsv (j) = mtc(j)*oidpsv(j)
  +    omoidpsv(j) = c1e3*((one-mtc(j))*oidpsv(j))
  +    yv1     (j) = ejf0v(j)/ejfv(j)*yv1(j)
  +    yv2     (j) = ejf0v(j)/ejfv(j)*yv2(j)
  +  end do
  +#else
  +  ! Copy particle data back and do path length stuff; check for absorption
  +  ! Add orbit offset back.
  +  do j=1,napx
  +    ! In order to get rid of numerical errors, just do the treatment for impacting particles
  +    if(part_hit_pos(j) == ie .and. part_hit_turn(j) == iturn) then
  +      ! For zero length element track back half collimator length
  +      if(stracki == zero) then
  +        if(iexact) then
  +          zpj    = sqrt(one-rcxp(j)**2-rcyp(j)**2)
             rcx(j) = rcx(j) - half*c_length*(rcxp(j)/zpj)
             rcy(j) = rcy(j) - half*c_length*(rcyp(j)/zpj)
  +        else
  +          rcx(j) = rcx(j) - half*c_length*rcxp(j)
  +          rcy(j) = rcy(j) - half*c_length*rcyp(j)
           end if
         end if
   
  -!++  Now copy data back to original verctor
  +      ! Now copy data back to original verctor
         xv1(j) = rcx(j)  * c1e3 + torbx(ie)
         yv1(j) = rcxp(j) * c1e3 + torbxp(ie)
         xv2(j) = rcy(j)  * c1e3 + torby(ie)
         yv2(j) = rcyp(j) * c1e3 + torbyp(ie)
         ejv(j) = rcp(j)  * c1e3
   
  -
  -! Update mtc and other arrays.
  -!      ejv(j)    = rcp(j)  * c1e3
  -!!!  write(lout,*) 'ejfv', ejv(j), nucm(j)
  -!      ejfv  (j) = sqrt((ejv(j)-nucm(j))*(ejv(j)+nucm(j)))   ! hisix: ion mass
  -      ejfv(j)=sqrt(ejv(j)**2-nucm(j)**2)
  -!!!  write(lout,*) 'ejfv ok', ejfv(j)
  -      rvv   (j) = (ejv(j)*e0f)/(e0*ejfv(j))                 ! hisix: remains unchanged
  -      dpsv  (j) = (ejfv(j)*(nucm0/nucm(j))-e0f)/e0f         ! hisix: new delta
  -      oidpsv(j) = one/(one+dpsv(j))
  -      dpsv1 (j) = (dpsv(j)*c1e3)*oidpsv(j)
  -      mtc     (j) = (nqq(j)*nucm0)/(qq0*nucm(j))            ! hisix: mass to charge
  -      moidpsv (j) = mtc(j)*oidpsv(j)                        ! hisix
  -      omoidpsv(j) = c1e3*((one-mtc(j))*oidpsv(j))           ! hisix
  -
  -!++  Energy update, as recommended by Frank
  -!      ejfv(j)=sqrt(ejv(j)**2-nucm(j)**2)
  -!      rvv(j)=(ejv(j)*e0f)/(e0*ejfv(j))
  -!      dpsv(j)=(ejfv(j)*(nucm0/nucm(j))-e0f)/e0f
  -!      oidpsv(j)=one/(one+dpsv(j))
  -!      moidpsv(j)=mtc(j)/(one+dpsv(j))
  -!      omoidpsv(j)=c1e3*((one-mtc(j))*oidpsv(j))
  -!      dpsv1(j)=(dpsv(j)*c1e3)*oidpsv(j)
  -      yv1(j)   = ejf0v(j)/ejfv(j)*yv1(j)
  -      yv2(j)   = ejf0v(j)/ejfv(j)*yv2(j)
  -!!!  write(lout,*) 'Coordinate loop end'
  -end do
  -#endif
  -
  -#ifndef G4COLLIMATION
  -!++  Copy particle data back and do path length stuff; check for absorption
  -!++  Add orbit offset back.
  -  do j = 1, napx
  -
  -!APRIL2005 IN ORDER TO GET RID OF NUMERICAL ERRORS, JUST DO THE TREATMENT FOR
  -!APRIL2005 IMPACTING PARTICLES...
  -    if(part_hit_pos(j) .eq.ie .and. part_hit_turn(j).eq.iturn) then
  -!++  For zero length element track back half collimator length
  -! DRIFT PART
  -      if (stracki.eq.0.) then
  -        if(iexact) then
  -          zpj=sqrt(one-rcxp(j)**2-rcyp(j)**2)
  -          rcx(j) = rcx(j) - half*c_length*(rcxp(j)/zpj)
  -          rcy(j) = rcy(j) - half*c_length*(rcyp(j)/zpj)
  -        else
  -          rcx(j)  = rcx(j) - half*c_length*rcxp(j)
  -          rcy(j)  = rcy(j) - half*c_length*rcyp(j)
  -        end if
  -      end if
  -
  -!++  Now copy data back to original verctor
  -      xv1(j) = rcx(j)  * c1e3 + torbx(ie)
  -      yv1(j) = rcxp(j) * c1e3 + torbxp(ie)
  -      xv2(j) = rcy(j)  * c1e3 + torby(ie)
  -      yv2(j) = rcyp(j) * c1e3 + torbyp(ie)
  -      ejv(j)  = rcp(j)  * c1e3
  -
  -!++  Energy update, as recommended by Frank
  -      ejfv(j)=sqrt(ejv(j)**2-nucm(j)**2)
  -      rvv(j)=(ejv(j)*e0f)/(e0*ejfv(j))
  -      dpsv(j)=(ejfv(j)*(nucm0/nucm(j))-e0f)/e0f
  -      oidpsv(j)=one/(one+dpsv(j))
  -      mtc(j) = (nqq(j)*nucm0)/(qq0*nucm(j))            ! hisix: mass to charge
  -      moidpsv(j)=mtc(j)/(one+dpsv(j))
  -      omoidpsv(j)=c1e3*((one-mtc(j))*oidpsv(j))
  -      dpsv1(j)=(dpsv(j)*c1e3)*oidpsv(j)
  -      yv1(j)   = ejf0v(j)/ejfv(j)*yv1(j)
  -      yv2(j)   = ejf0v(j)/ejfv(j)*yv2(j)
  -
  -!++   For absorbed particles set all coordinates to zero. Also
  -!++   include very large offsets, let's say above 100mm or
  -!++   100mrad.
  +      !  Energy update, as recommended by Frank
  +      ejfv(j)     = sqrt(ejv(j)**2-nucm(j)**2)
  +      rvv(j)      = (ejv(j)*e0f)/(e0*ejfv(j))
  +      dpsv(j)     = (ejfv(j)*(nucm0/nucm(j))-e0f)/e0f
  +      oidpsv(j)   = one/(one+dpsv(j))
  +      mtc(j)      = (nqq(j)*nucm0)/(qq0*nucm(j))
  +      moidpsv(j)  = mtc(j)/(one+dpsv(j))
  +      omoidpsv(j) = c1e3*((one-mtc(j))*oidpsv(j))
  +      dpsv1(j)    = (dpsv(j)*c1e3)*oidpsv(j)
  +      yv1(j)      = (ejf0v(j)/ejfv(j))*yv1(j)
  +      yv2(j)      = (ejf0v(j)/ejfv(j))*yv2(j)
  +
  +      ! For absorbed particles set all coordinates to zero. Also include very
  +      ! large offsets, let's say above 100mm or 100mrad.
         if((part_abs_pos(j) /= 0 .and. part_abs_turn(j) /= 0) .or. &
           xv1(j) > c1e2 .or. yv1(j) > c1e2 .or. xv2(j) > c1e2 .or. yv2(j) > c1e2) then
           xv1(j)           = zero
  
positive.code:             end if
  
  !++ 2D eff
              do ieffdpop = 1, numeffdpop
                if(counted2d(j,ieff,ieffdpop).eq.0 .and.abs((ejv(j)-c_enom)/c_enom).ge.dpopbins(ieffdpop)) then
                  neff2d(ieff,ieffdpop) = neff2d(ieff,ieffdpop)+one
                  counted2d(j,ieff,ieffdpop)=1
  
negative.code:   use mathlib_bouncer
  
    integer i, j, coll_mingap_id
    character(len=mNameLen) coll_mingap
    real(kind=fPrec) dummy,nsig_err,sig_offset,mingap,gap_h1,gap_h2,gap_h3,gap_h4,c_rmstilt,c_systilt
  
    ! Initialisation
    do i=1,numeff
  