#!/usr/bin/python

from ROOT import TFile, TH1D, TF1, gPad, TCanvas, TPad, TLatex, SetOwnership, TDirectory, gStyle, TH1
from ROOT import kFullCircle, TGraphErrors, TEllipse, TLine, Double, gROOT, TGaxis, kBlue, TArrow
from uncertainties import *
from uncertainties.umath import sqrt, exp
from physics_calcs import *
from math import pi
from os import path
import numpy as np

Asm = -.5706
Pol = .9913
saves = []

kWriteDelete = 4

def find_nearest(array,value):
    idx  = (np.abs(array-value)).argmin()
    idx2 = (np.abs(array-value)).argmax()
    #return [array[idx], array[idx2]]
    return [idx, idx2]

def FierzAnalysis() :
    bar = array('d', np.linspace(-0.03, 0.03, 1200))
    c = array('d', [])
    hname = "Physics_narrowcut_summed"
    #mean_e_loss = 198
    eloss_file = TFile("eloss_profile.root")
    eloss_profile = eloss_file.Get("up_pfx")
    
    for b in bar :
        comp_d = HCompare(dfilename = "asymmetry_D.root",
                          sfilename = "asymmetry_polBoth_M275_Be229_S295_emstandard_opt3_cfs1um_fr0.002_AB_gA0.923.root",
                          hname = hname)
        comp_b = HCompare(dfilename = "asymmetry_B.root",
                          sfilename = "asymmetry_polBoth_M275_Be229_S295_emstandard_opt3_cfs1um_fr0.002_AB_gA0.923.root",
                          hname = hname)
        comp_c = HCompare(dfilename = "asymmetry_C.root",
                          sfilename = "asymmetry_polBoth_M275_Be229_S295_emstandard_opt3_cfs1um_fr0.002_AB_gA0.923.root",
                          hname = hname)
        chi2 = 0
        #for comp in [comp_b, comp_c, comp_d] :
        for comp in [comp_d] :
            comp.AddFierzTerm(b, eloss_profile)
            comp.residuals = GetChi2Hists(comp.data, comp.simu)[2]
            chi2 = chi2 + comp.Chi2InRange(500, 4500)[0]
        
            
        c.append(chi2)
    g = r.TGraph(len(bar), bar, c)
    g.Draw("AL")
    gPad.Update()
    
    minimum = min(c)
    #minimum = c[599]
    yerr = minimum + 1.0
    foundlow = False
    foundhigh = False
    low = 0
    high = 0
    for i in range(len(bar)) :
        if (foundlow and foundhigh) :
            continue
        #print "%d\t%f" % (i, c[i])
        if (not foundlow) :
            if (c[i] < yerr) :
                low = bar[i]
                foundlow = True
                print "Foudn low when i = ", i
        else :
            if (c[i] > yerr) :
                high = bar[i-1]
                foundhigh = True
                print "Found high when i = ", i
    print yerr
    print low
    print high
    fitb = ufloat((high+low)/2., (high-low)/2.)
    print fitb
    raw_input("enter to end")

def ShiftTCutG(cut, xshift, yshift) :
    newcut = cut.Clone()

    for i in range(newcut.GetN()) :
        x = Double(0)
        y = Double(0)
        newcut.GetPoint(i, x, y)
        y = y + yshift
        x = x + xshift
        newcut.SetPoint(i, x, y)
    return newcut

def AsymmetryT_thesis() :
    df = "asymmetry_B.root"
    df = TFile(df)
    time = []
    #for dir in ["All_scint", "Scint*eMCP", "Scint*BB1", "Scint*BB1*eMCPcuts"] :
    for dir in ["All_scint", "Scint*BB1", "Scint*eMCP", "Scint*BB1*eMCPcuts"] :
        hup_plus = df.Get("Counts/" + dir + "/Upper_BB1_Sigma_plus_time")
        hup_mins = df.Get("Counts/" + dir + "/Upper_BB1_Sigma_minus_time")
        hlo_plus = df.Get("Counts/" + dir + "/Lower_BB1_Sigma_plus_time")
        hlo_mins = df.Get("Counts/" + dir + "/Lower_BB1_Sigma_minus_time")
        temp = PlotAsymE([hup_mins, hlo_plus], [hup_plus, hlo_mins], xax = "time")
        temp.Draw("HIST")
        gPad.Update()
        time.append(gPad.GetPrimitive("Physics"))
        time[len(time)-1].Draw("AP")
        
    c = TCanvas("c", "c", 1100, 600)
    c.cd()
    gPad.SetBottomMargin(0.6*gPad.GetBottomMargin())
    time[0].Draw("E1X0 0")
    time[0].GetFunction("Fit").SetBit(512)
    time[0].GetYaxis().SetRangeUser(-0.60, 0.05)
    time[0].GetXaxis().SetRangeUser(0, 4100)
    time[0].GetXaxis().SetTitleSize(0.06)
    time[0].GetXaxis().SetLabelSize(0.05)
    time[0].SetMarkerSize(0.7)
    time[0].SetTitle("All scintillator")
    time[1].SetTitle(" w/ eMCP coincidence")
    time[2].SetTitle(" w/ BB1  coincidence")
    time[3].SetTitle(" w/ eMCP and BB1 coincidence")
    for i in range(1, len(time)) :
        time[i].GetFunction("Fit").SetBit(512)
        time[i].SetLineColor(mycol[i])
        time[i].SetMarkerColor(mycol[i])
        time[i].SetMarkerStyle(mymark[i])
        time[i].SetMarkerSize(0.7)
        time[i].Draw("E1X0 SAME 0")
    gPad.Update()
    leg = gPad.BuildLegend(0.56, 0.31, 0.94, 0.68)
    leg.SetTextSize(0.045)
    leg.SetShadowColor(0)
    leg.SetTextFont(132)
    leg.Draw()
    gPad.Update()
    ch = raw_input("Enter \"SAVE\" to save as .eps for thesis\t")
    if (ch == "SAVE") :
        c.SaveAs("~/Documents/Reports/Images/asymmetry_T_overlay.eps")
        
def AsymmetryOverlayForThesis() :
    for s in ["B", "C", "D"] :
        
        df = "asymmetry_%s.root" % s
        sf = ""
        if (s == "B") :
            sf = "asymmetry_polBoth_M275_Be229_S295_emstandard_opt3_cfs1um_fr0.002_AB_gA0.923.root"
        else :
            sf = "asymmetry_polBoth_M275_Be229_S295_emstandard_opt3_cfs1um_fr0.002_CDE_gA0.923.root"

        c = HCompare(dfilename = df, sfilename = sf, hname = "Physics_narrowcut")
        c.Draw()
        result = c.Chi2InRange(400, 4500)
        text = "#chi^{2}/%d = %1.2f" % (result[1], result[0]/result[1])
        l = TLatex(3000, -0.46, text)
        l.SetTextSize(0.07)
        l.SetTextFont(132)
        l.Draw()
        text = "Set %s" % s
        l2 = TLatex(3000, -0.44, text)
        l2.SetTextSize(0.07)
        l2.SetTextFont(132)
        l2.Draw()
        gPad.Update()
        ch = raw_input("Enter \"SAVE\" to save as .eps for thesis\t") 
        if (ch == "SAVE") :
            outf = "~/Documents/Reports/Images/asymmetry_overlay_%s.eps" % s
            gPad.GetCanvas().SaveAs(outf)


