anchor.diff:        end if
   
         p_in(j) = (one + dpop) * p0
  -!     SR, 30-08-2005: add the initial position of the slice
         s_in(j) = sp + (real(j_slices,fPrec)-one) * c_length
  -
       else
         x_in(j) = x
         y_in(j) = z
  -    end if !if(x.lt.99.0d-3) then
  -
  +    end if
     end do ! End of loop over all particles
   
   end subroutine k2coll_collimate
  
positive.code:     lint(j)   = -one
      indiv(j)  = -one
  
      x   = x_in(j)
      xp  = xp_in(j)
      z   = y_in(j)
      zp  = yp_in(j)
      p   = p_in(j)
      sp   = zero
      dpop = (p - p0)/p0
      x_flk  = zero
      y_flk  = zero
      xp_flk = zero
      yp_flk = zero
  
  !++  Transform particle coordinates to get into collimator coordinate
  !++  system
  !
  !++  First check whether particle was lost before
  !        if (x.lt.99d-3 .and. z.lt.99d-3) then
  !++  First do rotation into collimator frame
      x  =  x_in(j)*cRot + sRot*y_in(j)
      z  =  y_in(j)*cRot - sRot*x_in(j)
      xp = xp_in(j)*cRot + sRot*yp_in(j)
      zp = yp_in(j)*cRot - sRot*xp_in(j)
  !
  !++  For one-sided collimators consider only positive X. For negative
  !++  X jump to the next particle
  
  
      if((onesided .and. x < zero) .and. ((icoll /= ipencil) .or. (iturn /= 1))) then
        ! RB: Adding exception from goto if it's
        cycle
      end if
  
  !++  Now mirror at the horizontal axis for negative X offset
      if(x.lt.zero) then
        mirror = -one
        tiltangle = -one*c_tilt(2)
      end if
  
      if(x.ge.zero) then
        mirror = one
        tiltangle = c_tilt(1)
      end if
  
      x  = mirror * x
      xp = mirror * xp
  
  !++  Shift with opening and offset
  !
      x  = (x - c_aperture/two) - mirror*c_offset
  !
  !++  Include collimator tilt
  !
      if(tiltangle.gt.zero) then
        xp = xp - tiltangle
      end if
  
      if(tiltangle.lt.zero) then
        x  = x + sin_mb(tiltangle) * c_length
        xp = xp - tiltangle
      end if
  
  !++  For selected collimator, first turn reset particle distribution
  !++  to simple pencil beam
  !
  ! -- TW why did I set this to 0, seems to be needed for getting
  !       right amplitude => no "tilt" of jaw for the first turn !!!!
  !          c_tilt(1) = 0d0
  !          c_tilt(2) = 0d0
  
      nprim = 3
  
      if(( (icoll.eq.ipencil .and. iturn.eq.1) .or. (iturn.eq.1.and. ipencil.eq.999 .and. icoll.le.nprim .and. &
   &    (j.ge.(icoll-1)*nev/nprim) .and. (j.le.(icoll)*nev/nprim))).and.(pencil_distr.ne.3)) then
  ! RB addition : don't go in this if-statement if pencil_distr=3. This distribution is generated in main loop instead
  
  ! -- TW why did I set this to 0, seems to be needed for getting
  !       right amplitude => no "tilt" of jaw for the first turn !!!!
        c_tilt(1) = zero
        c_tilt(2) = zero
  
  !AUGUST2006: Standard pencil beam as implemented by GRD ------- TW
        if(pencil_rmsx.eq.zero .and. pencil_rmsy.eq.zero) then
          x  = pencil_dx(icoll)
          xp = zero
          z  = zero
          zp = zero
        end if
  
  !AUGUST2006: Rectangular (pencil-beam) sheet-beam with  ------ TW
  !            pencil_offset is the rectangulars center
  !            pencil_rmsx defines spread of impact parameter
  !            pencil_rmsy defines spread parallel to jaw surface
  
        if(pencil_distr.eq.0 .and.(pencil_rmsx.ne.0..or.pencil_rmsy.ne.0.)) then
  ! how to assure that all generated particles are on the jaw ?!
          x  = pencil_dx(icoll)+pencil_rmsx*(real(rndm4(),fPrec)-half)
          xp = zero
          z  = pencil_rmsy*(real(rndm4(),fPrec)-half)
          zp = zero
        end if
  
  !AUGUST2006: Gaussian (pencil-beam) sheet-beam with ------- TW
  !            pencil_offset is the mean  gaussian distribution
  !            pencil_rmsx defines spread of impact parameter
  !            pencil_rmsy defines spread parallel to jaw surface
  
        if(pencil_distr.eq.1 .and.(pencil_rmsx.ne.zero.or.pencil_rmsy.ne.zero )) then
          x  = pencil_dx(icoll) + pencil_rmsx*ran_gauss(two)
  ! all generated particles are on the jaw now
          x  = sqrt(x**2)
          xp = zero
          z  = pencil_rmsy*ran_gauss(two)
          zp = zero
        end if
  
  !AUGUST2006: Gaussian (pencil-beam) sheet-beam with ------- TW
  !            pencil_offset is the mean  gaussian distribution
  !            pencil_rmsx defines spread of impact parameter
  !                        here pencil_rmsx is not gaussian!!!
  !            pencil_rmsy defines spread parallel to jaw surface
  
        if(pencil_distr.eq.2 .and.(pencil_rmsx.ne.zero.or.pencil_rmsy.ne.zero )) then
          x  = pencil_dx(icoll) + pencil_rmsx*(real(rndm4(),fPrec)-half)
  ! all generated particles are on the jaw now
          x  = sqrt(x**2)
          xp = zero
          z  = pencil_rmsy*ran_gauss(two)
          zp = zero
        end if
  
  !JULY2007: Selection of pos./neg. jaw  implemented by GRD ---- TW
  
  ! ensure that for onesided only particles on pos. jaw are created
        if(onesided) then
          mirror = one
        else
          if(rndm4().lt.half) then
            mirror = -one
          else
            mirror = one
          end if
        end if
  
  ! -- TW SEP07 if c_tilt is set to zero before entering pencil beam
  !             section the assigning of the tilt will result in
  !             assigning zeros
        if(mirror.lt.zero) then
  !!     tiltangle = -one*c_tilt(2)
          tiltangle = c_tilt(2)
        else
          tiltangle = c_tilt(1)
        end if
  !!!!--- commented this out since particle is tilted after leaving
  !!!!--- collimator -> remove this  code fragment in final verion
  !!             x  = mirror * x
  !!             xp = mirror * xp
  
  
  !++  Include collimator tilt
  ! this is propably not correct
  !
  !             xp =  (xp_pencil0*cRot+                         &
  !     &            sRot*yp_pencil0)
  !             if (tiltangle.gt.0.) then
  !                xp = xp - tiltangle
  !!             endif
  !!             elseif (tiltangle.lt.0.) then
  !             else
  !               x  = x + sin_mb(tiltangle) * c_length
  !               xp = xp - tiltangle
  !             endif
  !
        write(coll_pencilUnit,'(f10.8,(2x,f10.8),(2x,f10.8),(2x,f10.8),(2x,f10.8))') x, xp, z, zp, tiltangle
  
      end if !if(( (icoll.eq.ipencil .and. iturn.eq.1) .or. (itu
  
  !     SR, 18-08-2005: after finishing the coordinate transformation,
  !     or the coordinate manipulations in case of pencil beams,
  !     write down the initial coordinates of the impacting particles
      xinn  = x
      xpinn = xp
      yinn  = z
      ypinn = zp
  !
  !++  Possibility to slice here (RA,SR: 29-08-2005)
  !
  !++  particle passing above the jaw are discarded => take new event
  !++  entering by the face, shorten the length (zlm) and keep track of
  !++  entrance longitudinal coordinate (keeps) for histograms
  !
  !++  The definition is that the collimator jaw is at x>=0.
  !
  !++  1) Check whether particle hits the collimator
      hit   = .false.
      s     = zero
      keeps = zero
      zlm   = -one * length
  
      if(x.ge.zero) then
  
  !++  Particle hits collimator and we assume interaction length ZLM equal
  !++  to collimator length (what if it would leave collimator after
  !++  small length due to angle???)
        zlm = length
        impact(j) = x
        indiv(j) = xp
      else if(xp.le.zero) then
  
  !++  Particle does not hit collimator. Interaction length ZLM is zero.
        zlm = zero
      else
  !
  !++  Calculate s-coordinate of interaction point
        s = (-one*x) / xp
        if(s.le.0) then
          write(lout,*) 'S.LE.0 -> This should not happen'
          call prror
        end if
  
        if(s .lt. length) then
          zlm = length - s
          impact(j) = zero
          indiv(j) = xp
        else
          zlm = zero
        end if
      end if !if(x.ge.0.d0) then
  
  !++  First do the drift part
  ! DRIFT PART
      drift_length = length - zlm
      if(drift_length.gt.zero) then
        if(iexact) then
          zpj = sqrt(one-xp**2-zp**2)
          x = x + drift_length*(xp/zpj)
          z = z + drift_length*(zp/zpj)
          sp = sp + drift_length
        else
          x  = x + xp* drift_length
          z  = z + zp * drift_length
          sp = sp + drift_length
        end if
      end if
  !
  !++  Now do the scattering part
  !
      if (zlm.gt.zero) then
        if(.not.linside(j)) then
          ! first time particle hits collimator: entering jaw
          linside(j)=.true.
          if(dowrite_impact) then
            if ( tiltangle.gt.zero ) then
              x_Dump=(x+c_aperture/two+tiltangle*sp)*mirror+c_offset
            else
              x_Dump=(x+c_aperture/two+tiltangle*(sp-c_length))*mirror+c_offset
  
negative.code:   real(kind=fPrec) s_impact
    integer flagsec(npart)
  
    ! SR, 18-08-2005: add temporary variable to write in FirstImpacts
    ! the initial distribution of the impacting particles in the collimator frame.
    real(kind=fPrec) xinn,xpinn,yinn,ypinn
  
    ! SR, 29-08-2005: add the slice number to calculate the impact
    ! location within the collimator.
    ! j_slices = 1 for the a non sliced collimator!
    integer j_slices
  
    mat    = cdb_cMaterialID(icoll)
    length = c_length
    nev    = np
    p0     = enom
  
    ! Initialize scattering processes
    call k2coll_scatin(p0)
  
    nhit   = 0
    fracab = zero
    mirror = one
  
    ! Compute rotation factors for collimator rotation
    cRot   = cos_mb(c_rotation)
    sRot   = sin_mb(c_rotation)
    cRRot  = cos_mb(-c_rotation)
    sRRot  = sin_mb(-c_rotation)
  
    do j=1,nev
  
      if(part_abs_pos_local(j) /= 0 .and. part_abs_turn_local(j) /= 0) then
        ! Don't do scattering process for particles already absorbed
  