anchor.diff:      sigma   = zero
     end if
   
  -  if(cn_impact(icoll).gt.0) then
  +  if(cn_impact(icoll) > 0) then
       caverage(icoll) = csum(icoll)/cn_impact(icoll)
   
  -    if((caverage(icoll)**2).gt.(csqsum(icoll)/cn_impact(icoll))) then
  +    if(caverage(icoll)**2 > csqsum(icoll)/cn_impact(icoll)) then
         csigma(icoll) = 0
       else
         csigma(icoll) = sqrt(csqsum(icoll)/cn_impact(icoll) - caverage(icoll)**2)
  
positive.code:   integer j
    real(kind=fPrec) average,sigma,zpj,sum,sqsum
  
  !++  Output information:
  !++
  !++  PART_HIT_POS (MAX_NPART)  Hit flag for last hit
  !++  PART_HIT_TURN(MAX_NPART)  Hit flag for last hit
  !++  PART_ABS_POS (MAX_NPART)  Abs flag
  !++  PART_ABS_TURN(MAX_NPART)  Abs flag
  !++  PART_IMPACT  (MAX_NPART)  Impact parameter (0 for inner face)
  !++  PART_INDIV   (MAX_NPART)  Divergence of impacting particles
  !------------------------------------------------------------------------------
  !++  Calculate average impact parameter and save info for all
  !++  collimators. Copy information back and do negative drift.
    n_impact = 0
    n_absorbed = 0
    sum      = zero
    sqsum    = zero
  
  #ifdef G4COLLIMATION
    do j=1,napx
  
negative.code: 
    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(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
  
    if(dowrite_tracks) then
      call coll_writeTracks2(1)
    end if
  
    if(n_impact > 0) then
      ! Calculate statistical observables
      average = sum/n_impact
  
      if(sqsum/n_impact >= average**2) then
  