def DrawTOFcuts_thesis(hexcut = False) :
    gStyle.SetOptFit(0)
    f = 0
    if (hexcut) :
        f = TFile("asymmetry_B_hex5sig.root")
    else :
        f = TFile("asymmetry_B.root")
    print f.GetName()
    global saves
    cfname = ["broadcut.root", "narrowcut.root", "lincut.root", "verynarrowcut.root"]
    can = TCanvas("can", "can", 1100, 556)
    #can.Divide(2, 1)
    can.cd()
    can2 = TCanvas("can2", "can2", 1100, 556)
    gPad.SetBottomMargin(0.6*gPad.GetBottomMargin())
    gPad.SetLeftMargin(1.05*gPad.GetLeftMargin())
    hup = 0
    if (hexcut) :
        hup = f.Get("Counts/Scint*BB1*HEXpos/Upper_BB1_Sigma_minus_Energy_v_SOE_tof")
        hup.Add(f.Get("Counts/Scint*BB1*HEXpos/Upper_BB1_Sigma_plus_Energy_v_SOE_tof"))
    else :
        hup = f.Get("Counts/Scint*BB1*eMCP/Upper_BB1_Sigma_minus_Energy_v_SOE_tof")
        hup.Add(f.Get("Counts/Scint*BB1*eMCP/Upper_BB1_Sigma_plus_Energy_v_SOE_tof"))

    proj = hup.ProjectionY("proj", hup.GetXaxis().FindBin(400), hup.GetXaxis().FindBin(5200))
    title = "Counts / %1.1f ns" % proj.GetBinWidth(10)
    proj.GetYaxis().SetTitle(title)

    hup.Draw("COL")
    hup.SetDirectory(0)
    saves.append(hup)
    hup.GetXaxis().SetRangeUser(0., 5200.)
    hup.GetYaxis().SetRangeUser(2., 16.)

    i = 0
    tof = 9.867
    sig = 0.712
    for c in cfname :
        rf = TFile(c)
        cut = rf.Get("CUTG")
        cut.SetName(c)
        if (c == "narrowcut.root") :
            cut.SetLineWidth(4)
            cut.SetLineStyle(5)
        cut.Draw("SAME")
        
        cut.SetLineColor(mycol[i])
        i = i + 1

    i = 0
    alllines = []
    for i in range(1, 6) :
        min = tof - (i*sig)
        max = tof + (i*sig)
        for m in [min, max] :
            l = TLine(0, m, 5200, m)
            l.SetLineWidth(2)
            l.SetLineStyle((i%2)+1)
            l.Draw()
            alllines.append(l)
    gPad.Update()
    ch = raw_input("Enter \"SAVE\" to save as .eps for thesis")
    if (ch == "SAVE") :
        can.SaveAs("~/Documents/Reports/Images/soe_2d.eps")
    #can.cd(2)
    can.Update()
    can2.cd()
    gPad.SetLeftMargin(1.05*gPad.GetLeftMargin())

    proj.Draw("HIST")
    saves.append(proj)
    proj.SetDirectory(0)
    proj.GetXaxis().SetRangeUser(-10., 110.)
    gPad.SetLogy(1)

    for fitr in [[4., 15.], [19, 30], [35, 42]] :
        x0guess = (fitr[1]+fitr[0])/2.
        bmin = proj.FindBin(fitr[0])
        bmax = proj.FindBin(fitr[1])
        fit = 0

        if (fitr[0] < 32) :
            #fit = TF1("fit", "[0]*[1]/(((x-[2])*(x-[2]))+([1]*[1]))+expo(3)",
            #          fitr[0], fitr[1])
            #fit.SetParameters(0.5*proj.GetBinContent(proj.FindBin(x0guess)),
            #                  0.2, x0guess, 2.2, 0.19)
            fit = TF1("fit", "gaus(0)+expo(3)", fitr[0], fitr[1])
            fit.SetParameters(proj.GetBinContent(proj.FindBin(x0guess)),
                              x0guess, 0.7, 3., 0.16)
            fit.SetParLimits(4, 0., 10.)
            print fitr[0]
            if (fitr[0] < 10) :
                fit.FixParameter(3, 2.60*log(2.1)) # Set B
                fit.FixParameter(4, 0.187) # Set B
                #fit.FixParameter(3, 0.083)
                #fit.FixParameter(4, 0.248)
                #Fit.FixParameter(3, 0.72) # Set C
                #fit.FixParameter(4, 0.147) # Set C
                #fit.FixParameter(3, 3.04) # Set D
                #fit.FixParameter(4, 0.14) # Set D
            elif (fitr[0]< 32) :
                print "DOING IT"
                fit.FixParameter(3, 7.483)
                fit.FixParameter(4, -0.116)
        else :
            fit = TF1("fit", "gaus(0)+expo(3)", fitr[0], fitr[1])
            fit.SetParameters(proj.GetBinContent(proj.FindBin(x0guess)),
                              x0guess, 0.7, 2.2, -0.19)
            fit.SetParLimits(4, -10, 0)
        proj.Fit(fit, "R")
        #fit.Draw("LSAME")

        juste = TF1("juste" + str(fitr[0]), "expo(0)", fitr[0], fitr[1])
        juste.SetParameter(0, fit.GetParameter(3))
        juste.SetParameter(1, fit.GetParameter(4))
        juste.SetLineWidth(2)
        juste.SetLineStyle(2)
        juste.Draw("LSAME")
        saves.append(juste)
        gPad.Update()
        
    gPad.Update()
    print gPad.GetUymax()
    l1 = TLine(28.0, gPad.GetUymin(), 28.0, pow(10., gPad.GetUymax()))
    l1.Draw()
    saves.append(l1)
    l2 = TLine(35.0, gPad.GetUymin(), 35.0, pow(10., gPad.GetUymax()))
    l2.Draw()
    saves.append(l2)
    t1 = TLatex(33.2, 263.5, "Background region")
    t1.SetTextFont(132)
    t1.SetTextAngle(90.)
    t1.Draw()
    saves.append(t1)
    l3 = TLine(42, gPad.GetUymin(), 42, pow(10., gPad.GetUymax()))
    l3.Draw()
    saves.append(l3)
    l4 = TLine(80, gPad.GetUymin(), 80, pow(10., gPad.GetUymax()))
    l4.Draw()
    saves.append(l4)
    t2 = TLatex(42.5, 1000., "Background")
    t3 = TLatex(49.3, 400., "region")
    for t in [t2, t3] :
        t.SetTextFont(132)
        t.Draw()
        saves.append(t)

    simulated = r.TGraph("tof_spec_levinger_150vpcm_poisson.out")
    gint = 0.
    for i in range(simulated.GetN()) :
        gint = gint + simulated.GetY()[i]

    scale = proj.Integral(proj.FindBin(5), proj.FindBin(15.))
    print "Data integral: %f" % scale
    scale = scale*.012/gint/2.
    print "Scale factor: %f" % scale
    print "Bin width = %f" % proj.GetBinWidth(1)
    for i in range(simulated.GetN()) :
        simulated.GetY()[i] = simulated.GetY()[i]*scale
        simulated.GetX()[i] = simulated.GetX()[i]+1.5
    gint = 0.
    for i in range(simulated.GetN()) :
        gint = gint + simulated.GetY()[i]
    print "Simulated integral: %f" % gint
    simulated.SetMarkerSize(0)
    simulated.SetLineStyle(3)
    simulated.SetLineColor(kBlue)
    simulated.SetLineWidth(2)
    saves.append(simulated)
    simulated.Draw("L")

    # FROM HERE
    simulated_h4 = r.TGraph("tof_spec_levinger_150vpcm_poisson_H4_both.out")
    gint = 0.
    for i in range(simulated_h4.GetN()) :
        gint = gint + simulated_h4.GetY()[i]

    scale = proj.Integral(proj.FindBin(5), proj.FindBin(15.))
    print "Data integral: %f" % scale
    scale = scale*.00034/gint/2.
    print "Scale factor: %f" % scale
    print "Bin width = %f" % proj.GetBinWidth(1)
    for i in range(simulated_h4.GetN()) :
        simulated_h4.GetY()[i] = simulated_h4.GetY()[i]*scale
        simulated_h4.GetX()[i] = simulated_h4.GetX()[i]+1.5
    gint = 0.
    for i in range(simulated_h4.GetN()) :
        gint = gint + simulated_h4.GetY()[i]
    print "Simulated_H4 integral: %f" % gint
    simulated_h4.SetMarkerSize(0)
    simulated_h4.SetLineStyle(3)
    simulated_h4.SetLineColor(kBlue)
    simulated_h4.SetLineWidth(2)
    #simulated_h4.Draw("L")
    saves.append(simulated_h4)
    # TO HERE


    # FROM HERE
    simulated_h3 = r.TGraph("tof_spec_levinger_150vpcm_poisson_H3_fracp24.out")
    gint = 0.
    for i in range(simulated_h3.GetN()) :
        gint = gint + simulated_h3.GetY()[i]

    scale = proj.Integral(proj.FindBin(5), proj.FindBin(15.))
    print "Data integral: %f" % scale
    scale = scale*.00016/gint/2.
    print "Scale factor: %f" % scale
    print "Bin width = %f" % proj.GetBinWidth(1)
    for i in range(simulated_h3.GetN()) :
        simulated_h3.GetY()[i] = simulated_h3.GetY()[i]*scale
        simulated_h3.GetX()[i] = simulated_h3.GetX()[i]+1.5
    gint = 0.
    for i in range(simulated_h3.GetN()) :
        gint = gint + simulated_h3.GetY()[i]
    print "Simulated_H3 integral: %f" % gint
    simulated_h3.SetMarkerSize(0)
    simulated_h3.SetLineStyle(3)
    simulated_h3.SetLineColor(kBlue)
    simulated_h3.SetLineWidth(2)
    #simulated_h3.Draw("L")
    saves.append(simulated_h3)
    # TO HERE

    # FROM HERE
    simulated_h2 = r.TGraph("tof_spec_levinger_150vpcm_poisson_H2.out")
    gint = 0.
    for i in range(simulated_h2.GetN()) :
        gint = gint + simulated_h2.GetY()[i]

    scale = proj.Integral(proj.FindBin(5), proj.FindBin(15.))
    print "Data integral: %f" % scale
    scale = scale*.00016/gint/2.
    print "Scale factor: %f" % scale
    print "Bin width = %f" % proj.GetBinWidth(1)
    for i in range(simulated_h2.GetN()) :
        simulated_h2.GetY()[i] = simulated_h2.GetY()[i]*scale
        simulated_h2.GetX()[i] = simulated_h2.GetX()[i]+1.5
    gint = 0.
    for i in range(simulated_h2.GetN()) :
        gint = gint + simulated_h2.GetY()[i]
    print "Simulated_H2 integral: %f" % gint
    simulated_h2.SetMarkerSize(0)
    simulated_h2.SetLineStyle(3)
    simulated_h2.SetLineColor(kBlue)
    simulated_h2.SetLineWidth(2)
    #simulated_h2.Draw("L")
    saves.append(simulated_h2)
    # TO HERE

    hsum_x = array('d', [0]*simulated_h2.GetN())
    hsum_y = array('d', [0]*simulated_h2.GetN())
    for i in range(simulated_h2.GetN()) :
        hsum_x[i] = simulated_h2.GetX()[i]
        hsum_y[i] = simulated_h4.GetY()[i] + simulated_h3.GetY()[i] + simulated_h2.GetY()[i]
    hsum = TGraph(simulated_h2.GetN(), hsum_x, hsum_y)
    hsum.SetMarkerSize(0)
    hsum.SetLineStyle(3)
    hsum.SetLineColor(kBlue)
    hsum.SetLineWidth(2)
    hsum.Draw("L")
    saves.append(hsum)
    gPad.Update()
    
    ch = raw_input("Enter \"SAVE\" to save for thesis\t")
    if (ch == "SAVE") :
        can.SaveAs("~/Documents/Reports/Images/soe_cuts_2d1d.eps")

def DrawTOFcuts_thesis_A() :
    gStyle.SetOptFit(0)
    f = TFile("asymmetry_A.root")
    global saves
    cfname = ["broadcut.root", "narrowcut.root", "lincut.root", "verynarrowcut.root"]
    can = TCanvas("can", "can", 1100, 556)
    #can.Divide(2, 1)
    can.cd()

    gPad.SetBottomMargin(0.6*gPad.GetBottomMargin())
    gPad.SetLeftMargin(1.05*gPad.GetLeftMargin())
    hup = f.Get("Counts/Scint*BB1*eMCP/Upper_BB1_Sigma_minus_Energy_v_SOE_tof")
    hup.Add(f.Get("Counts/Scint*BB1*eMCP/Upper_BB1_Sigma_plus_Energy_v_SOE_tof"))

    proj = hup.ProjectionY("proj", hup.GetXaxis().FindBin(400), hup.GetXaxis().FindBin(5200))
    title = "Counts / %1.1f ns" % proj.GetBinWidth(10)
    proj.GetYaxis().SetTitle(title)

    saves.append(hup)
    hup.GetXaxis().SetRangeUser(0., 5200.)
    hup.GetYaxis().SetRangeUser(2., 16.)


    proj.Draw("HIST")
    saves.append(proj)
    proj.SetDirectory(0)
    proj.GetXaxis().SetRangeUser(-10., 110.)
    gPad.SetLogy(1)

    ef = TF1("ef", "expo(0)", 9., 18.)
    ef.SetParameters(2.41, 0.15)
    ef.SetLineWidth(2)
    ef.SetLineStyle(2)
    ef.Draw("L SAME")
    gPad.Update()
        
    print "Data: %f" % (proj.Integral(proj.FindBin(11.1), proj.FindBin(17.1)))
    print "Back: %f" % ef.Integral(11.1, 17.1)
    sn = (proj.Integral(proj.FindBin(11.1), proj.FindBin(17.1))-ef.Integral(11.1, 17.1))/ef.Integral(11.1, 17.1)
    sntext = "Signal-to-noise: %2.0f" % sn
    print sntext
    c = (1+sn)/sn
    ctext = "Correction: %1.4f" % c
    t1 = TLatex(25, 1500, "Set A")
    t2 = TLatex(25, 600, sntext)
    t3 = TLatex(25, 300, ctext)
    for t in [t1, t2, t3] :
        t.SetTextFont(132)
        t.Draw()
    ch = raw_input("Enter \"SAVE\" to save for thesis\t")
    if (ch == "SAVE") :
        can.SaveAs("~/Documents/Reports/Images/soe_cuts_setA.eps")
        
