makeCovariance.py

Go to the documentation of this file.

import ROOT,array,math,sys

#beam_config = "reference"
beam_config = "opt"


ROOT.gStyle.SetPadRightMargin(0.15);   
ROOT.gStyle.SetPadLeftMargin(0.15);   


plots = ["numu_neutrino_ND","numubar_neutrino_ND","nue_neutrino_ND","nuebar_neutrino_ND","numu_antineutrino_ND","numubar_antineutrino_ND","nue_antineutrino_ND","nuebar_antineutrino_ND","numu_neutrino_FD","numubar_neutrino_FD","nue_neutrino_FD","nuebar_neutrino_FD","numu_antineutrino_FD","numubar_antineutrino_FD","nue_antineutrino_FD","nuebar_antineutrino_FD"]

start_indices = [0,19,38,45,52,71,90,97,104,123,142,149,156,175,194,201]

# Helper function that makes a error histogram out of a bunch of universes
# by calculating rms of universes in each bin
def MakeErrorHistogram(universes,asFrac=True):
    n_universes = len(universes)
    avg_hist = universes[0].Clone(universes[0].GetName()+"_avg")
    error_hist = universes[0].Clone(universes[0].GetName()+"_error")
    avg_hist.Reset()
    error_hist.Reset()
    avg_hist.Reset()
    error_hist.Reset()
    for i in range(n_universes):
        avg_hist.Add(universes[i])
    avg_hist.Scale(1/float(n_universes))

    # Calculate Error = RMS of universes
    for j in range(error_hist.GetNbinsX()+2):
        error = 0
        for i in range(n_universes):
            error += 1/float(n_universes)*ROOT.TMath.Power(universes[i].GetBinContent(j)-avg_hist.GetBinContent(j),2)
        error_hist.SetBinContent(j,ROOT.TMath.Sqrt(error))
    if asFrac:
        error_hist.Divide(avg_hist)
    return error_hist

# Helper function that makes a covariance matrix out of a bunch of universes
# by calculating rms of universes in each bin
def MakeCovariance(universes,asFrac=True):
    n_universes = len(universes)
    avg_hist = universes[0].Clone(universes[0].GetName()+"_avg")
    avg_hist.Reset()
    nbins = universes[0].GetNbinsX()

    for i in range(n_universes):
        avg_hist.Add(universes[i])
    avg_hist.Scale(1/float(n_universes))

    covariance = ROOT.TH2D("Flux_Covariance","Flux_Covariance",nbins,0.5,nbins+0.5,nbins,0.5,nbins+0.5)
    # Calculate Error = RMS of universes
    for j in range(avg_hist.GetNbinsX()):
      for k in range(avg_hist.GetNbinsX()):
        cov = 0
        for i in range(n_universes):
            cov += 1/float(n_universes)*(universes[i].GetBinContent(j+1)-avg_hist.GetBinContent(j+1))*(universes[i].GetBinContent(k+1)-avg_hist.GetBinContent(k+1))
        if asFrac:
            if(avg_hist.GetBinContent(j+1)*avg_hist.GetBinContent(k+1) != 0):
                cov = cov / avg_hist.GetBinContent(j+1) / avg_hist.GetBinContent(k+1)
            else:
                cov = 0
        covariance.SetBinContent(j+1,k+1,cov)
    return covariance

def MergeHistograms(cvs,universes,name):
    nbins = 0
    for cv in cvs:
        nbins += cv.GetNbinsX()
    new_cv = ROOT.TH1D("Flux_CV","Flux_CV",nbins,0.5,nbins+0.5)
    new_universes = []
    for universe in range(len(universes[0])):
        new_universes.append(ROOT.TH1D("Flux_Universe_"+str(universe)+name,"Flux_Universe_"+str(universe)+name,nbins,0.5,nbins+0.5))
    global_bin_iter = 1
    for cv in range(len(cvs)):
        for local_bin_iter in range(cvs[cv].GetNbinsX()):

            new_cv.SetBinContent(global_bin_iter,cvs[cv].GetBinContent(local_bin_iter+1))
            for universe in range(len(universes[0])):
                new_universes[universe].SetBinContent(global_bin_iter,universes[cv][universe].GetBinContent(local_bin_iter+1))
            global_bin_iter = global_bin_iter+1
    return_array = [new_cv,new_universes]
    return return_array

