55

#!/usr/bin/python
# -*- coding: utf8 -*-
import sys
import ROOT

canvas = ROOT.TCanvas("canvas", "Canvas")
canvas.Divide(2,2)
canvas.cd(1)


g = ROOT.TGraphErrors("variant2.txt")

t = []
f = open("variant2.txt")
t = [tuple(map(float, i.split(' '))) for i in f]

func = ROOT.TF1("f", "[0]*([2]/2)^2/((x-[1])^2+([2]/2)^2)", 0, 1)

func.SetParameter(0, 400)
func.SetParameter(1, 0.77)
func.SetParameter(2, 0.005)

g.SetMarkerStyle(20)

#func.Draw()
#raw_input('--->')

g.Draw("AP")
raw_input('--->')

fit_result = g.Fit(func, "SQ")

# Get value A, M, G
print " "
print "A:"
print "VALUE: ", func.GetParameter(0)
print "ERROR: ", func.GetParError(0)
print " "

print "M:"
print "VALUE: ", func.GetParameter(1)
print "ERROR: ", func.GetParError(1)
print " "

print "G:"
print "VALUE: ", func.GetParameter(2)
print "ERROR: ", func.GetParError(2)
print " "

# Correlation between parameters.
print "Correlation(A,M): ", fit_result.Correlation(0, 1)
print "Correlation(A,G): ", fit_result.Correlation(0, 2)
print "Correlation(M,G): ", fit_result.Correlation(0, 1)
print " "

# hi**2
print "Coefficient hi", fit_result.Chi2()
print " "
canvas.Update()


# Define hist.
histA = ROOT.TH1F("histA", "A", 100, 160, 175)
histM = ROOT.TH1F("histM", "M", 100, 0.7818, 0.783)
histG = ROOT.TH1F("histG", "G", 100, -0.00950, -0.00880)

raw_input('--->')

sigma_E = 0.1 / 1000
for i in range(0, 100):
    delta_E = ROOT.gRandom.Gaus(0, sigma_E)
    g2 = ROOT.TGraphErrors()
    for j in range(0, len(t)):
        E  = t[j][0] + delta_E
        Si = t[j][1]
        g2.SetPoint(j, E, Si)
        g2.SetPointError(j, 0, t[j][3])
    fit_result2 = g2.Fit("f", "SQ")
    histA.Fill(func.GetParameter(0))
    histM.Fill(func.GetParameter(1))
    histG.Fill(func.GetParameter(2))

print "\nENERGY CONTRIBUTION"
print "A RMS:", histA.GetRMS()
print "M RMS:", histM.GetRMS()
print "G RMS:", histG.GetRMS()

canvas.cd(2)
histA.Draw()
canvas.cd(3)
histM.Draw()
canvas.cd(4)
histG.Draw()
raw_input('--->')

for i in range(0, 100):
    g2 = ROOT.TGraphErrors()
    for j in range(0, len(t)):
        sigma_epsilon = 0.01 * t[j][0] / 1.02
        delta_epsilon = ROOT.gRandom.Gaus(0, sigma_epsilon)
        E  = t[j][0]
        Si = t[j][1] * (1 + delta_epsilon)
        g2.SetPoint(j, E, Si)
        g2.SetPointError(j, 0, t[j][3])
    fit_result2 = g2.Fit("f", "SQ")
    histA.Fill(func.GetParameter(0))
    histM.Fill(func.GetParameter(1))
    histG.Fill(func.GetParameter(2))

print "\nEFFICIENCY CONTRIBUTION"
print "A RMS:", histA.GetRMS()
print "M RMS:", histM.GetRMS()
print "G RMS:", histG.GetRMS()

canvas.cd(2)
histA.Draw()
canvas.cd(3)
histM.Draw()
canvas.cd(4)
histG.Draw()
raw_input('--->')

sigma_E = 0.1 / 1000
for i in range(0, 100):
    delta_E = ROOT.gRandom.Gaus(0, sigma_E)
    g2 = ROOT.TGraphErrors()
    for j in range(0, len(t)):
        sigma_epsilon = 0.01 * t[j][0] / 1.02
        delta_epsilon = ROOT.gRandom.Gaus(0, sigma_epsilon)
        E  = t[j][0] + delta_E
        Si = t[j][1] * (1 + delta_epsilon)
        g2.SetPoint(j, E, Si)
        g2.SetPointError(j, 0, t[j][3])
    fit_result2 = g2.Fit("f", "SQ")
    histA.Fill(func.GetParameter(0))
    histM.Fill(func.GetParameter(1))
    histG.Fill(func.GetParameter(2))

print "\nTOTAL CONTRIBUTION"
print "A RMS:", histA.GetRMS()
print "M RMS:", histM.GetRMS()
print "G RMS:", histG.GetRMS()

canvas.cd(2)
histA.Draw()
canvas.cd(3)
histM.Draw()
canvas.cd(4)
histG.Draw()

raw_input('finish --->')