makeScanPlots.py
Go to the documentation of this file.
1 import ROOT, glob, sys, os, re
2 from array import array
3 import Optimizations, OptimizationUtils
4 
5 physics_list = "QGSP_BERT"
6 
7 # options: cp, mh, cp_fhcnumu, cp_fhcnumubar, cp_fhcnuenuebar, cp_rhcnumu, cp_rhcnumubar, cp_rhcnuenuebar
8 fitness_type = "cp"
9 
10 ##### PARSE INPUT ARGUMENT #####
11 
12 if not len(sys.argv) > 2:
13  print "Please specify an optimization name and configuration number"
14  sys.exit()
15 
16 optimization_name = sys.argv[1]
17 configuration = sys.argv[2]
18 
19 optimization = Optimizations.Optimization(optimization_name)
20 
21 #### merge configs if need ####
22 
23 n_scan_points = 30
24 
25 configs_already_merged = []
26 configs_merged_now = []
27 configs_not_finished = []
28 
29 for par_iter in range(len(optimization.parameter_names)):
30 
31  configs_already_merged.append([])
32  configs_merged_now.append([])
33  configs_not_finished.append([])
34 
35  #if "HornALength" not in optimization.parameter_names[par_iter]:
36  # continue
37 
38  lower_limit = optimization.parameter_lower_limits[par_iter]
39  upper_limit = optimization.parameter_upper_limits[par_iter]
40 
41  # Loop over values to scan for this parameter
42  for scan_iter in range(n_scan_points+1):
43  scan_value = lower_limit + scan_iter*(upper_limit-lower_limit)/n_scan_points
44  fhc_flux_dir = "/lbne/data/users/ljf26/fluxfiles/g4lbne/"+optimization.g4lbne_version+"/"+physics_list+"/Optimizations-"+optimization_name+"-"+configuration+"_Scan_"+optimization.parameter_names[par_iter]+"_"+str(scan_value);
45  rhc_flux_dir = "/lbne/data/users/ljf26/fluxfiles/g4lbne/"+optimization.g4lbne_version+"/"+physics_list+"/Optimizations-"+optimization_name+"-"+configuration+"_Scan_"+optimization.parameter_names[par_iter]+"_"+str(scan_value);
46  if(optimization.rhc_parameters_float_separate):
47  fhc_flux_dir = "/lbne/data/users/ljf26/fluxfiles/g4lbne/"+optimization.g4lbne_version+"/"+physics_list+"/Optimizations-"+optimization_name+"FHC-"+configuration+"_Scan_"+optimization.parameter_names[par_iter]+"_"+str(scan_value);
48  rhc_flux_dir = "/lbne/data/users/ljf26/fluxfiles/g4lbne/"+optimization.g4lbne_version+"/"+physics_list+"/Optimizations-"+optimization_name+"RHC-"+configuration+"_Scan_"+optimization.parameter_names[par_iter]+"_"+str(scan_value);
49 
50  merged_file = glob.glob(fhc_flux_dir+"/*/flux/histos_*_LBNEFD_fastmc.root")
51 
52  if not len(merged_file) == 2 :
53 
54  temp_fhc_dir = fhc_flux_dir.replace("/lbne/data/users/","/pnfs/lbne/scratch/users/")
55  temp_rhc_dir = rhc_flux_dir.replace("/lbne/data/users/","/pnfs/lbne/scratch/users/")
56 
57  fhc_files = glob.glob(temp_fhc_dir+"/neutrino/flux/histos/histos_*_neutrino_*_LBNEFD_fastmc.root")
58  rhc_files = glob.glob(temp_rhc_dir+"/antineutrino/flux/histos/histos_*_antineutrino_*_LBNEFD_fastmc.root")
59 
60  print "TEMP_FHC_FLUX_DIR",temp_fhc_dir+"/neutrino/flux/histos/histos_*_neutrino_*_LBNEFD_fastmc.root"
61 
62  if len(fhc_files) == 5 and len(rhc_files) == 5:
63 
64  print "Merging",optimization.parameter_names[par_iter],scan_value
65 
66  if not optimization.rhc_parameters_float_separate:
67  os.system("python ../../ProductionScripts/merge_histograms.py --mode neutrino --n 100000 -p "+physics_list+" -m Optimizations-"+optimization_name+"-"+configuration+"_Scan_"+optimization.parameter_names[par_iter]+"_"+str(scan_value))
68  os.system("python ../../ProductionScripts/merge_histograms.py --mode antineutrino --n 100000 -p "+physics_list+" -m Optimizations-"+optimization_name+"-"+configuration+"_Scan_"+optimization.parameter_names[par_iter]+"_"+str(scan_value))
69  else :
70  os.system("python ../../ProductionScripts/merge_histograms.py --mode neutrino --n 100000 -p "+physics_list+" -m Optimizations-"+optimization_name+"FHC-"+configuration+"_Scan_"+optimization.parameter_names[par_iter]+"_"+str(scan_value))
71  os.system("python ../../ProductionScripts/merge_histograms.py --mode antineutrino --n 100000 -p "+physics_list+" -m Optimizations-"+optimization_name+"RHC-"+configuration+"_Scan_"+optimization.parameter_names[par_iter]+"_"+str(scan_value))
72  configs_merged_now[par_iter].append(scan_value)
73  else:
74  print "Not all files finished for ",optimization.parameter_names[par_iter],scan_value
75  configs_not_finished[par_iter].append(scan_value)
76  else:
77  configs_already_merged[par_iter].append(scan_value)
78 
79 
80 for par_iter in range(len(optimization.parameter_names)):
81 
82  #if "HornALength" not in optimization.parameter_names[par_iter]:
83  # continue
84 
85 
86  par_name = optimization.parameter_names[par_iter]
87  lower_limit = optimization.parameter_lower_limits[par_iter]
88  upper_limit = optimization.parameter_upper_limits[par_iter]
89 
90  # Loop over values to scan for this parameter
91  par_values = []
92  fit_values = []
93  zero_errors = []
94  fit_errors = []
95 
96  for scan_iter in range(n_scan_points+1):
97  scan_value = lower_limit + scan_iter*(upper_limit-lower_limit)/n_scan_points
98 
99  macro_name_fhc = "Optimizations-"+optimization_name+"-"+configuration+"_Scan_"+optimization.parameter_names[par_iter]+"_"+str(scan_value)
100  macro_name_rhc = "Optimizations-"+optimization_name+"-"+configuration+"_Scan_"+optimization.parameter_names[par_iter]+"_"+str(scan_value)
101  if optimization.rhc_parameters_float_separate:
102  macro_name_fhc = "Optimizations-"+optimization_name+"FHC-"+configuration+"_Scan_"+optimization.parameter_names[par_iter]+"_"+str(scan_value)
103  macro_name_rhc = "Optimizations-"+optimization_name+"RHC-"+configuration+"_Scan_"+optimization.parameter_names[par_iter]+"_"+str(scan_value)
104 
105  nomfhcfile = "/lbne/data/users/ljf26/fluxfiles/g4lbne/v3r2p4/"+physics_list+"/Nominal/neutrino/flux/histos_g4lbne_v3r2p4_"+physics_list+"_Nominal_neutrino_LBNEFD_fastmc.root"
106 
107  nomrhcfile = "/lbne/data/users/ljf26/fluxfiles/g4lbne/v3r2p4/"+physics_list+"/Nominal/antineutrino/flux/histos_g4lbne_v3r2p4_"+physics_list+"_Nominal_antineutrino_LBNEFD_fastmc.root"
108 
109  varfhcfile = "/lbne/data/users/ljf26/fluxfiles/g4lbne/"+optimization.g4lbne_version+"/"+physics_list+"/"+macro_name_fhc+"/neutrino/flux/histos_g4lbne_"+optimization.g4lbne_version+"_"+physics_list+"_"+macro_name_fhc+"_neutrino_LBNEFD_fastmc.root"
110  varrhcfile = "/lbne/data/users/ljf26/fluxfiles/g4lbne/"+optimization.g4lbne_version+"/"+physics_list+"/"+macro_name_rhc+"/antineutrino/flux/histos_g4lbne_"+optimization.g4lbne_version+"_"+physics_list+"_"+macro_name_rhc+"_antineutrino_LBNEFD_fastmc.root"
111 
112  # look up power for this configuration
113  rhcpower = 0
114  fhcpower = 0
115  if(par_name != "ProtonEnergy" and
116  par_name != "FHCProtonEnergy" and
117  par_name != "RHCProtonEnergy"):
118  power_iter = optimization.getParameterIndex("ProtonEnergy")
119 
120  fhcpower = optimization.getEnergy(configuration,"FHC")
121  rhcpower = optimization.getEnergy(configuration,"RHC")
122  elif par_name == "ProtonEnergy":
123  fhcpower = scan_value
124  rhcpower = scan_value
125  elif par_name == "FHCProtonEnergy":
126  fhcpower = scan_value
127  rhcpower = optimization.getEnergy(configuration,"RHC")
128  elif par_name == "RHCProtonEnergy":
129  rhcpower = scan_value
130  fhcpower = optimization.getEnergy(configuration,"FHC")
131 
132  # look up antinu fraction for this configuration
133  antinufrac = optimization.getParameterValue("AntinuFraction",configuration)
134  if(antinufrac == -1):
135  antinufrac = 0.5
136  print antinufrac
137 
138  # Get whatever fitness the user asked for
139  [fitness_names,fitnesses,fiterrors] = OptimizationUtils.ComputeMean75PercentSensitivityAndError(varfhcfile,varrhcfile,fhcpower,rhcpower,antinufrac,"cp")
140  if(fitness_type=="mh"):
141  [fitness_names,fitnesses,fiterrors] = OptimizationUtils.ComputeMean75PercentSensitivityAndError(varfhcfile,varrhcfile,fhcpower,rhcpower,antinufrac,"mh")
142  if fitness_type in ["mh","cp"]:
143  t_fitness = fitnesses[0]
144  t_error = fiterrors[0]
145  if fitness_type =="cp_fhcnumu":
146  t_fitness = fitnesses[1]
147  t_error = fiterrors[1]
148  if fitness_type =="cp_fhcnumubar":
149  t_fitness = fitnesses[2]
150  t_error = fiterrors[2]
151  if fitness_type =="cp_fhcnuenuebar":
152  t_fitness = fitnesses[3]
153  t_error = fiterrors[3]
154  if fitness_type =="cp_rhcnumu":
155  t_fitness = fitnesses[4]
156  t_error = fiterrors[4]
157  if fitness_type =="cp_rhcnumubar":
158  t_fitness = fitnesses[5]
159  t_error = fiterrors[5]
160  if fitness_type =="cp_rhcnuenuebar":
161  t_fitness = fitnesses[6]
162  t_error = fiterrors[6]
163 
164  print t_fitness
165 
166 # if t_fitness >0:
167  if 1>0:
168  par_values.append(scan_value)
169  fit_values.append(t_fitness)
170  if(fiterrors[0]<10):
171  fit_errors.append(t_error)
172  else:
173  fit_errors.append(0);
174  zero_errors.append(0)
175 
176 # Make TGraphs
177  canv = ROOT.TCanvas("MyCanvas","MyCanvas")
178  print par_values
179  print zero_errors
180  print fit_values
181  print fit_errors
182  if len(par_values)==0:
183  continue
184  ROOT.gStyle.SetMarkerStyle(2);
185 # fitness_vs_par = ROOT.TGraphErrors(len(par_values),array('d',par_values),array('d',zero_errors),array('d',fit_values),array('d',fit_errors));
186  fitness_vs_par = ROOT.TGraphErrors(len(par_values),array('d',par_values),array('d',fit_values),array('d',zero_errors),array('d',fit_errors));
187  fitness_vs_par.GetXaxis().SetTitle(par_name);
188  fitness_vs_par.GetYaxis().SetTitle("Fitness");
189  fitness_vs_par.SetMinimum(1.5);
190  fitness_vs_par.SetMaximum(2.0);
191  if(fitness_type=="mh"):
192  fitness_vs_par.SetMinimum(5.0);
193  fitness_vs_par.SetMaximum(6.5);
194  if(fitness_type in ["cp_fhcnumu","cp_fhcnumubar","cp_fhcnuenuebar","cp_rhcnumu","cp_rhcnumubar","cp_rhcnuenuebar"]):
195  fitness_vs_par.SetMinimum(-0.1);
196  fitness_vs_par.SetMaximum(0.4);
197  ROOT.gPad.SetGrid(1);
198  fitness_vs_par.SetTitle("")
199  fitness_vs_par.Draw("AP")
200 
201  # Draw line at value chosen by optimization
202  optimized_value = optimization.getParameterValue(optimization.parameter_names[par_iter],configuration)
203 
204  my_line = ROOT.TLine(optimized_value,fitness_vs_par.GetMinimum(),optimized_value,fitness_vs_par.GetMaximum())
205  my_line.SetLineColor(2);
206  my_line.Draw();
207 
208  canv.Print("scan_plots/"+optimization_name+"_"+configuration+"_Scan_"+par_name+"_"+physics_list+"_"+fitness_type+".eps")
209  canv.Print("scan_plots/"+optimization_name+"_"+configuration+"_Scan_"+par_name+"_"+physics_list+"_"+fitness_type+".png")
210 
211 
212 """
213  if fitness_type=="mh":
214  fitness_vs_parameter = ROOT.TGraphErrors(len(configurations),array('d',parameter_values[i]),array('d',mhfitnesses),array('d',zero_errors),array('d',mhfitness_errors))
215  if fitness_type=="ih":
216  fitness_vs_parameter = ROOT.TGraphErrors(len(configurations),array('d',parameter_values[i]),array('d',ihfitnesses),array('d',zero_errors),array('d',ihfitness_errors))
217  if fitness_type=="nh":
218  fitness_vs_parameter = ROOT.TGraphErrors(len(configurations),array('d',parameter_values[i]),array('d',nhfitnesses),array('d',zero_errors),array('d',nhfitness_errors))
219  fitness_vs_parameter.GetYaxis().SetTitle(fitness_type+" Fitness")
220  x_axis_title = optimization.parameter_names[i]
221  if(optimization.parameter_units[i]!=""):
222  x_axis_title = x_axis_title+" ("+optimization.parameter_units[i]+")"
223  fitness_vs_parameter.GetXaxis().SetTitle(x_axis_title);
224  fitness_vs_parameter.GetXaxis().CenterTitle()
225  fitness_vs_parameter.GetYaxis().CenterTitle()
226  fitness_vs_parameter.Draw("AP");
227  canv.Print("plots/"+optimization_name+"_"+fitness_type+"_fitness_vs_"+parameter_name+".eps")
228  canv.Print("plots/"+optimization_name+"_"+fitness_type+"_fitness_vs_"+parameter_name+".png")
229 
230 
231  for j in range(i+1,n_parameters):
232  parameter_vs_parameter = ROOT.TGraphErrors(1000,array('d',parameter_values[i][-1000:]),array('d',parameter_values[j][-1000:]),array('d',zero_errors[-1000:]),array('d',zero_errors[-1000:]))
233  x_axis_title = optimization.parameter_names[i]
234  y_axis_title = optimization.parameter_names[j]
235  if(optimization.parameter_units[i]!=""):
236  x_axis_title = x_axis_title+" ("+optimization.parameter_units[i]+")"
237  if(optimization.parameter_units[j]!=""):
238  y_axis_title = y_axis_title+" ("+optimization.parameter_units[j]+")"
239  parameter_vs_parameter.GetXaxis().SetTitle(x_axis_title);
240  parameter_vs_parameter.GetYaxis().SetTitle(y_axis_title);
241  parameter_vs_parameter.GetXaxis().CenterTitle()
242  parameter_vs_parameter.GetYaxis().CenterTitle()
243  parameter_vs_parameter.Draw("AP");
244  canv.Print("plots/"+optimization_name+"_"+parameter_name+"_vs_"+optimization.parameter_names[j]+".eps")
245  canv.Print("plots/"+optimization_name+"_"+parameter_name+"_vs_"+optimization.parameter_names[j]+".png")
246 
247 """
OptimizationUtils.ComputeMean75PercentSensitivityAndError
def ComputeMean75PercentSensitivityAndError(varfhcfile, varrhcfile, fhcvarenergy=120.0, rhcvarenergy=120.0, antinufrac=0.5, quantity="cp")
Definition: OptimizationUtils.py:81
lar::dump::array
auto array(Array const &a)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:228
if
if(!yymsg) yymsg
str
static QCString str
Definition: fortrancode.cpp:27098
Generated by   doxygen 1.8.11