def GetCorrelation(cov_matrix):
    cor_matrix = cov_matrix.Clone()
    for i in range(cor_matrix.GetNbinsX()):
        for j in range(cor_matrix.GetNbinsX()):
            if cov_matrix.GetBinContent(i+1,i+1)!=0 and cov_matrix.GetBinContent(j+1,j+1)!=0:
                cor_matrix.SetBinContent(i+1,j+1,cov_matrix.GetBinContent(i+1,j+1)/ROOT.TMath.Sqrt(cov_matrix.GetBinContent(i+1,i+1))/ROOT.TMath.Sqrt(cov_matrix.GetBinContent(j+1,j+1)))
            else:
                if i==j: cor_matrix.SetBinContent(i+1,j+1,1)
                else: cor_matrix.SetBinContent(i+1,j+1,0)
    cor_matrix.SetMinimum(-1.0)
    cor_matrix.SetMaximum(1.0)
    return cor_matrix


# Open ppfx output file
in_files = []
if(beam_config=="reference"):
    in_files.append(ROOT.TFile("amit_histos_10August2016/histos_g4lbne_v3r4p3_QGSP_BERT_DUNECDR_Ref80GeV_neutrino_LBNFND_ppfx_500_rebinned.root"))
    in_files.append(ROOT.TFile("amit_histos_10August2016/histos_g4lbne_v3r4p3_QGSP_BERT_DUNECDR_Ref80GeV_antineutrino_LBNFND_ppfx_500_rebinned.root"))
    in_files.append(ROOT.TFile("amit_histos_10August2016/histos_g4lbne_v3r4p3_QGSP_BERT_DUNECDR_Ref80GeV_neutrino_LBNFFD_ppfx_500_rebinned.root"))
    in_files.append(ROOT.TFile("amit_histos_10August2016/histos_g4lbne_v3r4p3_QGSP_BERT_DUNECDR_Ref80GeV_antineutrino_LBNFFD_ppfx_500_rebinned.root"))

else:
    in_files.append(ROOT.TFile("amit_histos_10August2016/histos_g4lbne_v3r4p3_QGSP_BERT_CP_run15_12388_neutrino_LBNFND_500_rebinned.root"))
    in_files.append(ROOT.TFile("amit_histos_10August2016/histos_g4lbne_v3r4p3_QGSP_BERT_CP_run15_12388_antineutrino_LBNFND_500_rebinned.root"))
    in_files.append(ROOT.TFile("amit_histos_10August2016/histos_g4lbne_v3r4p3_QGSP_BERT_CP_run15_12388_neutrino_LBNFFD_500_rebinned.root"))
    in_files.append(ROOT.TFile("amit_histos_10August2016/histos_g4lbne_v3r4p3_QGSP_BERT_CP_run15_12388_antineutrino_LBNFFD_500_rebinned.root"))

# Get the central value flux out
cv_histos = []
for in_file in in_files:
    cv_histos.append(in_file.Get("nom/hcv_numu"))
    cv_histos.append(in_file.Get("nom/hcv_numubar"))
    cv_histos.append(in_file.Get("nom/hcv_nue"))
    cv_histos.append(in_file.Get("nom/hcv_nuebar"))

# Get the error universes out
n_universes = 100
error_universes = []
error_types = ["total","others","targetpcpi","targetpck","targetncpi","targetpcnu","targetmesonInc","targetnuA","tgtabsorption","totabsorption"]
error_types2 = ["total","others","targetpcpi","targetpcK","targetncpi","targetpcnu","targetmesonInc","targetnucleonA","targetabsorption","absorption"]
error_titles = ["Total","Other","pC#rightarrow#pi","pC#rightarrowK","nC#rightarrow#pi","pc#rightarrowN","Meson Inc","NucleonA","Target Absorption","Other Absorption"]
for error in error_types:
    error_universes.append([])
          
for in_file in in_files:
  for error_iter in range(len(error_types)):
    error_universes[error_iter].append([])
    for i in range(n_universes):
        error_universes[error_iter][-1].append(in_file.Get("numu_"+error_types[error_iter]+"/h"+error_types2[error_iter]+"_numu_"+str(i)))

    error_universes[error_iter].append([])
    for i in range(n_universes):
        error_universes[error_iter][-1].append(in_file.Get("numubar_"+error_types[error_iter]+"/h"+error_types2[error_iter]+"_numubar_"+str(i)))

    error_universes[error_iter].append([])
    for i in range(n_universes):
        error_universes[error_iter][-1].append(in_file.Get("nue_"+error_types[error_iter]+"/h"+error_types2[error_iter]+"_nue_"+str(i)))

    error_universes[error_iter].append([])
    for i in range(n_universes):
        error_universes[error_iter][-1].append(in_file.Get("nuebar_"+error_types[error_iter]+"/h"+error_types2[error_iter]+"_nuebar_"+str(i)))

