plotSensitivity.py

Go to the documentation of this file.

import ROOT, array, os

# Set information related to the beams you want to compare
x_title = "Beam Options"
macros = ["Nominal","CP_run5_9116_80GeV","CP_run15_12388","CP_run17_6432","CP_run18_3849"]
xs = [1,2,3,4,5]
energies = [120.0,80.0,62.4,111.4,108.8]
versions = ["v3r4p2","v3r4p2","v3r4p2","v3r4p2","v3r4p2"]
users = ["ljf26","ljf26","ljf26","ljf26","ljf26"]
locs = ["data","data","data","data","data"]
sens_min = 1.4
sens_max = 2.4
flux_min = 0e-12
flux_max = 150e-12
text_axis = True # puts words rather than numbers on the x axis
axis_labels = ["Nominal","2Horn","3HornNuMI","3HornSphere","3HornCyl"]

"""
x_title = "Horn A Downstream Endcap Thickness (mm)"
macros = ["CP_run15_12388_NoWater","CP_run15_12388_Horn1DE4mm_NoWater","CP_run15_12388_Horn1DE6mm_NoWater","CP_run15_12388_Horn1DE8mm_NoWater","CP_run15_12388_Horn1DE10mm_NoWater","CP_run15_12388_Horn1DE15mm_NoWater","CP_run15_12388_Horn1DE20mm_NoWater"]
xs = [2.0,4.0,6.0,8.0,10.0,15.0,20.0]
locs = ["data","data","data","data","data","data","data"]
energies = [62.4,62.4,62.4,62.4,62.4,62.4,62.4]
versions = ["v3r4p2","v3r4p2","v3r4p2","v3r4p2","v3r4p2","v3r4p2","v3r4p2"]
sens_min = 1.8
sens_max = 2.2
flux_min = 50e-12
flux_max = 100e-12
text_axis = False
"""

# Helper function for calculating 75th percentile of a set of numbers
def GetPercentile(y,X):
    n = len(y);

    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]


cp_75th_percentiles = []
sens_v_dcp = []
dcps = []
max = 0

# Loop over the beams to plot
for i in range(len(xs)):
    version = versions[i]
    user = users[i]
    macro = macros[i]

    # Look up the flux files
    fhc_file = "/dune/"+locs[i]+"/users/"+user+"/fluxfiles/g4lbne/"+version+"/QGSP_BERT/"+macros[i]+"/neutrino/flux/histos_g4lbne_"+version+"_QGSP_BERT_"+macros[i]+"_neutrino_LBNEFD_fastmc.root"
    rhc_file = "/dune/"+locs[i]+"/users/"+user+"/fluxfiles/g4lbne/"+version+"/QGSP_BERT/"+macros[i]+"/antineutrino/flux/histos_g4lbne_"+version+"_QGSP_BERT_"+macros[i]+"_antineutrino_LBNEFD_fastmc.root"
    fhc_tfile = ROOT.TFile(fhc_file)
    rhc_tfile = ROOT.TFile(rhc_file)


    # Get sensitivity
    sensitivity_file = open("/dune/data/users/"+os.getenv("USER")+"/sensitivities/"+macro+"_cpsens.dat")
    lines = sensitivity_file.readlines()
    sensitivity_file.close()
    cp_sens = []
    for line in lines:
        splitline = line.split()
        deltacp = float(splitline[0])
        delta_chi2 = float(splitline[1])
        if delta_chi2 <0: delta_chi2 = 0
        cp_sens.append(ROOT.TMath.Sqrt(delta_chi2))
        if cp_sens[-1]>max:
            max = cp_sens[-1]
        dcps.append(deltacp/ROOT.TMath.Pi())

    cp_75th_percentiles.append(GetPercentile(cp_sens,75))
    sens_v_dcp.append(cp_sens)

#    Plot 75% CP sensitivity
if(text_axis):
    ROOT.gStyle.SetPadBottomMargin(0.25);
c1 = ROOT.TCanvas()
ROOT.gPad.SetGrid(1)
g1 = ROOT.TGraph(len(xs),array.array('d',xs),array.array('d',cp_75th_percentiles))
g1.GetXaxis().SetTitle(x_title)
g1.GetYaxis().SetTitle("GLoBES 75% CP Sensitivity (NH)")
if(text_axis):
    g1.GetXaxis().SetTitle("")
    for i in range(len(xs)):
        g1.GetXaxis().SetBinLabel(g1.GetXaxis().FindBin(i+1),axis_labels[i])

g1.SetTitle("")
g1.SetLineStyle(2)
g1.SetMinimum(sens_min)
g1.SetMaximum(sens_max)
g1.SetMarkerStyle(21)
g1.Draw("ACP")

c1.Print(x_title.replace(" ","")+"_Sensitivity.eps")
c1.Print(x_title.replace(" ","")+"_Sensitivity.png")

# Plot CP sensitivity vs delta_CP w. different beams overlaid

ROOT.gStyle.SetTitleY(0.978);
ROOT.gStyle.SetOptStat(False);

height = 1.5*max
n = len(sens_v_dcp[0])

c = ROOT.TCanvas("c","c",500,500);
 
frame = ROOT.TH2D("frame","frame",100,-1,1,100,0,height);
frame.SetTitle("CP violation sensitivity");
frame.GetXaxis().SetTitle("#delta_{cp}/#pi");
frame.GetXaxis().SetTitleSize(0.05);
frame.GetXaxis().CenterTitle()
frame.GetXaxis().SetTitleOffset(1.1);
frame.GetXaxis().SetLabelSize(0.035);

frame.GetYaxis().SetTitle("#sigma = #sqrt{#Delta#chi^{2}}");
frame.GetYaxis().SetTitleSize(0.05);
frame.GetYaxis().CenterTitle()
frame.GetYaxis().SetTitleOffset(1.1);
frame.GetYaxis().SetLabelSize(0.035);


frame.Draw();
   
graphs = []
colors = [2,4,6,8,ROOT.kOrange]
legend = ROOT.TLegend(0.2,0.7,0.8,0.9);
legend.SetFillStyle(0);
legend.SetBorderSize(0);
for k in range(len(macros)):

    temp_graph = ROOT.TGraph(n,array.array('d',dcps),array.array('d',sens_v_dcp[k]))
    temp_graph.SetLineWidth(3);
    temp_graph.SetLineColor(colors[k])

    temp_graph.Draw("lsame")

    graphs.append(temp_graph)

    legend.AddEntry(temp_graph,axis_labels[k],"l");


legend.Draw();

c.SetTopMargin(0.9);
c.SetBottomMargin(0.15);
c.SetLeftMargin(0.15);
c.Update();

c.Print("CP_sens.eps");
c.Print("CP_sens.png");