drawSensitivitySummaries.py

Go to the documentation of this file.

import ROOT, os, array, sys

use_realistic_powers = True

def GetMedian(histo):
    n = histo.GetNbinsX();
    x = []
    y = []
    
    for i in range(0,n):
        x.append(histo.GetBinCenter(i+1));
        y.append(histo.GetBinContent(i+1));
    y_array = array.array('d',y)
    
    return ROOT.TMath.Median(n,y_array);

def GetPercentile(histo,X):
    n = histo.GetNbinsX();
    x = []
    y = []
    
    for i in range(0,n):
        x.append(histo.GetBinCenter(i+1));
        y.append(histo.GetBinContent(i+1));
    y_array = array.array('d',y)

    percentiles = array.array('d',[0.0]);
    probs = array.array('d',[1-float(X)/100.0])

    ROOT.TMath.Quantiles(n,1,y_array,percentiles,probs,False);
    return percentiles[0]

def GetPercentileCheck(histo,X):
    n = histo.GetNbinsX();
    x = []
    y = []
    
    for i in range(0,n):
        x.append(histo.GetBinCenter(i+1));
        y.append(histo.GetBinContent(i+1));
    y.sort()
    print y
    pos = float(100-X)/100.0 * n;
    if ROOT.TMath.Floor(pos)==pos:
        return y[pos]
    else:
        i = ROOT.TMath.FloorNint(pos)
        return (y[i]+y[i+1])/2

def GetFractionAtMoreThanXSigma(histo, X):
    n = histo.GetNbinsX();

    greaterThanX = 0.0;
    total = 0.0;

    for i in range(0,n):
        total = total + float(histo.GetBinWidth(i+1));
        if(histo.GetBinContent(i+1)>=X):
            greaterThanX = greaterThanX + float(histo.GetBinWidth(i+1));

    return greaterThanX/total;


def GetMinimum(histo):
    min = 99999
    for i in range(0,histo.GetNbinsX()):
        if histo.GetBinContent(i+1) < min:
            min = histo.GetBinContent(i+1)
    return min



sigma = 3

energies = [20,30,40,50,60,70,80,90,100,120,130]
energies_array = array.array('d',energies)

nh_cp_medians = [[],[],[]]
ih_cp_medians = [[],[],[]]
nh_cp_75thpercentiles = [[],[],[]]
ih_cp_75thpercentiles = [[],[],[]]
nh_cp_3sigmafracs = [[],[],[]]
ih_cp_3sigmafracs = [[],[],[]]
nh_cp_5sigmafracs = [[],[],[]]
ih_cp_5sigmafracs = [[],[],[]]
nh_mh_minimums = [[],[],[]]
ih_mh_minimums = [[],[],[]]