# Rebin
numu_bins = [0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,7 ,8,12,16,20,40,100] 
flavors = ["neutrino_near_numu","neutrino_near_numubar","neutrino_near_nue","neutrino_near_nuebar","antineutrino_near_numu","antineutrino_near_numubar","antineutrino_near_nue","antineutrino_near_nuebar","neutrino_far_numu","neutrino_far_numubar","neutrino_far_nue","neutrino_far_nuebar","antineutrino_far_numu","antineutrino_far_numubar","antineutrino_far_nue","antineutrino_far_nuebar"]
#for i in range(len(flavors)):
#    cv_histos[i] = cv_histos[i].Rebin(len(new_bins)-1,"cv_flux_"+flavors[i],array.array('d',new_bins))
#    for j in range(n_universes):
#        total_error_universes[i][j] = total_error_universes[i][j].Rebin(len(new_bins)-1,"universe_"+str(i)+"_"+flavors[i],array.array('d',new_bins))

# Merge into one histogram
merged_histograms = []
merged_universes = []
flux_uncertainties = []
for error_iter in range(len(error_types)):

    merged_histograms.append(MergeHistograms(cv_histos,error_universes[error_iter],error_types[error_iter]))
    merged_universes.append(merged_histograms[error_iter][1])

    flux_uncertainties.append(MakeErrorHistogram(merged_universes[error_iter]))


c1 = ROOT.TCanvas("c1","",1000,800)

flux_uncertainties[error_types.index("total")].Draw()
print "BLAH6",flux_uncertainties[error_types.index("total")].GetBinContent(start_indices[plots.index("numubar_antineutrino_FD")]+5)
c1.Print("HP_error_"+beam_config+".eps")
c1.Print("HP_error_"+beam_config+".png")

flux_covariance = MakeCovariance(merged_universes[error_types.index("total")])
# For some reason, the total doesn't equal the sum of the individual componenents
flux_covariance.Reset()
for i in range(1,len(error_types)):
    flux_covariance.Add(MakeCovariance(merged_universes[i]))

flux_covariance.SetMaximum(0.06)
flux_covariance.SetMinimum(-0.005)
flux_covariance.SetLabelSize(0.001,"X")
flux_covariance.SetLabelSize(0.001,"Y")
flux_covariance.SetTitle("Flux Hadron Production Covariance Matrix")
flux_covariance.Draw("COLZ")

latex0 = ROOT.TLatex( 0.5, 0.93, "Flux Hadron Production Fractional Covariance Matrix" );
latex0.SetTextSize(0.04)
latex0.SetTextAlign(22)
latex0.SetNDC()
latex0.SetTextFont(42)
latex0.Draw()

