chisquared.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 
3 import sys, os
4 from os.path import join
5 from math import sqrt
6 from HandyFuncs import VerticalRange
7 
8 if "PYTHONSTARTUP" in os.environ:
9  if sys.argv[0]:
10  execfile(os.environ["PYTHONSTARTUP"])
11 elif os.path.exists("rootlogon.py"):
12  execfile("rootlogon.py")
13 
14 
15 hpred = TH1D("hpred", "", 20, 0, 20)
16 hrand = TH1D("hrand", "", 20, 0, 20)
17 hpred.GetXaxis().CenterTitle()
18 hpred.GetYaxis().CenterTitle()
19 
20 hpred.SetLineColor(kRed)
21 hrand.SetLineColor(kBlack)
22 
23 rand = TRandom3(27565)
24 
25 for i in range(0, 20):
26  val = 100*TMath.Exp(-i/10.)
27  hpred.Fill(i, val)
28  hrand.Fill(i, rand.Poisson(val))
29  hpred.SetBinError(i+1, 0)
30  hrand.SetBinError(i+1, sqrt(hrand.GetBinContent(i+1)))
31 
32 hpred.SetMinimum(0)
33 
34 VerticalRange([hpred,hrand])
35 
36 c1 = TCanvas("c1","c1")
37 #hpred.SetLineWidth(3)
38 #hpred.SetFillColor(kBlue-7)
39 hpred.Draw("hist")
40 hrand.Draw("pesame")
41 
42 
43 def CalcChi2(hmc, hdata, alpha = 1.):
44  chi2 = 0
45  for i in range(1, hmc.GetNbinsX()+1):
46  ei = hmc.GetBinContent(i)*alpha
47  oi = hdata.GetBinContent(i)
48  sigma = sqrt(ei)
49  chi2 += (ei - oi)**2 / sigma**2
50  return chi2
51 
52 chi2 = CalcChi2(hpred, hrand)
53 Ndof = hpred.GetNbinsX()
54 p = TMath.Prob(chi2, Ndof)
55 
56 print "chi2 =", chi2
57 print "Ndof =", Ndof
58 print "p =", p
59 
60 c1.Print("chisq1.png")
61 
62 
63 halt = TH1D("halt", "", 20, 0, 20)
64 halt.SetLineColor(kBlack)
65 
66 for i in range(0, 20):
67  val = hpred.GetBinContent(i+1) + 7
68  halt.Fill(i, val)
69  halt.SetBinError(i+1, sqrt(val))
70 
71 c2 = TCanvas("c2","c2")
72 VerticalRange([hpred, halt])
73 hpred.Draw("hits")
74 halt.Draw("pesame")
75 
76 chi2 = CalcChi2(hpred, halt)
77 Ndof = hpred.GetNbinsX()
78 p = TMath.Prob(chi2, Ndof)
79 
80 print "chi2 =", chi2
81 print "Ndof =", Ndof
82 print "p =", p
83 
84 c2.Print("chisq2.png")
85 
86 
87 
88 hpred2 = hpred.Clone()
89 hpred2.Scale(2)
90 hpred3 = hpred.Clone()
91 hpred3.Scale(0.5)
92 
93 c3 = TCanvas("c3","c3")
94 
95 VerticalRange([hpred2, hpred3, hrand])
96 hpred2.Draw("hist")
97 hpred3.Draw("same")
98 hrand.Draw("pesame")
99 
100 c3.Print("chisq3.png")
101 
102 
103 print "chi2 2 =", CalcChi2(hpred2, hrand)
104 print "chi2 0.5 =", CalcChi2(hpred3, hrand)
105 
106 hpred.Draw("same")
107 c3.Print("chisq4.png")
108 
109 
110 hlike = TH1D("hlike", ";#theta;#chi^{2}", 400, 0.8, 1.2)
111 hlike.GetXaxis().CenterTitle()
112 hlike.GetYaxis().CenterTitle()
113 
114 min = 1000
115 minval = 0
116 
117 for bin in range(hlike.GetNbinsX()):
118  theta = hlike.GetBinLowEdge(bin)
119  chi2 = CalcChi2(hpred, hrand, theta)
120  hlike.SetBinContent(bin+1, chi2)
121  if chi2 < min:
122  min = chi2
123  minval = theta
124  print bin, theta, chi2
125 
126 c4 = TCanvas("c4", "c4")
127 hlike.Draw("hist")
128 
129 c4.Print("chisq5.png")
130 print min, minval
T sqrt(T number)
Definition: d0nt_math.hpp:156
def CalcChi2(hmc, hdata, alpha=1.)
Definition: chisquared.py:43
def VerticalRange(hists, xrange=(0, 0), ratio=False, forceOne=True, ignoreError=False, maxerr=0.25, absMax=-999, absMin=-999, buffer=0.05, verbose=0, log=False)
Definition: HandyFuncs.py:38