for energy in energies:
    nh_cp_file = ROOT.TFile("/lbne/data/users/lblpwg_tools/FastMC_Data/outputs/ljf26/Sensitivity_Plots/ProtonP"+str(energy)+"GeV_nh_cp_histos.root"); 
    if(use_realistic_powers and energy!=120):
        nh_cp_file = ROOT.TFile("/lbne/data/users/lblpwg_tools/FastMC_Data/outputs/ljf26/Sensitivity_Plots/ProtonP"+str(energy)+"GeV_RealisticPower_nh_cp_histos.root"); 

    nh_cp_medians[0].append(GetMedian(nh_cp_file.Get("h1")))
    nh_cp_medians[1].append(GetMedian(nh_cp_file.Get("h2")))
    nh_cp_medians[2].append(GetMedian(nh_cp_file.Get("h3")))

    nh_cp_75thpercentiles[0].append(GetPercentile(nh_cp_file.Get("h1"),75))
    nh_cp_75thpercentiles[1].append(GetPercentile(nh_cp_file.Get("h2"),75))
    nh_cp_75thpercentiles[2].append(GetPercentile(nh_cp_file.Get("h3"),75))

    nh_cp_3sigmafracs[0].append(GetFractionAtMoreThanXSigma(nh_cp_file.Get("h1"),3.0)*100)
    nh_cp_3sigmafracs[1].append(GetFractionAtMoreThanXSigma(nh_cp_file.Get("h2"),3.0)*100)
    nh_cp_3sigmafracs[2].append(GetFractionAtMoreThanXSigma(nh_cp_file.Get("h3"),3.0)*100)

    nh_cp_5sigmafracs[0].append(GetFractionAtMoreThanXSigma(nh_cp_file.Get("h1"),5.0)*100)
    nh_cp_5sigmafracs[1].append(GetFractionAtMoreThanXSigma(nh_cp_file.Get("h2"),5.0)*100)
    nh_cp_5sigmafracs[2].append(GetFractionAtMoreThanXSigma(nh_cp_file.Get("h3"),5.0)*100)

    ih_cp_file = ROOT.TFile("/lbne/data/users/lblpwg_tools/FastMC_Data/outputs/ljf26/Sensitivity_Plots/ProtonP"+str(energy)+"GeV_ih_cp_histos.root");
    if(use_realistic_powers and energy!=120):
        ih_cp_file = ROOT.TFile("/lbne/data/users/lblpwg_tools/FastMC_Data/outputs/ljf26/Sensitivity_Plots/ProtonP"+str(energy)+"GeV_RealisticPower_ih_cp_histos.root"); 
        
    ih_cp_medians[0].append(GetMedian(ih_cp_file.Get("h1")))
    ih_cp_medians[1].append(GetMedian(ih_cp_file.Get("h2")))
    ih_cp_medians[2].append(GetMedian(ih_cp_file.Get("h3")))

    ih_cp_75thpercentiles[0].append(GetPercentile(ih_cp_file.Get("h1"),75))
    ih_cp_75thpercentiles[1].append(GetPercentile(ih_cp_file.Get("h2"),75))
    ih_cp_75thpercentiles[2].append(GetPercentile(ih_cp_file.Get("h3"),75))

    ih_cp_3sigmafracs[0].append(GetFractionAtMoreThanXSigma(ih_cp_file.Get("h1"),3.0)*100)
    ih_cp_3sigmafracs[1].append(GetFractionAtMoreThanXSigma(ih_cp_file.Get("h2"),3.0)*100)
    ih_cp_3sigmafracs[2].append(GetFractionAtMoreThanXSigma(ih_cp_file.Get("h3"),3.0)*100)

    ih_cp_5sigmafracs[0].append(GetFractionAtMoreThanXSigma(ih_cp_file.Get("h1"),5.0)*100)
    ih_cp_5sigmafracs[1].append(GetFractionAtMoreThanXSigma(ih_cp_file.Get("h2"),5.0)*100)
    ih_cp_5sigmafracs[2].append(GetFractionAtMoreThanXSigma(ih_cp_file.Get("h3"),5.0)*100)

    nh_mh_file = ROOT.TFile("/lbne/data/users/lblpwg_tools/FastMC_Data/outputs/ljf26/Sensitivity_Plots/ProtonP"+str(energy)+"GeV_nh_mh_histos.root"); 
    if(use_realistic_powers and energy !=120):
        nh_mh_file = ROOT.TFile("/lbne/data/users/lblpwg_tools/FastMC_Data/outputs/ljf26/Sensitivity_Plots/ProtonP"+str(energy)+"GeV_RealisticPower_nh_mh_histos.root"); 

    nh_mh_minimums[0].append(GetMinimum(nh_mh_file.Get("h1")))
    nh_mh_minimums[1].append(GetMinimum(nh_mh_file.Get("h2")))
    nh_mh_minimums[2].append(GetMinimum(nh_mh_file.Get("h3")))

    ih_mh_file = ROOT.TFile("/lbne/data/users/lblpwg_tools/FastMC_Data/outputs/ljf26/Sensitivity_Plots/ProtonP"+str(energy)+"GeV_ih_mh_histos.root"); 
    if(use_realistic_powers and energy !=120):
        ih_mh_file = ROOT.TFile("/lbne/data/users/lblpwg_tools/FastMC_Data/outputs/ljf26/Sensitivity_Plots/ProtonP"+str(energy)+"GeV_RealisticPower_ih_mh_histos.root"); 

    ih_mh_minimums[0].append(GetMinimum(ih_mh_file.Get("h1")))
    ih_mh_minimums[1].append(GetMinimum(ih_mh_file.Get("h2")))
    ih_mh_minimums[2].append(GetMinimum(ih_mh_file.Get("h3")))



ROOT.gStyle.SetMarkerStyle(21);
ROOT.gStyle.SetMarkerStyle(21);

c_nh_cp_m = ROOT.TCanvas("nh_cp_m");
c_ih_cp_m = ROOT.TCanvas("ih_cp_m");
c_nh_cp_p = ROOT.TCanvas("nh_cp_p");
c_ih_cp_p = ROOT.TCanvas("ih_cp_p");
c_nh_cp_3 = ROOT.TCanvas("nh_cp_3");
c_ih_cp_3 = ROOT.TCanvas("ih_cp_3");
c_nh_cp_5 = ROOT.TCanvas("nh_cp_5");
c_ih_cp_5 = ROOT.TCanvas("ih_cp_5");