if(1>0):
    latex1 = ROOT.TLatex( 0.24, 0.035, "Near Detector" );
    latex1.SetTextSize(0.035);
    latex1.SetNDC()
    latex1.SetTextFont(42);     
    latex1.Draw()
    latex2 = latex1.Clone()
    latex2.SetText(0.61,0.035,"Far Detector")
    latex2.Draw()
    latex3 = latex1.Clone()
    latex3.SetText(0.05,0.24,"Near Detector")
    latex3.SetTextAngle(90)
    latex3.Draw()
    latex4 = latex1.Clone()
    latex4.SetText(0.05,0.62,"Far Detector")
    latex4.SetTextAngle(90)
    latex4.Draw()
    
    latex5 = ROOT.TLatex(0.2,0.08,"#nu Mode")
    latex5.SetTextSize(0.025);
    latex5.SetNDC()
    latex5.SetTextFont(42)
    latex5.Draw()
    latex6 = latex5.Clone()
    latex6.SetText(0.38,0.08,"#bar{#nu} Mode")
    latex6.Draw()
    latex7 = latex5.Clone()
    latex7.SetText(0.56,0.08,"#nu Mode")
    latex7.Draw()
    latex8 = latex5.Clone()
    latex8.SetText(0.74,0.08,"#bar{#nu} Mode")
    latex8.Draw()
    
    latex9 = ROOT.TLatex(0.085,0.21,"#nu Mode")
    latex9.SetTextSize(0.025);
    latex9.SetNDC()
    latex9.SetTextFont(42)
    latex9.SetTextAngle(90)
    latex9.Draw()
    latex10 = latex9.Clone()
    latex10.SetText(0.085,0.4,"#bar{#nu} Mode")
    latex10.Draw()
    latex11 = latex9.Clone()
    latex11.SetText(0.085,0.59,"#nu Mode")
    latex11.Draw()
    latex12 = latex9.Clone()
    latex12.SetText(0.085,0.77,"#bar{#nu} Mode")
    latex12.Draw()

    latex13 = ROOT.TLatex( 0.175, 0.124, "#nu_{#mu}" );
    latex13.SetTextSize(0.02);
    latex13.SetNDC()
    latex13.SetTextFont(42);     
    latex13.Draw()
    latex14 = ROOT.TLatex( 0.245, 0.124, "#bar{#nu}_{#mu}" );
    latex14.SetTextSize(0.02);
    latex14.SetNDC()
    latex14.SetTextFont(42);     
    latex14.Draw()
    latex15 = ROOT.TLatex( 0.285, 0.124, "#nu_{e}" );
    latex15.SetTextSize(0.02);
    latex15.SetNDC()
    latex15.SetTextFont(42);     
    latex15.Draw()
    latex16 = ROOT.TLatex( 0.31, 0.124, "#bar{#nu}_{e}" );
    latex16.SetTextSize(0.02);
    latex16.SetNDC()
    latex16.SetTextFont(42);     
    latex16.Draw()

    latex17 = ROOT.TLatex( 0.35, 0.124, "#nu_{#mu}" );
    latex17.SetTextSize(0.02);
    latex17.SetNDC()
    latex17.SetTextFont(42);     
    latex17.Draw()
    latex18 = ROOT.TLatex( 0.415, 0.124, "#bar{#nu}_{#mu}" );
    latex18.SetTextSize(0.02);
    latex18.SetNDC()
    latex18.SetTextFont(42);     
    latex18.Draw()
    latex19 = ROOT.TLatex( 0.46, 0.124, "#nu_{e}" );
    latex19.SetTextSize(0.02);
    latex19.SetNDC()
    latex19.SetTextFont(42);     
    latex19.Draw()
    latex20 = ROOT.TLatex( 0.485, 0.124, "#bar{#nu}_{e}" );
    latex20.SetTextSize(0.02);
    latex20.SetNDC()
    latex20.SetTextFont(42);     
    latex20.Draw()

    latex21 = ROOT.TLatex( 0.53, 0.124, "#nu_{#mu}" );
    latex21.SetTextSize(0.02);
    latex21.SetNDC()
    latex21.SetTextFont(42);     
    latex21.Draw()
    latex22 = ROOT.TLatex( 0.595, 0.124, "#bar{#nu}_{#mu}" );
    latex22.SetTextSize(0.02);
    latex22.SetNDC()
    latex22.SetTextFont(42);     
    latex22.Draw()
    latex23 = ROOT.TLatex( 0.631, 0.124, "#nu_{e}" );
    latex23.SetTextSize(0.02);
    latex23.SetNDC()
    latex23.SetTextFont(42);     
    latex23.Draw()
    latex24 = ROOT.TLatex( 0.654, 0.124, "#bar{#nu}_{e}" );
    latex24.SetTextSize(0.02);
    latex24.SetNDC()
    latex24.SetTextFont(42);     
    latex24.Draw()

    latex25 = ROOT.TLatex( 0.70, 0.124, "#nu_{#mu}" );
    latex25.SetTextSize(0.02);
    latex25.SetNDC()
    latex25.SetTextFont(42);     
    latex25.Draw()
    latex26 = ROOT.TLatex( 0.77, 0.124, "#bar{#nu}_{#mu}" );
    latex26.SetTextSize(0.02);
    latex26.SetNDC()
    latex26.SetTextFont(42);     
    latex26.Draw()
    latex27 = ROOT.TLatex( 0.81, 0.124, "#nu_{e}" );
    latex27.SetTextSize(0.02);
    latex27.SetNDC()
    latex27.SetTextFont(42);     
    latex27.Draw()
    latex28 = ROOT.TLatex( 0.835, 0.124, "#bar{#nu}_{e}" );
    latex28.SetTextSize(0.02);
    latex28.SetNDC()
    latex28.SetTextFont(42);     
    latex28.Draw()

    latex29 = ROOT.TLatex( 0.11,0.18, "#nu_{#mu}" );
    latex29.SetTextAngle(90)
    latex29.SetTextSize(0.02);
    latex29.SetNDC()
    latex29.SetTextFont(42);     
    latex29.Draw()
    latex30 = ROOT.TLatex( 0.11,0.25, "#bar{#nu}_{#mu}" );
    latex30.SetTextAngle(90)
    latex30.SetTextSize(0.02);
    latex30.SetNDC()
    latex30.SetTextFont(42);     
    latex30.Draw()
    latex31 = ROOT.TLatex( 0.11,0.295, "#nu_{e}" );
    latex31.SetTextAngle(90)
    latex31.SetTextSize(0.02);
    latex31.SetNDC()
    latex31.SetTextFont(42);     
    latex31.Draw()
    latex32 = ROOT.TLatex( 0.11,0.32, "#bar{#nu}_{e}" );
    latex32.SetTextAngle(90)
    latex32.SetTextSize(0.02);
    latex32.SetNDC()
    latex32.SetTextFont(42);     
    latex32.Draw()

    latex33 = ROOT.TLatex( 0.11,0.37, "#nu_{#mu}" );
    latex33.SetTextAngle(90)
    latex33.SetTextSize(0.02);
    latex33.SetNDC()
    latex33.SetTextFont(42);     
    latex33.Draw()
    latex34 = ROOT.TLatex( 0.11,0.435, "#bar{#nu}_{#mu}" );
    latex34.SetTextAngle(90)
    latex34.SetTextSize(0.02);
    latex34.SetNDC()
    latex34.SetTextFont(42);     
    latex34.Draw()
    latex35 = ROOT.TLatex( 0.11,0.48, "#nu_{e}" );
    latex35.SetTextAngle(90)
    latex35.SetTextSize(0.02);
    latex35.SetNDC()
    latex35.SetTextFont(42);     
    latex35.Draw()
    latex36 = ROOT.TLatex( 0.11,0.505, "#bar{#nu}_{e}" );
    latex36.SetTextAngle(90)
    latex36.SetTextSize(0.02);
    latex36.SetNDC()
    latex36.SetTextFont(42);     
    latex36.Draw()

    latex37 = ROOT.TLatex( 0.11,0.56, "#nu_{#mu}" );
    latex37.SetTextAngle(90)
    latex37.SetTextSize(0.02);
    latex37.SetNDC()
    latex37.SetTextFont(42);     
    latex37.Draw()
    latex38 = ROOT.TLatex( 0.11,0.62, "#bar{#nu}_{#mu}" );
    latex38.SetTextAngle(90)
    latex38.SetTextSize(0.02);
    latex38.SetNDC()
    latex38.SetTextFont(42);     
    latex38.Draw()
    latex39 = ROOT.TLatex( 0.11,0.67, "#nu_{e}" );
    latex39.SetTextAngle(90)
    latex39.SetTextSize(0.02);
    latex39.SetNDC()
    latex39.SetTextFont(42);     
    latex39.Draw()
    latex40 = ROOT.TLatex( 0.11,0.70, "#bar{#nu}_{e}" );
    latex40.SetTextAngle(90)
    latex40.SetTextSize(0.02);
    latex40.SetNDC()
    latex40.SetTextFont(42);     
    latex40.Draw()

    latex41 = ROOT.TLatex( 0.11,0.75, "#nu_{#mu}" );
    latex41.SetTextAngle(90)
    latex41.SetTextSize(0.02);
    latex41.SetNDC()
    latex41.SetTextFont(42);     
    latex41.Draw()
    latex42 = ROOT.TLatex( 0.11,0.81, "#bar{#nu}_{#mu}" );
    latex42.SetTextAngle(90)
    latex42.SetTextSize(0.02);
    latex42.SetNDC()
    latex42.SetTextFont(42);     
    latex42.Draw()
    latex43 = ROOT.TLatex( 0.11,0.855, "#nu_{e}" );
    latex43.SetTextAngle(90)
    latex43.SetTextSize(0.02);
    latex43.SetNDC()
    latex43.SetTextFont(42);     
    latex43.Draw()
    latex44 = ROOT.TLatex( 0.11,0.88, "#bar{#nu}_{e}" );
    latex44.SetTextAngle(90)
    latex44.SetTextSize(0.02);
    latex44.SetNDC()
    latex44.SetTextFont(42);     
    latex44.Draw()
    
    line_1 = ROOT.TLine(19.5,0.5,19.5,208)
    line_2 = ROOT.TLine(38.5,0.5,38.5,208)
    line_3 = ROOT.TLine(45.5,0.5,45.5,208)
    line_4 = ROOT.TLine(52.5,0.5,52.5,208)
    line_5 = ROOT.TLine(71.5,0.5,71.5,208)
    line_6 = ROOT.TLine(90.5,0.5,90.5,208)
    line_7 = ROOT.TLine(97.5,0.5,97.5,208)
    line_8 = ROOT.TLine(104.5,0.5,104.5,208)
    line_9 = ROOT.TLine(123.5,0.5,123.5,208)
    line_10 = ROOT.TLine(142.5,0.5,142.5,208)
    line_11 = ROOT.TLine(149.5,0.5,149.5,208)
    line_12 = ROOT.TLine(156.5,0.5,156.5,208)
    line_13 = ROOT.TLine(175.5,0.5,175.5,208)
    line_14 = ROOT.TLine(194.5,0.5,194.5,208)
    line_15 = ROOT.TLine(201.5,0.5,201.5,208)
    lineb_1 = ROOT.TLine(0.5,19.5,208,19.5)
    lineb_2 = ROOT.TLine(0.5,38.5,208,38.5)
    lineb_3 = ROOT.TLine(0.5,45.5,208,45.5)
    lineb_4 = ROOT.TLine(0.5,52.5,208,52.5)
    lineb_5 = ROOT.TLine(0.5,71.5,208,71.5)
    lineb_6 = ROOT.TLine(0.5,90.5,208,90.5)
    lineb_7 = ROOT.TLine(0.5,97.5,208,97.5)
    lineb_8 = ROOT.TLine(0.5,104.5,208,104.5)
    lineb_9 = ROOT.TLine(0.5,123.5,208,123.5)
    lineb_10 = ROOT.TLine(0.5,142.5,208,142.5)
    lineb_11 = ROOT.TLine(0.5,149.5,208,149.5)
    lineb_12 = ROOT.TLine(0.5,156.5,208,156.5)
    lineb_13 = ROOT.TLine(0.5,175.5,208,175.5)
    lineb_14 = ROOT.TLine(0.5,194.5,208,194.5)
    lineb_15 = ROOT.TLine(0.5,201.5,208,201.5)
    
    line_8.SetLineWidth(4)
    lineb_8.SetLineWidth(4)
    
    line_4.SetLineWidth(2)
    lineb_4.SetLineWidth(2)
    line_12.SetLineWidth(2)
    lineb_12.SetLineWidth(2)

    for line in [line_1,line_2,line_3,line_4,line_5,line_6,line_7,line_8,line_9,line_10,line_11,line_12,line_13,line_14,line_15,lineb_1,lineb_2,lineb_3,lineb_4,lineb_5,lineb_6,lineb_7,lineb_8,lineb_9,lineb_10,lineb_11,lineb_12,lineb_13,lineb_14,lineb_15]:
        line.SetLineColor(13)
        line.Draw()
    

