anchor.diff:    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.
  +  ! 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!
  +  ! 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
  -
  -  save
   
     mat    = cdb_cMaterialID(icoll)
     length = c_length
     nev    = np
     p0     = enom
   
  -!++  Initialize scattering processes
  +  ! Initialize scattering processes
     call k2coll_scatin(p0)
   
  -! EVENT LOOP,  initial distribution is here a flat distribution with
  -! xmin=x-, xmax=x+, etc. from the input file
  -
  -  nhit    = 0
  -  fracab  = zero
  -  mirror  = one
  -  cRot    = cos_mb(c_rotation)
  -  sRot    = sin_mb(c_rotation)
  -  cRRot   = cos_mb(-c_rotation)
  -  sRRot   = sin_mb(-c_rotation)
  -
  -  do j = 1, nev
  +  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
  
positive.code:       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 do ! End of loop over all particles
  
  end subroutine k2coll_collimate
  
negative.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 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
        cycle
      end if
  
      ! Now mirror at the horizontal axis for negative X offset
      if(x < zero) then
        mirror    = -one
        tiltangle = -one*c_tilt(2)
      end if
      if(x >= 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 > zero) then
        xp = xp - tiltangle
      end if
      if(tiltangle < 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
      nprim = 3
      if(((icoll == ipencil .and. iturn == 1) .or. (iturn == 1 .and. ipencil == 999 .and. icoll <= nprim .and. &
         (j >= (icoll-1)*nev/nprim) .and. (j <= (icoll)*nev/nprim))) .and. (pencil_distr /= 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
  
        ! Standard pencil beam as implemented by GRD ------- TW
        if(pencil_rmsx == zero .and. pencil_rmsy == zero) then
          x  = pencil_dx(icoll)
          xp = zero
          z  = zero
          zp = zero
        end if
  
        ! 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 == 0 .and.(pencil_rmsx /= zero .or. pencil_rmsy /= zero)) 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
  
        ! 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 == 1 .and. (pencil_rmsx /= zero .or. pencil_rmsy /= 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
  
        ! 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 == 2 .and. (pencil_rmsx /= zero .or. pencil_rmsy /= 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
  
        ! 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() < half) then
            mirror = -one
          else
            mirror = one
          end if
        end if
  
        ! TW if c_tilt is set to zero before entering pencil beam
        !    section the assigning of the tilt will result in assigning zeros
        if(mirror < zero) then
          tiltangle = c_tilt(2)
        else
          tiltangle = c_tilt(1)
        end if
  !
        write(coll_pencilUnit,"(f10.8,4(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
  
      ! 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 >= 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 <= 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 <= zero) then
          write(lerr,"(a)") "COLLK2> ERROR S <= zero. This should not happen!"
          call prror
        end if
        if(s < length) then
          zlm       = length - s
          impact(j) = zero
          indiv(j)  = xp
        else
          zlm = zero
        end if
      end if
  
      ! First do the drift part
      ! DRIFT PART
      drift_length = length - zlm
      if(drift_length > 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 > zero) then
        if(.not.linside(j)) then
          ! first time particle hits collimator: entering jaw
          linside(j)=.true.
          if(dowrite_impact) then
            if(tiltangle > 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
  