c_nh_mh = ROOT.TCanvas("nh_mh");
c_ih_mh = ROOT.TCanvas("ih_mh");

nh_cp_mgraphs = []
ih_cp_mgraphs = []
nh_cp_pgraphs = []
ih_cp_pgraphs = []
nh_cp_3graphs = []
ih_cp_3graphs = []
nh_cp_5graphs = []
ih_cp_5graphs = []
nh_mh_mgraphs = []
ih_mh_mgraphs = []

line_styles = [1,2,3]
line_colors = [1,2,4]

for i in 0,1,2:
    temp = ROOT.TGraph(len(energies),energies_array,array.array('d',nh_cp_medians[i]));
    nh_cp_mgraphs.append(temp)
    temp = ROOT.TGraph(len(energies),energies_array,array.array('d',ih_cp_medians[i]));
    ih_cp_mgraphs.append(temp)

    temp = ROOT.TGraph(len(energies),energies_array,array.array('d',nh_cp_75thpercentiles[i]));
    nh_cp_pgraphs.append(temp)
    temp = ROOT.TGraph(len(energies),energies_array,array.array('d',ih_cp_75thpercentiles[i]));
    ih_cp_pgraphs.append(temp)

    temp = ROOT.TGraph(len(energies),energies_array,array.array('d',nh_cp_3sigmafracs[i]));
    nh_cp_3graphs.append(temp)
    temp = ROOT.TGraph(len(energies),energies_array,array.array('d',ih_cp_3sigmafracs[i]));
    ih_cp_3graphs.append(temp)
    temp = ROOT.TGraph(len(energies),energies_array,array.array('d',nh_cp_5sigmafracs[i]));
    nh_cp_5graphs.append(temp)
    temp = ROOT.TGraph(len(energies),energies_array,array.array('d',ih_cp_5sigmafracs[i]));
    ih_cp_5graphs.append(temp)
    temp = ROOT.TGraph(len(energies),energies_array,array.array('d',nh_mh_minimums[i]));
    nh_mh_mgraphs.append(temp)
    temp = ROOT.TGraph(len(energies),energies_array,array.array('d',ih_mh_minimums[i]));
    ih_mh_mgraphs.append(temp)

    nh_cp_mgraphs[i].SetMinimum(1);
    ih_cp_mgraphs[i].SetMinimum(1);
    nh_cp_pgraphs[i].SetMinimum(1);
    ih_cp_pgraphs[i].SetMinimum(1);
    nh_cp_3graphs[i].SetMinimum(30);
    ih_cp_3graphs[i].SetMinimum(30);
    nh_cp_5graphs[i].SetMinimum(30);
    ih_cp_5graphs[i].SetMinimum(30);
    nh_cp_3graphs[i].SetMaximum(80);
    ih_cp_3graphs[i].SetMaximum(80);
    nh_cp_5graphs[i].SetMaximum(100);
    ih_cp_5graphs[i].SetMaximum(100);
    nh_mh_mgraphs[i].SetMinimum(2.5);
    ih_mh_mgraphs[i].SetMinimum(2.5);

    nh_cp_mgraphs[i].GetXaxis().SetTitle("Proton Energy (GeV)");
    ih_cp_mgraphs[i].GetXaxis().SetTitle("Proton Energy (GeV)");
    nh_cp_pgraphs[i].GetXaxis().SetTitle("Proton Energy (GeV)");
    ih_cp_pgraphs[i].GetXaxis().SetTitle("Proton Energy (GeV)");
    nh_cp_3graphs[i].GetXaxis().SetTitle("Proton Energy (GeV)");
    ih_cp_3graphs[i].GetXaxis().SetTitle("Proton Energy (GeV)");
    nh_cp_5graphs[i].GetXaxis().SetTitle("Proton Energy (GeV)");
    ih_cp_5graphs[i].GetXaxis().SetTitle("Proton Energy (GeV)");
    nh_mh_mgraphs[i].GetXaxis().SetTitle("Proton Energy (GeV)");
    ih_mh_mgraphs[i].GetXaxis().SetTitle("Proton Energy (GeV)");

    nh_cp_mgraphs[i].GetYaxis().SetTitle("Median CP Sensitivity (#sigma)");
    ih_cp_mgraphs[i].GetYaxis().SetTitle("Median CP Sensitivity (#sigma)");
    nh_cp_pgraphs[i].GetYaxis().SetTitle("75% CP Sensitivity) (#sigma)");
    ih_cp_pgraphs[i].GetYaxis().SetTitle("75% CP Sensitivity) (#sigma)");
    nh_cp_3graphs[i].GetYaxis().SetTitle("3 Sigma Coverage (%)");
    ih_cp_3graphs[i].GetYaxis().SetTitle("3 Sigma Coverage (%)");
    nh_cp_5graphs[i].GetYaxis().SetTitle("5 Sigma Coverage (%)");
    ih_cp_5graphs[i].GetYaxis().SetTitle("5 Sigma Coverage (%)");
    nh_mh_mgraphs[i].GetYaxis().SetTitle("Minimum MH Sensitivity (#sigma)");
    ih_mh_mgraphs[i].GetYaxis().SetTitle("Minimum MH Sensitivity (#sigma)");

    nh_cp_mgraphs[i].SetTitle("");
    ih_cp_mgraphs[i].SetTitle("");
    nh_cp_pgraphs[i].SetTitle("");
    ih_cp_pgraphs[i].SetTitle("");
    nh_cp_3graphs[i].SetTitle("");
    ih_cp_3graphs[i].SetTitle("");
    nh_cp_5graphs[i].SetTitle("");
    ih_cp_5graphs[i].SetTitle("");
    nh_mh_mgraphs[i].SetTitle("");
    ih_mh_mgraphs[i].SetTitle("");
    
    #nh_cp_mgraphs[i].SetLineStyle(line_styles[i])
    #ih_cp_mgraphs[i].SetLineStyle(line_styles[i])
    #nh_mh_mgraphs[i].SetLineStyle(line_styles[i])
    #ih_mh_mgraphs[i].SetLineStyle(line_styles[i])
    
    nh_cp_mgraphs[i].SetLineColor(line_colors[i])
    ih_cp_mgraphs[i].SetLineColor(line_colors[i])
    nh_cp_pgraphs[i].SetLineColor(line_colors[i])
    ih_cp_pgraphs[i].SetLineColor(line_colors[i])
    nh_cp_3graphs[i].SetLineColor(line_colors[i])
    ih_cp_3graphs[i].SetLineColor(line_colors[i])
    nh_cp_5graphs[i].SetLineColor(line_colors[i])
    ih_cp_5graphs[i].SetLineColor(line_colors[i])
    nh_mh_mgraphs[i].SetLineColor(line_colors[i])
    ih_mh_mgraphs[i].SetLineColor(line_colors[i])

    nh_cp_mgraphs[i].SetMarkerColor(line_colors[i])
    ih_cp_mgraphs[i].SetMarkerColor(line_colors[i])
    nh_cp_pgraphs[i].SetMarkerColor(line_colors[i])
    ih_cp_pgraphs[i].SetMarkerColor(line_colors[i])
    nh_cp_3graphs[i].SetMarkerColor(line_colors[i])
    ih_cp_3graphs[i].SetMarkerColor(line_colors[i])
    nh_cp_5graphs[i].SetMarkerColor(line_colors[i])
    ih_cp_5graphs[i].SetMarkerColor(line_colors[i])
    nh_mh_mgraphs[i].SetMarkerColor(line_colors[i])
    ih_mh_mgraphs[i].SetMarkerColor(line_colors[i])

    c_nh_cp_m.cd()
    if i == 0:
        nh_cp_mgraphs[i].Draw("ACP");
    else:
        nh_cp_mgraphs[i].Draw("SAME CP");

    c_ih_cp_m.cd()
    if i == 0:
        ih_cp_mgraphs[i].Draw("ACP");
    else:
        ih_cp_mgraphs[i].Draw("SAME CP");

    c_nh_cp_p.cd()
    if i == 0:
        nh_cp_pgraphs[i].Draw("ACP");
    else:
        nh_cp_pgraphs[i].Draw("SAME CP");

    c_ih_cp_p.cd()
    if i == 0:
        ih_cp_pgraphs[i].Draw("ACP");
    else:
        ih_cp_pgraphs[i].Draw("SAME CP");

    c_nh_cp_3.cd()
    if i == 0:
        nh_cp_3graphs[i].Draw("ACP");
    else:
        nh_cp_3graphs[i].Draw("SAME CP");

    c_ih_cp_3.cd()
    if i == 0:
        ih_cp_3graphs[i].Draw("ACP");
    else:
        ih_cp_3graphs[i].Draw("SAME CP");

    c_nh_cp_5.cd()
    if i == 0:
        nh_cp_5graphs[i].Draw("ACP");
    else:
        nh_cp_5graphs[i].Draw("SAME CP");

    c_ih_cp_5.cd()
    if i == 0:
        ih_cp_5graphs[i].Draw("ACP");
    else:
        ih_cp_5graphs[i].Draw("SAME CP");

    c_nh_mh.cd()
    if i == 0:
        nh_mh_mgraphs[i].Draw("ACP");
    else:
        nh_mh_mgraphs[i].Draw("SAME CP");

    c_ih_mh.cd()
    if i == 0:
        ih_mh_mgraphs[i].Draw("ACP");
    else:
        ih_mh_mgraphs[i].Draw("SAME CP");

