OptimizationUtils.py
Go to the documentation of this file.
1 import ROOT, random, array, sys, os, math
2 
3 def ConvertGraphToHisto(new_bins,pGraph):
4 
5  name = "h_"+pGraph.GetName();
6  histo = ROOT.TH1D(name,name,len(new_bins)-1,array.array('d',new_bins));
7 
8  pGraph.Sort();
9 
10  if pGraph.GetN() != len(new_bins)-1:
11  print "OptimizationUtils.ConvertGraphToHisto: WARNING converting graph with",pGraph.GetN(),"points to histogram with",len(new_bins)-1,"bins"
12 
13  for i in range(0,len(new_bins)):
14  x = ROOT.Double(0.0)
15  y = ROOT.Double(0.0)
16  pGraph.GetPoint(i,x,y);
17  histo.SetBinContent(i+1,y)
18  histo.SetBinError(i+1,0)
19 
20  return histo
21 
22 def GetPercentile(histo,X):
23  n = histo.GetNbinsX();
24  x = []
25  y = []
26 
27  for i in range(0,n):
28  x.append(histo.GetBinCenter(i+1));
29  y.append(histo.GetBinContent(i+1));
30  y_array = array.array('d',y)
31 
32  percentiles = array.array('d',[0.0]);
33  probs = array.array('d',[1-float(X)/100.0])
34 
35  ROOT.TMath.Quantiles(n,1,y_array,percentiles,probs,False);
36  return percentiles[0]
37 
38 def natural_sort(l):
39  convert = lambda text: int(text) if text.isdigit() else text.lower()
40  alphanum_key = lambda key: [ convert(c) for c in re.split('([0-9]+)', key) ]
41  return sorted(l, key = alphanum_key)
42 
43 
44 
45 def GetFractionAtMoreThanXSigma(histo, X):
46  n = histo.GetNbinsX();
47 
48  greaterThanX = 0.0;
49  total = 0.0;
50 
51  for i in range(0,n):
52  total = total + float(histo.GetBinWidth(i+1));
53  if(histo.GetBinContent(i+1)>=X):
54  greaterThanX = greaterThanX + float(histo.GetBinWidth(i+1));
55 
56 
57  return greaterThanX/total;
58 
59 def GetPowerPOTScaleFactor(energy):
60  if(energy >= 120):
61  return 120/energy
62  if(energy < 60):
63  return 1.71e21/1e21; # limited to 0.7 s cycle time
64  # linear interpolate
65  energies = [60,70,80,90,100,110,120]
66  potperyear = [1.71e21,1.5e21,1.33e21,1.23e21,1.14e21,1.07e21,1e21]
67  for i in range(0,len(energies)):
68  if energy < energies[i]:
69  interp_high = i
70  interp_low = i-1
71  break
72 
73  scale_factor = (potperyear[interp_high]+(potperyear[interp_low]-potperyear[interp_high])/(energies[interp_high]-energies[interp_low])*(energies[interp_high]-energy))/1e21
74  return scale_factor
75 
76 
77 def MakeSumw2(histo):
78  if histo.GetSumw2N() == 0:
79  histos.Sumw2();
80 
81 def ComputeMean75PercentSensitivityAndError(varfhcfile,varrhcfile, fhcvarenergy=120.0,rhcvarenergy = 120.0, antinufrac=0.5,quantity="cp"):
82 
83  nomfhcfile = "/dune/data/users/ljf26/fluxfiles/g4lbne/v3r2p4/QGSP_BERT/Nominal/200kA/flux/histos_g4lbne_v3r2p4_QGSP_BERT_Nominal_200kA_LBNEFD_fastmc.root"
84 
85  nomrhcfile = "/dune/data/users/ljf26/fluxfiles/g4lbne/v3r2p4/QGSP_BERT/Nominal/-200kA/flux/histos_g4lbne_v3r2p4_QGSP_BERT_Nominal_-200kA_LBNEFD_fastmc.root"
86 
87  # change sensitivity type string to get file names right
88  if quantity == "cp":
89  quantity = "cp_75thpercentile"
90  if quantity == "mh":
91  quantity = "mh_minimum"
92 
93  # Have to use same bins that were used for the FMC bin-by-bin variation study
94  new_bins = [0.0,0.5,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,15.0,20.0,120.0]
95 
96  null_return = [[-1,-1,-1,-1,-1,-1,-1],[-1,-1,-1,-1,-1,-1,-1],[-1,-1,-1,-1,-1,-1,-1]]
97 
98  # Open files with flux histograms and make sure they are good
99  nomfhctfile = ROOT.TFile(nomfhcfile)
100  nomrhctfile = ROOT.TFile(nomrhcfile)
101  varfhctfile = ROOT.TFile(varfhcfile)
102  varrhctfile = ROOT.TFile(varrhcfile)
103 
104  fhcnumu_nh_weightfile = ROOT.TFile("sensitivity_plots/fhcnumu_nh_"+quantity+"_fixedScale.root");
105  fhcnumu_ih_weightfile = ROOT.TFile("sensitivity_plots/fhcnumu_ih_"+quantity+"_fixedScale.root");
106  rhcnumu_nh_weightfile = ROOT.TFile("sensitivity_plots/rhcnumu_nh_"+quantity+"_fixedScale.root");
107  rhcnumu_ih_weightfile = ROOT.TFile("sensitivity_plots/rhcnumu_ih_"+quantity+"_fixedScale.root");
108  fhcnumubar_nh_weightfile = ROOT.TFile("sensitivity_plots/fhcnumubar_nh_"+quantity+"_fixedScale.root");
109  fhcnumubar_ih_weightfile = ROOT.TFile("sensitivity_plots/fhcnumubar_ih_"+quantity+"_fixedScale.root");
110  rhcnumubar_nh_weightfile = ROOT.TFile("sensitivity_plots/rhcnumubar_nh_"+quantity+"_fixedScale.root");
111  rhcnumubar_ih_weightfile = ROOT.TFile("sensitivity_plots/rhcnumubar_ih_"+quantity+"_fixedScale.root");
112  fhcnuenuebar_nh_weightfile = ROOT.TFile("sensitivity_plots/fhcnuenuebar_nh_"+quantity+"_fixedScale.root");
113  fhcnuenuebar_ih_weightfile = ROOT.TFile("sensitivity_plots/fhcnuenuebar_ih_"+quantity+"_fixedScale.root");
114  rhcnuenuebar_nh_weightfile = ROOT.TFile("sensitivity_plots/rhcnuenuebar_nh_"+quantity+"_fixedScale.root");
115  rhcnuenuebar_ih_weightfile = ROOT.TFile("sensitivity_plots/rhcnuenuebar_ih_"+quantity+"_fixedScale.root");
116 
117  files_to_check = [nomfhctfile,nomrhctfile,varfhctfile,varrhctfile]
118  for a_file in files_to_check:
119  if(a_file.IsZombie()):
120  print "Optimizationutils::ComputeMeanPercentCPSensitivity: Could not open",a_file.GetName()
121  return null_return
122 
123  # Get Flux histograms for nominal case and variation
124  nom_fhcnumu = nomfhctfile.Get("numu_flux").Rebin(len(new_bins)-1,"h_nomfhcnumu_flux",array.array('d',new_bins))
125  nom_fhcnumubar = nomfhctfile.Get("numubar_flux").Rebin(len(new_bins)-1,"h_nomfhcnumubar_flux",array.array('d',new_bins))
126  nom_fhcnue = nomfhctfile.Get("nue_flux").Rebin(len(new_bins)-1,"h_nomfhcnue_flux",array.array('d',new_bins))
127  nom_fhcnuebar = nomfhctfile.Get("nuebar_flux").Rebin(len(new_bins)-1,"h_nomfhcnuebar_flux",array.array('d',new_bins))
128 
129  nom_rhcnumu = nomrhctfile.Get("numu_flux").Rebin(len(new_bins)-1,"h_nomrhcnumu_flux",array.array('d',new_bins))
130  nom_rhcnumubar = nomrhctfile.Get("numubar_flux").Rebin(len(new_bins)-1,"h_nomrhcnumubar_flux",array.array('d',new_bins))
131  nom_rhcnue = nomrhctfile.Get("nue_flux").Rebin(len(new_bins)-1,"h_nomrhcnue_flux",array.array('d',new_bins))
132  nom_rhcnuebar = nomrhctfile.Get("nuebar_flux").Rebin(len(new_bins)-1,"h_nomrhcnuebar_flux",array.array('d',new_bins))
133 
134  var_fhcnumu = varfhctfile.Get("numu_flux").Rebin(len(new_bins)-1,"h_fhcnumu_flux",array.array('d',new_bins))
135  var_fhcnumubar = varfhctfile.Get("numubar_flux").Rebin(len(new_bins)-1,"h_fhcnumubar_flux",array.array('d',new_bins))
136  var_fhcnue = varfhctfile.Get("nue_flux").Rebin(len(new_bins)-1,"h_fhcnue_flux",array.array('d',new_bins))
137  var_fhcnuebar = varfhctfile.Get("nuebar_flux").Rebin(len(new_bins)-1,"h_fhcnuebar_flux",array.array('d',new_bins))
138 
139  var_rhcnumu = varrhctfile.Get("numu_flux").Rebin(len(new_bins)-1,"h_rhcnumu_flux",array.array('d',new_bins))
140  var_rhcnumubar = varrhctfile.Get("numubar_flux").Rebin(len(new_bins)-1,"h_rhcnumubar_flux",array.array('d',new_bins))
141  var_rhcnue = varrhctfile.Get("nue_flux").Rebin(len(new_bins)-1,"h_rhcnue_flux",array.array('d',new_bins))
142  var_rhcnuebar = varrhctfile.Get("nuebar_flux").Rebin(len(new_bins)-1,"h_rhcnuebar_flux",array.array('d',new_bins))
143 
144  nom_histos = [nom_fhcnumu, nom_fhcnumubar, nom_fhcnue, nom_fhcnuebar, nom_rhcnumu, nom_rhcnumubar, nom_rhcnue, nom_rhcnuebar]
145  var_histos = [var_fhcnumu, var_fhcnumubar, var_fhcnue, var_fhcnuebar, var_rhcnumu, var_rhcnumubar, var_rhcnue, var_rhcnuebar]
146 
147  # Rebin makes the last bin screwy sometimes, so fix these bins
148  # Fluxes should be small (<10^-10) postive numbers
149  for histo in nom_histos+var_histos:
150  for i in range(0,histo.GetNbinsX()):
151  if(histo.GetBinContent(i+1) != histo.GetBinContent(i+1) or
152  histo.GetBinContent(i+1)<-10000 or histo.GetBinContent(i+1)>1):
153  print "WARNING: the",i+1,"th bin in",histo.GetName(),"is",histo.GetBinContent(i+1),"... Resetting it to zero"
154  histo.SetBinContent(i+1,0)
155  histo.SetBinError(i+1,0)
156 
157 
158  # Make sure errors get propagated
159  for histo in nom_histos+var_histos:
160  MakeSumw2(histo)
161 
162  # Scale variations for powers different than 120
163  fhc_scale_factor = GetPowerPOTScaleFactor(fhcvarenergy)
164  rhc_scale_factor = GetPowerPOTScaleFactor(rhcvarenergy)
165 
166  # Scale variations for antinu fraction other than 0.5
167  antinu_scale = antinufrac/0.5*rhc_scale_factor;
168  nu_scale = (1.0-antinufrac)/0.5*fhc_scale_factor;
169 
170  var_fhcnumu.Scale(nu_scale);
171  var_fhcnumubar.Scale(nu_scale);
172  var_fhcnue.Scale(nu_scale);
173  var_fhcnuebar.Scale(nu_scale);
174 
175  var_rhcnumu.Scale(antinu_scale);
176  var_rhcnumubar.Scale(antinu_scale);
177  var_rhcnue.Scale(antinu_scale);
178  var_rhcnuebar.Scale(antinu_scale);
179 
180  # Compute change in sensitivity
181  fhcnumu = var_fhcnumu.Clone()
182  fhcnumu.Add(nom_fhcnumu,-1)
183  fhcnumu.Divide(nom_fhcnumu)
184  #fhcnumu.Scale(1/0.1)
185  fhcnumu.Scale(100)
186  fhcnumu_ih = fhcnumu.Clone()
187  fhcnumu_nh = fhcnumu.Clone()
188 
189  fhcnumubar = var_fhcnumubar.Clone()
190  fhcnumubar.Add(nom_fhcnumubar,-1)
191  fhcnumubar.Divide(nom_fhcnumubar)
192  fhcnumubar.Scale(100)
193  fhcnumubar_ih = fhcnumubar.Clone()
194  fhcnumubar_nh = fhcnumubar.Clone()
195 
196  fhcnuenuebar = var_fhcnue.Clone()
197  fhcnuenuebar.Add(nom_fhcnue,-1)
198  fhcnuenuebar.Divide(nom_fhcnue)
199  fhcnuenuebar.Scale(100)
200  fhcnuenuebar_ih = fhcnuenuebar.Clone()
201  fhcnuenuebar_nh = fhcnuenuebar.Clone()
202 
203  rhcnumu = var_rhcnumu.Clone()
204  rhcnumu.Add(nom_rhcnumu,-1)
205  rhcnumu.Divide(nom_rhcnumu)
206  rhcnumu.Scale(100)
207  rhcnumu_ih = rhcnumu.Clone()
208  rhcnumu_nh = rhcnumu.Clone()
209 
210  rhcnumubar = var_rhcnumubar.Clone()
211  rhcnumubar.Add(nom_rhcnumubar,-1)
212  rhcnumubar.Divide(nom_rhcnumubar)
213  rhcnumubar.Scale(100)
214  rhcnumubar_ih = rhcnumubar.Clone()
215  rhcnumubar_nh = rhcnumubar.Clone()
216 
217  rhcnuenuebar = var_rhcnuebar.Clone()
218  rhcnuenuebar.Add(nom_rhcnuebar,-1)
219  rhcnuenuebar.Divide(nom_rhcnuebar)
220  rhcnuenuebar.Scale(100)
221  rhcnuenuebar_ih = rhcnuenuebar.Clone()
222  rhcnuenuebar_nh = rhcnuenuebar.Clone()
223 
224  tot_fhcnumu_ih = 0;
225  tot_fhcnumu_nh = 0;
226  tot_fhcnumubar_ih = 0;
227  tot_fhcnumubar_nh = 0;
228  tot_fhcnuenuebar_ih = 0;
229  tot_fhcnuenuebar_nh = 0;
230  tot_rhcnumu_ih = 0;
231  tot_rhcnumu_nh = 0;
232  tot_rhcnumubar_ih = 0;
233  tot_rhcnumubar_nh = 0;
234  tot_rhcnuenuebar_ih = 0;
235  tot_rhcnuenuebar_nh = 0;
236 
237  err_fhcnumu_ih = 0;
238  err_fhcnumu_nh = 0;
239  err_fhcnumubar_ih = 0;
240  err_fhcnumubar_nh = 0;
241  err_fhcnuenuebar_ih = 0;
242  err_fhcnuenuebar_nh = 0;
243  err_rhcnumu_ih = 0;
244  err_rhcnumu_nh = 0;
245  err_rhcnumubar_ih = 0;
246  err_rhcnumubar_nh = 0;
247  err_rhcnuenuebar_ih = 0;
248  err_rhcnuenuebar_nh = 0;
249 
250  for i in range(fhcnumu.GetNbinsX()):
251 
252  fhcnumu_ih.SetBinContent(i+1,fhcnumu_ih_weightfile.Get("Graph;"+str(i+1)).Eval(fhcnumu.GetBinContent(i+1)+100));
253  if(fhcnumu.GetBinContent(i+1) != 0):
254  fhcnumu_ih.SetBinError(i+1,fhcnumu.GetBinError(i+1)/fhcnumu.GetBinContent(i+1)*fhcnumu_ih.GetBinContent(i+1));
255  rhcnumu_ih.SetBinContent(i+1,rhcnumu_ih_weightfile.Get("Graph;"+str(i+1)).Eval(rhcnumu.GetBinContent(i+1)+100));
256  if(rhcnumu.GetBinContent(i+1) != 0):
257  rhcnumu_ih.SetBinError(i+1,rhcnumu.GetBinError(i+1)/rhcnumu.GetBinContent(i+1)*rhcnumu_ih.GetBinContent(i+1));
258 
259  fhcnumubar_ih.SetBinContent(i+1,fhcnumubar_ih_weightfile.Get("Graph;"+str(i+1)).Eval(fhcnumubar.GetBinContent(i+1)+100));
260  if(fhcnumubar.GetBinContent(i+1)!=0):
261  fhcnumubar_ih.SetBinError(i+1,fhcnumubar.GetBinError(i+1)/fhcnumubar.GetBinContent(i+1)*fhcnumubar_ih.GetBinContent(i+1));
262  rhcnumubar_ih.SetBinContent(i+1,rhcnumubar_ih_weightfile.Get("Graph;"+str(i+1)).Eval(rhcnumubar.GetBinContent(i+1)+100));
263  if(rhcnumubar.GetBinContent(i+1)!=0):
264  rhcnumubar_ih.SetBinError(i+1,rhcnumubar.GetBinError(i+1)/rhcnumubar.GetBinContent(i+1)*rhcnumubar_ih.GetBinContent(i+1));
265 
266  fhcnuenuebar_ih.SetBinContent(i+1,fhcnuenuebar_ih_weightfile.Get("Graph;"+str(i+1)).Eval(fhcnuenuebar.GetBinContent(i+1)+100));
267  if(fhcnuenuebar.GetBinContent(i+1)!=0):
268  fhcnuenuebar_ih.SetBinError(i+1,fhcnuenuebar.GetBinError(i+1)/fhcnuenuebar.GetBinContent(i+1)*fhcnuenuebar_ih.GetBinContent(i+1));
269  rhcnuenuebar_ih.SetBinContent(i+1,rhcnuenuebar_ih_weightfile.Get("Graph;"+str(i+1)).Eval(rhcnuenuebar.GetBinContent(i+1)+100));
270  if(rhcnuenuebar.GetBinContent(i+1)!=0):
271  rhcnuenuebar_ih.SetBinError(i+1,rhcnuenuebar.GetBinError(i+1)/rhcnuenuebar.GetBinContent(i+1)*rhcnuenuebar_ih.GetBinContent(i+1));
272 
273  fhcnumu_nh.SetBinContent(i+1,fhcnumu_nh_weightfile.Get("Graph;"+str(i+1)).Eval(fhcnumu.GetBinContent(i+1)+100));
274  if(fhcnumu.GetBinContent(i+1)!=0):
275  fhcnumu_nh.SetBinError(i+1,fhcnumu.GetBinError(i+1)/fhcnumu.GetBinContent(i+1)*fhcnumu_nh.GetBinContent(i+1));
276  rhcnumu_nh.SetBinContent(i+1,rhcnumu_nh_weightfile.Get("Graph;"+str(i+1)).Eval(rhcnumu.GetBinContent(i+1)+100));
277  if(rhcnumu.GetBinContent(i+1)!=0):
278  rhcnumu_nh.SetBinError(i+1,rhcnumu.GetBinError(i+1)/rhcnumu.GetBinContent(i+1)*rhcnumu_nh.GetBinContent(i+1));
279 
280  fhcnumubar_nh.SetBinContent(i+1,fhcnumubar_nh_weightfile.Get("Graph;"+str(i+1)).Eval(fhcnumubar.GetBinContent(i+1)+100));
281  if(fhcnumubar.GetBinContent(i+1)!=0):
282  fhcnumubar_nh.SetBinError(i+1,fhcnumubar.GetBinError(i+1)/fhcnumubar.GetBinContent(i+1)*fhcnumubar_nh.GetBinContent(i+1));
283  rhcnumubar_nh.SetBinContent(i+1,rhcnumubar_nh_weightfile.Get("Graph;"+str(i+1)).Eval(rhcnumubar.GetBinContent(i+1)+100));
284  if(rhcnumubar.GetBinContent(i+1)!=0):
285  rhcnumubar_nh.SetBinError(i+1,rhcnumubar.GetBinError(i+1)/rhcnumubar.GetBinContent(i+1)*rhcnumubar_nh.GetBinContent(i+1));
286 
287  fhcnuenuebar_nh.SetBinContent(i+1,fhcnuenuebar_nh_weightfile.Get("Graph;"+str(i+1)).Eval(fhcnuenuebar.GetBinContent(i+1)+100));
288  if(fhcnuenuebar.GetBinContent(i+1)!=0):
289  fhcnuenuebar_nh.SetBinError(i+1,fhcnuenuebar.GetBinError(i+1)/fhcnuenuebar.GetBinContent(i+1)*fhcnuenuebar_nh.GetBinContent(i+1));
290  rhcnuenuebar_nh.SetBinContent(i+1,rhcnuenuebar_nh_weightfile.Get("Graph;"+str(i+1)).Eval(rhcnuenuebar.GetBinContent(i+1)+100));
291  if(rhcnuenuebar.GetBinContent(i+1)!=0):
292  rhcnuenuebar_nh.SetBinError(i+1,rhcnuenuebar.GetBinError(i+1)/rhcnuenuebar.GetBinContent(i+1)*rhcnuenuebar_nh.GetBinContent(i+1));
293 
294  tot_fhcnumu_ih += fhcnumu_ih.GetBinContent(i+1)
295  tot_fhcnumu_nh += fhcnumu_nh.GetBinContent(i+1)
296  tot_fhcnumubar_ih += fhcnumubar_ih.GetBinContent(i+1)
297  tot_fhcnumubar_nh += fhcnumubar_nh.GetBinContent(i+1)
298  tot_fhcnuenuebar_ih += fhcnuenuebar_ih.GetBinContent(i+1)
299  tot_fhcnuenuebar_nh += fhcnuenuebar_nh.GetBinContent(i+1)
300 
301  tot_rhcnumu_ih += rhcnumu_ih.GetBinContent(i+1)
302  tot_rhcnumu_nh += rhcnumu_nh.GetBinContent(i+1)
303  tot_rhcnumubar_ih += rhcnumubar_ih.GetBinContent(i+1)
304  tot_rhcnumubar_nh += rhcnumubar_nh.GetBinContent(i+1)
305  tot_rhcnuenuebar_ih += rhcnuenuebar_ih.GetBinContent(i+1)
306  tot_rhcnuenuebar_nh += rhcnuenuebar_nh.GetBinContent(i+1)
307 
308  err_fhcnumu_ih +=math.pow(fhcnumu_ih.GetBinError(i+1),2)
309  err_fhcnumu_nh +=math.pow(fhcnumu_nh.GetBinError(i+1),2)
310  err_fhcnumubar_ih +=math.pow(fhcnumubar_ih.GetBinError(i+1),2)
311  err_fhcnumubar_nh +=math.pow(fhcnumubar_nh.GetBinError(i+1),2)
312  err_fhcnuenuebar_ih +=math.pow(fhcnuenuebar_ih.GetBinError(i+1),2)
313  err_fhcnuenuebar_nh +=math.pow(fhcnuenuebar_nh.GetBinError(i+1),2)
314 
315  err_rhcnumu_ih +=math.pow(rhcnumu_ih.GetBinError(i+1),2)
316  err_rhcnumu_nh +=math.pow(rhcnumu_nh.GetBinError(i+1),2)
317  err_rhcnumubar_ih +=math.pow(rhcnumubar_ih.GetBinError(i+1),2)
318  err_rhcnumubar_nh +=math.pow(rhcnumubar_nh.GetBinError(i+1),2)
319  err_rhcnuenuebar_ih +=math.pow(rhcnuenuebar_ih.GetBinError(i+1),2)
320  err_rhcnuenuebar_nh +=math.pow(rhcnuenuebar_nh.GetBinError(i+1),2)
321 
322  err_fhcnumu_ih = math.sqrt(err_fhcnumu_ih)
323  err_fhcnumu_nh = math.sqrt(err_fhcnumu_nh)
324  err_fhcnumubar_ih = math.sqrt(err_fhcnumubar_ih)
325  err_fhcnumubar_nh = math.sqrt(err_fhcnumubar_nh)
326  err_fhcnuenuebar_ih = math.sqrt(err_fhcnuenuebar_ih)
327  err_fhcnuenuebar_nh = math.sqrt(err_fhcnuenuebar_nh)
328 
329  err_rhcnumu_ih = math.sqrt(err_rhcnumu_ih)
330  err_rhcnumu_nh = math.sqrt(err_rhcnumu_nh)
331  err_rhcnumubar_ih = math.sqrt(err_rhcnumubar_ih)
332  err_rhcnumubar_nh = math.sqrt(err_rhcnumubar_nh)
333  err_rhcnuenuebar_ih = math.sqrt(err_rhcnuenuebar_ih)
334  err_rhcnuenuebar_nh = math.sqrt(err_rhcnuenuebar_nh)
335 
336 
337  # Get nominal 75% CP sensitivities
338  if(quantity=="cp_75thpercentile"):
339  quantity = "cp"
340  if(quantity=="mh_minimum"):
341  quantity = "mh"
342  FMC_SENSIT = "/dune/data/users/lblpwg_tools/FastMC_Data/outputs/ljf26/Sensitivity_Plots"
343  nh_cp_file = ROOT.TFile(FMC_SENSIT+"/ProtonP120GeV_nh_"+quantity+"_histos.root")
344  ih_cp_file = ROOT.TFile(FMC_SENSIT+"/ProtonP120GeV_ih_"+quantity+"_histos.root")
345 
346  nh_sensitivity = GetPercentile(nh_cp_file.Get("h2"),75)
347  ih_sensitivity = GetPercentile(ih_cp_file.Get("h2"),75)
348 
349  # Calculate estimated sensitivity given changes in flux
350  fitness_types = ["standard","fhcnumu","fhcnumubar","fhcnuenuebar","rhcnumu","rhcnumubar","rhcnuenuebar","standard_nh","standard_ih"]
351  nh_fitness = [0,0,0,0,0,0,0]
352  ih_fitness = [0,0,0,0,0,0,0]
353  nh_fiterror = [0,0,0,0,0,0,0]
354  ih_fiterror = [0,0,0,0,0,0,0]
355  # standard
356  nh_fitness[0] = nh_sensitivity+ tot_fhcnumu_nh + tot_fhcnumubar_nh + tot_rhcnumu_nh + tot_rhcnumubar_nh + tot_fhcnuenuebar_nh + tot_rhcnuenuebar_nh;
357  ih_fitness[0] = ih_sensitivity+ tot_fhcnumu_ih + tot_fhcnumubar_ih + tot_rhcnumu_ih + tot_rhcnumubar_ih + tot_fhcnuenuebar_ih + tot_rhcnuenuebar_ih;
358  nh_fiterror[0] = ROOT.TMath.Sqrt(ROOT.TMath.Power(err_fhcnumu_nh,2)+ROOT.TMath.Power(err_fhcnumubar_nh,2)+ROOT.TMath.Power(err_rhcnumu_nh,2)+ROOT.TMath.Power(err_rhcnumubar_nh,2)+ROOT.TMath.Power(err_fhcnuenuebar_nh,2)+ROOT.TMath.Power(err_rhcnuenuebar_nh,2))
359  ih_fiterror[0] = ROOT.TMath.Sqrt(ROOT.TMath.Power(err_fhcnumu_ih,2)+ROOT.TMath.Power(err_fhcnumubar_ih,2)+ROOT.TMath.Power(err_rhcnumu_ih,2)+ROOT.TMath.Power(err_rhcnumubar_ih,2)+ROOT.TMath.Power(err_fhcnuenuebar_ih,2)+ROOT.TMath.Power(err_rhcnuenuebar_ih,2))
360 
361  nh_fitness[1] = tot_fhcnumu_nh
362  ih_fitness[1] = tot_fhcnumu_ih
363  nh_fiterror[1] = err_fhcnumu_nh
364  ih_fiterror[1] = err_fhcnumu_ih
365 
366  nh_fitness[2] = tot_fhcnumubar_nh
367  ih_fitness[2] = tot_fhcnumubar_ih
368  nh_fiterror[2] = err_fhcnumubar_nh
369  ih_fiterror[2] = err_fhcnumubar_ih
370 
371  nh_fitness[3] = tot_fhcnuenuebar_nh
372  ih_fitness[3] = tot_fhcnuenuebar_ih
373  nh_fiterror[3] = err_fhcnuenuebar_nh
374  ih_fiterror[3] = err_fhcnuenuebar_ih
375 
376  nh_fitness[4] = tot_rhcnumu_nh
377  ih_fitness[4] = tot_rhcnumu_ih
378  nh_fiterror[4] = err_rhcnumu_nh
379  ih_fiterror[4] = err_rhcnumu_ih
380 
381  nh_fitness[5] = tot_rhcnumubar_nh
382  ih_fitness[5] = tot_rhcnumubar_ih
383  nh_fiterror[5] = err_rhcnumubar_nh
384  ih_fiterror[5] = err_rhcnumubar_ih
385 
386  nh_fitness[6] = tot_rhcnuenuebar_nh
387  ih_fitness[6] = tot_rhcnuenuebar_ih
388  nh_fiterror[6] = err_rhcnuenuebar_nh
389  ih_fiterror[6] = err_rhcnuenuebar_ih
390 
391  fitness = []
392  fiterror = []
393  for i in range(len(nh_fitness)):
394  fitness.append((nh_fitness[i]+ih_fitness[i])/2.0)
395  fiterror.append(ROOT.TMath.Sqrt(nh_fiterror[i]*nh_fiterror[i]+ih_fiterror[i]*ih_fiterror[i])/2.0)
396 
397  fitness.append(nh_fitness[0])
398  fitness.append(ih_fitness[0])
399  fiterror.append(nh_fiterror[0])
400  fiterror.append(ih_fiterror[0])
401 
402  fhcnumu_nh_weightfile.Close()
403  fhcnumubar_nh_weightfile.Close()
404  fhcnuenuebar_nh_weightfile.Close()
405  rhcnumu_nh_weightfile.Close()
406  rhcnumubar_nh_weightfile.Close()
407  rhcnuenuebar_nh_weightfile.Close()
408  rhcnumu_nh_weightfile.Close()
409  rhcnumubar_nh_weightfile.Close()
410  rhcnuenuebar_nh_weightfile.Close()
411  rhcnumu_nh_weightfile.Close()
412  rhcnumubar_nh_weightfile.Close()
413  rhcnuenuebar_nh_weightfile.Close()
414 
415  nomfhctfile.Close()
416  nomrhctfile.Close()
417  varfhctfile.Close()
418  varrhctfile.Close()
419 
420  return [fitness_types,fitness,fiterror]
421 
422 def AddPlotLabel(label,x,y,size=0.05,color=1,font=62,align=22,angle=0):
423 
424  latex = ROOT.TLatex( x, y, label );
425  latex.SetNDC();
426  latex.SetTextSize(4);
427  latex.SetTextColor(color);
428  latex.SetTextFont(font);
429  latex.SetTextAlign(align);
430  latex.SetTextAngle(angle);
431  latex.DrawClone();
432 
433 
434 def MakeMacro(optimization,macro_template,newmac_path,new_values,parameter_prefix=""):
435 
436  oldmac = open(macro_template)
437  newmac = open(newmac_path, 'w')
438 
439  # Some preliminar calculations for
440  # Optimizations with three configurable horns
441 
442  if optimization.optname in ["CP_run11","CP_run15","CP_run16","CP_run17","CP_run18","CP_run19","CP_run20","CP_run21","CP_run22","CP_run23","CP_run24","CP_run25","CP_run26","CP_run27","CP_run28","CP_run29","CP_run30","CP_run31","LE_run2"]:
443  HornBF1 = new_values[optimization.getParameterIndex(parameter_prefix+"HornBF1")];
444  HornBF2 = new_values[optimization.getParameterIndex(parameter_prefix+"HornBF2")];
445  HornBF3 = new_values[optimization.getParameterIndex(parameter_prefix+"HornBF3")];
446  HornBF4 = new_values[optimization.getParameterIndex(parameter_prefix+"HornBF4")];
447  HornBF5 = new_values[optimization.getParameterIndex(parameter_prefix+"HornBF5")];
448 
449  Btot = HornBF1+HornBF2+HornBF3+HornBF4+HornBF5
450  HornBF1 = HornBF1 / Btot
451  HornBF2 = HornBF2 / Btot
452  HornBF3 = HornBF3 / Btot
453  HornBF4 = HornBF4 / Btot
454  HornBF5 = HornBF5 / Btot
455 
456  new_values[optimization.getParameterIndex(parameter_prefix+"HornBF1")] = HornBF1
457  new_values[optimization.getParameterIndex(parameter_prefix+"HornBF2")] = HornBF2
458  new_values[optimization.getParameterIndex(parameter_prefix+"HornBF3")] = HornBF3
459  new_values[optimization.getParameterIndex(parameter_prefix+"HornBF4")] = HornBF4
460  new_values[optimization.getParameterIndex(parameter_prefix+"HornBF5")] = HornBF5
461 
462  HornCF1 = new_values[optimization.getParameterIndex(parameter_prefix+"HornCF1")];
463  HornCF2 = new_values[optimization.getParameterIndex(parameter_prefix+"HornCF2")];
464  HornCF3 = new_values[optimization.getParameterIndex(parameter_prefix+"HornCF3")];
465  HornCF4 = new_values[optimization.getParameterIndex(parameter_prefix+"HornCF4")];
466  HornCF5 = new_values[optimization.getParameterIndex(parameter_prefix+"HornCF5")];
467 
468  Ctot = HornCF1+HornCF2+HornCF3+HornCF4+HornCF5
469  HornCF1 = HornCF1 / Ctot
470  HornCF2 = HornCF2 / Ctot
471  HornCF3 = HornCF3 / Ctot
472  HornCF4 = HornCF4 / Ctot
473  HornCF5 = HornCF5 / Ctot
474 
475  new_values[optimization.getParameterIndex(parameter_prefix+"HornCF1")] = HornCF1
476  new_values[optimization.getParameterIndex(parameter_prefix+"HornCF2")] = HornCF2
477  new_values[optimization.getParameterIndex(parameter_prefix+"HornCF3")] = HornCF3
478  new_values[optimization.getParameterIndex(parameter_prefix+"HornCF4")] = HornCF4
479  new_values[optimization.getParameterIndex(parameter_prefix+"HornCF5")] = HornCF5
480 
481 
482  # Start looping over old macro lines
483  for s in oldmac.xreadlines():
484 
485  if s.find("/LBNE/det/seHornCurrent")>=0:
486  if not ("HornCurrent" in optimization.parameter_names or
487  "FHCHornCurrent" in optimization.parameter_names):
488  newmac.write(s)
489  for j in range(0,len(optimization.parameter_names)):
490  if optimization.rhc_parameters_float_separate:
491  if (optimization.parameter_names[j].startswith("FHC") and macro_num==1):
492  continue
493  if (optimization.parameter_names[j].startswith("RHC") and macro_num==0):
494  continue
495 
496  parameter = optimization.parameter_names[j]
497  macro_set_stage = optimization.macro_set_stages[j]
498 
499  if macro_set_stage == "PreInit":
500 
501  value = new_values[j]
502 
503  parameter = optimization.parameter_names[j]
504 
505 # if("MultiSphereTargetRadius" in parameter):
506 # newmac.write("/LBNE/det/UseMultiSphereTarget True\n")
507 # newmac.write("/LBNE/det/TargetMaterial Beryllium\n")
508  if("SimpleTargetLength" in parameter):
509  newmac.write("/LBNE/det/UseSimpleCylindricalTarget True\n")
510 
511  newmac.write(optimization.macro_commands[j]+" "+str(value)+" "+optimization.parameter_units[j]+"\n")
512 
513  # setup special horn shape
514 
515 
516  if optimization.optname in ["CP_run5","CP_run7","CP_run8","CP_run13"]:
517  Horn1OCRadius = new_values[optimization.getParameterIndex(parameter_prefix+"Horn1OCRadius")];
518  Horn1ICRadius1 = new_values[optimization.getParameterIndex(parameter_prefix+"Horn1ICRadius1")];
519  Horn1ICRadius2 = new_values[optimization.getParameterIndex(parameter_prefix+"Horn1ICRadius2")];
520  Horn1ICRadius3 = new_values[optimization.getParameterIndex(parameter_prefix+"Horn1ICRadius3")];
521  Horn1ICRadius4 = new_values[optimization.getParameterIndex(parameter_prefix+"Horn1ICRadius4")];
522  Horn1ICLength1 = new_values[optimization.getParameterIndex(parameter_prefix+"Horn1ICLength1")];
523  Horn1ICLength2 = new_values[optimization.getParameterIndex(parameter_prefix+"Horn1ICLength2")];
524  Horn1ICLength3 = new_values[optimization.getParameterIndex(parameter_prefix+"Horn1ICLength3")];
525  Horn1ICLength4 = new_values[optimization.getParameterIndex(parameter_prefix+"Horn1ICLength4")];
526  Horn1ICLength5 = new_values[optimization.getParameterIndex(parameter_prefix+"Horn1ICLength5")];
527  Horn1ICLength6 = new_values[optimization.getParameterIndex(parameter_prefix+"Horn1ICLength6")];
528  Horn1ICLength7 = new_values[optimization.getParameterIndex(parameter_prefix+"Horn1ICLength7")];
529  newmac.write("/LBNE/det/UseHorn1Polycone True \n")
530  newmac.write("/LBNE/det/NumInnerPtsHorn1Polycone 10 \n")
531  newmac.write("/LBNE/det/Horn1PolyPt0RinThickZ "+str(Horn1OCRadius)+" 2.000 0.00 mm\n");
532  newmac.write("/LBNE/det/Horn1PolyPt1RinThickZ "+str(Horn1ICRadius1)+" 2.000 1.00 mm\n");
533  newmac.write("/LBNE/det/Horn1PolyPt2RinThickZ "+str(Horn1ICRadius1)+" 2.000 "+str(1.00+Horn1ICLength1)+" mm\n");
534  newmac.write("/LBNE/det/Horn1PolyPt3RinThickZ "+str(Horn1ICRadius2)+" 2.000 "+str(1.00+Horn1ICLength1+Horn1ICLength2)+" mm\n")
535  newmac.write("/LBNE/det/Horn1PolyPt4RinThickZ "+str(Horn1ICRadius2)+" 2.000 "+str(1.00+Horn1ICLength1+Horn1ICLength2+Horn1ICLength3)+" mm\n")
536  newmac.write("/LBNE/det/Horn1PolyPt5RinThickZ "+str(Horn1ICRadius3)+" 2.000 "+str(1.00+Horn1ICLength1+Horn1ICLength2+Horn1ICLength3+Horn1ICLength4)+" mm\n")
537  newmac.write("/LBNE/det/Horn1PolyPt6RinThickZ "+str(Horn1ICRadius3)+" 2.000 "+str(1.00+Horn1ICLength1+Horn1ICLength2+Horn1ICLength3+Horn1ICLength4+Horn1ICLength5)+" mm\n")
538  newmac.write("/LBNE/det/Horn1PolyPt6RinThickZ "+str(Horn1ICRadius4)+" 2.000 "+str(1.00+Horn1ICLength1+Horn1ICLength2+Horn1ICLength3+Horn1ICLength4+Horn1ICLength5+Horn1ICLength6)+" mm\n")
539  newmac.write("/LBNE/det/Horn1PolyPt8RinThickZ "+str(Horn1ICRadius4)+" 2.000 "+str(1.00+Horn1ICLength1+Horn1ICLength2+Horn1ICLength3+Horn1ICLength4+Horn1ICLength5+Horn1ICLength6+Horn1ICLength7)+" mm\n")
540  newmac.write("/LBNE/det/Horn1PolyPt9RinThickZ "+str(Horn1OCRadius)+" 2.000 "+str(2.00+Horn1ICLength1+Horn1ICLength2+Horn1ICLength3+Horn1ICLength4+Horn1ICLength5+Horn1ICLength6+Horn1ICLength7)+" mm\n")
541  newmac.write("/LBNE/det/Horn1PolyOuterRadius "+str(Horn1OCRadius)+" mm\n")
542 
543 
544  elif optimization.optname in ["CP_run9","CP_run10","CP_run11","CP_run14","CP_run15","CP_run16","CP_run17","CP_run18","CP_run19","CP_run20","CP_run21","CP_run22","CP_run23","CP_run24","CP_run25","CP_run26","CP_run27","CP_run28","CP_run29","CP_run30","CP_run31","LE_run2"]:
545  HornAThick1 = "2.0"
546  HornAThick2 = "2.0"
547  HornBThick1 = "2.0"
548  HornBThick2 = "2.0"
549  HornBThick3 = "2.0"
550  HornCThick1 = "4.0"
551  HornCThick2 = "4.0"
552  HornCThick3 = "4.0"
553  HornEndcapThick = "2.0"
554 
555  if optimization.optname in ["CP_run22","CP_run25"]:
556 
557  HornAThick1 = "2.5"
558  HornAThick2 = "2.0"
559  HornBThick1 = "3.0"
560  HornBThick2 = "4.0"
561  HornBThick3 = "3.0"
562  HornCThick1 = "4.0"
563  HornCThick2 = "4.0"
564  HornCThick3 = "4.0"
565  HornEndcapThick = "8.0"
566 
567  HornARadiusOC = new_values[optimization.getParameterIndex(parameter_prefix+"HornARadiusOC")];
568  HornARadius1 = new_values[optimization.getParameterIndex(parameter_prefix+"HornARadius1")];
569  if(optimization.optname not in ["CP_run27","CP_run29","CP_run30"]):
570  HornARadius2 = new_values[optimization.getParameterIndex(parameter_prefix+"HornARadius2")];
571  HornALength = new_values[optimization.getParameterIndex(parameter_prefix+"HornALength")];
572  if(optimization.optname not in ["CP_run27","CP_run29","CP_run30","CP_run31"]):
573  HornAF1 = new_values[optimization.getParameterIndex(parameter_prefix+"HornAF1")];
574  HornBRadiusOC = new_values[optimization.getParameterIndex(parameter_prefix+"HornBRadiusOC")];
575  HornBRadius1 = new_values[optimization.getParameterIndex(parameter_prefix+"HornBRadius1")];
576  HornBRadius2 = new_values[optimization.getParameterIndex(parameter_prefix+"HornBRadius2")];
577  HornBRadius3 = new_values[optimization.getParameterIndex(parameter_prefix+"HornBRadius3")];
578  HornBLength = new_values[optimization.getParameterIndex(parameter_prefix+"HornBLength")];
579  HornBF1 = new_values[optimization.getParameterIndex(parameter_prefix+"HornBF1")];
580  HornBF2 = new_values[optimization.getParameterIndex(parameter_prefix+"HornBF2")];
581  HornBF3 = new_values[optimization.getParameterIndex(parameter_prefix+"HornBF3")];
582  HornBF4 = new_values[optimization.getParameterIndex(parameter_prefix+"HornBF4")];
583  HornBF5 = new_values[optimization.getParameterIndex(parameter_prefix+"HornBF5")];
584 
585 
586  HornCRadiusOC = new_values[optimization.getParameterIndex(parameter_prefix+"HornCRadiusOC")];
587  HornCRadius1 = new_values[optimization.getParameterIndex(parameter_prefix+"HornCRadius1")];
588  HornCRadius2 = new_values[optimization.getParameterIndex(parameter_prefix+"HornCRadius2")];
589  HornCRadius3 = new_values[optimization.getParameterIndex(parameter_prefix+"HornCRadius3")];
590  HornCLength = new_values[optimization.getParameterIndex(parameter_prefix+"HornCLength")];
591  HornCF1 = new_values[optimization.getParameterIndex(parameter_prefix+"HornCF1")];
592  HornCF2 = new_values[optimization.getParameterIndex(parameter_prefix+"HornCF2")];
593  HornCF3 = new_values[optimization.getParameterIndex(parameter_prefix+"HornCF3")];
594  HornCF4 = new_values[optimization.getParameterIndex(parameter_prefix+"HornCF4")];
595  HornCF5 = new_values[optimization.getParameterIndex(parameter_prefix+"HornCF5")];
596 
597  HornBLongPosition = new_values[optimization.getParameterIndex(parameter_prefix+"HornBLongPosition")]
598  HornCLongPosition = new_values[optimization.getParameterIndex(parameter_prefix+"HornCLongPosition")]
599 
600  # CP_run22 has varying conductor thicknesses, so
601  # it has longer transition segments than
602  # than other optimizations, and has to be
603  # configured separately
604  if optimization.optname in ["CP_run22","CP_run25"]:
605 
606  newmac.write("/LBNE/det/NumberOfHornsPolycone 3 \n")
607  newmac.write("/LBNE/det/Horn2PolyZStartPos "+str(HornBLongPosition)+" \n")
608  newmac.write("/LBNE/det/Horn3PolyZStartPos "+str(HornCLongPosition)+" \n")
609  newmac.write("/LBNE/det/NumInnerPtsHorn2Polycone 12 \n")
610  newmac.write("/LBNE/det/NumInnerPtsHorn3Polycone 12 \n")
611 
612  newmac.write("/LBNE/det/Horn1PolyPt0RinThickZ "+str(HornARadiusOC-15)+" "+HornEndcapThick+" 0.00 mm\n");
613  newmac.write("/LBNE/det/Horn1PolyPt1RinThickZ "+str(HornARadius1)+" "+HornEndcapThick+" 20.00 mm\n");
614  newmac.write("/LBNE/det/Horn1PolyPt2RinThickZ "+str(HornARadius1)+" "+HornAThick1+" 40.00 mm\n");
615  newmac.write("/LBNE/det/Horn1PolyPt3RinThickZ "+str(HornARadius1)+" "+HornAThick1+" "+str(40.00+HornALength*HornAF1)+" mm\n");
616  newmac.write("/LBNE/det/Horn1PolyPt4RinThickZ "+str(HornARadius1)+" "+HornAThick2+" "+str(41.00+HornALength*HornAF1)+" mm\n");
617  newmac.write("/LBNE/det/Horn1PolyPt5RinThickZ "+str(HornARadius2)+" "+HornAThick2+" "+str(41.00+HornALength)+" mm\n")
618  newmac.write("/LBNE/det/Horn1PolyPt6RinThickZ "+str(HornARadius2)+" "+HornEndcapThick+" "+str(61.00+HornALength)+" mm\n")
619  newmac.write("/LBNE/det/Horn1PolyPt7RinThickZ "+str(HornARadiusOC-15)+" "+HornEndcapThick+" "+str(81.00+HornALength)+" mm\n")
620  newmac.write("/LBNE/det/Horn1PolyOuterRadius "+str(HornARadiusOC)+" mm\n")
621 
622  newmac.write("/LBNE/det/Horn2PolyPt0RinThickZ "+str(HornBRadiusOC-15)+" "+HornEndcapThick+" 0.00 mm\n");
623  newmac.write("/LBNE/det/Horn2PolyPt1RinThickZ "+str(HornBRadius1)+" "+HornEndcapThick+" 20.00 mm\n");
624  newmac.write("/LBNE/det/Horn2PolyPt2RinThickZ "+str(HornBRadius1)+" "+HornBThick1+" 40.00 mm\n");
625  newmac.write("/LBNE/det/Horn2PolyPt3RinThickZ "+str(HornBRadius1)+" "+HornBThick1+" "+str(40.00+HornBLength*HornBF1)+" mm\n");
626  newmac.write("/LBNE/det/Horn2PolyPt4RinThickZ "+str(HornBRadius2)+" "+HornBThick1+" "+str(40.00+HornBLength*(HornBF1+HornBF2))+" mm\n")
627  newmac.write("/LBNE/det/Horn2PolyPt5RinThickZ "+str(HornBRadius2)+" "+HornBThick2+" "+str(41.00+HornBLength*(HornBF1+HornBF2))+" mm\n")
628  newmac.write("/LBNE/det/Horn2PolyPt6RinThickZ "+str(HornBRadius2)+" "+HornBThick2+" "+str(41.00+HornBLength*(HornBF1+HornBF2+HornBF3))+" mm\n")
629  newmac.write("/LBNE/det/Horn2PolyPt7RinThickZ "+str(HornBRadius2)+" "+HornBThick3+" "+str(42.00+HornBLength*(HornBF1+HornBF2+HornBF3))+" mm\n")
630  newmac.write("/LBNE/det/Horn2PolyPt8RinThickZ "+str(HornBRadius3)+" "+HornBThick3+" "+str(42.00+HornBLength*(HornBF1+HornBF2+HornBF3+HornBF4))+" mm\n")
631  newmac.write("/LBNE/det/Horn2PolyPt9RinThickZ "+str(HornBRadius3)+" "+HornBThick3+" "+str(42.00+HornBLength)+" mm\n")
632  newmac.write("/LBNE/det/Horn2PolyPt10RinThickZ "+str(HornBRadius3)+" "+HornEndcapThick+" "+str(62.00+HornBLength)+" mm\n")
633  newmac.write("/LBNE/det/Horn2PolyPt11RinThickZ "+str(HornBRadiusOC-15)+" "+HornEndcapThick+" "+str(82.00+HornBLength)+" mm\n")
634  newmac.write("/LBNE/det/Horn2PolyOuterRadius "+str(HornBRadiusOC)+" mm\n")
635 
636  newmac.write("/LBNE/det/Horn3PolyPt0RinThickZ "+str(HornCRadiusOC-15)+" "+HornEndcapThick+" 0.00 mm\n");
637  newmac.write("/LBNE/det/Horn3PolyPt1RinThickZ "+str(HornCRadius1)+" "+HornEndcapThick+" 20.00 mm\n");
638  newmac.write("/LBNE/det/Horn3PolyPt2RinThickZ "+str(HornCRadius1)+" "+HornCThick1+" 40.00 mm\n");
639  newmac.write("/LBNE/det/Horn3PolyPt3RinThickZ "+str(HornCRadius1)+" "+HornCThick1+" "+str(40.00+HornCLength*HornCF1)+" mm\n");
640  newmac.write("/LBNE/det/Horn3PolyPt4RinThickZ "+str(HornCRadius2)+" "+HornCThick1+" "+str(40.00+HornCLength*(HornCF1+HornCF2))+" mm\n")
641  newmac.write("/LBNE/det/Horn3PolyPt5RinThickZ "+str(HornCRadius2)+" "+HornCThick2+" "+str(41.00+HornCLength*(HornCF1+HornCF2))+" mm\n")
642  newmac.write("/LBNE/det/Horn3PolyPt6RinThickZ "+str(HornCRadius2)+" "+HornCThick2+" "+str(41.00+HornCLength*(HornCF1+HornCF2+HornCF3))+" mm\n")
643  newmac.write("/LBNE/det/Horn3PolyPt7RinThickZ "+str(HornCRadius2)+" "+HornCThick3+" "+str(42.00+HornCLength*(HornCF1+HornCF2+HornCF3))+" mm\n")
644  newmac.write("/LBNE/det/Horn3PolyPt8RinThickZ "+str(HornCRadius3)+" "+HornCThick3+" "+str(42.00+HornCLength*(HornCF1+HornCF2+HornCF3+HornCF4))+" mm\n")
645  newmac.write("/LBNE/det/Horn3PolyPt9RinThickZ "+str(HornCRadius3)+" "+HornCThick3+" "+str(42.00+HornCLength)+" mm\n")
646  newmac.write("/LBNE/det/Horn3PolyPt10RinThickZ "+str(HornCRadius3)+" "+HornEndcapThick+" "+str(62.00+HornCLength)+" mm\n")
647  newmac.write("/LBNE/det/Horn3PolyPt11RinThickZ "+str(HornCRadiusOC-15)+" "+HornEndcapThick+" "+str(82.00+HornCLength)+" mm\n")
648  newmac.write("/LBNE/det/Horn3PolyOuterRadius "+str(HornCRadiusOC)+" mm\n")
649  else:
650  newmac.write("/LBNE/det/NumberOfHornsPolycone 3 \n")
651  newmac.write("/LBNE/det/Horn2PolyZStartPos "+str(HornBLongPosition)+" \n")
652  newmac.write("/LBNE/det/Horn3PolyZStartPos "+str(HornCLongPosition)+" \n")
653  if(optimization.optname not in ["CP_run27","CP_run29","CP_run30","CP_run31"]):
654  newmac.write("/LBNE/det/NumInnerPtsHorn1Polycone 5 \n")
655  else:
656  newmac.write("/LBNE/det/NumInnerPtsHorn1Polycone 4 \n")
657  newmac.write("/LBNE/det/NumInnerPtsHorn2Polycone 8 \n")
658  newmac.write("/LBNE/det/NumInnerPtsHorn3Polycone 8 \n")
659  if(optimization.optname not in ["cp_run27","CP_run29","CP_run30","CP_run31"]):
660  newmac.write("/LBNE/det/Horn1PolyPt0RinThickZ "+str(HornARadiusOC-15)+" 2.000 0.00 mm\n");
661  newmac.write("/LBNE/det/Horn1PolyPt1RinThickZ "+str(HornARadius1)+" 2.000 1.00 mm\n");
662  newmac.write("/LBNE/det/Horn1PolyPt2RinThickZ "+str(HornARadius1)+" 2.000 "+str(1.00+HornALength*HornAF1)+" mm\n");
663  newmac.write("/LBNE/det/Horn1PolyPt3RinThickZ "+str(HornARadius2)+" 2.000 "+str(1.00+HornALength)+" mm\n")
664  newmac.write("/LBNE/det/Horn1PolyPt4RinThickZ "+str(HornARadiusOC-15)+" 2.000 "+str(2.00+HornALength)+" mm\n")
665  elif optimization.optname not in ["CP_run31"]:
666  newmac.write("/LBNE/det/Horn1PolyPt0RinThickZ "+str(HornARadiusOC-15)+" 2.000 0.00 mm\n");
667  newmac.write("/LBNE/det/Horn1PolyPt1RinThickZ "+str(HornARadius1)+" 2.000 1.00 mm\n");
668  newmac.write("/LBNE/det/Horn1PolyPt2RinThickZ "+str(HornARadius1)+" 2.000 "+str(1.00+HornALength)+" mm\n");
669  newmac.write("/LBNE/det/Horn1PolyPt3RinThickZ "+str(HornARadiusOC-15)+" 2.000 "+str(2.00+HornALength)+" mm\n")
670  else:
671  newmac.write("/LBNE/det/Horn1PolyPt0RinThickZ "+str(HornARadiusOC-15)+" 2.000 0.00 mm\n");
672  newmac.write("/LBNE/det/Horn1PolyPt1RinThickZ "+str(HornARadius1)+" 2.000 1.00 mm\n");
673  newmac.write("/LBNE/det/Horn1PolyPt2RinThickZ "+str(HornARadius2)+" 2.000 "+str(1.00+HornALength)+" mm\n");
674  newmac.write("/LBNE/det/Horn1PolyPt3RinThickZ "+str(HornARadiusOC-15)+" 2.000 "+str(2.00+HornALength)+" mm\n")
675 
676  newmac.write("/LBNE/det/Horn1PolyOuterRadius "+str(HornARadiusOC)+" mm\n")
677 
678  newmac.write("/LBNE/det/Horn2PolyPt0RinThickZ "+str(HornBRadiusOC-15)+" 2.000 0.00 mm\n");
679  newmac.write("/LBNE/det/Horn2PolyPt1RinThickZ "+str(HornBRadius1)+" 2.000 1.00 mm\n");
680  newmac.write("/LBNE/det/Horn2PolyPt2RinThickZ "+str(HornBRadius1)+" 2.000 "+str(1.00+HornBLength*HornBF1)+" mm\n");
681  newmac.write("/LBNE/det/Horn2PolyPt3RinThickZ "+str(HornBRadius2)+" 2.000 "+str(1.00+HornBLength*(HornBF1+HornBF2))+" mm\n")
682  newmac.write("/LBNE/det/Horn2PolyPt4RinThickZ "+str(HornBRadius2)+" 2.000 "+str(1.00+HornBLength*(HornBF1+HornBF2+HornBF3))+" mm\n")
683  newmac.write("/LBNE/det/Horn2PolyPt5RinThickZ "+str(HornBRadius3)+" 2.000 "+str(1.00+HornBLength*(HornBF1+HornBF2+HornBF3+HornBF4))+" mm\n")
684  newmac.write("/LBNE/det/Horn2PolyPt6RinThickZ "+str(HornBRadius3)+" 2.000 "+str(1.00+HornBLength)+" mm\n")
685  newmac.write("/LBNE/det/Horn2PolyPt7RinThickZ "+str(HornBRadiusOC-15)+" 2.000 "+str(2.00+HornBLength)+" mm\n")
686  newmac.write("/LBNE/det/Horn2PolyOuterRadius "+str(HornBRadiusOC)+" mm\n")
687 
688  newmac.write("/LBNE/det/Horn3PolyPt0RinThickZ "+str(HornCRadiusOC-15)+" 4.000 0.00 mm\n");
689  newmac.write("/LBNE/det/Horn3PolyPt1RinThickZ "+str(HornCRadius1)+" 4.000 1.00 mm\n");
690  newmac.write("/LBNE/det/Horn3PolyPt2RinThickZ "+str(HornCRadius1)+" 4.000 "+str(1.00+HornCLength*HornCF1)+" mm\n");
691  newmac.write("/LBNE/det/Horn3PolyPt3RinThickZ "+str(HornCRadius2)+" 4.000 "+str(1.00+HornCLength*(HornCF1+HornCF2))+" mm\n")
692  newmac.write("/LBNE/det/Horn3PolyPt4RinThickZ "+str(HornCRadius2)+" 4.000 "+str(1.00+HornCLength*(HornCF1+HornCF2+HornCF3))+" mm\n")
693  newmac.write("/LBNE/det/Horn3PolyPt5RinThickZ "+str(HornCRadius3)+" 4.000 "+str(1.00+HornCLength*(HornCF1+HornCF2+HornCF3+HornCF4))+" mm\n")
694  newmac.write("/LBNE/det/Horn3PolyPt6RinThickZ "+str(HornCRadius3)+" 4.000 "+str(1.00+HornCLength)+" mm\n")
695  newmac.write("/LBNE/det/Horn3PolyPt7RinThickZ "+str(HornCRadiusOC-15)+" 4.000 "+str(2.00+HornCLength)+" mm\n")
696  newmac.write("/LBNE/det/Horn3PolyOuterRadius "+str(HornCRadiusOC)+" mm\n")
697 
698  if optimization.optname=="CP_run19":
699  newmac.write("/LBNE/det/SetPolyconeHornParabolic 1 \n")
700  newmac.write("/LBNE/det/SetPolyconeHornParabolic 2 \n")
701  newmac.write("/LBNE/det/SetPolyconeHornParabolic 3 \n")
702  newmac.write("/LBNE/det/TargetLengthOutsideHorn 0 m\n")
703 
704 
705  elif s.find("/run/initialize")>=0:
706  newmac.write(s)
707  for j in range(0,len(optimization.parameter_names)):
708  parameter = optimization.parameter_names[j]
709 
710  if optimization.rhc_parameters_float_separate:
711  if (optimization.parameter_names[j].startswith("FHC") and
712  macro_num==1):
713  continue
714  if (optimization.parameter_names[j].startswith("RHC") and
715  macro_num==0):
716  continue
717 
718  macro_set_stage = optimization.macro_set_stages[j]
719  if(macro_set_stage == "Idle"):
720  value = new_values[j]
721  newmac.write(optimization.macro_commands[j]+" "+str(value)+" "+optimization.parameter_units[j]+"\n")
722  # Also scale beam width to match target fin size
723  if("GraphiteTargetFinWidth" in parameter):
724  # Scale nominal size (1.7 mm) based on
725  # change in fin width from nominal (.10 m)
726  if optimization.optname in ["CP_run5","CP_run7","CP_run8","CP_run13"]:
727  new_beam_size = 1.7 * new_values[j]/10
728  newmac.write("/LBNE/generator/beamSigmaX "+str(new_beam_size)+" mm \n");
729  newmac.write("/LBNE/generator/beamSigmaY "+str(new_beam_size)+" mm \n");
730 
731  if parameter == "BeamSigma":
732  temp = optimization.macro_commands[j].replace("SigmaX","SigmaY")
733  newmac.write(temp+" "+str(value)+" "+optimization.parameter_units[j]+"\n")
734 
735 
736  elif s.find("/LBNE/run/NEvents")>=0: # set to small number of events for now (because we're gonna test it interactively)
737  newmac.write("/LBNE/run/NEvents 10 \n")
738  else:
739  newmac.write(s);
740 
741  newmac.close();
742  oldmac.close();
743 
744 def ParametersOkay(optimization,i):
745 
746  if(optimization.optname=="CP_run5" or optimization.optname=="CP_run6" or optimization.optname=="CP_run7" or optimization.optname=="CP_run8" or optimization.optname=="CP_run13"):
747  horn1length = optimization.getParameterValue(parameter_prefix+"Horn1ICLength1",i)+optimization.getParameterValue(parameter_prefix+"Horn1ICLength2",i)+optimization.getParameterValue(parameter_prefix+"Horn1ICLength3",i)+optimization.getParameterValue(parameter_prefix+"Horn1ICLength4",i)+optimization.getParameterValue(parameter_prefix+"Horn1ICLength5",i)+optimization.getParameterValue(parameter_prefix+"Horn1ICLength6",i)+optimization.getParameterValue(parameter_prefix+"Horn1ICLength7",i)
748  horn2pos = optimization.getParameterValue(parameter_prefix+"Horn2LongPosition",i)*1000 # convert to mm
749  horn2length = 3.63 * optimization.getParameterValue(parameter_prefix+"Horn2LongRescale",i) * 1000 # convert to mm
750  if(horn2pos < horn1length):
751  print "Warning: Configuration not valid because horn2 position is less than horn1 length.. rethrowing"
752  return False
753 
754  if(optimization.optname=="CP_run9" or optimization.optname=="CP_run10"):
755  HornALength = optimization.getParameterValue("HornALength",i)
756  HornBLength = optimization.getParameterValue("HornBLength",i)
757  HornCLength = optimization.getParameterValue("HornCLength",i)
758  HornBLongPosition = optimization.getParameterValue("HornBLongPosition",i)
759  HornCLongPosition = optimization.getParameterValue("HornCLongPosition",i)
760  HornBF1 = optimization.getParameterValue("HornBF1",i);
761  HornBF2 = optimization.getParameterValue("HornBF2",i);
762  HornBF3 = optimization.getParameterValue("HornBF3",i);
763  HornBF4 = optimization.getParameterValue("HornBF4",i);
764  HornCF1 = optimization.getParameterValue("HornCF1",i);
765  HornCF2 = optimization.getParameterValue("HornCF2",i);
766  HornCF3 = optimization.getParameterValue("HornCF3",i);
767  HornCF4 = optimization.getParameterValue("HornCF4",i);
768 
769  if HornBF1+HornBF2+HornBF3+HornBF4>0.99:
770  print "Failing BF"
771  print HornBF1,HornBF2,HornBF3,HornBF4
772  return False
773  if HornCF1+HornCF2+HornCF3+HornCF4>0.99:
774  print HornCF1,HornCF2,HornCF3,HornCF4
775  print "Failing CF"
776  return False
777  if HornBLongPosition < HornALength:
778  print "Failing HornB Start"
779  return False
780  if HornCLongPosition < HornBLength + HornBLongPosition:
781  print "Failing HornC Start"
782  return False
783  if HornCLongPosition + HornCLength > 21000:
784  print "Failing HornC End"
785  return False
786 
787  if(optimization.optname in ["CP_run11","CP_run15","CP_run16","CP_run17","CP_run18","CP_run19","CP_run21","CP_run22","CP_run23","CP_run24","CP_run25","CP_run26","CP_run27","CP_run28","CP_run29","CP_run30","CP_run31","LE_run2"]):
788  HornALength = optimization.getParameterValue("HornALength",i)
789  HornBLength = optimization.getParameterValue("HornBLength",i)
790  HornCLength = optimization.getParameterValue("HornCLength",i)
791  HornBLongPosition = optimization.getParameterValue("HornBLongPosition",i)
792  HornCLongPosition = optimization.getParameterValue("HornCLongPosition",i)
793  if HornBLongPosition < HornALength:
794  print "Failing HornB Start"
795  return False
796  if HornCLongPosition < HornBLength + HornBLongPosition:
797  print "Failing HornC Start"
798  return False
799  if HornCLongPosition + HornCLength > 21000:
800  print "Failing HornC End"
801  return False
802  if(optimization.optname in ["CP_run20"]):
803  HornALength = optimization.getParameterValue("HornALength",i)
804  HornBLength = optimization.getParameterValue("HornBLength",i)
805  HornCLength = optimization.getParameterValue("HornCLength",i)
806  HornBLongPosition = optimization.getParameterValue("HornBLongPosition",i)
807  HornCLongPosition = optimization.getParameterValue("HornCLongPosition",i)
808  if HornBLongPosition < HornALength:
809  print "Failing HornB Start"
810  return False
811  if HornCLongPosition < HornBLength + HornBLongPosition:
812  print "Failing HornC Start"
813  return False
814 
815  if(optimization.optname=="CP_run12"):
816  HornALength = optimization.getParameterValue("HornAICLength1",i)+optimization.getParameterValue("HornAICLength2",i)
817  HornBLength = optimization.getParameterValue("HornBICLength1",i)+optimization.getParameterValue("HornBICLength2",i)+optimization.getParameterValue("HornBICLength3",i)+optimization.getParameterValue("HornBICLength4",i)+optimization.getParameterValue("HornBICLength5",i)
818  HornCLength = optimization.getParameterValue("HornCICLength1",i)+optimization.getParameterValue("HornCICLength2",i)+optimization.getParameterValue("HornCICLength3",i)+optimization.getParameterValue("HornCICLength4",i)+optimization.getParameterValue("HornCICLength5",i)
819  HornBLongPosition = optimization.getParameterValue("HornBLongPosition",i)
820  HornCLongPosition = optimization.getParameterValue("HornCLongPosition",i)
821  if HornBLongPosition < HornALength:
822  print "Failing HornB Start"
823  return False
824  if HornCLongPosition < HornBLength + HornBLongPosition:
825  print "Failing HornC Start"
826  return False
827  if HornCLongPosition + HornCLength > 21000:
828  print "Failing HornC End"
829  return False
830  return True
831 
OptimizationUtils.GetFractionAtMoreThanXSigma
def GetFractionAtMoreThanXSigma(histo, X)
Definition: OptimizationUtils.py:45
open
int open(const char *, int)
Opens a file descriptor.
OptimizationUtils.MakeMacro
def MakeMacro(optimization, macro_template, newmac_path, new_values, parameter_prefix="")
Definition: OptimizationUtils.py:434
OptimizationUtils.ComputeMean75PercentSensitivityAndError
def ComputeMean75PercentSensitivityAndError(varfhcfile, varrhcfile, fhcvarenergy=120.0, rhcvarenergy=120.0, antinufrac=0.5, quantity="cp")
Definition: OptimizationUtils.py:81
OptimizationUtils.ParametersOkay
def ParametersOkay(optimization, i)
Definition: OptimizationUtils.py:744
OptimizationUtils.GetPowerPOTScaleFactor
def GetPowerPOTScaleFactor(energy)
Definition: OptimizationUtils.py:59
keras_to_tensorflow.int
int
Definition: keras_to_tensorflow.py:69
OptimizationUtils.GetPercentile
def GetPercentile(histo, X)
Definition: OptimizationUtils.py:22
if
if(!yymsg) yymsg
OptimizationUtils.AddPlotLabel
def AddPlotLabel(label, x, y, size=0.05, color=1, font=62, align=22, angle=0)
Definition: OptimizationUtils.py:422
str
static QCString str
Definition: fortrancode.cpp:27098
OptimizationUtils.ConvertGraphToHisto
def ConvertGraphToHisto(new_bins, pGraph)
Definition: OptimizationUtils.py:3
Generated by   doxygen 1.8.11