anchor.diff:  
     do j=1,napx
       if(part_hit_pos(j) == ie .and. part_hit_turn(j) == iturn) then
  -      ! Particle has impacted
  +
  +      ! Calculate impact observables, fill histograms, save collimator info
  +      n_impact         = n_impact + 1
  +      sum              = sum   + part_impact(j)
  +      sqsum            = sqsum + part_impact(j)**2
  +
  +      cn_impact(icoll) = cn_impact(icoll) + 1
  +      csum(icoll)      = csum(icoll)   + part_impact(j)
  +      csqsum(icoll)    = csqsum(icoll) + part_impact(j)**2
  +
         if(part_abs_pos(j) /= 0 .and. part_abs_turn(j) /= 0) then
  -        continue
  -      else if(part_abs_pos (j) == 0 .and. part_abs_turn(j) == 0) then
  +        ! If the interacting particle was lost, add-up counters for absorption
  +        ! Note: a particle with x/y >= 99. never hits anything any more in the logic of this program.
  +        ! Be careful to always fulfill this!
  +        n_absorbed         = n_absorbed + 1
  +        cn_absorbed(icoll) = cn_absorbed(icoll) + 1
  +        n_tot_absorbed     = n_tot_absorbed + 1
  +        iturn_last_hit     = part_hit_before_turn(j)
  +        iturn_absorbed     = part_hit_turn(j)
  +        if(iturn_last_hit == 0) then
  +          iturn_last_hit = iturn_absorbed
  +        end if
  +
  +      elseif(part_abs_pos (j) == 0 .and. part_abs_turn(j) == 0) then
           nhit_type(j) = ior(nhit_type(j),cdb_cType(icoll)) ! Record the hit type
  +
         else
           write(lerr,"(a)")          "COLL> ERROR Particle cannot be both absorbed and not absorbed"
           write(lerr,"(a,2(1x,i0))") "COLL>      ",part_abs_pos (j),part_abs_turn(j)
           call prror
         end if
  +
       end if
     end do
  -
  -  do j=1,napx
  -    if(part_hit_pos(j) == ie .and. part_hit_turn(j) == iturn) then
  -!++  Calculate impact observables, fill histograms, save collimator info, ...
  -      n_impact = n_impact + 1
  -      sum = sum + part_impact(j)
  -      sqsum = sqsum + part_impact(j)**2
  -      cn_impact(icoll) = cn_impact(icoll) + 1
  -      csum(icoll) = csum(icoll) + part_impact(j)
  -      csqsum(icoll) = csqsum(icoll) + part_impact(j)**2
  -
  -!++  If the interacting particle was lost, add-up counters for absorption
  -!++  Note: a particle with x/y >= 99. never hits anything any more in
  -!++        the logic of this program. Be careful to always fulfill this!
  -      if(part_abs_pos(j).ne.0 .and. part_abs_turn(j).ne.0) then
  -        n_absorbed = n_absorbed + 1
  -        cn_absorbed(icoll) = cn_absorbed(icoll) + 1
  -        n_tot_absorbed = n_tot_absorbed + 1
  -        iturn_last_hit = part_hit_before_turn(j)
  -        iturn_absorbed = part_hit_turn(j)
  -        if(iturn_last_hit == 0) then
  -          iturn_last_hit = iturn_absorbed
  -        end if
  -      end if
  -
  -!++  End of check for hit this turn and element
  -    end if
  -  end do ! end do j = 1, napx
   
     if(dowrite_tracks) then
       call coll_writeTracks2(1)
     end if
   
  -!++  Calculate statistical observables and save into files...
  -  if (n_impact.gt.0) then
  +  if(n_impact > 0) then
  +    ! Calculate statistical observables
       average = sum/n_impact
   
       if(sqsum/n_impact >= average**2) then
  
positive.code:     sigma   = zero
    end if
  
    if(cn_impact(icoll).gt.0) then
      caverage(icoll) = csum(icoll)/cn_impact(icoll)
  
      if((caverage(icoll)**2).gt.(csqsum(icoll)/cn_impact(icoll))) then
        csigma(icoll) = 0
      else
        csigma(icoll) = sqrt(csqsum(icoll)/cn_impact(icoll) - caverage(icoll)**2)
  
negative.code:   integer j
    real(kind=fPrec) average,sigma,zpj,sum,sqsum
  
    n_impact   = 0
    n_absorbed = 0
    sum        = zero
    sqsum      = zero
  
  #ifdef G4COLLIMATION
    do j=1,napx
  