anchor.diff:    nsig     = cdb_cNSig(icoll)
     c_length = zero
   
  -  ! SR, 23-11-2005: To avoid binary entries in 'amplitude.dat'
     if(firstrun) then
  -
  -    gammax = (one + talphax(ie)**2)/tbetax(ie)
  -    gammay = (one + talphay(ie)**2)/tbetay(ie)
  -
  -    if(iturn == 1) then
  -      sum_ax(ie) = zero
  -      sum_ay(ie) = zero
  -    end if
  -
  -    do j=1,napx
  -      xj  = (xv1(j)-torbx(ie))/c1e3
  -      xpj = (yv1(j)-torbxp(ie))/c1e3
  -      yj  = (xv2(j)-torby(ie))/c1e3
  -      ypj = (yv2(j)-torbyp(ie))/c1e3
  -      pj  = ejv(j)/c1e3
  -
  -      ! DRIFT PART
  -      if(stracki == zero) then
  -        if(iexact) then
  -          zpj = sqrt(one-xpj**2-ypj**2)
  -          xj  = xj + half*c_length*(xpj/zpj)
  -          yj  = yj + half*c_length*(ypj/zpj)
  -        else
  -          xj  = xj + half*c_length*xpj
  -          yj  = yj + half*c_length*ypj
  -        end if
  -      end if
  -
  -      if(part_abs_pos(j) == 0 .and. part_abs_turn(j) == 0) then
  -        nspx = sqrt(abs(gammax*xj**2 + two*talphax(ie)*xj*xpj + tbetax(ie)*xpj**2)/c_emitx0_collgap)
  -        nspy = sqrt(abs(gammay*yj**2 + two*talphay(ie)*yj*ypj + tbetay(ie)*ypj**2)/c_emity0_collgap)
  -        sum_ax(ie)   = sum_ax(ie) + nspx
  -        sqsum_ax(ie) = sqsum_ax(ie) + nspx**2
  -        sum_ay(ie)   = sum_ay(ie) + nspy
  -        sqsum_ay(ie) = sqsum_ay(ie) + nspy**2
  -        nampl(ie)    = nampl(ie) + 1
  -      else
  -        nspx = zero
  -        nspy = zero
  -      end if
  -    end do
  +    call coll_computeStats
     end if
   
   end subroutine collimate_start_collimator
  
positive.code: #endif
  
    integer j
    real(kind=fPrec) gammax,gammay,xdisp,nspxd,xj,xpj,yj,ypj,pj
  
  #ifdef HDF5
    integer hdfturn,hdfpid,hdftyp
  
negative.code: #endif
  
    if(firstrun) then
      call coll_computeStats
    end if
  
    ! Note: Not in firstrun
    if(dowritetracks) then
      gammax = (one + talphax(ie)**2)/tbetax(ie)
      gammay = (one + talphay(ie)**2)/tbetay(ie)
  
      do j=1,napx
        xj  = (xv1(j)-torbx(ie)) /c1e3
        xpj = (yv1(j)-torbxp(ie))/c1e3
  