def DrawTOFcuts() :
    f = TFile("asymmetry_B.root")

    cfname = ["broadcut.root", "narrowcut.root", "lincut.root", "verynarrowcut.root"]
    can = TCanvas()
    can.Divide(1, 2)
    can.cd(1)
    hup = f.Get("Counts/Scint*BB1*eMCP/Upper_BB1_Sigma_minus_Energy_v_SOE_tof")
    hup.Add(f.Get("Counts/Scint*BB1*eMCP/Upper_BB1_Sigma_plus_Energy_v_SOE_tof"))
    
    hup.Draw("COL")
    
    i = 0
    tof = 9.867
    sig = 0.712
    for c in cfname :
        rf = TFile(c)
        cut = rf.Get("CUTG")
        cut.SetName(c)
        cut.Draw("SAME")
        cut.SetLineColor(mycol[i])
        i = i + 1

    i = 0
    alllines = []
    for i in range(1, 6) :
        min = tof - (i*sig)
        max = tof + (i*sig)
        for m in [min, max] :
            l = TLine(0, m, 5200, m)
            l.SetLineWidth(2)
            l.SetLineStyle((i%2)+1)
            l.Draw()
            alllines.append(l)
    hlo = f.Get("Counts/Scint*BB1*eMCP/Lower_BB1_Sigma_minus_Energy_v_SOE_tof")
    hlo.Add(f.Get("Counts/Scint*BB1*eMCP/Lower_BB1_Sigma_plus_Energy_v_SOE_tof"))
    
    can.cd(2)
    hlo.Draw("COL")
    
    i = 0
    tof = 10.163
    sig = 0.810
    for c in cfname :
        rf = TFile(c)
        cut = rf.Get("CUTG")
        cut.SetName(c)
        cut.Draw("SAME")
        cut.SetLineColor(mycol[i])
        i = i + 1

    i = 0

    for i in range(1, 6) :
        min = tof - (i*sig)
        max = tof + (i*sig)
        for m in [min, max] :
            l = TLine(0, m, 5200, m)
            l.SetLineWidth(2)
            l.SetLineStyle((i%2)+1)
            l.Draw()
            alllines.append(l)
        

    gPad.Update()
    raw_input("Enter to end")

def AsymmetryTOFmap(filename) :
    f = TFile(filename)

    hup_mins_2d = f.Get("Counts/Scint*BB1*eMCP/Upper_BB1_Sigma_minus_Energy_v_SOE_tof")
    hlo_plus_2d = f.Get("Counts/Scint*BB1*eMCP/Lower_BB1_Sigma_plus_Energy_v_SOE_tof")
    hup_plus_2d = f.Get("Counts/Scint*BB1*eMCP/Upper_BB1_Sigma_plus_Energy_v_SOE_tof")
    hlo_mins_2d = f.Get("Counts/Scint*BB1*eMCP/Lower_BB1_Sigma_minus_Energy_v_SOE_tof")

    asym = TH2D("asym", "asym", 250, 0, 10000, 120, -10, 50)
    
    for x in range(asym.GetNbinsX()) :
        for y in range(asym.GetNbinsY()) :
            xval = asym.GetXaxis().GetBinCenter(x)
            yval = asym.GetYaxis().GetBinCenter(y)

            um = hup_mins_2d.GetBinContent(hup_mins_2d.FindBin(xval, yval))
            um = ufloat(um, sqrt(um))
            
            up = hup_plus_2d.GetBinContent(hup_plus_2d.FindBin(xval, yval))
            up = ufloat(up, sqrt(up))
    
            lm = hlo_mins_2d.GetBinContent(hlo_mins_2d.FindBin(xval, yval))
            lm = ufloat(lm, sqrt(lm))
            
            lp = hlo_plus_2d.GetBinContent(hlo_plus_2d.FindBin(xval, yval))
            lp = ufloat(lp, sqrt(lp))

            if (up.nominal_value == 0 or lm.nominal_value == 0) :
                asym.SetBinContent(x, y, 0.)
                asym.SetBinError(x, y, 0.)
            else :
                s = sqrt(um*lp/up/lm)
                A = (1.-s)/(1.+s)
                asym.SetBinContent(x, y, A.nominal_value)
                asym.SetBinError(x, y, A.std_dev)
    asym.Draw("COLZ")
    asym.GetZaxis().SetRangeUser(-0.75, 0)
    raw_input("Enter to end")

def GetScalingFactor(modfname) :
    canonical_fname = "asymmetry_polBoth_M275_Be229_S295_emstandard_opt3_cfs1um_fr0.002_AB_gA0.907.root"
    hname = "Physics_scint_energy_rmax15.5"
    
    Abeta  = -0.56800
    dAbeta =  0.00135
    canonical_file = TFile(canonical_fname)
    modfile = TFile(modfname)
    
    chist = canonical_file.Get(hname)
    mhist = modfile.Get(hname)

    d = chist.Clone()
    d.SetDirectory(0)
    d.Divide(mhist)
    d.Draw("E1X0")

    f = TF1("f", "pol0(0)", 400, 4500)
    d.Fit(f, "R")
    
    #f.Draw("L SAME")
    d.GetYaxis().SetRangeUser(0.99, 1.01)
    scale = f.GetParameter(0)
    print "scale = %f" % scale

    dA = abs(Abeta - Abeta/scale)
    print "Orig  A_beta = %1.5f" % Abeta
    print "New   A_beta = %1.5f" % (Abeta/scale)
    print "Delta A_beta = %1.5f" % dA
    if (abs(1.0-scale) < f.GetParError(0)) :
        print "NOTE: consistent with no change"
    raw_input("Enter to end")
    canonical_file.Close()
    modfile.Close()

def AsymmetryWithTCut(hfilename, cfilename, htype = "Energy") :
    hf = TFile(hfilename, "UPDATE")
    cf = TFile(cfilename)

    cut = cf.Get("CUTG")
    cname = "soecut%s" % cfilename[:-5].lower()
    print cname
    cut.SetName(cname)
    opt = "[%s]" % cname
    hup_mins_2d = hf.Get("Counts/Scint*BB1*HEXpos/Upper_BB1_Sigma_minus_" + htype  + "_v_SOE_tof")
    hlo_plus_2d = hf.Get("Counts/Scint*BB1*HEXpos/Lower_BB1_Sigma_plus_"  + htype  + "_v_SOE_tof")
    hup_plus_2d = hf.Get("Counts/Scint*BB1*HEXpos/Upper_BB1_Sigma_plus_"  + htype  + "_v_SOE_tof")
    hlo_mins_2d = hf.Get("Counts/Scint*BB1*HEXpos/Lower_BB1_Sigma_minus_" + htype  + "_v_SOE_tof")

    hup_mins_2d.Draw("COL")
    hup_mins_2d.GetYaxis().SetRangeUser(-10, 90)
    cut.Draw("SAME")
    gPad.Update()

    hup_mins = hup_mins_2d.ProjectionX("hup_mins", 0, -1, opt)
    hup_plus = hup_plus_2d.ProjectionX("hup_plus", 0, -1, opt)
    hlo_mins = hlo_mins_2d.ProjectionX("hlo_mins", 0, -1, opt)
    hlo_plus = hlo_plus_2d.ProjectionX("hlo_plus", 0, -1, opt)


    print hup_mins.GetEntries()
    print hup_plus.GetEntries()
    print hlo_mins.GetEntries()
    print hlo_plus.GetEntries()
    
    cphysics = PlotAsymE([hup_mins, hlo_plus], [hup_plus, hlo_mins],
                         xax = "scint_energy")
    cphysics.SetName("Physics_%s" % cfilename[:-5])
    phys_hist = cphysics.GetPrimitive("pup").GetPrimitive("Physics")
    phys_hist.SetName("Physics_%s" % cfilename[:-5])
    
    hf.cd()
    phys_hist.Write(phys_hist.GetName(), kWriteDelete)
    
def SNandAsymForAllPeak(fname) :
    f = TFile(fname)

    det = "Upper"
    s = "Counts/Scint*BB1*HEXpos/%s_BB1_Sigma_minus_Energy_v_SOE_tof" % det
    h2d = f.Get(s)
    s = "Counts/Scint*BB1*HEXpos/%s_BB1_Sigma_plus_Energy_v_SOE_tof" % det
    h2d.Add(f.Get(s))
    
    det = "Lower"
    s = "Counts/Scint*BB1*HEXpos/%s_BB1_Sigma_minus_Energy_v_SOE_tof" % det
    h2d_2 = f.Get(s)
    s = "Counts/Scint*BB1*HEXpos/%s_BB1_Sigma_plus_Energy_v_SOE_tof" % det
    h2d_2.Add(f.Get(s))
    
    htemp = h2d.ProjectionY("up", h2d.GetXaxis().FindBin(400),
                            h2d.GetXaxis().FindBin(5100))
    htemp_2 = h2d_2.ProjectionY("lo", h2d_2.GetXaxis().FindBin(400),
                                h2d_2.GetXaxis().FindBin(5100))
    htemp.Add(htemp_2)
    htemp.Draw("HIST")
    htemp.SetDirectory(0)
    saves.append(htemp)
    for fitr in [[4., 15.], [19, 30], [35, 42]] :
        x0guess = (fitr[1]+fitr[0])/2.
        bmin = htemp.FindBin(fitr[0])
        bmax = htemp.FindBin(fitr[1])
        fit = 0

        if (fitr[0] < 19) :
            #fit = TF1("fit", "[0]*[1]/(((x-[2])*(x-[2]))+([1]*[1]))+expo(3)",
            #          fitr[0], fitr[1])
            #fit.SetParameters(0.5*htemp.GetBinContent(htemp.FindBin(x0guess)),
            #                  0.2, x0guess, 2.2, 0.19)
            fit = TF1("fit", "gaus(0)+expo(3)", fitr[0], fitr[1])
            fit.SetParameters(htemp.GetBinContent(htemp.FindBin(x0guess)),
                              x0guess, 0.7, 3., 0.16)
            fit.SetParLimits(4, 0., 10.)

            if (fitr[0] < 10) :
                #fit.FixParameter(3, 2.60) # Set B (no pos)
                #fit.FixParameter(4, 0.187) # Set B (no pos)
                #fit.FixParameter(3, 2.60*.5) # Set B (pos cut)
                #fit.FixParameter(4, 0.187) # Set B (pos cut)

                #fit.FixParameter(3, 0.72) # Set C (no pos)
                #fit.FixParameter(4, 0.147) # Set C (no pos) 
                #fit.FixParameter(3, 0.72*.1) # Set C (pos cut)
                #fit.FixParameter(4, 0.147) # Set C (pos cut) 

                #fit.FixParameter(3, 3.04) # Set D (no pos) 
                #fit.FixParameter(4, 0.14) # Set D (no pos)
                fit.FixParameter(3, 3.04*.6) # Set D (pos cut)
                fit.FixParameter(4, 0.14) # Set D (pos cut)
        else :
            fit = TF1("fit", "gaus(0)+expo(3)", fitr[0], fitr[1])
            fit.SetParameters(htemp.GetBinContent(htemp.FindBin(x0guess)),
                              x0guess, 0.7, 2.2, -0.19)
            fit.SetParLimits(4, -10, 0)
        htemp.Fit(fit, "R")
        fit.Draw("LSAME")
        gPad.Update()

        A = ufloat(fit.GetParameter(3), fit.GetParError(3))
        B = ufloat(fit.GetParameter(4), fit.GetParError(4))

        if (fitr[0] == 4.) :
            x0 = 10.15
            mu = 0.75
            
            #A = ufloat(fit.GetParameter(3), 0.18) # Set B
            #B = ufloat(fit.GetParameter(4), 0.013) # Set B
            #A = ufloat(fit.GetParameter(3), 0.30) # Set C
            #B = ufloat(fit.GetParameter(4), 0.023) # Set C
            A = ufloat(fit.GetParameter(3), 0.13) # Set D
            B = ufloat(fit.GetParameter(4), 0.010) # Set D

            for nsig in [0.50, 0.75, 1.0, 2.0, 3.0, 4.0, 5.0] :
                range_min = (x0 - (nsig*mu))
                range_max = (x0 + (nsig*mu))

                #range_min = 8.
                #range_max = 12.

                noise_int = exp(A)/B*(exp(B*range_max)-exp(B*range_min))

                sig_int = htemp.Integral(htemp.FindBin(range_min), htemp.FindBin(range_max))
                sig_int = ufloat(sig_int, sqrt(sig_int))

                print sig_int, "\t", noise_int
                sn = (sig_int-noise_int)/noise_int
                print "nsig: ", nsig, "s/n = ", sn
            raw_input("Enter to continue")
        else :
            noise_int = exp(A)/B*(exp(B*fitr[1])-exp(B*fitr[0]))
            print noise_int
            
            sighist = ufloat(htemp.Integral(bmin, bmax), sqrt(htemp.Integral(bmin, bmax)))        
            print sighist
            
            sn = (sighist-noise_int)/noise_int
            print sn
            
            raw_input("Enter to continue")
    