c1.Print("HP_covariance_"+beam_config+".eps")
c1.Print("HP_covariance_"+beam_config+".png")

# Draw total matrix (with focusing matrix included)
f_foc = ROOT.TFile("focusing_covariance_"+beam_config+".root")
covariance_focusing = f_foc.Get("Focusing_Covariance")

print "BLAH5",ROOT.TMath.Sqrt(flux_covariance.GetBinContent(start_indices[plots.index("numubar_antineutrino_FD")]+5,start_indices[plots.index("numubar_antineutrino_FD")]+5))


total_covariance = flux_covariance.Clone("total_covariance")
total_covariance.Add(covariance_focusing)

total_covariance.SetMaximum(0.06)
total_covariance.SetMinimum(-0.005)
total_covariance.SetLabelSize(0.001,"X")
total_covariance.SetLabelSize(0.001,"Y")
total_covariance.SetTitle("Fractional Flux Covariance Matrix")
total_covariance.Draw("COLZ")

latex0b = ROOT.TLatex( 0.5, 0.93, "Fractional Flux Covariance Matrix" );
latex0b.SetTextSize(0.04)
latex0b.SetTextAlign(22)
latex0b.SetNDC()
latex0b.SetTextFont(42)
latex0b.Draw()

for line in [line_1,line_2,line_3,line_4,line_5,line_6,line_7,line_8,line_9,line_10,line_11,line_12,line_13,line_14,line_15,lineb_1,lineb_2,lineb_3,lineb_4,lineb_5,lineb_6,lineb_7,lineb_8,lineb_9,lineb_10,lineb_11,lineb_12,lineb_13,lineb_14,lineb_15]:
        line.SetLineColor(13)
        line.Draw()

