inverse_fit.py

Go to the documentation of this file.

from ROOT import *

filein = TFile("mymccheck_uniform_50GeV.root", "read")
tree = filein.Get("detsimana/tree")
tree.SetBranchStatus("*",0)
tree.SetBranchStatus("nu_eng_mcn",1)

h_energy = TH1D("heng", "heng", 500, 0., 50)
for entry in tree:
  h_energy.Fill(entry.nu_eng_mcn)

hinv = h_energy.Clone()
hloginv = h_energy.Clone()
for i in range(h_energy.GetNbinsX()):
  if h_energy.GetBinContent(i+1) != 0:
    hinv.SetBinContent(i+1, h_energy.Integral()/h_energy.GetBinContent(i+1))
    hloginv.SetBinContent(i+1, TMath.Log(h_energy.Integral()/h_energy.GetBinContent(i+1)))

def fit_log_fcn(x, par):
  energy = x[0]
  #  return p[0] - p[1]*tmath.power(energy, p[2])+p[3]*tmath.power(energy, -p[4])
  def profile(y, p): return p[0]+p[1]*TMath.Power(y,-p[2])+p[3]*TMath.Power(y, -p[4])+p[5]*y
  threshold = 0.1
  if energy > threshold:
      return profile(energy, par)
  else:
      return profile(100., par)
      #  return fthreshold*TMath.Exp(energy-threshold)
      #  return fthreshold-dthreshold*(threshold-energy)

def fit_fcn(x, par):
    return TMath.Exp(fit_log_fcn(x, par))

fit_log = TF1("fit_log", fit_log_fcn, 0, 50, 6)
fit_log.SetLineColor(kRed)
hinv.Fit("fit_log")

#  fit_inv = TF1("fit_inv", fit_fcn, 0.01, 50, 4)
#  fit_inv.SetLineColor(kRed)
#  for i in range(fit_log.GetNumberFreeParameters()):
#      pari = fit_log.GetParameter(i)
#      print pari
#      fit_inv.FixParameter(i, pari)
#  fit_inv.SetNpx(100000000)
fit_log.SetNpx(200000)

#  print "-----------"
hflat = TH1D("hflat", "hflat", 500, 0., 50)
gRandom.SetSeed(0)
#  for i in range(100):
#      print "---------"
#  hflat.FillRandom("fit_log", tree.GetEntries())
for j in range(tree.GetEntries()):
    hflat.Fill(fit_log.GetRandom())
    
#  hflat2 = hflat.Clone()
#  for i in range(hflat.GetNbinsX()):
#      hflat2.SetBinContent(i+1, hflat.GetBinContent(i+1)*h_energy.GetBinContent(i+2))
#  hflat.Multiply(h_energy)

#  hflat.Scale(hinv.Integral()/hflat.Integral())
hflat.Scale(hinv.Integral()/(2500*hflat.Integral()))
hflat.SetLineColor(kRed)
hflat.Multiply(h_energy)
c = TCanvas()
#  hinv.Draw("hist same")
#  hinv.GetXaxis().SetRangeUser(0, 20)
#  fit_log.Draw("same")
hflat.Rebin(5)
hflat.GetYaxis().SetRangeUser(0, 2500)
hflat.GetXaxis().SetRangeUser(0, 50)
hflat.Draw("hist same")
#  hinv.Draw("hist")
#  fit_inv.Draw("same")
#  c.SetLogy()
c.Print("test_fit.pdf")

filein.Close()