inverse_fit.py
Go to the documentation of this file.
1 from ROOT import *
2 
3 filein = TFile("mymccheck_uniform_50GeV.root", "read")
4 tree = filein.Get("detsimana/tree")
5 tree.SetBranchStatus("*",0)
6 tree.SetBranchStatus("nu_eng_mcn",1)
7 
8 h_energy = TH1D("heng", "heng", 500, 0., 50)
9 for entry in tree:
10  h_energy.Fill(entry.nu_eng_mcn)
11 
12 hinv = h_energy.Clone()
13 hloginv = h_energy.Clone()
14 for i in range(h_energy.GetNbinsX()):
15  if h_energy.GetBinContent(i+1) != 0:
16  hinv.SetBinContent(i+1, h_energy.Integral()/h_energy.GetBinContent(i+1))
17  hloginv.SetBinContent(i+1, TMath.Log(h_energy.Integral()/h_energy.GetBinContent(i+1)))
18 
19 def fit_log_fcn(x, par):
20  energy = x[0]
21  # return p[0] - p[1]*tmath.power(energy, p[2])+p[3]*tmath.power(energy, -p[4])
22  def profile(y, p): return p[0]+p[1]*TMath.Power(y,-p[2])+p[3]*TMath.Power(y, -p[4])+p[5]*y
23  threshold = 0.1
24  if energy > threshold:
25  return profile(energy, par)
26  else:
27  return profile(100., par)
28  # return fthreshold*TMath.Exp(energy-threshold)
29  # return fthreshold-dthreshold*(threshold-energy)
30 
31 def fit_fcn(x, par):
32  return TMath.Exp(fit_log_fcn(x, par))
33 
34 fit_log = TF1("fit_log", fit_log_fcn, 0, 50, 6)
35 fit_log.SetLineColor(kRed)
36 hinv.Fit("fit_log")
37 
38 # fit_inv = TF1("fit_inv", fit_fcn, 0.01, 50, 4)
39 # fit_inv.SetLineColor(kRed)
40 # for i in range(fit_log.GetNumberFreeParameters()):
41 # pari = fit_log.GetParameter(i)
42 # print pari
43 # fit_inv.FixParameter(i, pari)
44 # fit_inv.SetNpx(100000000)
45 fit_log.SetNpx(200000)
46 
47 # print "-----------"
48 hflat = TH1D("hflat", "hflat", 500, 0., 50)
49 gRandom.SetSeed(0)
50 # for i in range(100):
51 # print "---------"
52 # hflat.FillRandom("fit_log", tree.GetEntries())
53 for j in range(tree.GetEntries()):
54  hflat.Fill(fit_log.GetRandom())
55 
56 # hflat2 = hflat.Clone()
57 # for i in range(hflat.GetNbinsX()):
58 # hflat2.SetBinContent(i+1, hflat.GetBinContent(i+1)*h_energy.GetBinContent(i+2))
59 # hflat.Multiply(h_energy)
60 
61 # hflat.Scale(hinv.Integral()/hflat.Integral())
62 hflat.Scale(hinv.Integral()/(2500*hflat.Integral()))
63 hflat.SetLineColor(kRed)
64 hflat.Multiply(h_energy)
65 c = TCanvas()
66 # hinv.Draw("hist same")
67 # hinv.GetXaxis().SetRangeUser(0, 20)
68 # fit_log.Draw("same")
69 hflat.Rebin(5)
70 hflat.GetYaxis().SetRangeUser(0, 2500)
71 hflat.GetXaxis().SetRangeUser(0, 50)
72 hflat.Draw("hist same")
73 # hinv.Draw("hist")
74 # fit_inv.Draw("same")
75 # c.SetLogy()
76 c.Print("test_fit.pdf")
77 
78 filein.Close()
def fit_log_fcn(x, par)
Definition: inverse_fit.py:19
def fit_fcn(x, par)
Definition: inverse_fit.py:31