def CompareSOE(f1, f2) :
    htocomp = []

    for f in [f1, f2] :
        if (f == "") :
            continue
        rfile = TFile(f)

        det = "Upper"
        s = "Counts/Scint*BB1*HEXpos/%s_BB1_Sigma_minus_Energy_v_SOE_tof" % det
        h2d = rfile.Get(s)
        s = "Counts/Scint*BB1*HEXpos/%s_BB1_Sigma_plus_Energy_v_SOE_tof" % det
        h2d.Add(rfile.Get(s))

        det = "Lower"
        s = "Counts/Scint*BB1*HEXpos/%s_BB1_Sigma_minus_Energy_v_SOE_tof" % det
        h2d_2 = rfile.Get(s)
        s = "Counts/Scint*BB1*HEXpos/%s_BB1_Sigma_plus_Energy_v_SOE_tof" % det
        h2d_2.Add(rfile.Get(s))

        
        htemp = h2d.ProjectionY(f1, h2d.GetXaxis().FindBin(400),
                                h2d.GetXaxis().FindBin(5100))
        htemp_2 = h2d_2.ProjectionY(f1+"L", h2d_2.GetXaxis().FindBin(400),
                                    h2d_2.GetXaxis().FindBin(5100))
        htemp.SetDirectory(0)
        htemp_2.SetDirectory(0)
        htocomp.append(htemp)
        htocomp.append(htemp_2)

    htocomp[0].Draw("HIST")
    htocomp[1].SetLineColor(kRed)
    htocomp[1].Draw("SAME")
    if (len(f) > 0) :
        htocomp[2].Draw("SAME")
        htocomp[3].SetLineColor(kRed)
        htocomp[3].Draw("SAME")
    gPad.Update()


def CheckConsistency(fn1, fn2) :
    f1 = TFile(fn1)
    f2 = TFile(fn2)

    h1 = f1.Get("Counts/Scint*BB1*eMCPcuts/Upper_BB1_Sigma_minus_scint_energy")
    h2 = f2.Get("Counts/Scint*BB1*eMCPcuts/Upper_BB1_Sigma_minus_scint_energy")
    if (not h1.GetEntries() == h2.GetEntries()) :
        print "Entries: ", h1.GetEntries(), " \t", h2.GetEntries()
        return False

    suma = 0
    for i in range(h1.GetNbinsX()+1) :
        a = h1.GetBinContent(i) - h2.GetBinContent(i)
        suma = suma + a
        
        if (abs(a) > 0) :
            print "Bin %d: %s\t%s" % i, h1.GetBinContent(i), h2.GetBinContent(i)

    h1 = f1.Get("Counts/Scint*BB1*eMCPcuts/Lower_BB1_Sigma_plus_scint_energy")
    h2 = f2.Get("Counts/Scint*BB1*eMCPcuts/Lower_BB1_Sigma_plus_scint_energy")

    if (not h1.GetEntries() == h2.GetEntries()) :
        print "Entries: ", h1.GetEntries(), " \t", h2.GetEntries()
        return False

    for i in range(h1.GetNbinsX()+1) :
        a = h1.GetBinContent(i) - h2.GetBinContent(i)
        suma = suma + a
        
        if (abs(a) > 0) :
            print "Bin %d: %s\t%s" % i, h1.GetBinContent(i), h2.GetBinContent(i)

    f1.Close()
    f2.Close()
    return suma == 0
    

def GetStatisticsForUncertainty(Aset, goal_uncert) :
    print "Goal uncert: dA/A = %5.3f%%" % goal_uncert

    r = (1.-Aset)/(1.+Aset)
    this_uncert = 100
    N = 1
    while (this_uncert > goal_uncert) :
        N = N*1.5
        f1 = ufloat(N*r, sqrt(N*r))
        f2 = ufloat(N*r, sqrt(N*r))
        f3 = ufloat(N, sqrt(N))
        f4 = ufloat(N, sqrt(N))

        s = sqrt(f1*f2/f3/f4)
        A = (1.-s)/(1.+s)
        this_uncert = 100.*abs((A.std_dev/A.nominal_value))

        print "%e\t%f" % (N, this_uncert)
    print N*r
    print "Requires N1 = N2 = %e" % (N*r)
    print "         N3 = N4 = %e" % N
    print "  total          = %e" % (N*r+N*r+N+N)

def PlotInPairs(fname, dir = "Scint*BB1*eMCPcuts", rbin = 16) :
    f = TFile(fname)
    
    #c1 = r.TCanvas()
    #c1.Divide(2, 2)
    d = ["Upper", "Lower"]
    p = ["plus", "minus"]
    htail = ["evr_scint_energy", "evr_summed_energy", "bb1_energy_v_rad"]
    allh = []
    allc = []
    for a in range(len(htail)) :
        c = TCanvas()
        allc.append(c)
        c.Divide(2, 2)
        for i in range(4) :
            c.cd(i+1)    
            hname = "Counts/%s/%s_BB1_Sigma_%s_%s" % (dir, d[i/2], p[i%2], htail[a])
            print hname
            scint_2d = f.Get(hname)

            scint_up_in  = scint_2d.ProjectionX(hname + "_in",
                                                scint_2d.GetYaxis().FindBin(0.),
                                                scint_2d.GetYaxis().FindBin(11.)-1)
            scint_up_out =  scint_2d.ProjectionX(hname + "_out",
                                                 scint_2d.GetYaxis().FindBin(11.),
                                                 scint_2d.GetYaxis().FindBin(15.5)-1)
            allh.append(scint_up_in)
            allh.append(scint_up_out)
            scint_up_in.SetLineColor(r.kRed)
            if (a < 2) :
                scint_up_in.Rebin(rbin)
                scint_up_out.Rebin(rbin)
                scint_up_in.GetXaxis().SetRangeUser(0, 5200)
                scint_up_out.GetXaxis().SetRangeUser(0, 5200)
            else :
                scint_up_in.GetXaxis().SetRangeUser(0, 350)
                scint_up_out.GetXaxis().SetRangeUser(0, 350)
                
            scint_up_in.Draw("HIST")
            scint_up_out.Draw("SAME")

            r.gPad.Update()
    
    cpos = r.TCanvas()
    cpos.Divide(2, 2)
    for i in range(4) :
        cpos.cd(i+1)    
        hname = "Counts/%s/%s_BB1_Sigma_%s_bb1_position" % (dir, d[i/2], p[i%2])
        print hname 
        h2d = f.Get(hname)
        allh.append(h2d)
        h2d.Draw("COL")

        e1 = TEllipse(0., 0., 11., 11.)
        e1.SetFillStyle(0)
        e1.SetLineWidth(3)
        e1.Draw()
        allh.append(e1)
        r.gPad.Update()
    raw_input("Enter to end")

def ProduceResultFile(gA_list, rmax = 15.5, Emin = 1000, Emax = 4500,
                      dirname = "Counts/Scint*BB1*eMCPcuts/",
                      hname = "Physics_scint_energy",
                      modifier = "",
                      setlist = ["A", "B", "C", "D"],
                      rmin = 0.,
                      projname = "evr",
                      mu_min = 0.965, mu_max = 1.000,
                      use_cut = False) :

    modifier_sim = modifier
    modifier_dat = modifier
    if (modifier == "_calU-1" or modifier == "_calU+1" or
        modifier == "_calL-1" or modifier == "_calL+1" or
        modifier == "_calU+1L-1" or modifier == "_calU-1L+1" or
        modifier == "_tof1sig" or modifier == "_tof2sig" or
        modifier == "_tof4sig" or modifier == "_tof5sig" or
        modifier == "_pos19sig" or modifier == "_tof.75sig" or
        modifier == "_tof.50sig" or  modifier == "_main_echo1"  or
        modifier == "_main_echo2" or modifier == "_hex5sig" or
        modifier == "_hex19sig") :
        modifier_sim = ""
    if (use_cut) :
        modifier_sim = ""
    flist_AB = []
    flist_CD = []
    for gA in gA_list :
        fname = "asymmetry_polBoth_M275_Be229_S295_emstandard_opt3_cfs1um_fr0.002_AB_gA%1.3f%s.root" % (gA, modifier_sim)
        if (not path.isfile(fname)) :
            print "%s not found. aborting" % fname
            return 0
        
        flist_AB.append(fname)
        fname = "asymmetry_polBoth_M275_Be229_S295_emstandard_opt3_cfs1um_fr0.002_CDE_gA%1.3f%s.root" % (gA, modifier_sim)
        if (not path.isfile(fname)) :
            print "%s not found. aborting" % fname
            return 0
        flist_CD.append(fname)

    thish = ""
    if (projname == "evr") :
        thish = "%s_rmax%3.1f" % (hname, rmax)
        if (rmin > 0.) :
            thish = "%s_rmin%3.1f" % (thish, rmin)
    elif (projname == "evmu") :
        thish = "%s_mumax%4.1f_mumin%4.1f" % (hname, mu_min, mu_max)
    else :
        thish = "Physics_%s" % projname

    print "Looking for hist %s" % thish
    outf = open("gafit.dat", "w")
    for i in range(len(gA_list)) :
        gA = gA_list[i]
        name_AB = flist_AB[i]
        name_CD = flist_CD[i]

        
        chi2 = 0
        
        print "Comparing to %s" % name_AB
        if "A" in setlist :
            chi2 = chi2 + CompareAll(name_AB, ["A"], Emin, Emax, thish,
                                     datsuffix = modifier_dat)[0]
        if "B" in setlist :
            chi2 = chi2 + CompareAll(name_AB, ["B"], Emin, Emax, thish,
                                     datsuffix = modifier_dat)[0]
        print "Comparing to %s" % name_CD
        if "C" in setlist :
            chi2 = chi2 + CompareAll(name_CD, ["C"], Emin, Emax, thish,
                                     datsuffix = modifier_dat)[0]
        if "D" in setlist :
            chi2 = chi2 + CompareAll(name_CD, ["D"], Emin, Emax, thish,
                                     datsuffix = modifier_dat)[0]

        print "%f\t%f" % (gA, chi2)
        outf.write("%f\t%f\n" % (gA, chi2))
    outf.close()

    GetResult("gafit.dat")
        

