optimization_test.py

Go to the documentation of this file.

import ROOT, os, array, sys, OptimizationUtils,Optimizations

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 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"]

print "Blah"
optimization = Optimizations.Optimization("CP_run3")
print "Blah2"

for hierarchy in hierarchies:
    for plotvar in plotvars:
        amounts_varied = [177,3544,4177,6177,6911,7177,7862,8177,8186]

        # 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_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=="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")))

        medians = [[],[],[]]

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

                # get varied medians
                var = "Optimizations_CP_run3_"+str(amount_varied)
                file_prefix = var+"_"+hierarchy+"_"+mixing_param
                
                t_file = ROOT.TFile("/lbne/data/users/lblpwg_tools/FastMC_Data/outputs/ljf26/Sensitivity_Plots/"+file_prefix+"_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))
                    medians[1].append(GetPercentile(t_file.Get("h2"),75))
                    print amount_varied,GetPercentile(t_file.Get("h2"),75),t_file.Get("h2").Integral()
                    medians[2].append(GetPercentile(t_file.Get("h3"),75))

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

        # 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"
            for var in range(len(amounts_varied_good_files)):
                print var
                varfhcfile = "/lbne/data/users/ljf26/fluxfiles/g4lbne/v3r3p2/QGSP_BERT/Optimizations-CP_run3-"+str(amounts_varied_good_files[var])+"/200kA/flux/histos_g4lbne_"+optimization.g4lbne_version+"_QGSP_BERT_Optimizations-CP_run3-"+str(amounts_varied_good_files[var])+"_200kA_LBNEFD_fastmc.root"
                varrhcfile = "/lbne/data/users/ljf26/fluxfiles/g4lbne/v3r3p2/QGSP_BERT/Optimizations-CP_run3-"+str(amounts_varied_good_files[var])+"/-200kA/flux/histos_g4lbne_"+optimization.g4lbne_version+"_QGSP_BERT_Optimizations-CP_run3-"+str(amounts_varied_good_files[var])+"_-200kA_LBNEFD_fastmc.root"
        
                energy = optimization.getEnergy(float(amounts_varied_good_files[var]),"")
                antinufrac = 0.5

                type = "cp"
                if plotvar.startswith("mh"): 
                    type = "mh"

                [fitness_names,fitnesses,fiterrors] = OptimizationUtils.ComputeMean75PercentSensitivityAndError(varfhcfile,varrhcfile,energy,energy,antinufrac,type)

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

                    
        # 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)
            
            t_graphs.append(temp)
                              
            if(plotvar == "cp_75thpercentile"):
                t_graphs[i].GetYaxis().SetTitle("75% CP Sensitivity (#sigma)");
            elif(plotvar == "mh_minimum"):
                t_graphs[i].GetYaxis().SetTitle("Change in Minimum MH Sensitivity (#sigma)");
   
            t_graphs[i].GetXaxis().SetTitle("Optimization Configuration");

            t_graphs[i].SetLineColor(line_colors[i])
            t_graphs[i].SetMarkerColor(line_colors[i])
                    
            c1.cd()

            if(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(3.5);
                t_graphs[i].SetMinimum(0);
                t_graphs[i].Draw("ACP")
                m_graph.Draw("SAME CP")

                leg = ROOT.TLegend(0.4,0.2,0.6,0.4);
                #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/optimization_test_"+hierarchy+"_"+plotvar

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