pro mean1,n,ss

;checks influence of systematics on arithmetic mean
;in comparison with "mean_with_systematics" 

;difference to previous code: now, each measurement has a systematic
;component, but individually different

;n:  number of samples for building arithmetic mean: 10,100,1000,...
;ss: number of arithmetic means, to be investigated w.r.t. expectation value
;    and variance  
  
  meanx=5.
  meansys=0.1
;  meansys=0.
  varx=2.^2
  varsys=0.5^2

  distmean1=fltarr(ss)
  distmean2=distmean1
; calculation of ss arithmetic means for n samples
  for i=0,ss-1 do begin
    distx1=randomn(seed,n)*sqrt(varx)+meanx; x-sample
; normally distributed systematics, different for each measurement
    dist_sys=randomn(seed,n)*sqrt(varsys)+meansys

    distx2=distx1+dist_sys; x-sample + individual systematics 
    distmean1(i)=total(distx1)/n; mean for x
    distmean2(i)=total(distx2)/n; mean for x+sys
  endfor

  print,'individual'
  print,'x  :',meanx,'+/-',sqrt(varx)
  print,'sys:',meansys,'+/-',sqrt(varsys)
  print,'mean'
  print,'x    :',  meanx,'+/-',sqrt(varx/n)
  print,'x+sys:',meanx+meansys,'+/-',sqrt(varx/n+varsys/n)
; predicted: E(mean) = E(x)+E(sys)  approx = mean(arithmetic means)
; predicted: Var(mean) =  Var(x)/n + Var(sys)/n 

  print  
  print,'without systematics, should decrease with sqrt(n)'
  print,mean(distmean1),stdev(distmean1)
  print,'with individual systematics, should also decrease with sqrt(n)'
  print,mean(distmean2),stdev(distmean2)

  return
  end