def GetResult(fname) :
    global saves
    gStyle.SetOptFit(0)

    c = TCanvas("c", "c", 1100, 550)
    c.SetBottomMargin(0.6*c.GetBottomMargin())
    c.SetLeftMargin(0.7*c.GetLeftMargin())
    c.SetTopMargin(3.*c.GetTopMargin())
    c.cd()


    gPad.Update()
    gPad.SetTickx(0)    
    chi2g = TGraph(fname)
    chi2g.GetXaxis().SetTitle("g_{A}")
    chi2g.GetYaxis().SetTitle("#chi^{2}")
    chi2g.GetYaxis().SetTitleOffset(0.7)
    #chi2g.Draw("AP")
    chi2g.GetXaxis().SetRangeUser(0.82, 0.96)
    poly = TF1("chi2_parabola", "pol2(0)", 0.82, 0.96)
    poly.SetParameters(18000, -40000, 22000)

    chi2g.Fit(poly, "WN")
    poly.GetYaxis().SetRangeUser(80, 210)
    poly.Draw("L")
    chi2g.Draw("P SAME")
    gPad.Update()
    sbresult = ufloat(-0.563, 0.009)
    sbgamin = ((sbresult.nominal_value-sbresult.std_dev)+0.3721)/-0.2146
    sbgamax = ((sbresult.nominal_value+sbresult.std_dev)+0.3721)/-0.2146
    print "SB gamin: %f" % sbgamin
    sbl1 = TLine(sbgamin, gPad.GetUymin(), sbgamin, gPad.GetUymax())
    sbl2 = TLine(sbgamax, gPad.GetUymin(), sbgamax, gPad.GetUymax())
    for s in [sbl1, sbl2] :
        s.SetLineColor(kBlue)
        s.SetLineWidth(2)
        s.SetLineStyle(2)
        s.Draw()

    ptext = TLatex(0.846, 140.44, "Previous work")
    mtext = TLatex(0.915, 140.44, "This work")
    for t in [ptext, mtext] :
        t.SetTextFont(132)
        t.SetTextAngle(90)
        t.Draw()
    #sb = AddHatchedRegion(sbgamax, sbgamin, gPad.GetUymin(), gPad.GetUymax())
    #sb.SetFillStyle(3007)
    #sb.Draw("F")

    gPad.Update()
    saves.append(chi2g)

    minx = poly.GetMinimumX()
    miny = poly.Eval(minx)

    print ""
    text = "Minimum at g_{A} = %5.3f has #chi^{2}_{fit}/123 = %3.2f" % (minx, miny/123.)
    print text
    t1 = TLatex(0.88, 145, text)
    t1.SetTextFont(132)
    t1.Draw()

    saves.append(t1)

    yplus1 = miny + 1.0
    xlo = poly.GetX(yplus1, 0,    minx)
    xup = poly.GetX(yplus1, minx, 2.)

    up_err = xup - minx
    lo_err = minx - xlo
    avg_err = 0.5*(up_err+lo_err)

    result = ufloat(minx, avg_err)
    print "Symmetric errors: ", result
    #text = "g_{A} =" + " {:L}".format(result)
    text = "g_{A} = %1.3f #pm %1.3f" % (result.nominal_value, result.std_dev)
    print text
    t2 = TLatex(0.885, 138, text)
    t2.SetTextFont(132)
    t2.Draw()
    saves.append(t2)

    print "#### NOTE: This (quick) script just calculates Abeta from a linear"
    print "           approximation for Abeta(gA)."


    Abeta = -0.3721 - (0.2146*result)
    print "Abeta(0)           = {:.2u}".format(Abeta)
    #print "Asymmetryic errors: %5.3f + %5.3f - %5.3f" % (minx, up_err, lo_err)
    correction = 1.0013
    Abeta_corr = Abeta.nominal_value*correction
    dAbeta_corr = sqrt(fabs(Abeta_corr-Abeta.nominal_value)**2+Abeta.std_dev**2)
    print "Abeta(0) corrected = %1.4f+/-%1.4f" % (Abeta_corr, dAbeta_corr)
    print "NOTE: If parabola is unbalanced, I will run more simulation to"
    print "      balance it out. This may change Abeta"

    print ""
    print "Recommended next points: %5.3f, %5.3f, %5.3f" % (minx-avg_err, minx, minx+avg_err)
    
    print result.nominal_value, "\t", result.std_dev
    gPad.Update()
    l = (-0.3721 - (0.2146*gPad.GetUxmin()))
    r = (-0.3721 - (0.2146*gPad.GetUxmax()))
    print gPad.GetUxmin()
    print gPad.GetUxmax()
    print gPad.GetUymax()
    print l
    print r
    a1 = TGaxis(gPad.GetUxmax(), gPad.GetUymax(), gPad.GetUxmin(), gPad.GetUymax(), r, l, 510, "+")
    a1.SetTitle("A_{#beta} (E_{e} = 0 keV)")
    a1.SetTitleFont(132)
    a1.SetLabelFont(132)
    a1.SetLabelSize(chi2g.GetXaxis().GetLabelSize())
    a1.SetTitleSize(chi2g.GetXaxis().GetTitleSize())
    a1.SetLabelOffset(10.*a1.GetLabelOffset())
    a1.Draw("SAME")
    gPad.Update()

    print result
    bf = AddHatchedRegion(result.nominal_value-result.std_dev,
                          result.nominal_value+result.std_dev,
                          gPad.GetUymin(), gPad.GetUymax())
    bf.SetFillColor(kRed)
    bf.SetFillStyle(3006)
    bf.Draw("F SAME")
    gPad.Update()

    l1 = TLine(result.nominal_value, gPad.GetUymin(), result.nominal_value, miny)
    l2 = TLine(result.nominal_value, miny, gPad.GetUxmin(), miny)
    l3 = TLine(result.nominal_value-result.std_dev, gPad.GetUymin(),
               result.nominal_value-result.std_dev, miny+1)
    l3.SetLineStyle(2)
    l4 = TLine(result.nominal_value+result.std_dev, gPad.GetUymin(),
               result.nominal_value+result.std_dev, miny+1)
    l4.SetLineStyle(2)
    l5 = TLine(result.nominal_value+result.std_dev, miny+1,
               gPad.GetUxmin(), miny+1)
    l5.SetLineStyle(2)
    

    for l in [l1, l2] :##, l3, l4, l5] :
        l.SetLineWidth(2)
        l.Draw()
    gPad.Update()

    ch = raw_input("Enter \"SAVE\" to save as .eps for thesis\t")
    if (ch == "SAVE") :
        c.SaveAs("~/Documents/Reports/Images/chi2_parabola.eps")
        
def CompareAll(simfilename, setlist,  xmin, xmax,
               hname = "Physics_scint_energy_rmax15.5",
               datsuffix = "") :
#    setlist = ["A", "B", "C", "D"]
#    setlist = ["C", "D"]
        
    sumc = 0
    sumn = 0
    for s in setlist :
        fname = "asymmetry_%s%s.root" % (s, datsuffix)
        print fname
        print "********"
        print fname
        print simfilename
        print hname
        print "********"
        comp = HCompare(dfilename = fname,
                        sfilename = simfilename,
                        hname = hname)
        result = comp.Chi2InRange(xmin, xmax)
        sumc = sumc + result[0]
        sumn = sumn + result[1]
        #sumc = sumc + comp.chi2
        #sumn = sumn + comp.ndf
        print "%s\t%f\t%f" % (s, result[0], result[1])
    return [sumc, sumn]
    
def AnalyzeAngle(fname, dirname = "Counts/Scint*BB1*eMCPcuts/",
                 suffix = "scint_energy", nbins = 10., rmax = 19.) :    
    har = []
    det = ["Upper", "Lower"]
    sig = ["plus" , "minus"]
    i = 0
    f = TFile(fname)
    for i in range(4) :
        hname = "%s%s_BB1_Sigma_%s_evmu_%s" % (dirname, det[i/2], sig[i%2], suffix)
        print hname
        h2d = f.Get(hname)
        har.append(h2d)
        
    hup_plus = har[0]
    hup_mins = har[1]
    hlo_plus = har[2]
    hlo_mins = har[3]

    A = pi*rmax*rmax/(nbins)

    r = [0]
    rbar = array('d', [])
    dr = array('d', [])
    asym = array('d', [])
    dasym = array('d', [])
    i = 0
    dn = (1.0-rmax)/nbins
    n = rmax
    while(n <= 1.) :
#         n = sqrt((A/pi) + r[i]**2)
#         if (n > rmax) :
#             break
#         n = round((2*n))/2
        print n
        r.append(n)
        #thisrbar = (n + r[i])/2.

        #costheta = 103.5 / sqrt(thisrbar**2 + 103.5**2)
        costheta = n
        thisrbar = n
        rbar.append(thisrbar)
        print "%f\t%f" % (n, costheta)
        firstxbin = hup_plus.GetXaxis().FindBin(2000)
        lastxbin  = hup_plus.GetXaxis().FindBin(4500)
        firstybin = hup_plus.GetYaxis().FindBin(r[i])
        lastybin  = hup_plus.GetYaxis().FindBin(r[i+1])-1
        x0 = hup_mins.Integral(firstxbin, lastxbin, firstybin, lastybin)
        x2 = hup_plus.Integral(firstxbin, lastxbin, firstybin, lastybin)

        firstxbin = hlo_plus.GetXaxis().FindBin(2000)
        lastxbin  = hlo_plus.GetXaxis().FindBin(4500)
        firstybin = hlo_plus.GetYaxis().FindBin(r[i])
        lastybin  = hlo_plus.GetYaxis().FindBin(r[i+1])-1
        x1 = hlo_plus.Integral(firstxbin, lastxbin, firstybin, lastybin)
        x3 = hlo_mins.Integral(firstxbin, lastxbin, firstybin, lastybin)
        #print "%f\t%f\t%f\t%f\t" % (x0, x1, x2, x3)
        n0 = ufloat(x0, sqrt(x0))
        n1 = ufloat(x1, sqrt(x1))
        n2 = ufloat(x2, sqrt(x2))
        n3 = ufloat(x3, sqrt(x3))

        s = sqrt(n0*n1/n2/n3)
        thisasym = (1.-s)/(1.+s)/costheta

        snom = s.nominal_value
        man_uncert = (snom/(1.+snom)**2)*sqrt((1./x0)+(1./x1)+(1./x2)+(1./x3))/costheta

        asym.append(thisasym.nominal_value)
        dasym.append(thisasym.std_dev)
        dr.append(n-thisrbar)
        #print "%f\t%f\t%f\t%f\t%f\t%f\t%f" % (firstxbin, lastxbin, firstybin, lastybin, thisasym.nominal_value, thisasym.std_dev, man_uncert)

        i = i + 1
        n = n + dn

    zero = array('d', [0]*len(rbar))
    g = TGraphErrors(len(rbar), rbar, asym, dr, dasym)
    g.GetYaxis().SetTitle("Asymmetry / cos(#theta)")
    g.GetXaxis().SetTitle("BB1 radius [mm]")
    return g