for latex in [latex1,latex2,latex3,latex4,latex5,latex6,latex7,latex8,latex9,latex10,latex11,latex12,latex13,latex14,latex15,latex16,latex17,latex18,latex19,latex20,latex21,latex22,latex23,latex24,latex25,latex26,latex27,latex28,latex29,latex30,latex31,latex32,latex33,latex34,latex35,latex35,latex36,latex37,latex38,latex39,latex40,latex41,latex42,latex43,latex44]:
        latex.Draw();

c1.Print("Total_covariance_"+beam_config+".eps")
c1.Print("Total_covariance_"+beam_config+".png")

f_out = ROOT.TFile("total_covariance_DUNE_"+beam_config+".root","RECREATE")
total_covariance.Write()
f_out.Close()

total_error = ROOT.TH1D("Total Error","Total Error",208,0.5,208.5);
for i in range(208):
    total_error.SetBinContent(i+1,ROOT.TMath.Sqrt(total_covariance.GetBinContent(i+1,i+1)))
total_error.SetLabelSize(0.001,"X")
total_error.SetLabelSize(0.05,"Y")
total_error.GetYaxis().SetTitleOffset(1.3)
total_error.GetYaxis().SetTitle("Fractional Flux Uncertainty")

total_error.SetMaximum(0.4)
total_error.SetMinimum(0)
total_error.Draw()


