bestFitnesses.py
Go to the documentation of this file.
1 import ROOT, glob, sys, os
2 from array import array
3 import Optimizations, OptimizationUtils
4 
5 ##### PARSE INPUT ARGUMENT #####
6 
7 if not len(sys.argv) > 1:
8  print "Please specify an optimization name"
9  sys.exit()
10 
11 optimization_name = sys.argv[1]
12 
13 max_configurations = 3200
14 
15 if len(sys.argv)>2:
16  max_configurations = int(sys.argv[2])
17 
18 ##### Find processed configurations and calculate fitnesses #####
19 
20 configurations = []
21 fitnesses = []
22 fitness_errors = []
23 proton_energies = []
24 zero_errors = []
25 
26 n_var_plots = 1;
27 
28 optimization = Optimizations.Optimization(optimization_name)
29 
30 fluxpath = optimization.output_location+"users/ljf26/fluxfiles/g4lbne/"+optimization.g4lbne_version+"/"+optimization.physics_list+"/"
31 fluxdirs = glob.glob(fluxpath+"*")
32 if optimization.rhc_parameters_float_separate:
33  fluxdirs=glob.glob(fluxpath+"*FHC*")
34 
35 for dir in fluxdirs:
36 
37  if "Scan" in dir:
38  continue
39 
40  if os.path.basename(dir).startswith("Optimizations-"+optimization_name):
41  configuration = dir.split("/")[9].split("-")[2]
42 
43  if float(configuration) >max_configurations: continue
44  [fitness_types,tfitnesses,fiterrors] = optimization.getFitness(float(configuration))
45  fitness = tfitnesses[0]
46  fitness_error = fiterrors[0]
47  print "AAAA",configuration, fitness, fitness_error
48 
49  if(fitness>0):
50  configurations.append(float(configuration))
51  fitnesses.append(fitness)
52  fitness_errors.append(fitness_error)
53  zero_errors.append(0.0)
54 
55 # Git list of list indices ordered by fitness
56 ordered_indices = [i[0] for i in sorted(enumerate(fitnesses), key=lambda x:x[1],reverse=True)]
57 
58 knobturns = optimization.getKnobturns()
59 
60 for i in range(0,4):
61  config_number = configurations[ordered_indices[i]]
62 
63  print "config_number",config_number
64 
65  knob_index = -1
66  for j in range(0,len(knobturns)):
67  if float(knobturns[j][0])==float(config_number):
68  knob_index = j
69  break
70  if knob_index == -1:
71  print "ERROR: Could not find knobs for configuration",config_number
72  knob_names = optimization.parameter_names
73 
74  print "\n\nFitness rank:",i+1
75  print "Configuration number:",int(config_number)
76  all_fitnesses = optimization.getFitness(int(config_number))
77  print "Fitnesses:"
78  print " ",all_fitnesses
79  print "MH Fitnesses"
80  print " ",optimization.getFitness(int(config_number),"mh")
81  #print fitnesses[ordered_indices[i]]
82  for j in range(0,len(knobturns[knob_index][1])):
83  print " ",optimization.parameter_names[j],knobturns[knob_index][1][j]
84 
85 # Get Nominal Flux
86 nom_flux_file_fhc = ROOT.TFile("/lbne/data/users/ljf26/fluxfiles/g4lbne/v3r2p4/QGSP_BERT/Nominal/200kA/flux/histos_g4lbne_v3r2p4_QGSP_BERT_Nominal_200kA_LBNEFD_fastmc.root")
87 nom_flux_file_rhc = ROOT.TFile("/lbne/data/users/ljf26/fluxfiles/g4lbne/v3r2p4/QGSP_BERT/Nominal/-200kA/flux/histos_g4lbne_v3r2p4_QGSP_BERT_Nominal_-200kA_LBNEFD_fastmc.root")
88 
89 nom_flux = []
90 nom_flux.append(nom_flux_file_fhc.Get("numu_flux"))
91 nom_flux.append(nom_flux_file_fhc.Get("numubar_flux"))
92 nom_flux.append(nom_flux_file_fhc.Get("nue_flux"))
93 nom_flux.append(nom_flux_file_fhc.Get("nuebar_flux"))
94 nom_flux.append(nom_flux_file_rhc.Get("numu_flux"))
95 nom_flux.append(nom_flux_file_rhc.Get("numubar_flux"))
96 nom_flux.append(nom_flux_file_rhc.Get("nue_flux"))
97 nom_flux.append(nom_flux_file_rhc.Get("nuebar_flux"))
98 
99 # Get Optimized Fluxes
100 n_fluxes_to_plot = 4
101 var_flux_files_fhc = []
102 var_flux_files_rhc = []
103 var_fluxes = []
104 new_bins = [0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8]
105 for i in range(0,n_fluxes_to_plot):
106 
107  configuration = configurations[ordered_indices[i]]
108  print "configuration",configuration
109  if optimization.rhc_parameters_float_separate:
110  var_flux_files_fhc.append(ROOT.TFile(optimization.output_location+"users/ljf26/fluxfiles/g4lbne/"+optimization.g4lbne_version+"/"+optimization.physics_list+"/Optimizations-"+optimization_name+"FHC-"+str(int(configuration))+"/200kA/flux/histos_g4lbne_"+optimization.g4lbne_version+"_"+optimization.physics_list+"_Optimizations-"+optimization_name+"FHC-"+str(int(configuration))+"_200kA_"+optimization.detector_location_name+"_fastmc.root"))
111  var_flux_files_rhc.append(ROOT.TFile(optimization.output_location+"users/ljf26/fluxfiles/g4lbne/"+optimization.g4lbne_version+"/"+optimization.physics_list+"/Optimizations-"+optimization_name+"RHC-"+str(int(configuration))+"/-200kA/flux/histos_g4lbne_"+optimization.g4lbne_version+"_"+optimization.physics_list+"_Optimizations-"+optimization_name+"RHC-"+str(int(configuration))+"_-200kA_"+optimization.detector_location_name+"_fastmc.root"))
112  else:
113  var_flux_files_fhc.append(ROOT.TFile(optimization.output_location+"users/ljf26/fluxfiles/g4lbne/"+optimization.g4lbne_version+"/"+optimization.physics_list+"/Optimizations-"+optimization_name+"-"+str(int(configuration))+"/neutrino/flux/histos_g4lbne_"+optimization.g4lbne_version+"_"+optimization.physics_list+"_Optimizations-"+optimization_name+"-"+str(int(configuration))+"_neutrino_"+optimization.detector_location_name+"_fastmc.root"))
114  var_flux_files_rhc.append(ROOT.TFile(optimization.output_location+"users/ljf26/fluxfiles/g4lbne/"+optimization.g4lbne_version+"/"+optimization.physics_list+"/Optimizations-"+optimization_name+"-"+str(int(configuration))+"/antineutrino/flux/histos_g4lbne_"+optimization.g4lbne_version+"_"+optimization.physics_list+"_Optimizations-"+optimization_name+"-"+str(int(configuration))+"_antineutrino_"+optimization.detector_location_name+"_fastmc.root"))
115 
116  var_fluxes.append([])
117  var_fluxes[i].append(var_flux_files_fhc[i].Get("numu_flux"))
118  var_fluxes[i].append(var_flux_files_fhc[i].Get("numubar_flux"))
119  var_fluxes[i].append(var_flux_files_fhc[i].Get("nue_flux"))
120  var_fluxes[i].append(var_flux_files_fhc[i].Get("nuebar_flux"))
121  var_fluxes[i].append(var_flux_files_rhc[i].Get("numu_flux"))
122  var_fluxes[i].append(var_flux_files_rhc[i].Get("numubar_flux"))
123  var_fluxes[i].append(var_flux_files_rhc[i].Get("nue_flux"))
124  var_fluxes[i].append(var_flux_files_rhc[i].Get("nuebar_flux"))
125  fhcenergy = optimization.getEnergy(configuration,"FHC")
126  rhcenergy = optimization.getEnergy(configuration,"RHC")
127  print "Energies",fhcenergy,rhcenergy
128 
129  # Assume 7.5e13 protons / cycle and 1.2 s cycle time @120 GeV
130  # Scale other energies as a function of power
131 
132  scale_factor_fhc = 7.5e13/1.2*OptimizationUtils.GetPowerPOTScaleFactor(fhcenergy)
133  scale_factor_rhc = 7.5e13/1.2*OptimizationUtils.GetPowerPOTScaleFactor(rhcenergy)
134  print scale_factor_fhc
135  for j in range(0,4):
136  var_fluxes[i][j] = var_fluxes[i][j].Rebin(len(new_bins)-1,"var_flux_"+str(j),array('d',new_bins))
137  var_fluxes[i][j].Scale(scale_factor_fhc)
138  var_fluxes[i][j].Scale(1.0,"width") # per GeV
139  for j in range(4,8):
140  var_fluxes[i][j] = var_fluxes[i][j].Rebin(len(new_bins)-1,"var_flux_"+str(j),array('d',new_bins))
141  var_fluxes[i][j].Scale(scale_factor_rhc)
142  var_fluxes[i][j].Scale(1.0,"width") # per GeV
143 
144 
145 for j in range(0,8):
146  nom_flux[j] = nom_flux[j].Rebin(len(new_bins)-1,"nom_flux",array('d',new_bins))
147  nom_flux[j].Scale(7.5e13/1.2,"width") # per GeV
148 
149 # Plot all fluxes
150  canv = ROOT.TCanvas("MyCanvas","MyCanvas")
151  maxes = [nom_flux[j].GetMaximum()]
152 #for var_flux in var_fluxes:
153 # maxes.append[var_flux.GetMaximum()]
154  plot_max = max(maxes)*1.7
155 
156  colors = [2,4,8,6]
157  nom_flux[j].SetMaximum(plot_max)
158  nom_flux[j].SetTitle("")
159  nom_flux[j].Draw()
160 
161  frame = ROOT.TH2D("frame","frame",100,0,8,100,0,plot_max);
162  frame.GetXaxis().SetTitle(nom_flux[j].GetXaxis().GetTitle())
163  frame.GetYaxis().SetTitle(nom_flux[j].GetYaxis().GetTitle().replace("POT","s"))
164  frame.GetXaxis().CenterTitle()
165  frame.GetYaxis().CenterTitle()
166  frame.SetTitle("")
167  frame.Draw()
168 
169  nom_flux[j].Draw("HSAME")
170 
171 # for i in range(0,len(var_fluxes)):
172  for i in range(0,n_var_plots):
173  var_fluxes[i][j].SetLineColor(colors[i])
174  if(i==0):
175  var_fluxes[i][j].SetLineWidth(4);
176  else:
177  var_fluxes[i][j].SetLineWidth(2);
178  var_fluxes[i][j].Draw("HSAME")
179 
180  nom_flux[j].Draw("HSAME")
181 
182  leg = ROOT.TLegend(0.6,0.6,0.8,0.8)
183  leg.SetBorderSize(0);
184  leg.SetFillStyle(0);
185  leg.AddEntry(nom_flux[j],"Nominal","l");
186  labels = []
187  labels.append("Optimized (#"+str(int(configurations[ordered_indices[0]]))+")")
188  labels.append("2nd Best (#"+str(int(configurations[ordered_indices[1]]))+")")
189  labels.append("3rd Best (#"+str(int(configurations[ordered_indices[2]]))+")")
190  labels.append("4th Best (#"+str(int(configurations[ordered_indices[3]]))+")")
191 
192  for i in range(0,n_var_plots):
193  leg.AddEntry(var_fluxes[i][j],labels[i],"l")
194  leg.Draw()
195 
196  flavor_strings = ["numu_fhc","numubar_fhc","nue_fhc","nuebar_fhc","numu_rhc","numubar_rhc","nue_rhc","nuebar_rhc"]
197  canv.Print("plots/"+optimization_name+"_"+flavor_strings[j]+"_best_fluxes.eps");
198  canv.Print("plots/"+optimization_name+"_"+flavor_strings[j]+"_best_fluxes.png");
199 
200 
actions::Scale
Scale(size_t pos, T factor) -> Scale< T >
util::enumerate
auto enumerate(Iterables &&...iterables)
Range-for loop helper tracking the number of iteration.
Definition: enumerate.h:69
lar::dump::array
auto array(Array const &a)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:228
OptimizationUtils.GetPowerPOTScaleFactor
def GetPowerPOTScaleFactor(energy)
Definition: OptimizationUtils.py:59
max
static int max(int a, int b)
Definition: fortranscanner.cpp:60366
keras_to_tensorflow.int
int
Definition: keras_to_tensorflow.py:69
cet::split
void split(std::string const &s, char c, OutIter dest)
Definition: split.h:35
if
if(!yymsg) yymsg
str
static QCString str
Definition: fortrancode.cpp:27098
Generated by   doxygen 1.8.11