def CompareCounts(rfile, rfile2, basepath="/", showonlybad = False) :
    for key in rfile.GetListOfKeys() :
        kname = key.GetName()
        if key.IsFolder() :
            CompareCounts(rfile.Get(kname), rfile2, basepath+kname+"/",
                          showonlybad)
        else :
            h = 0
            h2 = 0
            if (basepath == "/") :
                h = rfile.Get(kname)
                h2 = rfile2.Get(kname)
            else :
                h = rfile.Get(basepath + kname)
                h2 = rfile2.Get(basepath + kname)
            if (not h or not h2) :
                print "Problem with %s" % (basepath+kname)
            elif (not h.InheritsFrom("TH1") or not h2.InheritsFrom("TH1")) :
                print "Not inheriting from TH1: %s" % (basepath+kname)
            else :
                r = float(h.GetEntries())/float(h2.GetEntries())
                if (showonlybad and r < 10. or not showonlybad) :
                    print "%-70s\t%10d\t%10d\t%2.1f" % (basepath + h.GetName(),
                                                        h.GetEntries(),
                                                        h2.GetEntries(),
                                                        r)
                

def Abeta_SM(x, par) :
    vc = ElectronVoverC(x[0])
    Abeta = par[0]
    
    return (vc*Abeta)

def GetMOTpolarization(fname) :
    f = TFile(fname)
    h = f.Get("SuperRatioT")
    h.SetDirectory(0)


    tmin = 0    
    tmax = int(raw_input("Enter time when OP starts: [330]") or "330")

    h.Draw("E1X0")

    l = TF1("line", "[0]", tmin, tmax)
    h.Fit(l, "R")
    A = ufloat(l.GetParameter(0), l.GetParError(0))
    P = A / Asm
    print "A_MOT = ", A
    print "P_MOT = ", P
    return P

def UnevenRebin2DX(horig, htofill) :

    for ybin in range(0, horig.GetNbinsY() + 1) :
        y = horig.GetYaxis().GetBinCenter(ybin)
        for i in range(0, horig.GetNbinsX()+1) :
            xlo = horig.GetXaxis().GetBinLowEdge(i)
            xup = xlo + horig.GetXaxis().GetBinWidth(i)
            
            cont = horig.GetBinContent(i, ybin)
            
            if (cont == 0) :
                continue

            if (htofill.GetXaxis().FindBin(xlo) == htofill.GetXaxis().FindBin(xup)) :
                b = htofill.FindBin(xlo, y)
                htofill.SetBinContent(b, htofill.GetBinContent(b)+cont)
            else :
                b1 = htofill.FindBin(xlo, y)
                b2 = htofill.FindBin(xup, y)
                edge = htofill.GetXaxis().GetBinLowEdge(htofill.GetXaxis().FindBin(xup))
                frac1 = (edge - xlo) / (xup - xlo)
                frac2 = 1.0 - frac1

                htofill.SetBinContent(b1, htofill.GetBinContent(b1) + cont*frac1)
                htofill.SetBinContent(b2, htofill.GetBinContent(b2) + cont*frac2)

def UnevenRebin(horig, htofill) :
    debug = 0

    if (debug == 1) :
        print "htofill: nbins = %d, low = %f, high = %f" % (htofill.GetNbinsX(),
                                                            htofill.GetBinLowEdge(1),
                                                            htofill.GetBinLowEdge(htofill.GetNbinsX())+htofill.GetBinWidth(htofill.GetNbinsX()))

    bmax = horig.GetNbinsX()
    #    for i in range(0, horig.GetNbinsX()+1) :
    for i in range(0, bmax+1) :
        xlo = horig.GetBinLowEdge(i)
        xup = xlo + horig.GetBinWidth(i)

        cont = horig.GetBinContent(i)
        if (cont == 0) :
            continue

        if (debug > 0 and cont > 0) :
            print "Source bin %d: lo edge = %f\thi edge = %f\tcontent = %f" % (i, xlo, xup, cont)

        if (htofill.FindBin(xlo) == htofill.FindBin(xup)) : # all in one bin            
            b = htofill.FindBin(xlo)
            htofill.SetBinContent(b, htofill.GetBinContent(b) + cont)
            if (debug > 0) :
                print "Dest bin %d incremented by %f" % (b, cont)
        else :                  # split it up (linearly)
            b1 = htofill.FindBin(xlo)
            b2 = htofill.FindBin(xup)
            if (not (b2 == b1 + 1)) :
                print "Uhoh, unsupported UnevenRebin!!"
                return 0
            edge = htofill.GetBinLowEdge(b2)
            frac1 = (edge - xlo) / (xup - xlo)
            frac2 = 1.0 - frac1

            htofill.SetBinContent(b1, htofill.GetBinContent(b1) + cont*frac1)
            htofill.SetBinContent(b2, htofill.GetBinContent(b2) + cont*frac2)
            if (debug > 0) :
                print "Dest bin %d incremented by %f" % (b1, cont*frac1)
                print "Dest bin %d incremented by %f" % (b2, cont*frac2)
                print "Total increment: %f" % (cont*(frac1+frac2))
                raw_input("e")
    if (abs(horig.Integral(0, bmax) - htofill.Integral(0, bmax)) > 0.1) :

        print "Entries not conserved!"
        print "Source:      ", horig.Integral(0, bmax)
        print "Destination: ", htofill.Integral(0, bmax)
    else :
        print "Entries conserved"
        print "Integral = %f" % htofill.Integral(0, bmax)
        #for i in range(100) :
        #print "%f\t%f" % (htofill.GetBinError(i), sqrt(htofill.GetBinContent(i)))
def GetSuperRatio(h0, h1, h2, h3, nbins, xlo, xup) :
    
    h0_mod = TH1D("h0_mod", "h0_mod", nbins, xlo, xup)
    h1_mod = TH1D("h1_mod", "h1_mod", nbins, xlo, xup)
    h2_mod = TH1D("h2_mod", "h2_mod", nbins, xlo, xup)
    h3_mod = TH1D("h3_mod", "h3_mod", nbins, xlo, xup)
    
    UnevenRebin(htofill = h0_mod, horig = h0)
    UnevenRebin(htofill = h1_mod, horig = h1)
    UnevenRebin(htofill = h2_mod, horig = h2)
    UnevenRebin(htofill = h3_mod, horig = h3)
    
    sr = TH1D("super_ratio", "super_ratio", nbins, xlo, xup)
    for i in range(nbins + 1) :
        none = ufloat(h0_mod.GetBinContent(i),
                      sqrt(h0_mod.GetBinContent(i)))
        ntwo = ufloat(h1_mod.GetBinContent(i),
                      sqrt(h1_mod.GetBinContent(i)))
        nthr = ufloat(h2_mod.GetBinContent(i),
                      sqrt(h2_mod.GetBinContent(i)))
        nfou = ufloat(h3_mod.GetBinContent(i),
                      sqrt(h3_mod.GetBinContent(i)))
        if (nthr == 0 or nfou == 0) :
            sr.SetBinContent(i, 0.)
            sr.SetBinError(i, 0.) 
        else :
            s = sqrt((none*ntwo)/(nthr*nfou))
            A = (1. - s) / (1. + s)
            sr.SetBinContent(i, A.nominal_value)
            sr.SetBinError(  i, A.std_dev)
        # print "%f\t%f\t%f\t%f\t%f\t%f\t%f" % (h0_mod.GetBinCenter(i),
        #                                       h0_mod.GetBinContent(i),
        #                                       h1_mod.GetBinContent(i),
        #                                       h2_mod.GetBinContent(i),
        #                                       h3_mod.GetBinContent(i),
        #                                       sr.GetBinContent(i),
        #                                       sr.GetBinError(i))
    return sr

def PlotPhysicsSuperRatioTOF(fname, tof_min, tof_max) :
    f = TFile(fname)

    har = []
    det = ["Upper", "Lower"]
    sig = ["plus" , "minus"]
    i = 0
    for i in range(4) :

        hname = "Counts/Scint*BB1*eMCP/%s_BB1_Sigma_%s_Energy_v_SOE_tof" % (det[i/2], sig[i%2])
        
        print hname
        h2d = f.Get(hname)
        bmin = 0
        bmax = 0
        bmin = h2d.GetYaxis().FindBin(tof_min[i/2])
        bmax = h2d.GetYaxis().FindBin(tof_max[i/2]) - 1
        har.append(h2d.ProjectionX(hname+"_px", bmin, bmax))
        
    hup_plus = har[0]
    hup_mins = har[1]
    hlo_plus = har[2]
    hlo_mins = har[3]
    print "Entries: %d\t%d\t%d\t%d" % (har[0].GetEntries(),
                                       har[1].GetEntries(),
                                       har[2].GetEntries(),
                                       har[3].GetEntries())
    print "Computing physics super ratio"
    cphysics = PlotAsymE([hup_mins, hlo_plus], [hup_plus, hlo_mins],
                         xax = "scint_energy")



def PlotAllSuperRatiosV2(fname, dir = "Counts/Scint*BB1*eMCPcuts/",
                         suffix = "scint_energy", bb1_rmax = 15.5,
                         bb1_rmin = 0., projname = "evr", mu_min = 0.965,
                         mu_max = 1.000) :
