drawSensitivitySummaries_simplevars_improved.py

Go to the documentation of this file.

import ROOT, os, array, sys, math

use_realistic_powers = True

use_fixed_scale = True

tempgraphs = []
for i in range(20):
    tempgraphs.append(0)
def myfunc(x, par):
    return tempgraphs[int(par[0])].Eval(x[0]);

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 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 GetMinimum(histo):
    min = 99999
    for i in range(0,histo.GetNbinsX()):
        if histo.GetBinContent(i+1) < min:
            min = histo.GetBinContent(i+1)
    return min

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

amounts_varied = [0,50,90,110,200,500]

modes = ["fhc","rhc"]
nus = ["numu","numubar","nuenuebar"]
#nus = ["numu","numubar"]
vars = ["0to0p5GeVPlus10per","0p5to1GeVPlus10per","1to2GeVPlus10per","2to3GeVPlus10per","3to4GeVPlus10per","4to5GeVPlus10per","5to6GeVPlus10per","6to7GeVPlus10per","7to8GeVPlus10per","8to9GeVPlus10per","9to10GeVPlus10per","10to15GeVPlus10per","15to20GeVPlus10per","20to120GeVPlus10per"]
bin_centers = [0.25,0.75,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,12.5,17.5,70.0]
#bin_centers = [0.25,0.75,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,12.5,17.5]
bin_centers_array = array.array('d',bin_centers)
hierarchies = ["nh","ih"]
#plotvars = ["cp_median","cp_75thpercentile","cp_3sigma","cp_5sigma","mh_minimum"]
plotvars = ["cp_75thpercentile","mh_minimum"]