leg = ROOT.TLegend(0.45,0.2,0.85,0.45);
leg.SetHeader("Sig/Bkgd Uncertainties");
leg.AddEntry(nh_cp_mgraphs[0],"1%/5%","lp");
leg.AddEntry(nh_cp_mgraphs[1],"2%/5%","lp");
leg.AddEntry(nh_cp_mgraphs[2],"5%/10%","lp");
leg.SetFillStyle(0);
leg.SetBorderSize(0);

c_nh_cp_m.cd();
leg.Draw();
c_ih_cp_m.cd();
leg.Draw();
c_nh_cp_p.cd();
leg.Draw();
c_ih_cp_p.cd();
leg.Draw();
c_nh_cp_3.cd();
leg.Draw();
c_ih_cp_3.cd();
leg.Draw();
c_nh_cp_5.cd();
leg.Draw();
c_ih_cp_5.cd();
leg.Draw();
c_nh_mh.cd();
leg.Draw();
c_ih_mh.cd();
leg.Draw();


if(use_realistic_powers):

    c_nh_cp_m.Print("summary_rp_nh_cp_median.eps");
    c_nh_cp_m.Print("summary_rp_nh_cp_median.png");
    
    c_ih_cp_m.Print("summary_rp_ih_cp_median.eps");
    c_ih_cp_m.Print("summary_rp_ih_cp_median.png");

    c_nh_cp_p.Print("summary_rp_nh_cp_75thpercentile.eps");
    c_nh_cp_p.Print("summary_rp_nh_cp_75thpercentile.png");
    
    c_ih_cp_p.Print("summary_rp_ih_cp_75thpercentile.eps");
    c_ih_cp_p.Print("summary_rp_ih_cp_75thpercentile.png");

    c_nh_cp_3.Print("summary_rp_nh_cp_3sigmafrac.eps");
    c_nh_cp_3.Print("summary_rp_nh_cp_3sigmafrac.png");
    
    c_ih_cp_3.Print("summary_rp_ih_cp_3sigmafrac.eps");
    c_ih_cp_3.Print("summary_rp_ih_cp_3sigmafrac.png");

    c_nh_cp_5.Print("summary_rp_nh_cp_5sigmafrac.eps");
    c_nh_cp_5.Print("summary_rp_nh_cp_5sigmafrac.png");
    
    c_ih_cp_5.Print("summary_rp_ih_cp_5sigmafrac.eps");
    c_ih_cp_5.Print("summary_rp_ih_cp_5sigma.png");
    
    c_nh_mh.Print("summary_rp_nh_mh.eps");
    c_nh_mh.Print("summary_rp_nh_mh.png");
    
    c_ih_mh.Print("summary_rp_ih_mh.eps");
    c_ih_mh.Print("summary_rp_ih_mh.png");