#
# Generate physics, luminosity, and detector asymmetry plots
# file specified by fname is MODIFIED and existing histograms are overwritten
# The bb1_rmax is a strict less than not less than or equal to
# The bb1_rmin is a less that OR equal to 
#
#Counts/Scint*BB1*eMCPcuts/Upper_BB1_Sigma_plus_evr_scint
#Counts/Scint*BB1*eMCPcuts/Upper_BB1_Sigma_minus_evr_scint
#Counts/Scint*BB1*eMCPcuts/Lower_BB1_Sigma_plus_evr_scint
#Counts/Scint*BB1*eMCPcuts/Lower_BB1_Sigma_minus_evr_scint
    if (projname == "evr") :
        print "NOTICE: Using radial bins, ignoring mu_min and mu_max parameters"
    if (projname == "evmu") :
        print "NOTICE: Using cos(theta) bins, ignoring bb1_rmin and bb1_rmax",
        print " parameters"

    if (bb1_rmin != int(bb1_rmin) and bb1_rmin - 0.5 != int(bb1_rmin)) :
        print "WARNING: bb1_rmin bust be integer multipe of 0.5"
        print "   this is to required to land on a bin boundary  in r"
        bb1_rmin = floor(bb1_rmin*2)/2.
        print "   SET bb1_rmin to %f" % bb1_rmin
        raw_input("Enter to continue")

    if (bb1_rmax != int(bb1_rmax) and bb1_rmax - 0.5 != int(bb1_rmax)) :
        print "WARNING: bb1_rmax bust be integer multipe of 0.5"
        print "   this is to required to land on a bin boundary  in r"
        bb1_rmax = ceil(bb1_rmax*2)/2.
        print "   SET bb1_rmax to %f" % bb1_rmax
        raw_input("Enter to continue")
    f = TFile(fname, "UPDATE")

    har = []
    det = ["Upper", "Lower"]
    sig = ["plus" , "minus"]
    i = 0
    for i in range(4) :
        hname = "%s%s_BB1_Sigma_%s_%s_%s" % (dir, det[i/2], sig[i%2], projname,
                                             suffix)
        print hname
        h2d = f.Get(hname)
        bmin = 0
        bmax = 0
        if (projname == "evr") :
            bmin = h2d.GetYaxis().FindBin(bb1_rmin)
            bmax = h2d.GetYaxis().FindBin(bb1_rmax) - 1
        elif (projname == "evmu") :
            bmin = h2d.GetYaxis().FindBin(mu_min)
            bmax = h2d.GetYaxis().FindBin(mu_max) - 1
        har.append(h2d.ProjectionX("_px", bmin, bmax))
        
    hup_plus = har[0]
    hup_mins = har[1]
    hlo_plus = har[2]
    hlo_mins = har[3]
    print "Entries: %d\t%d\t%d\t%d" % (har[0].GetEntries(),
                                       har[1].GetEntries(),
                                       har[2].GetEntries(),
                                       har[3].GetEntries())
    # hup_plus = f.Get(dir + "Upper_BB1_Sigma_plus_"  + suffix)
    # hup_mins = f.Get(dir + "Upper_BB1_Sigma_minus_" + suffix)
    
    # hlo_plus = f.Get(dir + "Lower_BB1_Sigma_plus_" + suffix)
    # hlo_mins = f.Get(dir + "Lower_BB1_Sigma_minus_" + suffix)
    
    print "Computing physics super ratio"
    cphysics = PlotAsymE([hup_mins, hlo_plus], [hup_plus, hlo_mins],
                         xax = suffix)
    cphysics.SetName("Physics_" + suffix)
    phys_hist = cphysics.GetPrimitive("pup").GetPrimitive("Physics")
    hname = ""
    if (projname == "evr") :
        hname = "%s_%s_rmax%3.1f" % (phys_hist.GetName(), suffix, bb1_rmax)
        if (bb1_rmin > 0) :
            hname = "%s_rmin%3.1f" % (hname, bb1_rmin)
    elif (projname == "evmu") :
        hname = "%s_%s_mumax%1.4f" % (phys_hist.GetName(), suffix, mu_max)
        print hname
        if (mu_min > 0.965) :
            hname = "%s_mumin%1.4f" % (hname, mu_min)
            print hname

    #phys_hist.SetName(phys_hist.GetName() + "_" + suffix + "_rmax" + str(bb1_rmax))
    phys_hist.SetName(hname)

    print "Computing luminositry super ratio"
    cluminosity = PlotAsymE([hup_mins, hlo_mins], [hup_plus, hlo_plus],
                            prediction = 0., tfix = "Luminosity", xax = suffix)
    cluminosity.SetName("Luminosity_" + suffix)
    lumi_hist = cluminosity.GetPrimitive("pup").GetPrimitive("Luminosity")
    if (projname == "evr") :
        hname = "%s_%s_rmax%3.1f" % (lumi_hist.GetName(), suffix, bb1_rmax)
        if (bb1_rmin > 0) :
            hname = "%s_rmin%3.1f" % (hname, bb1_rmin)
    elif (projname == "evmu") :
        hname = "%s_%s_mumax%1.4f" % (lumi_hist.GetName(), suffix, mu_max)
        print hname
        if (mu_min > 0.965) :
            hname = "%s_mumin%1.4f" % (hname, mu_min)
            print hname

    #lumi_hist.SetName(lumi_hist.GetName() + "_" + suffix + "_rmax" + str(bb1_rmax))
    lumi_hist.SetName(hname)
    
    print "Computing physics super ratio"
    cdetector = PlotAsymE([hup_mins, hup_plus], [hlo_mins, hlo_plus],
                          prediction = 0., tfix = "Detector", xax = suffix)
    cdetector.SetName("Detector_" + suffix)
    dete_hist = cdetector.GetPrimitive("pup").GetPrimitive("Detector")
    if (projname == "evr") :
        hname = "%s_%s_rmax%3.1f" % (dete_hist.GetName(), suffix, bb1_rmax)
        if (bb1_rmin > 0) :
            hname = "%s_rmin%3.1f" % (hname, bb1_rmin)
    elif (projname == "evmu") :
        hname = "%s_%s_mumax%1.4f" % (dete_hist.GetName(), suffix, mu_max)
        print hname
        if (mu_min > 0.965) :
            hname = "%s_mumin%1.4f" % (hname, mu_min)
            print hname

    #dete_hist.SetName(dete_hist.GetName() + "_" + suffix + "_rmax" + str(bb1_rmax))
    dete_hist.SetName(hname)
    

    f.cd()
    phys_hist.Write(phys_hist.GetName(), kWriteDelete)
    lumi_hist.Write(lumi_hist.GetName(), kWriteDelete)
    dete_hist.Write(dete_hist.GetName(), kWriteDelete)

    if (projname == "evr") :
        print "NOTICE: Using radial bins, ignoring mu_min and mu_max parameters"
    if (projname == "evmu") :
        print "NOTICE: Using cos(theta) bins, ignoring bb1_rmin and bb1_rmax",
        print " parameters"

    # outf = open("v2.dat", "w")
    # for i in range(hlo_plus.GetNbinsX()+1) :
    #     outf.write("%s\t%f\t%f\n" % (hlo_plus.GetBinCenter(i),
    #                                  hlo_plus.GetBinContent(i),
    #                                  hlo_plus.GetBinError(i)))
    # outf.close()
    f.Close()
    

def PlotAllSuperRatios(fname, dir = "Counts/Scint*BB1*eMCPcuts/",
                       suffix = "scint_energy") :
    print "USE VERSION TWO INSTEAD"
    return 0
    f = TFile(fname, "UPDATE")

    hup_plus = f.Get(dir + "Upper_BB1_Sigma_plus_"  + suffix)
    hup_mins = f.Get(dir + "Upper_BB1_Sigma_minus_" + suffix)
    
    hlo_plus = f.Get(dir + "Lower_BB1_Sigma_plus_" + suffix)
    hlo_mins = f.Get(dir + "Lower_BB1_Sigma_minus_" + suffix)
    
    print "Computing physics super ratio"
    cphysics = PlotAsymE([hup_mins, hlo_plus], [hup_plus, hlo_mins],
                         xax = suffix)
    cphysics.SetName("Physics_" + suffix)
    phys_hist = cphysics.GetPrimitive("pup").GetPrimitive("Physics")
    phys_hist.SetName(phys_hist.GetName() + "_" + suffix)

    print "Computing luminositry super ratio"
    cluminosity = PlotAsymE([hup_mins, hlo_mins], [hup_plus, hlo_plus],
                            prediction = 0., tfix = "Luminosity", xax = suffix)
    cluminosity.SetName("Luminosity_" + suffix)
    lumi_hist = cluminosity.GetPrimitive("pup").GetPrimitive("Luminosity")
    lumi_hist.SetName(lumi_hist.GetName() + "_" + suffix)
    
    print "Computing physics super ratio"
    cdetector = PlotAsymE([hup_mins, hup_plus], [hlo_mins, hlo_plus],
                          prediction = 0., tfix = "Detector", xax = suffix)
    cdetector.SetName("Detector_" + suffix)
    dete_hist = cdetector.GetPrimitive("pup").GetPrimitive("Detector")
    dete_hist.SetName(dete_hist.GetName() + "_" + suffix)

    phys_hist.Write(phys_hist.GetName(), kWriteDelete)
    lumi_hist.Write(lumi_hist.GetName(), kWriteDelete)
    dete_hist.Write(dete_hist.GetName(), kWriteDelete)

    # outf = open("v1.dat", "w")
    # for i in range(hlo_plus.GetNbinsX()+1) :
    #     outf.write("%s\t%f\t%f\n" % (hlo_plus.GetBinCenter(i),
    #                                  hlo_plus.GetBinContent(i),
    #                                  hlo_plus.GetBinError(i)))
    # outf.close()

    f.Close()
    
    
def PlotAsymE(hg1, hg2, prediction = Asm, tfix = "Physics",
              xax = "scint_energy") :
    global saves

    if (len(hg1) != 2 or len(hg2) != 2) :
        print "Must have correct length"
        return

    
    c = TCanvas("asym", "asym", 850, 500)
    gStyle.SetOptFit(0)
    c.SetLeftMargin(1.2*c.GetLeftMargin())
    
    p1 = TPad("pup", "pup", 0., 0.25, 1., 1.)
#    saves['p1'] = p1
    saves.append(p1)
    p2 = TPad("plo", "plo", 0.,   0., 1., 0.25)
#    saves['ps'] = p2
    saves.append(p2)
    p1.SetBottomMargin(1E-5)
    p2.SetTopMargin(1E-5)
    p2.SetBottomMargin(2*p2.GetBottomMargin())
    p1.Draw()
    p2.Draw()
    p1.cd()
    
    xlo = 0
    xup = 6000
    xfitmax = 5200
    xfitmin = 350
    if (xax == "time") :
        xup = 5000
        xfitmax = 1900
        xfitmin = 432
        
    A = GetSuperRatio(hg1[0], hg1[1], hg2[0], hg2[1], 60, xlo, xup)
    A.SetName(tfix)