for amount_varied in amounts_varied:

  outfile = ROOT.TFile("/lbne/data/users/ljf26/fluxfiles/sensitivity_comps_"+str(amount_varied)+"_improved.root","RECREATE");

  for mode in modes:
    for nu in nus:
        for hierarchy in hierarchies:
            for plotvar in plotvars:
                medians = [[],[],[]]
                # get baseline median vailes
                baseline_medians = []
                mixing_param = plotvar[0:2]
                baseline_file = ROOT.TFile("/lbne/data/users/lblpwg_tools/FastMC_Data/outputs/ljf26/Sensitivity_Plots/ProtonP120GeV_"+hierarchy+"_"+mixing_param+"_histos.root");

                if(plotvar=="cp_median"):
                    baseline_medians.append(GetMedian(baseline_file.Get("h1")))
                    baseline_medians.append(GetMedian(baseline_file.Get("h2")))
                    baseline_medians.append(GetMedian(baseline_file.Get("h3")))
                if(plotvar=="cp_75thpercentile"):
                    baseline_medians.append(GetPercentile(baseline_file.Get("h1"),75))
                    baseline_medians.append(GetPercentile(baseline_file.Get("h2"),75))
                    baseline_medians.append(GetPercentile(baseline_file.Get("h3"),75))
                elif(plotvar=="cp_3sigma"):
                    baseline_medians.append(GetFractionAtMoreThanXSigma(baseline_file.Get("h1"),3.0))
                    baseline_medians.append(GetFractionAtMoreThanXSigma(baseline_file.Get("h2"),3.0))
                    baseline_medians.append(GetFractionAtMoreThanXSigma(baseline_file.Get("h3"),3.0))
                elif(plotvar=="cp_5sigma"):
                    baseline_medians.append(GetFractionAtMoreThanXSigma(baseline_file.Get("h1"),5.0))
                    baseline_medians.append(GetFractionAtMoreThanXSigma(baseline_file.Get("h2"),5.0))
                    baseline_medians.append(GetFractionAtMoreThanXSigma(baseline_file.Get("h3"),5.0))
                elif(plotvar=="mh_minimum"):
                    baseline_medians.append(GetMinimum(baseline_file.Get("h1")))
                    baseline_medians.append(GetMinimum(baseline_file.Get("h2")))
                    baseline_medians.append(GetMinimum(baseline_file.Get("h3")))

                for i in range(0,len(vars)):

                    var = vars[i]+str(amount_varied)
                    if(amount_varied==110): # special case has no suffix
                        var = vars[i]
                    file_prefix = mode+nu+var+"_"+hierarchy+"_"+mixing_param
                
                    t_file = ROOT.TFile("/lbne/data/users/lblpwg_tools/FastMC_Data/outputs/ljf26/Sensitivity_Plots/"+file_prefix+"_histos.root");
                    t_file.Print()
                    if(plotvar=="cp_median"):
                        medians[0].append(math.copysign(math.sqrt(math.fabs(math.pow(GetMedian(t_file.Get("h1")),2)-math.pow(baseline_medians[0],2))),math.pow(GetMedian(t_file.Get("h1")),2)-math.pow(baseline_medians[0],2)))
                        medians[1].append(math.copysign(math.sqrt(math.fabs(math.pow(GetMedian(t_file.Get("h2")),2)-math.pow(baseline_medians[1],2))),math.pow(GetMedian(t_file.Get("h2")),2)-math.pow(baseline_medians[1],2)))
                        medians[2].append(math.copysign(math.sqrt(math.fabs(math.pow(GetMedian(t_file.Get("h3")),2)-math.pow(baseline_medians[2],2))),math.pow(GetMedian(t_file.Get("h3")),2)-math.pow(baseline_medians[2],2)))

                    if(plotvar=="cp_75thpercentile"):
                        medians[0].append(math.copysign(math.sqrt(math.fabs(math.pow(GetPercentile(t_file.Get("h1"),75),2)-math.pow(baseline_medians[0],2))),math.pow(GetPercentile(t_file.Get("h1"),75),2)-math.pow(baseline_medians[0],2)))
                        medians[1].append(math.copysign(math.sqrt(math.fabs(math.pow(GetPercentile(t_file.Get("h2"),75),2)-math.pow(baseline_medians[1],2))),math.pow(GetPercentile(t_file.Get("h2"),75),2)-math.pow(baseline_medians[1],2)))
                        medians[2].append(math.copysign(math.sqrt(math.fabs(math.pow(GetPercentile(t_file.Get("h3"),75),2)-math.pow(baseline_medians[2],2))),math.pow(GetPercentile(t_file.Get("h3"),75),2)-math.pow(baseline_medians[2],2)))

                    elif(plotvar=="mh_minimum"):
                        medians[0].append(math.copysign(math.sqrt(math.fabs(math.pow(GetMinimum(t_file.Get("h1")),2)-math.pow(baseline_medians[0],2))),math.pow(GetMinimum(t_file.Get("h1")),2)-math.pow(baseline_medians[0],2)))
                        medians[1].append(math.copysign(math.sqrt(math.fabs(math.pow(GetMinimum(t_file.Get("h2")),2)-math.pow(baseline_medians[1],2))),math.pow(GetMinimum(t_file.Get("h2")),2)-math.pow(baseline_medians[1],2)))
                        medians[2].append(math.copysign(math.sqrt(math.fabs(math.pow(GetMinimum(t_file.Get("h3")),2)-math.pow(baseline_medians[2],2))),math.pow(GetMinimum(t_file.Get("h3")),2)-math.pow(baseline_medians[2],2)))
                
                # if the file is empty, just set to -999 so it can be caught later
                for i in range(3):
                    for var in range(len(vars)):
                        if medians[i][var]==-1.0*baseline_medians[i]:
                            medians[i][var]=-999

                ROOT.gStyle.SetMarkerStyle(21);
                ROOT.gStyle.SetMarkerStyle(21);
                
                c1 = ROOT.TCanvas(file_prefix);
                
                t_graphs = []
                
                line_styles = [1,2,3]
                line_colors = [1,2,4]

                for i in 0,1,2:

                    temp = ROOT.TGraph(len(bin_centers_array),bin_centers_array,array.array('d',medians[i]));
                    temp.SetName(plotvar+"_"+mode+"_"+nu+"_"+hierarchy)
                    nu_string = "#nu_{#mu}"
                    if(nu == "numubar"):
                        nu_string = "#bar{#nu}_{#mu}"
                    temp.SetTitle("Effect of increasing "+nu_string+" "+mode+" flux by 10%")
                    t_graphs.append(temp)
                    
                    if(plotvar == "cp_median"):
                        #t_graphs[i].SetMinimum(4);
                        t_graphs[i].GetYaxis().SetTitle("Change in Median CP Sensitivity (#sigma)");
                    elif(plotvar == "cp_75thpercentile"):
                        #t_graphs[i].SetMinimum(4);
                        t_graphs[i].GetYaxis().SetTitle("75% CP Sensitivity (#sigma)");
                    elif(plotvar == "cp_3sigma"):
                        #t_graphs[i].SetMinimum(4);
                        t_graphs[i].GetYaxis().SetTitle("Fraction > 3#sigma Coverage (%)");
                    elif(plotvar == "cp_5sigma"):
                        #t_graphs[i].SetMinimum(4);
                        t_graphs[i].GetYaxis().SetTitle("Fraction > 5#sigma Coverage (%)");
                    elif(plotvar == "mh_minimum"):
                        #t_graphs[i].SetMinimum(2.5);
                        t_graphs[i].GetYaxis().SetTitle("Change in Minimum MH Sensitivity (#sigma)");
                        
                        
                    t_graphs[i].GetXaxis().SetTitle("Neutrino Energy (GeV)");


                    t_graphs[i].SetLineColor(line_colors[i])
                    t_graphs[i].SetMarkerColor(line_colors[i])
                    
                    c1.cd()
                    if i == 0:
                        t_graphs[i].Draw("ACP");
                    else:
                        t_graphs[i].Draw("SAME CP");


                leg = ROOT.TLegend(0.45,0.3,0.85,0.55);
                leg.SetHeader("Sig/Bkgd Uncertainties");
                leg.AddEntry(t_graphs[0],"1%/5%","lp");
                leg.AddEntry(t_graphs[1],"2%/5%","lp");
                leg.AddEntry(t_graphs[2],"5%/10%","lp");
                leg.SetFillStyle(0);
                leg.SetBorderSize(0);
                
                leg.Draw();
                file_prefix = "sensitivity_plots/"+mode+str(amount_varied)+nu+"_"+hierarchy+"_"+plotvar

                outfile.cd()
                t_graphs[1].Write();

                c1.Print(file_prefix+"_improved.eps")
                c1.Print(file_prefix+"_improved.png")
        

  outfile.Close()


