Ticket #7383: FitLowStatsAsymChallenge.py

File FitLowStatsAsymChallenge.py, 2.6 KB (added by Anders Markvardsen, 7 years ago)
Line 
1import math
2# arrays to fill
3resx=numpy.zeros(200)
4resy=numpy.zeros(200)
5rese=numpy.zeros(200)
6resChi=numpy.zeros(200)
7res0=numpy.zeros(200)
8res1=numpy.ones(200)
9resNZ=numpy.zeros(200)
10resIC=numpy.zeros(200)
11resICE=numpy.zeros(200)
12a=0.25 # asymmetry to simulate
13NB=1000 # number of (raw) bins
14x12arr=numpy.zeros(2*NB)
15y12arr=numpy.zeros(2*NB)
16e12arr=numpy.zeros(2*NB)
17ws=CreateWorkspace(x12arr,y12arr,e12arr,2,OutputWorkspace="Hello")
18for x in range(200):
19        lam=math.exp((x-75.0)/10.0) # counts per bin, in absence of any signal. Log scale to show detail.
20        lam1=lam*(1+a)
21        lam2=lam*(1-a) # counts per bin in forward and backward banks
22        Xarr=range(NB) # "time"
23        Y1arr=numpy.random.poisson(lam1,NB) # measured forward counts
24        E1arr=numpy.sqrt(Y1arr) # and their errors by the "standard" formula, Fit() will treat error=0 specially
25        Y2arr=numpy.random.poisson(lam2,NB)
26        E2arr=numpy.sqrt(Y2arr)
27        ws.dataX(0)[:]=Xarr
28        ws.dataY(0)[:]=Y1arr
29        ws.dataE(0)[:]=E1arr
30        ws.dataX(1)[:]=Xarr
31        ws.dataY(1)[:]=Y2arr
32        ws.dataE(1)[:]=E2arr
33        AsymmetryCalc(InputWorkspace="Hello",OutputWorkspace="Asym",ForwardSpectra="0",BackwardSpectra="1",alpha=1.0) # re-generated asymmetry
34        (stat,chisq,Covar,params,curves)=Fit(Function="name=FlatBackground,A0=0.1",InputWorkspace="Asym",Output="Asym")
35        print lam," -> ",params.column(1)[0]," +- ",params.column(2)[0]," chisq=",chisq," st=",stat
36        resx[x]=lam
37        resy[x]=params.column(1)[0]
38        rese[x]=params.column(2)[0]
39        resNZ[x]=(len(Y1arr)-numpy.count_nonzero(Y1arr)+len(Y2arr)-numpy.count_nonzero(Y2arr))/(len(Y1arr)+len(Y2arr)+0.0)
40        resChi[x]=chisq
41        DeleteWorkspace("Asym")
42
43        YS1=numpy.sum(Y1arr,dtype=numpy.float) # integral asymmetry of same data for comparison
44        YS2=numpy.sum(Y2arr,dtype=numpy.float)
45        if(YS1+YS2>0):
46                resIC[x]=(YS1-YS2)/(YS1+YS2)
47                resICE[x]=2.0*math.sqrt(YS1*YS2)*(YS1+YS2)**(-1.5)
48        else:
49                resIC[x]=float('NaN')
50                resICE[x]=float('Inf') # no counts, asymmetry completely uncertain!
51
52DeleteWorkspace("Hello")
53CreateWorkspace(resx,resy,rese,1,OutputWorkspace="Summary") # what the fit thought the asymmetry was
54CreateWorkspace(resx,resNZ,res0,1,OutputWorkspace="FractionOfZeros")
55DiffOverErr=(resy-a)/rese
56CreateWorkspace(resx,DiffOverErr,res1,1,OutputWorkspace="NormalisedDifference") # if outside +-1 then systematic errors are significant
57CreateWorkspace(resx,resChi,res0,1,OutputWorkspace="ChiSquared")
58
59CreateWorkspace(resx,resIC,resICE,1,OutputWorkspace="IntegralCounted")
60
61p=plotSpectrum("Summary",0,True)
62l=p.activeLayer()
63l.setAxisScale(Layer.Bottom,0.001,1000000.0,Layer.Log10)
64
65p2=plotSpectrum("NormalisedDifference",0,True)
66l2=p2.activeLayer()
67l2.setAxisScale(Layer.Bottom,0.001,1000000.0,Layer.Log10)