#    A = GetSuperRatio(hlo_plus, hup_mins, hup_plus, hlo_mins, 60, 0, 6000)    

    saves.append(A)
    A.SetDirectory(0)

    A.SetTitle(tfix + " Super Ratio")
    A.GetXaxis().SetTitle("Scintillator Energy [keV]")
    if (xax == "time") :
        A.GetXaxis().SetTitle("Time since MOT off [#mus]")
    A.GetXaxis().SetRangeUser(xlo, xup)
    if (tfix == "Physics") :
        A.GetYaxis().SetRangeUser(-0.63, -0.42)
    else :
        A.GetYaxis().SetRangeUser(-0.12, 0.12)
        
    Atitle = tfix + " Super Ratio"
    f0 = 0
    if (tfix == "Physics" and xax == "scint_energy") :
        Atitle = Atitle# + " #pm 0.004"
        f0 = TF1("SM", Abeta_SM, xfitmin, xfitmax, 1)
    else :
        f0 = TF1("Fit", "pol0(0)", xfitmin, xfitmax)
        f0.SetTitle("Fit")
        f0.SetParameter(0, 0.)
        print xfitmax
                
    #f0.SetRange(xfitmin, xfitmax)
    A.GetYaxis().SetTitle(Atitle)
    FormatAxis(A, xtsize = 0, xlsize = 0)
    A.Draw("E1X0")
    gPad.Update()

    saves.append(f0)
    pred = prediction*Pol
    f0.SetParameter(0, pred)
    A.Fit(f0, "R")
    fit_result = ufloat(f0.GetParameter(0), f0.GetParError(0))
    fit_text = "A = {:}".format(fit_result)
    print fit_text
    f0.Draw("L SAME")
    leg = p1.BuildLegend(0.5, 0.64, 0.80, 0.88)
    leg.SetShadowColor(0)
    leg.SetTextFont(132)
    leg.SetName(tfix + " Legend")
    leg.Draw()
    saves.append(leg)
    p1.Update()
    p2.cd()
    p2.SetGridx(1)
    p2.SetGridy(1)
    res_arr = GetResiduals(f0, A, xfitmin, xfitmax)
    res = res_arr[0]
    saves.append(res)
    res.Draw("P")
    res.SetMarkerStyle(kFullCircle)
    res.SetMarkerSize(0.5)
    res.GetXaxis().SetRangeUser(p1.GetUxmin(), p1.GetUxmax())
    res.GetYaxis().SetRangeUser(-5, 5)
    res.GetYaxis().SetTitle("Residuals")
    res.GetXaxis().SetTitle("Scintillator energy [keV]")
    print xax
    if (xax == "time") :
        res.GetXaxis().SetTitle("Time since MOT off [#mus]")
    res.GetYaxis().SetNdivisions(405)
    FormatAxis(res, ytsize = 0.2, ytoffset = 0.25, ytfont = 132,
               ylsize = 0.17, ylfont = 132, xtsize = 0.2, xtoffset = 0.9,
               xtfont = 132, xlsize = 0.15)

    p2.Update()

    p1.cd()
    chi2txt = "#chi^{2}/%d = %1.1f" % (res_arr[2], res_arr[1])
    tl = 0
    tr = 0
    if (tfix == "Physics") :
        tl = TLatex(980, -0.43, chi2txt)
        tr = TLatex(980, -0.47, fit_text)
    else :
        tl = TLatex(980, 0.06, chi2txt)
        tr = TLatex(980, -0.06, fit_text)

    saves.append(tl)
    saves.append(tr)
    for t in [tl, tr] :
        t.SetTextFont(132)
        t.SetTextSize(0.08)
        t.Draw()
    p1.Update()
    
    saves.append(c)
    return c

def CheckSuperRatios() :
    gStyle.SetOptFit(0)

    bf = TFile("asymmetry_B.root")
    df = TFile("asymmetry_D.root")

    clum = TCanvas("lum", "lum", 550*2, 556)
    clum.Divide(2, 1)
    clum.cd(1)
    gPad.SetBottomMargin(0.6*gPad.GetBottomMargin())
    blum = bf.Get("Luminosity_scint_energy_rmax15.5")
    dlum = df.Get("Luminosity_scint_energy_rmax15.5")


    FormatAxis(blum, ytsize = 0.05, ytoffset = 1.2, ytfont = 132,
               ylsize = 0.04, ylfont = 132, xtsize = 0.05, xtoffset = 0.9,
               xtfont = 132, xlsize = 0.04)
    FormatAxis(dlum, ytsize = 0.05, ytoffset = 1.2, ytfont = 132,
               ylsize = 0.04, ylfont = 132, xtsize = 0.05, xtoffset = 0.9,
               xtfont = 132, xlsize = 0.04)
    blum.Draw("P")
    blum.GetXaxis().SetRangeUser(0., 5200)
    blum.GetYaxis().SetRangeUser(-0.12, 0.12)
    
    
    bfit = blum.GetFunction("Fit")
    bres = GetUfloatFromFit(bfit, 0)
    x = "A_{lum} =" + '{:L}'.format(bres)
    blat = TLatex(250, 0.07, x)
    blat.Draw()
    gPad.Update()
    clum.cd(2)
    gPad.SetBottomMargin(0.6*gPad.GetBottomMargin())
    dlum.Draw("P")
    dlum.GetXaxis().SetRangeUser(0., 5200)
    dlum.GetYaxis().SetRangeUser(-0.12, 0.12)
    dfit = dlum.GetFunction("Fit")
    dres = GetUfloatFromFit(dfit, 0)
    y = "A_{lum} =" + '{:L}'.format(dres)
    dlat = TLatex(250, 0.07, y)
    dlat.Draw()
    gPad.Update()
    ch = raw_input("Enter \"SAVE\" to save for .eps for theis\t")
    if (ch == "SAVE") :
        clum.SaveAs("~/Documents/Reports/Images/luminosity_sr.eps")
    
    cdet = TCanvas("det", "det", 550*2, 556)
    cdet.Divide(2, 1)
    cdet.cd(1)
    gPad.SetBottomMargin(0.6*gPad.GetBottomMargin())
    bdet = bf.Get("Detector_scint_energy_rmax15.5")
    ddet = df.Get("Detector_scint_energy_rmax15.5")


    FormatAxis(bdet, ytsize = 0.05, ytoffset = 1.2, ytfont = 132,
               ylsize = 0.04, ylfont = 132, xtsize = 0.05, xtoffset = 0.9,
               xtfont = 132, xlsize = 0.04)
    FormatAxis(ddet, ytsize = 0.05, ytoffset = 1.2, ytfont = 132,
               ylsize = 0.04, ylfont = 132, xtsize = 0.05, xtoffset = 0.9,
               xtfont = 132, xlsize = 0.04)
    bdet.Draw("P")
    bdet.GetXaxis().SetRangeUser(0., 5200)
    bdet.GetYaxis().SetRangeUser(-0.12, 0.12)
    
    
    bfit = bdet.GetFunction("Fit")
    bres = GetUfloatFromFit(bfit, 0)
    x = "A_{det} =" + '{:L}'.format(bres)
    blat = TLatex(250, 0.07, x)
    blat.Draw()
    gPad.Update()
    cdet.cd(2)
    gPad.SetBottomMargin(0.6*gPad.GetBottomMargin())
    ddet.Draw("P")
    ddet.GetXaxis().SetRangeUser(0., 5200)
    ddet.GetYaxis().SetRangeUser(-0.12, 0.12)
    dfit = ddet.GetFunction("Fit")
    dres = GetUfloatFromFit(dfit, 0)
    y = "A_{det} =" + '{:.2uL}'.format(dres)
    dlat = TLatex(250, 0.07, y)
    dlat.Draw()
    gPad.Update()
    ch = raw_input("Enter \"SAVE\" to save for .eps for theis\t")
    if (ch == "SAVE") :
        cdet.SaveAs("~/Documents/Reports/Images/detector_sr.eps")
    
    
    
def OverlaysForThesis(fname) :

    opt = ""
    hf = TFile(fname)

    c = TCanvas("ocan", "ocan", 550*2, 556)

    c.Divide(2, 1)
    c.cd(1)
    gPad.SetBottomMargin(0.6*gPad.GetBottomMargin())
    gPad.SetLeftMargin(1.5*gPad.GetLeftMargin())
    hup_plus = hf.Get("Counts/Scint*BB1*eMCP/Upper_BB1_Sigma_plus_Energy_v_SOE_tof")
    hup_mins = hf.Get("Counts/Scint*BB1*eMCP/Upper_BB1_Sigma_minus_Energy_v_SOE_tof")
    hlo_plus = hf.Get("Counts/Scint*BB1*eMCP/Lower_BB1_Sigma_plus_Energy_v_SOE_tof")
    hlo_mins = hf.Get("Counts/Scint*BB1*eMCP/Lower_BB1_Sigma_minus_Energy_v_SOE_tof")

    pup_plus = hup_plus.ProjectionX("pup_plus", 0, -1, opt)
    pup_mins = hup_mins.ProjectionX("pup_mins", 0, -1, opt)
    plo_plus = hlo_plus.ProjectionX("plo_plus", 0, -1, opt)
    plo_mins = hlo_mins.ProjectionX("plo_mins", 0, -1, opt)

    pup_mins.SetLineColor(kRed)
    pup_mins.SetMarkerColor(kRed)
    pup_mins.SetLineStyle(2)
    plo_mins.SetLineColor(kRed)
    plo_mins.SetMarkerColor(kRed)
    plo_mins.SetLineStyle(2)
    pup_plus.SetTitle("#sigma^{+}")
    pup_mins.SetTitle("#sigma^{-}")
    

    pup_mins.Draw("HIST")
    pup_plus.Draw("SAME")
    bw = pup_mins.GetBinWidth(10)
    yax = "Counts / %2.0f keV" % bw
    pup_mins.GetXaxis().SetRangeUser(0, 5500)
    pup_mins.GetYaxis().SetTitle(yax)
    pup_mins.GetYaxis().SetTitleOffset(1.5)
    

    gPad.Update()
    ymax = gPad.GetUymax()

    leg = gPad.BuildLegend(0.63, 0.77, 0.94, 0.90)
    leg.SetShadowColor(0)
    leg.Draw()

    tup = TLatex(380, 1160, "Upper")
    tup.SetTextFont(132)
    tup.Draw()
    t1 = TLine(400, 0, 400, 1000)
    a1 = TArrow(400, 1000, 1150, 1000)
    a1.SetArrowSize(0.02)
    t3 = TLine(4500, 0, 4500, 1000)
    a3 = TArrow(4500, 1000, 3750, 1000)
    a3.SetArrowSize(0.02)

    for t in [a1, t1, a3, t3] :
        t.SetLineWidth(2)
        t.SetLineStyle(2)
        t.Draw()

    c.cd(2)
    gPad.SetBottomMargin(0.6*gPad.GetBottomMargin())
    gPad.SetLeftMargin(1.5*gPad.GetLeftMargin())
    plo_plus.Draw("HIST")
    plo_mins.Draw("SAME")
    bw = plo_plus.GetBinWidth(10)
    yax = "Counts / %2.0f keV" % bw
    plo_plus.GetXaxis().SetRangeUser(0, 5500)
    plo_plus.GetYaxis().SetRangeUser(0, ymax)
    plo_plus.GetYaxis().SetTitle(yax)
    plo_plus.GetYaxis().SetTitleOffset(1.5)
    tlo = TLatex(380, 1160, "Lower")
    tlo.SetTextFont(132)
    tlo.Draw()

    t2 = TLine(4500, 0, 4500, 1000)
    a2 = TArrow(4500, 1000, 3750, 1000)
    a2.SetArrowSize(0.02)
    t4 = TLine(400, 0, 400, 1000)
    a4 = TArrow(400, 1000, 1150, 1000)
    a4.SetArrowSize(0.02)

    for t in [a2, t2, a4, t4] :
        t.SetLineWidth(2)
        t.SetLineStyle(2)
        t.Draw()

    gPad.Update()
    ch = raw_input("Enter \"SAVE\" to save as .eps for thesis\t")
    if (ch == "SAVE") :
        c.SaveAs("~/Documents/Reports/Images/scintillator_overlay.eps")
                                    