linec_1 = ROOT.TLine(19.5,0,19.5,total_error.GetMaximum())
linec_2 = ROOT.TLine(38.5,0,38.5,total_error.GetMaximum())
linec_3 = ROOT.TLine(45.5,0,45.5,total_error.GetMaximum())
linec_4 = ROOT.TLine(52.5,0,52.5,total_error.GetMaximum())
linec_5 = ROOT.TLine(71.5,0,71.5,total_error.GetMaximum())
linec_6 = ROOT.TLine(90.5,0,90.5,total_error.GetMaximum())
linec_7 = ROOT.TLine(97.5,0,97.5,total_error.GetMaximum())
linec_8 = ROOT.TLine(104.5,0,104.5,total_error.GetMaximum())
linec_9 = ROOT.TLine(123.5,0,123.5,total_error.GetMaximum())
linec_10 = ROOT.TLine(142.5,0,142.5,total_error.GetMaximum())
linec_11 = ROOT.TLine(149.5,0,149.5,total_error.GetMaximum())
linec_12 = ROOT.TLine(156.5,0,156.5,total_error.GetMaximum())
linec_13 = ROOT.TLine(175.5,0,175.5,total_error.GetMaximum())
linec_14 = ROOT.TLine(194.5,0,194.5,total_error.GetMaximum())
linec_15 = ROOT.TLine(201.5,0,201.5,total_error.GetMaximum())

for line in [linec_1,linec_2,linec_3,linec_4,linec_5,linec_6,linec_7,linec_8,linec_9,linec_10,linec_11,linec_12,linec_13,linec_14,linec_15]:
        line.SetLineColor(2)
        line.SetLineStyle(2)
        line.Draw()

for latex in [latex1,latex2,latex5,latex6,latex7,latex8,latex13,latex14,latex15,latex16,latex17,latex18,latex19,latex20,latex21,latex22,latex23,latex24,latex25,latex26,latex27,latex28]:
        latex.Draw();

c1.Print("Total_error_"+beam_config+".eps")
c1.Print("Total_error_"+beam_config+".png")

"""
const int NRGBs = 3, NCont = 999;
  gStyle.xSetNumberContours(NCont);
  Double_t stops[NRGBs] = { 0.00, 0.50, 1.00};
  Double_t red[NRGBs]   = { 0.00, 1.00, 1.00};
  Double_t green[NRGBs] = { 0.00, 1.00, 0.00};
  Double_t blue[NRGBs]  = { 1.00, 1.00, 0.00};
  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
"""
NRGBs = 3;
NCont = 999;
ROOT.gStyle.SetNumberContours(NCont);
stops = [ 0.00, 0.50, 1.00];
red   = [ 0.00, 1.00, 1.00];
green = [ 0.00, 1.00, 0.00];
blue  = [ 1.00, 1.00, 0.00];
ROOT.TColor.CreateGradientColorTable(NRGBs, array.array('d',stops), array.array('d',red), array.array('d',green), array.array('d',blue), NCont);

correlation = GetCorrelation(total_covariance)
correlation.Draw("COLZ")

for line in [line_1,line_2,line_3,line_4,line_5,line_6,line_7,line_8,line_9,line_10,line_11,line_12,line_13,line_14,line_15,lineb_1,lineb_2,lineb_3,lineb_4,lineb_5,lineb_6,lineb_7,lineb_8,lineb_9,lineb_10,lineb_11,lineb_12,lineb_13,lineb_14,lineb_15]:
        line.SetLineColor(13)
        line.Draw()
for latex in [latex1,latex2,latex3,latex4,latex5,latex6,latex7,latex8,latex9,latex10,latex11,latex12,latex13,latex14,latex15,latex16,latex17,latex18,latex19,latex20,latex21,latex22,latex23,latex24,latex25,latex26,latex27,latex28,latex29,latex30,latex31,latex32,latex33,latex34,latex35,latex36,latex37,latex38,latex39,latex40,latex41,latex42,latex43,latex44]:
        latex.Draw();

c1.Print("total_correlation_"+beam_config+".eps")
c1.Print("total_correlation_"+beam_config+".png")

# Draw pretty error histograms

