energyTest.py

Go to the documentation of this file.

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

optimization_name = "power_test"

##### Find processed configurations and calculate fitnesses #####

iterations = []
fitnesses = []
fitness_errors = []
proton_energies = []
zero_errors = []

max_iterations = 1e7

optimization = Optimizations.Optimization(optimization_name)

fluxpath = "/lbne/data/users/ljf26/fluxfiles/g4lbne/"+optimization.g4lbne_version+"/QGSP_BERT/"
fluxdirs = glob.glob(fluxpath+"*")



for dir in fluxdirs:

    if os.path.basename(dir).startswith("Optimizations-"+optimization_name):
        iteration = dir.split("/")[9].split("-")[2]

        varfhcfile = dir+"/200kA/flux/histos_g4lbne_"+optimization.g4lbne_version+"_QGSP_BERT_"+os.path.basename(dir)+"_200kA_LBNEFD_fastmc.root"
        varrhcfile = dir+"/-200kA/flux/histos_g4lbne_"+optimization.g4lbne_version+"_QGSP_BERT_"+os.path.basename(dir)+"_-200kA_LBNEFD_fastmc.root"
        
        #[tfitness_names,tfitnesses,tfiterrors] = optimization.getFitness(float(iteration))
        energy = optimization.getEnergy(float(iteration),"FHC")
        antinufrac = 0.5
        print "energy",energy
        [tfitness_names,tfitnesses,tfiterrors] = OptimizationUtils.ComputeMean75PercentSensitivityAndError(varfhcfile,varrhcfile,energy,energy,antinufrac,"cp")
        fitness = tfitnesses[0]
        print "fitness",fitness
        fitness_error = tfiterrors[0]

        if(fitness>0):
            iterations.append(float(iteration))
            fitnesses.append(fitness)
            fitness_errors.append(fitness_error)
            zero_errors.append(0.0)
    
# Read in knob turns from macros
iterations2 = []
parameter_values = []

parameter_histos = []

n_parameters = len(optimization.parameter_names);
for i in range(0,n_parameters):
    parameter_name = optimization.parameter_names[i]
    parameter_values.append([])
    parameter_histos.append(ROOT.TH1F("h_parameters_"+str(i),parameter_name,10,optimization.parameter_lower_limits[i],optimization.parameter_upper_limits[i]));

knobturns = optimization.getKnobturns()

for iteration in iterations:
    if float(iteration)>max_iterations: continue
    iteration_found = False
    for j in range(0,len(knobturns)):
        if knobturns[j][0]==iteration:
            iteration_found = True
            iterations2.append(knobturns[j][0])
            for i in range(0,n_parameters):
                parameter_histos[i].Fill(knobturns[j][1][i])
                parameter_values[i].append(knobturns[j][1][i])
    if not iteration_found:
        print "ERROR: iteration:",iteration," found in flux directories but can't find macros to look up knobturns... quitting."
        sys.exit()

# Make TGraphs
canv = ROOT.TCanvas("MyCanvas","MyCanvas")

ROOT.gStyle.SetMarkerStyle(2);

for i in range(0,n_parameters):
    
    fitness_vs_parameter = ROOT.TGraphErrors(len(iterations),array('d',parameter_values[i]),array('d',fitnesses),array('d',zero_errors),array('d',fitness_errors))
    fitness_vs_parameter.GetYaxis().SetTitle("Fitness")
    x_axis_title = optimization.parameter_names[i]
    if(optimization.parameter_units[i]!=""):
        x_axis_title = x_axis_title+" ("+optimization.parameter_units[i]+")"
    fitness_vs_parameter.GetXaxis().SetTitle(x_axis_title);
    fitness_vs_parameter.GetXaxis().CenterTitle()
    fitness_vs_parameter.GetYaxis().CenterTitle()
    fitness_vs_parameter.SetTitle("")
    fitness_vs_parameter.Draw("AP");

    # Cross check from FMC
    FMC_SENSIT = "/lbne/data/users/lblpwg_tools/FastMC_Data/outputs/ljf26/Sensitivity_Plots"

    energies = [20,30,50,60,70,80,90,110,120,130]
    mean_sensitivities = []
    for energy in energies:

        if energy != 120:
            nh_cp_file = ROOT.TFile(FMC_SENSIT+"/ProtonP"+str(int(energy))+"GeV_RealisticPower_nh_cp_histos.root") 
            ih_cp_file = ROOT.TFile(FMC_SENSIT+"/ProtonP"+str(int(energy))+"GeV_RealisticPower_ih_cp_histos.root") 
        else:
            nh_cp_file = ROOT.TFile(FMC_SENSIT+"/ProtonP"+str(int(energy))+"GeV_nh_cp_histos.root") 
            ih_cp_file = ROOT.TFile(FMC_SENSIT+"/ProtonP"+str(int(energy))+"GeV_ih_cp_histos.root") 

        nh_sensitivity = OptimizationUtils.GetPercentile(nh_cp_file.Get("h2"),75)
        ih_sensitivity = OptimizationUtils.GetPercentile(ih_cp_file.Get("h2"),75)
        
        mean_sensitivities.append((nh_sensitivity+ih_sensitivity)/2)

    blah = ROOT.TGraph(len(energies),array('d',energies),array('d',mean_sensitivities));
    blah.SetMarkerColor(2);
    blah.SetLineColor(2);
    blah.SetMarkerStyle(21);
    blah.GetXaxis().SetTitle("Proton Energy (GeV)")
    blah.GetYaxis().SetTitle("Average 75% #delta_{CP} Coverage (#sigma)")
    blah.SetTitle("")
    blah.Draw("ACP");
    fitness_vs_parameter.Draw("SAMEP");

    canv.Print("plots/"+optimization_name+"_fitness_vs_"+parameter_name+".eps")
    canv.Print("plots/"+optimization_name+"_fitness_vs_"+parameter_name+".png")