computeNormScales.py

Go to the documentation of this file.

import ROOT, os, array, sys, OptimizationUtils

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

hierarchies = ["nh","ih"]
plotvars = ["cp_75thpercentile","mh_minimum"]
#tweaks = ["fhc_numu"]
tweaks = ["fhc_numu","fhc_numubar","fhc_nue","rhc_numu","rhc_numu","rhc_numubar","rhc_nuebar"]
#tweaks = ["antinufrac"]

for tweak in tweaks:
 for hierarchy in hierarchies:
    for plotvar in plotvars:
        if tweak=="antinufrac":
            amounts_varied = [0,10,25,50,75,90,100]
        else:
            amounts_varied = [0,50,90,100,110,200]

        # get baseline median values
        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_3sigma"):
            baseline_medians.append(OptimizationUtils.GetFractionAtMoreThanXSigma(baseline_file.Get("h1"),3.0))
            baseline_medians.append(OptimizationUtils.GetFractionAtMoreThanXSigma(baseline_file.Get("h2"),3.0))
            baseline_medians.append(OptimizationUtils.GetFractionAtMoreThanXSigma(baseline_file.Get("h3"),3.0))
        if(plotvar=="cp_75thpercentile"):
            baseline_medians.append(GetPercentile(baseline_file.Get("h1"),75))
            baseline_medians.append(GetPercentile(baseline_file.Get("h2"),75))
            print "Baseline",GetPercentile(baseline_file.Get("h2"),75)
            baseline_medians.append(GetPercentile(baseline_file.Get("h3"),75))
        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")))
        elif(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")))
 
        medians = [[],[],[]]

        histos = [[],[],[]]
        print amounts_varied
        amounts_varied_good_files = []
        for amount_varied in amounts_varied:

                # get varied medians
                var = tweak+str(amount_varied)
                file_prefix = var+"_"+hierarchy+"_"+mixing_param
                
                if (tweak == "antinufrac" or amount_varied != 100):
                    t_file = ROOT.TFile("/lbne/data/users/lblpwg_tools/FastMC_Data/outputs/ljf26/Sensitivity_Plots/"+file_prefix+"_histos.root");
                else:
                    t_file = ROOT.TFile("/lbne/data/users/lblpwg_tools/FastMC_Data/outputs/ljf26/Sensitivity_Plots/ProtonP120GeV_"+hierarchy+"_"+mixing_param+"_histos.root");            
                
                ROOT.gROOT.cd()
                if(t_file.IsZombie() or t_file.Get("h2").Integral()==0):
                    continue
                amounts_varied_good_files.append(amount_varied)

                histos[0].append(t_file.Get("h1").Clone("h1_"+str(amount_varied)))
                histos[1].append(t_file.Get("h2").Clone("h2_"+str(amount_varied)))
                histos[2].append(t_file.Get("h3").Clone("h3_"+str(amount_varied)))
                
                if(plotvar=="cp_75thpercentile"):
                    medians[0].append(GetPercentile(t_file.Get("h1"),75)-baseline_medians[0])
                    medians[1].append(GetPercentile(t_file.Get("h2"),75)-baseline_medians[1])
                    medians[2].append(GetPercentile(t_file.Get("h3"),75)-baseline_medians[2])

                elif(plotvar=="mh_minimum"):
                    medians[0].append(GetMinimum(t_file.Get("h1"))-baseline_medians[0])
                    medians[1].append(GetMinimum(t_file.Get("h2"))-baseline_medians[1])
                    medians[2].append(GetMinimum(t_file.Get("h3"))-baseline_medians[2])

                elif(plotvar=="cp_median"):
                    medians[0].append(GetMedian(t_file.Get("h1"))-baseline_medians[0])
                    medians[1].append(GetMedian(t_file.Get("h2"))-baseline_medians[1])
                    medians[2].append(GetMedian(t_file.Get("h3"))-baseline_medians[2])

                elif(plotvar=="cp_3sigma"):
                    medians[0].append(OptimizationUtils.GetFractionAtMoreThanXSigma(t_file.Get("h1"),3.0)-baseline_medians[0])
                    medians[1].append(OptimizationUtils.GetFractionAtMoreThanXSigma(t_file.Get("h2"),3.0)-baseline_medians[1])
                    medians[2].append(OptimizationUtils.GetFractionAtMoreThanXSigma(t_file.Get("h3"),3.0)-baseline_medians[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(amounts_varied_good_files)):
        #        if medians[i][var]==0:
        #            medians[i][var]=-999

 
        if(plotvar == "cp_75thpercentile"):
            metrics = [[],[],[]]
                  
            # Get fitness metric, for comparison
            nomfhcfile = "/lbne/data/users/ljf26/fluxfiles/g4lbne/v3r2p4/QGSP_BERT/Nominal/200kA/flux/histos_g4lbne_v3r2p4_QGSP_BERT_Nominal_200kA_LBNEFD_fastmc.root"        
            nomrhcfile = "/lbne/data/users/ljf26/fluxfiles/g4lbne/v3r2p4/QGSP_BERT/Nominal/-200kA/flux/histos_g4lbne_v3r2p4_QGSP_BERT_Nominal_-200kA_LBNEFD_fastmc.root"
            if(tweak=="antinufrac"):
                varfhcfile = "/lbne/data/users/ljf26/fluxfiles/g4lbne/v3r2p4/QGSP_BERT/Nominal/200kA/flux/histos_g4lbne_v3r2p4_QGSP_BERT_Nominal_200kA_LBNEFD_fastmc.root"        
                varrhcfile = "/lbne/data/users/ljf26/fluxfiles/g4lbne/v3r2p4/QGSP_BERT/Nominal/-200kA/flux/histos_g4lbne_v3r2p4_QGSP_BERT_Nominal_-200kA_LBNEFD_fastmc.root"

            power = 120
            antinufrac = 0.5
            for var in range(len(amounts_varied_good_files)):

                if tweak != "antinufrac":
                    varfhcfile = "temp/fhc_"+tweak+str(amounts_varied_good_files[var])+".root"
                    varrhcfile = "temp/rhc_"+tweak+str(amounts_varied_good_files[var])+".root"

                else:
                    antinufrac = amounts_varied_good_files[var]/100.0
                
                [fitness_names,fitnesses,fiterrors] = OptimizationUtils.ComputeMean75PercentSensitivityAndError(varfhcfile,varrhcfile,power,antinufrac,"cp")

                print "metric",fitnesses[0],fitnesses[7],fitnesses[8]
                if(hierarchy=="nh"):
                    metrics[1].append(fitnesses[7]-baseline_medians[1])
                if(hierarchy=="ih"):
                    metrics[1].append(fitnesses[8]-baseline_medians[1])

                    
        # Make plots
        ROOT.gStyle.SetMarkerStyle(21);
        ROOT.gStyle.SetMarkerStyle(21);
                
        c1 = ROOT.TCanvas(file_prefix);
                
        t_graphs = []
        mgraph = 0
        
        line_styles = [1,2,3]
        line_colors = [1,2,4]
    
        for i in 0,1,2:

            temp = ROOT.TGraph(len(amounts_varied_good_files),array.array('d',amounts_varied_good_files),array.array('d',medians[i]));
            temp.SetName(plotvar+"_"+hierarchy)
#            if tweak == "antinufrac":
#                temp.SetTitle("CP sensitivity versus fraction of antineutrino running")
#            else:
#                temp.SetTitle("CP sensitivity versus change in "+tweak+" flux")
            t_graphs.append(temp)
                              
            if(plotvar == "cp_75thpercentile"):
                        #t_graphs[i].SetMinimum(4);
                t_graphs[i].GetYaxis().SetTitle("75% CP Sensitivity (#sigma)");
            elif(plotvar == "mh_minimum"):
                        #t_graphs[i].SetMinimum(2.5);
                t_graphs[i].GetYaxis().SetTitle("Change in Minimum MH Sensitivity (#sigma)");
            elif(plotvar == "cp_median"):
                        #t_graphs[i].SetMinimum(2.5);
                t_graphs[i].GetYaxis().SetTitle("Median CP Sensitivity (#sigma)");
   
            if tweak=="antinufrac":
                t_graphs[i].GetXaxis().SetTitle("Antineutrino Running Fraction of Total");
            else:
                t_graphs[i].GetXaxis().SetTitle("Normalization of "+tweak+" flux to nominal configuration (%)");    

            t_graphs[i].SetLineColor(line_colors[i])
            t_graphs[i].SetMarkerColor(line_colors[i])
                    
            c1.cd()
            if(plotvar!="cp_75thpercentile"):
                t_graphs[i].SetMinimum(-1);
                t_graphs[i].SetMaximum(1);
                if i == 0:
                    t_graphs[i].Draw("ACP");
                else:
                    t_graphs[i].Draw("SAME CP");
                if i ==2:
                    leg = ROOT.TLegend(0.35,0.2,0.65,0.5);
                    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();


            if(plotvar == "cp_75thpercentile" and i==1):
                m_graph = ROOT.TGraph(len(amounts_varied_good_files),array.array('d',amounts_varied_good_files),array.array('d',metrics[i]));
                m_graph.SetLineColor(line_colors[i])
                m_graph.SetMarkerColor(line_colors[i])
                m_graph.SetLineStyle(2)
                t_graphs[i].SetMaximum(1);
                t_graphs[i].SetMinimum(-1);
                t_graphs[i].Draw("ACP")
                m_graph.Draw("SAME CP")

                leg = ROOT.TLegend(0.4,0.7,0.6,0.9);
                #leg.SetHeader("Sig/Bkgd Uncertainties");
                leg.AddEntry(t_graphs[1],"Fast MC","l");
                leg.AddEntry(m_graph,"Metric","l");
                leg.SetFillStyle(0);
                leg.SetBorderSize(0);
                leg.Draw();
        
        
        file_prefix = "sensitivity_plots/normscales_"+tweak+"_"+hierarchy+"_"+plotvar

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

        
        ofile = ROOT.TFile(file_prefix+".root","RECREATE")
        ofile.cd()
        m_graph.Write()
        t_graphs[1].Write()
        ofile.Close()

        """
        c1.Clear()
        colors = [2,4,6,1,3,5,7]
        for i in range(len(histos[1])):
            histos[1][i].SetLineColor(colors[i])
            histos[1][i].SetMaximum(10)
            if(i==1):
                histos[1][i].Draw()
            else:
                histos[1][i].Draw("same")

        leg = ROOT.TLegend(0.45,0.5,0.85,0.85);
        if(tweak=="antinufrac"):
            leg.AddEntry(histos[1][0],"10/0");
            leg.AddEntry(histos[1][1],"9/1");
            leg.AddEntry(histos[1][2],"7.5/2.5");
            leg.AddEntry(histos[1][3],"5/5");
            leg.AddEntry(histos[1][4],"2.5/7.5");
            leg.AddEntry(histos[1][5],"1/9");
            leg.AddEntry(histos[1][6],"0/10");
        else:
            for i in range(len(histos[1])):
                leg.AddEntry(histos[1][i],str(amounts_varied[i]))
 

        leg.SetFillStyle(0);
        leg.SetBorderSize(0);
        leg.Draw();

        c1.Print(plotvar+"_"+hierarchy+"blah.eps")
        """