# read in the above and plot change in sensitivity as a function of change in flux
#unicorns (for searching)

modes = ["fhc","rhc"]
nus = ["numu","numubar","nuenuebar"]
vars = ["0to0p5GeVPlus10per","0p5to1GeVPlus10per","1to2GeVPlus10per","2to3GeVPlus10per","3to4GeVPlus10per","4to5GeVPlus10per","5to6GeVPlus10per","6to7GeVPlus10per","7to8GeVPlus10per","8to9GeVPlus10per","9to10GeVPlus10per","10to15GeVPlus10per","15to20GeVPlus10per","20to120GeVPlus10per"]
bin_centers_array = array.array('d',bin_centers)
hierarchies = ["nh","ih"]
#plotvars = ["cp_median","cp_75thpercentile","cp_3sigma","cp_5sigma","mh_minimum"]
plotvars = ["cp_75thpercentile","mh_minimum"]

amounts_varied = [0, 50,90,100,110,200,500,1000]

for mode in modes:
    for nu in nus:
        for hierarchy in hierarchies:
            for plotvar in plotvars:
                file_prefix = "sensitivity_plots/"+mode+nu+"_"+hierarchy+"_"+plotvar
                c1 = ROOT.TCanvas(file_prefix,file_prefix,1800,1000)
                c1.Divide(5,3)

                # declare and initialize list of sensitivity changes
                sensitivity_changes = []
                flux_changes = []
                for bin_iter in range(0,len(vars)):
                    sensitivity_changes.append([])
                    flux_changes.append([])

                for variation_iter in range(0,len(amounts_varied)):
                    variation = amounts_varied[variation_iter]
                    if(variation==100):
                        for bin_iter in range(len(vars)):
                            sensitivity_changes[bin_iter].append(0)
                            flux_changes[bin_iter].append(100)
                    else:
                        temp_file = ROOT.TFile("/lbne/data/users/ljf26/fluxfiles/sensitivity_comps_"+str(variation)+"_improved.root");                    
                        mygraph = temp_file.Get(plotvar+"_"+mode+"_"+nu+"_"+hierarchy)
                        for bin_iter in range(mygraph.GetN()):
                            if(mygraph.GetY()[bin_iter] != -999):
                                sensitivity_changes[bin_iter].append(mygraph.GetY()[bin_iter])
                                flux_changes[bin_iter].append(amounts_varied[variation_iter])

                temps = []
                lines = []
                tempfuncs = []
                for bin_iter in range(0,len(vars)):
                    c1.cd(bin_iter+1)
                    temps.append(ROOT.TGraph(len(flux_changes[bin_iter]),array.array('d',flux_changes[bin_iter]),array.array('d',sensitivity_changes[bin_iter])));
                    temps[bin_iter].Draw("AP")
                    xlow = -20
                    xhigh = 1000
                    temps[bin_iter].GetXaxis().SetLimits(xlow,xhigh);
                    if(use_fixed_scale):
                        maximum = 0.6*5
                        minimum = -0.6*5
                    else:
                        maximum = temps[bin_iter].GetHistogram().GetMaximum()
                        minimum = temps[bin_iter].GetHistogram().GetMinimum()
                    myrange = ROOT.TMath.Max(ROOT.TMath.Abs(maximum),ROOT.TMath.Abs(minimum));
                    temps[bin_iter].GetHistogram().SetMinimum(-1.0*myrange)
                    temps[bin_iter].GetHistogram().SetMaximum(myrange)
                    temps[bin_iter].SetMarkerStyle(21);
                    temps[bin_iter].SetMarkerSize(2);
                    temps[bin_iter].SetMarkerColor(4);
                    temps[bin_iter].GetXaxis().SetNdivisions(4);
                    temps[bin_iter].GetYaxis().SetNdivisions(4);
                    temps[bin_iter].Draw("AP")
                    c1.Update()
                    temps[bin_iter].GetXaxis().SetLabelSize(0.08)
                    ROOT.gStyle.SetOptTitle(1);
                    temps[bin_iter].GetXaxis().SetTitle("Change in Flux (%)")
                    temps[bin_iter].GetYaxis().SetTitle("Change in CP Sens. (#sigma)")
                    temps[bin_iter].GetXaxis().SetTitleSize(0.07)
                    temps[bin_iter].GetYaxis().SetTitleSize(0.07)
                    temps[bin_iter].GetXaxis().SetTitleOffset(0.9)
                    temps[bin_iter].GetYaxis().SetTitleOffset(1.05)
                    temps[bin_iter].GetYaxis().SetLabelSize(0.07)
                    temps[bin_iter].SetTitle(vars[bin_iter].split("Plus")[0])
                    ROOT.gPad.Update()
                    title = ROOT.gPad.GetPrimitive("title");
                    title.SetBorderSize(0);
                    ROOT.gStyle.SetTitleY(0.95)
                    title.SetTextSize(0.08)
                    title.SetTextColor(8)
                    #temps[bin_iter].SetTitle(vars[bin_iter])
                    #add_plot_label(Form("%1.1f < E < %1.1f GeV",h_nof[0]->GetBinLowEdge(k+1),h_nof[0]->GetBinLowEdge(k+1)+h_nof[0]->GetBinWidth(k+1)),0.55,0.96,0.09,38,62,22);

                    # draw a line through the points at 100,110% 
                    # this was the first estimator I used -- want to see how much the other points deviate from it
                    iter_100 = -1
                    iter_110 = -1
                    for i in range(len(flux_changes[bin_iter])):
                        if flux_changes[bin_iter][i]  == 100:
                            iter_100 = i
                        if flux_changes[bin_iter][i] == 110:
                            iter_110 = i

                    if(iter_100 == -1 or iter_110==-1):
                        continue # Don't draw a line
                    slope = (sensitivity_changes[bin_iter][iter_110]-sensitivity_changes[bin_iter][iter_100])/(flux_changes[bin_iter][iter_110]-flux_changes[bin_iter][iter_100]);
                    y_int = sensitivity_changes[bin_iter][iter_110]-slope*flux_changes[bin_iter][iter_110];
                    ylow_line = slope*xlow+y_int
                    yhigh_line = slope*xhigh+y_int
                    xlow_line = xlow
                    xhigh_line = xhigh
                    if ylow_line < -1.0*myrange:
                        ylow_line = -1.0*myrange
                        xlow_line = (ylow_line-y_int)/slope
                    if yhigh_line > 1.0*myrange:
                        yhigh_line = 1.0*myrange
                        xhigh_line = (yhigh_line-y_int)/slope
                    lines.append(ROOT.TLine(xlow_line,ylow_line,xhigh_line,yhigh_line))
                    lines[bin_iter].SetLineColor(2)
                    lines[bin_iter].Draw("SAME")
                                 

                    # Now draw a linear interpolation
                    tempgraphs[bin_iter] = temps[bin_iter]
                    tempfuncs.append(ROOT.TF1("f"+str(bin_iter),myfunc,xlow,xhigh,1))
                    tempfuncs[bin_iter].SetParameter(0,bin_iter)
                    tempfuncs[bin_iter].SetLineColor(38)
                    tempfuncs[bin_iter].Draw("SAME")

                scale_string = "fixedScale"
                if not use_fixed_scale:
                    scale_string = "variableScale"
                c1.Print(file_prefix+"_"+scale_string+"_improved.eps")
                c1.Print(file_prefix+"_"+scale_string+"_improved.png")
                outfile = ROOT.TFile(file_prefix+"_"+scale_string+"_improved.root","RECREATE")
                for graph in temps:
                    graph.Write()
                outfile.Close()