chisquared.py

Go to the documentation of this file.

#!/usr/bin/env python

import sys, os
from os.path import join
from math import sqrt
from HandyFuncs import VerticalRange

if "PYTHONSTARTUP" in os.environ:
    if sys.argv[0]:
        execfile(os.environ["PYTHONSTARTUP"])
elif os.path.exists("rootlogon.py"):
    execfile("rootlogon.py")


hpred = TH1D("hpred", "", 20, 0, 20)
hrand = TH1D("hrand", "", 20, 0, 20)
hpred.GetXaxis().CenterTitle()
hpred.GetYaxis().CenterTitle()

hpred.SetLineColor(kRed)
hrand.SetLineColor(kBlack)

rand = TRandom3(27565)

for i in range(0, 20):
    val = 100*TMath.Exp(-i/10.)
    hpred.Fill(i, val)
    hrand.Fill(i, rand.Poisson(val))
    hpred.SetBinError(i+1, 0)
    hrand.SetBinError(i+1, sqrt(hrand.GetBinContent(i+1)))

hpred.SetMinimum(0)

VerticalRange([hpred,hrand])

c1 = TCanvas("c1","c1")
#hpred.SetLineWidth(3)
#hpred.SetFillColor(kBlue-7)
hpred.Draw("hist")
hrand.Draw("pesame")


def CalcChi2(hmc, hdata, alpha = 1.):
    chi2 = 0
    for i in range(1, hmc.GetNbinsX()+1):
        ei = hmc.GetBinContent(i)*alpha
        oi = hdata.GetBinContent(i)
        sigma = sqrt(ei)
        chi2 += (ei - oi)**2 / sigma**2
    return chi2

chi2 = CalcChi2(hpred, hrand)
Ndof = hpred.GetNbinsX()
p = TMath.Prob(chi2, Ndof)

print "chi2 =", chi2
print "Ndof =", Ndof
print "p =", p

c1.Print("chisq1.png")


halt = TH1D("halt", "", 20, 0, 20)
halt.SetLineColor(kBlack)

for i in range(0, 20):
    val = hpred.GetBinContent(i+1) + 7
    halt.Fill(i, val)
    halt.SetBinError(i+1, sqrt(val))

c2 = TCanvas("c2","c2")
VerticalRange([hpred, halt])
hpred.Draw("hits")
halt.Draw("pesame")

chi2 = CalcChi2(hpred, halt)
Ndof = hpred.GetNbinsX()
p = TMath.Prob(chi2, Ndof)

print "chi2 =", chi2
print "Ndof =", Ndof
print "p =", p

c2.Print("chisq2.png")



hpred2 = hpred.Clone()
hpred2.Scale(2)
hpred3 = hpred.Clone()
hpred3.Scale(0.5)

c3 = TCanvas("c3","c3")

VerticalRange([hpred2, hpred3, hrand])
hpred2.Draw("hist")
hpred3.Draw("same")
hrand.Draw("pesame")

c3.Print("chisq3.png")


print "chi2 2 =", CalcChi2(hpred2, hrand)
print "chi2 0.5 =", CalcChi2(hpred3, hrand)

hpred.Draw("same")
c3.Print("chisq4.png")


hlike = TH1D("hlike", ";#theta;#chi^{2}", 400, 0.8, 1.2)
hlike.GetXaxis().CenterTitle()
hlike.GetYaxis().CenterTitle()

min = 1000
minval = 0

for bin in range(hlike.GetNbinsX()):
    theta = hlike.GetBinLowEdge(bin)
    chi2 = CalcChi2(hpred, hrand, theta)
    hlike.SetBinContent(bin+1, chi2)
    if chi2 < min:
        min = chi2
        minval = theta
    print bin, theta, chi2

c4 = TCanvas("c4", "c4")
hlike.Draw("hist")

c4.Print("chisq5.png")
print min, minval