Source Code for Module reflectometry.reduction.examples.polcorerr

 1  #!/usr/bin/env python 
 2   
 3  """ 
 4  Program to show the monte carlo error estimate of the polarization correction 
 5  operation, assuming a value and uncertainty for polarizer and flipper 
 6  efficiencies. 
 7  """ 
 8  n = 100000 
 9  err = 5 
10  Ic = 100000 
11  doplot = False 
12   
13   
14  import numpy,math,pylab 
15  import reflectometry.reduction as reflred 
16   
17  eff = reflred.PolarizationEfficiency() 
18   
19  # Seed polarizer/flipper efficiencies from a gaussian distribution 
20  eff.ff = numpy.random.normal(0.95,0.01*err,n) 
21  eff.fp = numpy.random.normal(0.90,0.01*err,n) 
22  eff.rf = numpy.random.normal(0.95,0.01*err,n) 
23  eff.rp = numpy.random.normal(0.90,0.01*err,n) 
24  eff.Ic = numpy.random.normal(Ic,numpy.sqrt(Ic),n) 
25   
26  data = reflred.PolarizedData() 
27  for V,v in [(data.pp,Ic), (data.pm,Ic/5), (data.mp,Ic/5), (data.mm,Ic)]: 
28      V.v = numpy.ones(n)*v 
29      V.variance = V.v   # Variance is poisson variance 
30      V.v = numpy.random.normal(V.v,V.dv)  # Randomize inputs 
31   
32  eff(data)  # Apply polarization efficiency correction to data 
33   
34  for plt,d,label,E in [(221,data.pp,'++',Ic), 
35                        (222,data.pm,'+-',Ic/5), 
36                        (223,data.mp,'-+',Ic/5), 
37                        (224,data.mm,'--',Ic)]: 
38      if doplot: 
39          pylab.subplot(plt) 
40          pylab.hist(d.v) 
41          pylab.legend(['%s %0.2f (%0.2f)'%(label,pylab.mean(d.v),pylab.std(d.v))]) 
42      print "%s measurement uncertainty %.2f, corrected uncertainty %.3f, value %.3f"\ 
43          %(label,math.sqrt(E),pylab.std(d.v),numpy.mean(d.v)) 
44  if doplot: pylab.show() 
45