for plot in plots:
    error_bins = [0,2,4,6,8,10]
    if plot.startswith("numu"):
        error_bins = [0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,7 ,8,12,20,40] 
    error_hists = []
    total_uncertainties = [0]*(len(error_bins)-1)
    for i in range(1,len(error_types)):

        temp_hist = ROOT.TH1D(error_types[i],error_types[i],len(error_bins)-1,array.array('d',error_bins))
        for j in range(len(error_bins)-1):
                temp_hist.SetBinContent(j+1,abs(flux_uncertainties[i].GetBinContent(start_indices[plots.index(plot)]+j+1)))
                total_uncertainties[j] = total_uncertainties[j]+math.pow(flux_uncertainties[i].GetBinContent(start_indices[plots.index(plot)]+j+1),2)
        error_hists.append(temp_hist)

    temp_hist = ROOT.TH1D("Total","Total",len(error_bins)-1,array.array('d',error_bins))
    for j in range(len(error_bins)-1):
        temp_hist.SetBinContent(j+1,math.sqrt(total_uncertainties[j]))
    error_hists.append(temp_hist)
                                     
    error_hists[-1].SetMinimum(0);
    error_hists[-1].SetMaximum(0.2);
    error_hists[-1].SetLineWidth(3)
    error_hists[-1].Draw("hist")
    error_hists[-1].GetXaxis().SetTitle("Neutrino Energy (GeV)")
    error_hists[-1].GetYaxis().SetTitle("Fractional Uncertainty")
    error_hists[-1].GetYaxis().SetTitleOffset(1.3)
    for i in range(0,len(error_hists)-1):
            colors = [2,ROOT.kOrange+3,6,ROOT.kOrange+7,ROOT.kGreen+2,3,ROOT.kAzure+7,ROOT.kBlue-4,7,ROOT.kViolet+1,13,15,17,19]
            error_hists[i].SetLineColor(colors[i])
            error_hists[i].SetLineWidth(3)
            error_hists[i].Draw("histsame")
    
    leg = ROOT.TLegend(0.25,0.6,0.5,0.88)
    leg.AddEntry(error_hists[-1],"Total","l")
    for i in range(0,len(error_hists)-1):
        leg.AddEntry(error_hists[i],error_titles[i+1],"l")
    leg.SetFillStyle(0)
    leg.Draw()

    c1.Print("HP_error_overlay_"+plot+"_"+beam_config+".eps")        
    c1.Print("HP_error_overlay_"+plot+"_"+beam_config+".png")

    # Draw plots with focusing and hp errors overlaid
    focusing_unc = ROOT.TH1D("Focusing","Focusing",len(error_bins)-1,array.array('d',error_bins))
    total_unc = ROOT.TH1D("Total","Total",len(error_bins)-1,array.array('d',error_bins))
    for j in range(len(error_bins)-1):
        focusing_unc.SetBinContent(j+1,ROOT.TMath.Sqrt(covariance_focusing.GetBinContent(start_indices[plots.index(plot)]+j+1,start_indices[plots.index(plot)]+j+1)))
        total_unc.SetBinContent(j+1,ROOT.TMath.Sqrt(total_covariance.GetBinContent(start_indices[plots.index(plot)]+j+1,start_indices[plots.index(plot)]+j+1)))

    total_unc.SetMinimum(0)
    total_unc.SetMaximum(0.2)
    total_unc.SetLineWidth(3)
    focusing_unc.SetLineWidth(3)
    total_unc.GetXaxis().SetTitle("Neutrino Energy (GeV)")
    total_unc.GetYaxis().SetTitle("Fractional Uncertainty")
    total_unc.GetYaxis().SetTitleOffset(1.3)
    total_unc.Draw("hist")
    focusing_unc.SetLineColor(2)
    error_hists[-1].SetLineColor(4)
    focusing_unc.Draw("histsame")
    error_hists[-1].Draw("histsame")

    if plot=="numubar_antineutrino_FD":
        print "blah3",total_unc.GetBinContent(5),focusing_unc.GetBinContent(5),error_hists[-1].GetBinContent(5)
        print "BLAH4",ROOT.TMath.Sqrt(total_covariance.GetBinContent(start_indices[plots.index("numubar_antineutrino_FD")]+5,start_indices[plots.index("numubar_antineutrino_FD")]+5))
        total_unc.Print("ALL")
        print plots.index("numubar_antineutrino_FD"),plots.index(plot)


    leg2 = ROOT.TLegend(0.25,0.7,0.6,0.88)
    leg2.AddEntry(total_unc,"Total","l")
    leg2.AddEntry(error_hists[-1],"Hadron Production","l")
    leg2.AddEntry(focusing_unc,"Focusing","l")
    leg2.SetFillStyle(0)
    leg2.Draw()

    c1.Print("error_overlay_"+plot+"_"+beam_config+".eps")        
    c1.Print("error_overlay_"+plot+"_"+beam_config+".png")