else:

    c_nh_cp_m.Print("summary_nh_cp_median.eps");
    c_nh_cp_m.Print("summary_nh_cp_median.png");
    
    c_ih_cp_m.Print("summary_ih_cp_median.eps");
    c_ih_cp_m.Print("summary_ih_cp_median.png");

    c_nh_cp_m.Print("summary_nh_cp_75thpercentile.eps");
    c_nh_cp_m.Print("summary_nh_cp_75thperceentile.png");
    
    c_ih_cp_m.Print("summary_ih_cp_75thpercentile.eps");
    c_ih_cp_m.Print("summary_ih_cp_75thpercentile.png");

    c_nh_cp_3.Print("summary_nh_cp_3sigmafrac.eps");
    c_nh_cp_3.Print("summary_nh_cp_3sigmafrac.png");
    
    c_ih_cp_3.Print("summary_ih_cp_3sigmafrac.eps");
    c_ih_cp_3.Print("summary_rp_ih_cp_3sigmafrac.png");

    c_nh_cp_5.Print("summary_nh_cp_5sigmafrac.eps");
    c_nh_cp_5.Print("summary_nh_cp_5sigmafrac.png");
    
    c_ih_cp_5.Print("summary_ih_cp_5sigmafrac.eps");
    c_ih_cp_5.Print("summary_ih_cp_5sigma.png");

    c_nh_mh.Print("summary_nh_mh.eps");
    c_nh_mh.Print("summary_nh_mh.png");
    
    c_ih_mh.Print("summary_ih_mh.eps");
    c_ih_mh.Print("summary_ih_mh.png");