task5Igor

#!/usr/bin/python
# -*- coding: utf8 -*-

import sys
import ROOT

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

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


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

#for i in xrange(0, t.GetEntries()):
#  g.SetPoint( i, t.GetV1()[i], t.GetV2()[i] )
#  g.SetPointError( i, 0, t.GetV3()[i] )



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

func.SetParameter(0, 1400)
func.SetParameter(1, 0.8)
func.SetParameter(2, 0.01)

g.SetMarkerStyle(20)
g.Draw("AP")
sys.stdin.readline()

ROOT.gStyle.SetOptFit(1111)

fit_result = g.Fit("f", "S")

sys.stdin.readline()

# 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", 50, 1420, 1540)
histM = ROOT.TH1F("histM", "M", 50, 0.783, 0.784)
histG = ROOT.TH1F("histG", "G", 50, 0.0085, 0.0092)

sigma_E = 0.1 / 1000
for i in range(0, 100):
    delta_E = ROOT.gRandom.Gaus(0, sigma_E)
    g2 = ROOT.TGraphErrors("variant2.root")
    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][2])
    fit_result2 = g2.Fit("f", "QS")
    histA.Fill(func.GetParameter(0))
    histM.Fill(func.GetParameter(1))
    histG.Fill(func.GetParameter(2))

print ""
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()
sys.stdin.readline("------>")
sys.stdin.readline("------>")