makeFocusingCovariance.py
Go to the documentation of this file.
1 import sys,ROOT,math,array
2 
3 ROOT.gStyle.SetOptStat(0)
4 
5 small_bins = False
6 
7 #beam_config = "reference"
8 beam_config = "opt"
9 
10 def GetCorrelation(cov_matrix):
11  cor_matrix = cov_matrix.Clone("Focusing_Correlation")
12  for i in range(cor_matrix.GetNbinsX()):
13  for j in range(cor_matrix.GetNbinsX()):
14  if cov_matrix.GetBinContent(i+1,i+1)!=0 and cov_matrix.GetBinContent(j+1,j+1)!=0:
15  cor_matrix.SetBinContent(i+1,j+1,cov_matrix.GetBinContent(i+1,j+1)/ROOT.TMath.Sqrt(cov_matrix.GetBinContent(i+1,i+1))/ROOT.TMath.Sqrt(cov_matrix.GetBinContent(j+1,j+1)))
16  else:
17  if i==j: cor_matrix.SetBinContent(i+1,j+1,1)
18  else: cor_matrix.SetBinContent(i+1,j+1,0)
19  cor_matrix.SetMinimum(-1.0)
20  cor_matrix.SetMaximum(1.0)
21  return cor_matrix
22 
23 
24 # Parse input files
25 
26 shifts = []
27 cv = []
28 
29 
30 if beam_config == "opt":
31 
32  uncertainties = ["HornCurrent","WaterThickness","DecayPipeRadius","Horn1XShift","Horn1XTilt"]
33 
34  uncertainty_titles = ["Horn Current","Water Layer","Decay Pipe Radius","Horn 1 Trans. Offset","Horn 1 Tilt"]
35 
36  indir = "/dune/app/users/pmadigan/beam_sim/ToleranceStudy/"
37 
38  # I'm assuming that uncertainties are 100% correlated for nu mode and antinumode
39  files = ["table_ND_numu.csv","table_ND_numubar.csv","table_ND_nue.csv","table_ND_nuebar.csv","table_ND_numubar.csv","table_ND_numu.csv","table_ND_nuebar.csv","table_ND_nue.csv","table_FD_numu.csv","table_FD_numubar.csv","table_FD_nue.csv","table_FD_nuebar.csv","table_FD_numubar.csv","table_FD_numu.csv","table_FD_nuebar.csv","table_FD_nue.csv"]
40 
41  for uncertainty in uncertainties:
42  shifts.append([])
43 
44  for afile in files:
45  f_temp = open(indir+afile)
46  lines = f_temp.readlines()
47  f_temp.close()
48  for i in range(len(lines)):
49  if lines[i].split(",")[0]=="Nominal value":
50  temp = lines[i].split(",")[7:]
51  for i in range(len(temp)):
52  temp[i] = float(temp[i])
53  cv = cv + temp
54 
55  for unc_iter in range(len(uncertainties)):
56  line_iter = -1;
57  for i in range(len(lines)):
58  if lines[i].split(",")[0]==uncertainties[unc_iter]:
59  line_iter = i
60  if line_iter==-1:
61  print "ERROR: Could not find uncertainty",uncertainties[unc_iter],"in file",afile
62  continue
63  temp = lines[line_iter].split(",")[7:]
64  for i in range(len(temp)):
65  temp[i] = float(temp[i].strip())
66 
67  # Peter's plots are in 0.5 GeV bins from 0 to 10 GeV for numu
68  # or 1 GeV bins from 0 to 10 GeV for nue,numubar
69  # But for NDTF, we need bins of
70  # numu/numubar: [0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,7 ,8,12,16,20,40,100]
71  # nue/nuebar:[0,2,4,6,8,10,20,100]
72  scale = 1.0
73  # Dave Pushka thinks that the radius uncertainty is more like 4 cm
74  # instead of the 10 cm difference simulated, so scale that unc accordingly
75  if uncertainties[unc_iter] == "DecayPipeRadius": scale = 0.4
76 
77  temp_2 = []
78 
79  nbins = 7
80  if small_bins:
81  nbins = 19
82  if "numu" in afile:
83  nbins = 19
84 
85  for i in range(nbins):
86  if "bar" in afile:
87  temp_2 = nbins*[0]
88  elif "numu" in afile:
89 
90  if i<12:
91  temp_2.append(scale*temp[i]/100)
92  elif i==12:
93  new_shift = (temp[12]/100*cv[12]+temp[13]/100*cv[13])/(cv[12]+cv[13])
94 
95  temp_2.append(scale*new_shift)
96  elif i==13:
97  new_shift = (temp[14]/100*cv[15]+temp[14]/100*cv[15])/(cv[14]+cv[15])
98  temp_2.append(scale*new_shift)
99  elif i==14:
100  new_shift = (temp[16]/100*cv[16]+temp[17]/100*cv[17]+temp[18]/100*cv[18]+temp[19]/100*cv[19])/(cv[16]+cv[17]+cv[18]+cv[19])
101  temp_2.append(scale*new_shift)
102  else:
103  temp_2.append(0)
104  else:
105  if i<5:
106  new_shift = (temp[2*i]/100*cv[2*i]+temp[2*i+1]/100*cv[4*i+1])/(cv[2*i+0]+cv[2*i+1])
107  temp_2.append(scale*new_shift)
108  else:
109  temp_2.append(0)
110 
111  shifts[unc_iter] = shifts[unc_iter]+temp_2
112 
113  # Add a few things peter didn't estimate
114  temp_uncertainties = ["Target","BeamSigma","BeamOffsetX","TargetDensity"]
115  uncertainties = uncertainties + temp_uncertainties
116  uncertainty_titles += ["Target Long. Offset","Beam Size","Beam Position","Target Density"]
117  sigmas = ["0p5mm","1p7mm","0p5mm","1p8156gcm3"]
118 
119  indir = "/dune/data/users/ljf26/fluxfiles/g4lbne/v3r4p3/QGSP_BERT/"
120 
121  for uncertainty in temp_uncertainties:
122  unc_iter2 = temp_uncertainties.index(uncertainty)
123  shifts.append([])
124  detectors = ["ND","FD"]
125  modes = ["neutrino","antineutrino"]
126  flavors = ["numu","numubar","nue","nuebar"]
127  for detector in detectors:
128  for mode in modes:
129  for flavor in flavors:
130 
131  varied_file = ROOT.TFile("/dune/data/users/ljf26/fluxfiles/g4lbne/v3r4p3/QGSP_BERT/CP_run15_12388_80GeV_"+uncertainty+sigmas[unc_iter2]+"/"+mode+"/flux/histos_g4lbne_v3r4p3_QGSP_BERT_CP_run15_12388_80GeV_"+uncertainty+sigmas[unc_iter2]+"_"+mode+"_LBNE"+detector+"_fastmc.root")
132  varied_histo = varied_file.Get(flavor+"_flux")
133  varied_histo.SetDirectory(0);
134  varied_file.Close()
135 
136  nominal_file = ROOT.TFile("/dune/data/users/ljf26/fluxfiles/g4lbne/v3r4p2/QGSP_BERT/CP_run15_12388_80GeV/"+mode+"/flux/histos_g4lbne_v3r4p2_QGSP_BERT_CP_run15_12388_80GeV_"+mode+"_LBNE"+detector+"_fastmc.root")
137  nominal_histo = nominal_file.Get(flavor+"_flux")
138  nominal_histo.SetDirectory(0);
139  nominal_file.Close()
140 
141  varied_bins = [0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,7 ,8,12,16,20,40,100]
142  if flavor in ["nue","nuebar"]:
143  if not small_bins:
144  varied_bins = [0,2,4,6,8,10,20,100]
145  varied_histo = varied_histo.Rebin(len(varied_bins)-1,uncertainty+flavor,array.array('d',varied_bins))
146  nominal_histo = nominal_histo.Rebin(len(varied_bins)-1,uncertainty+flavor,array.array('d',varied_bins))
147 
148  varied_histo.Add(nominal_histo,-1.0)
149  varied_histo.Divide(nominal_histo)
150 
151  if(uncertainty=="TargetDensity" or (mode=="neutrino" and "bar" not in flavor) or (mode=="antineutrino" and "bar" in flavor)):
152  for bin_iter in range(0,varied_histo.GetNbinsX()):
153  shifts[unc_iter+1+unc_iter2].append(varied_histo.GetBinContent(bin_iter+1))
154  else:
155  # assume these uncertainties are zero
156  # for unfocused flavors
157  for bin_iter in range(0,varied_histo.GetNbinsX()):
158  shifts[unc_iter+1+unc_iter2].append(0)
159 
160  # Add flat POT accounting uncertainty
161  shifts.append([])
162  uncertainties.append("POT Counting")
163  uncertainty_titles.append("POT Counting")
164  for i in range(208):
165  shifts[unc_iter+1+len(temp_uncertainties)].append(0.02)
166 
167  # Add baffle scraping uncertainty
168  shifts.append([])
169  uncertainties.append("Baffle Scraping")
170  uncertainty_titles.append("Baffle Scraping")
171  detectors = ["ND","FD"]
172  modes = ["neutrino","antineutrino"]
173  flavors = ["numu","numubar","nue","nuebar"]
174 
175  for detector in detectors:
176  for mode in modes:
177  for flavor in flavors:
178 
179  baffle_scraping_file = ROOT.TFile("/dune/data/users/ljf26/fluxfiles/g4lbne/v3r4p3/QGSP_BERT/CP_run15_12388_80GeV_BaffleScraping/"+mode+"/flux/histos_g4lbne_v3r4p3_QGSP_BERT_CP_run15_12388_80GeV_BaffleScraping_"+mode+"_LBNE"+detector+"_fastmc.root")
180  baffle_scraping_histo = baffle_scraping_file.Get(flavor+"_flux")
181  baffle_scraping_histo.SetDirectory(0);
182  baffle_scraping_file.Close()
183 
184  nominal_file = ROOT.TFile("/dune/data/users/ljf26/fluxfiles/g4lbne/v3r4p2/QGSP_BERT/CP_run15_12388_80GeV/"+mode+"/flux/histos_g4lbne_v3r4p2_QGSP_BERT_CP_run15_12388_80GeV_"+mode+"_LBNE"+detector+"_fastmc.root")
185  nominal_histo = nominal_file.Get(flavor+"_flux")
186  nominal_histo.SetDirectory(0);
187  nominal_file.Close()
188 
189  baffle_bins = [0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,7 ,8,12,16,20,40,100]
190  if flavor in ["nue","nuebar"]:
191  if not small_bins:
192  baffle_bins = [0,2,4,6,8,10,20,100]
193  baffle_scraping_histo = baffle_scraping_histo.Rebin(len(baffle_bins)-1,"baffle_scraping_histo",array.array('d',baffle_bins))
194  nominal_histo = nominal_histo.Rebin(len(baffle_bins)-1,"nominal_histo",array.array('d',baffle_bins))
195 
196  baffle_scraping_histo.Scale(0.0025) # 0.25% baffle scraping unc
197  for bin_iter in range(0,baffle_scraping_histo.GetNbinsX()):
198  if nominal_histo.GetBinContent(bin_iter+1):
199  shifts[unc_iter+2+len(temp_uncertainties)].append(baffle_scraping_histo.GetBinContent(bin_iter+1)/nominal_histo.GetBinContent(bin_iter+1))
200  else:
201  shifts[unc_iter+2+len(temp_uncertainties)].append(0)
202 
203 
204 
205 
206 if beam_config == "reference":
207 
208  uncertainties = ["HornCurrent","WaterLayer","DecayPipeRadius","Horn1X","Horn1XTilt","Target","BeamSigma","BeamOffsetX","TargetDensity"]
209  uncertainty_titles = ["Horn Current","Water Layer","Decay Pipe Radius","Horn 1 Trans. Offset","Horn 1 Tilt", "Target Long. Offset","Beam Size","Beam Position","Target Density"]
210 
211  sigmas = ["228","0p5mm","1p9m","0p5mm","0p5mm","45p5cm","1p8mm","0p5mm","1p8156gcm3"]
212 
213  indir = "/dune/data/users/ljf26/fluxfiles/g4lbne/v3r4p3/QGSP_BERT/"
214 
215  for uncertainty in uncertainties:
216  unc_iter = uncertainties.index(uncertainty)
217  shifts.append([])
218  detectors = ["ND","FD"]
219  modes = ["neutrino","antineutrino"]
220  flavors = ["numu","numubar","nue","nuebar"]
221  for detector in detectors:
222  for mode in modes:
223  for flavor in flavors:
224 
225  varied_file = ROOT.TFile("/dune/data/users/ljf26/fluxfiles/g4lbne/v3r4p3/QGSP_BERT/Nominal_80GeV_"+uncertainty+sigmas[unc_iter]+"/"+mode+"/flux/histos_g4lbne_v3r4p3_QGSP_BERT_Nominal_80GeV_"+uncertainty+sigmas[unc_iter]+"_"+mode+"_LBNE"+detector+"_fastmc.root")
226  varied_histo = varied_file.Get(flavor+"_flux")
227  varied_histo.SetDirectory(0);
228  varied_file.Close()
229 
230  nominal_file = ROOT.TFile("/dune/data/users/ljf26/fluxfiles/g4lbne/v3r4p3/QGSP_BERT/Nominal_80GeV_AlignmentStudy/"+mode+"/flux/histos_g4lbne_v3r4p3_QGSP_BERT_Nominal_80GeV_AlignmentStudy_"+mode+"_LBNE"+detector+"_fastmc.root")
231  nominal_histo = nominal_file.Get(flavor+"_flux")
232  nominal_histo.SetDirectory(0);
233  nominal_file.Close()
234 
235  varied_bins = [0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,7 ,8,12,16,20,40,100]
236  if flavor in ["nue","nuebar"]:
237  if not small_bins:
238  varied_bins = [0,2,4,6,8,10,20,100]
239  varied_histo = varied_histo.Rebin(len(varied_bins)-1,uncertainty+flavor,array.array('d',varied_bins))
240  nominal_histo = nominal_histo.Rebin(len(varied_bins)-1,uncertainty+flavor,array.array('d',varied_bins))
241 
242  varied_histo.Add(nominal_histo,-1.0)
243  varied_histo.Divide(nominal_histo)
244 
245  scale = 1.0
246  if(uncertainty=="DecayPipeRadius"): scale = 0.4
247  if((mode=="neutrino" and "bar" not in flavor) or (mode=="antineutrino" and "bar" in flavor)):
248  for bin_iter in range(0,varied_histo.GetNbinsX()):
249  shifts[unc_iter].append(scale*varied_histo.GetBinContent(bin_iter+1))
250  else:
251  # assume these uncertainties are zero
252  # for unfocused flavors
253  for bin_iter in range(0,varied_histo.GetNbinsX()):
254  shifts[unc_iter].append(0)
255 
256 
257  # Add flat POT accounting uncertainty
258  shifts.append([])
259  uncertainties.append("POT Counting")
260  uncertainty_titles.append("POT Counting")
261  for i in range(208):
262  shifts[unc_iter+1].append(0.02)
263 
264  # Add baffle scraping uncertainty
265 
266  shifts.append([])
267  uncertainties.append("Baffle Scraping")
268  uncertainty_titles.append("Baffle Scraping")
269  detectors = ["ND","FD"]
270  modes = ["neutrino","antineutrino"]
271  flavors = ["numu","numubar","nue","nuebar"]
272  for detector in detectors:
273  for mode in modes:
274  for flavor in flavors:
275 
276  baffle_scraping_file = ROOT.TFile("/dune/data/users/ljf26/fluxfiles/g4lbne/v3r4p3/QGSP_BERT/Nominal_80GeV_BaffleScraping/"+mode+"/flux/histos_g4lbne_v3r4p3_QGSP_BERT_Nominal_80GeV_BaffleScraping_"+mode+"_LBNE"+detector+"_fastmc.root")
277  baffle_scraping_histo = baffle_scraping_file.Get(flavor+"_flux")
278  baffle_scraping_histo.SetDirectory(0);
279  baffle_scraping_file.Close()
280 
281  nominal_file = ROOT.TFile("/dune/data/users/ljf26/fluxfiles/g4lbne/v3r4p3/QGSP_BERT/Nominal_80GeV_AlignmentStudy/"+mode+"/flux/histos_g4lbne_v3r4p3_QGSP_BERT_Nominal_80GeV_AlignmentStudy_"+mode+"_LBNE"+detector+"_fastmc.root")
282  nominal_histo = nominal_file.Get(flavor+"_flux")
283  nominal_histo.SetDirectory(0);
284  nominal_file.Close()
285 
286  baffle_bins = [0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,7 ,8,12,16,20,40,100]
287  if flavor in ["nue","nuebar"]:
288  if not small_bins:
289  baffle_bins = [0,2,4,6,8,10,20,100]
290  baffle_scraping_histo = baffle_scraping_histo.Rebin(len(baffle_bins)-1,"baffle_scraping_histo",array.array('d',baffle_bins))
291  nominal_histo = nominal_histo.Rebin(len(baffle_bins)-1,"nominal_histo",array.array('d',baffle_bins))
292 
293  baffle_scraping_histo.Scale(0.0025) # 0.25% baffle scraping unc
294 
295  for bin_iter in range(0,baffle_scraping_histo.GetNbinsX()):
296  if nominal_histo.GetBinContent(bin_iter+1):
297  shifts[unc_iter+2].append(baffle_scraping_histo.GetBinContent(bin_iter+1)/nominal_histo.GetBinContent(bin_iter+1))
298  else:
299  shifts[unc_iter+2].append(0)
300 
301 nbins = 208
302 if small_bins: nbins = 304
303 error = ROOT.TH1D("Focusing_Error","Focusing_Error",nbins,0.5,nbins+0.5)
304 covariance = ROOT.TH2D("Focusing_Covariance","Focusing_Covariance",nbins,0.5,nbins+0.5,nbins,0.5,nbins+0.5)
305 for i in range(nbins):
306  for j in range(nbins):
307  for k in range(len(uncertainties)):
308  covariance.SetBinContent(i+1,j+1,covariance.GetBinContent(i+1,j+1)+shifts[k][i]*shifts[k][j])
309 
310 for i in range(nbins):
311  error.SetBinContent(i+1,ROOT.TMath.Sqrt(covariance.GetBinContent(i+1,i+1)))
312 
313 c1 = ROOT.TCanvas("c1","",1000,800)
314 error.Draw("")
315 c1.Print("focusing_error_"+beam_config+".eps")
316 c1.Print("focusing_error_"+beam_config+".png")
317 
318 ROOT.gStyle.SetPadRightMargin(0.15);
319 ROOT.gStyle.SetPadLeftMargin(0.15);
320 c1 = ROOT.TCanvas("c1","",1000,800)
321 covariance.SetMaximum(0.005)
322 covariance.SetMinimum(-.00030)
323 #covariance.SetMaximum(0.06)
324 #covariance.SetMinimum(-.005)
325 covariance.SetLabelSize(0.001,"X")
326 covariance.SetLabelSize(0.001,"Y")
327 covariance.Draw("COLZ")
328 
329 latex0 = ROOT.TLatex( 0.5, 0.93, "Flux Alignment Fractional Covariance Matrix" );
330 latex0.SetTextSize(0.04)
331 latex0.SetTextAlign(22)
332 latex0.SetNDC()
333 latex0.SetTextFont(42)
334 latex0.Draw()
335 
336 if(1>0):
337  latex1 = ROOT.TLatex( 0.24, 0.035, "Near Detector" );
338  latex1.SetTextSize(0.035);
339  latex1.SetNDC()
340  latex1.SetTextFont(42);
341  latex1.Draw()
342  latex2 = latex1.Clone()
343  latex2.SetText(0.61,0.035,"Far Detector")
344  latex2.Draw()
345  latex3 = latex1.Clone()
346  latex3.SetText(0.05,0.24,"Near Detector")
347  latex3.SetTextAngle(90)
348  latex3.Draw()
349  latex4 = latex1.Clone()
350  latex4.SetText(0.05,0.62,"Far Detector")
351  latex4.SetTextAngle(90)
352  latex4.Draw()
353 
354  latex5 = ROOT.TLatex(0.2,0.08,"#nu Mode")
355  latex5.SetTextSize(0.025);
356  latex5.SetNDC()
357  latex5.SetTextFont(42)
358  latex5.Draw()
359  latex6 = latex5.Clone()
360  latex6.SetText(0.38,0.08,"#bar{#nu} Mode")
361  latex6.Draw()
362  latex7 = latex5.Clone()
363  latex7.SetText(0.56,0.08,"#nu Mode")
364  latex7.Draw()
365  latex8 = latex5.Clone()
366  latex8.SetText(0.74,0.08,"#bar{#nu} Mode")
367  latex8.Draw()
368 
369  latex9 = ROOT.TLatex(0.085,0.21,"#nu Mode")
370  latex9.SetTextSize(0.025);
371  latex9.SetNDC()
372  latex9.SetTextFont(42)
373  latex9.SetTextAngle(90)
374  latex9.Draw()
375  latex10 = latex9.Clone()
376  latex10.SetText(0.085,0.4,"#bar{#nu} Mode")
377  latex10.Draw()
378  latex11 = latex9.Clone()
379  latex11.SetText(0.085,0.59,"#nu Mode")
380  latex11.Draw()
381  latex12 = latex9.Clone()
382  latex12.SetText(0.085,0.77,"#bar{#nu} Mode")
383  latex12.Draw()
384 
385  latex13 = ROOT.TLatex( 0.16, 0.124, "#nu_{#mu}" );
386  latex13.SetTextSize(0.02);
387  latex13.SetNDC()
388  latex13.SetTextFont(42);
389  latex13.Draw()
390  latex14 = ROOT.TLatex( 0.225, 0.124, "#bar{#nu}_{#mu}" );
391  latex14.SetTextSize(0.02);
392  latex14.SetNDC()
393  latex14.SetTextFont(42);
394  latex14.Draw()
395  latex15 = ROOT.TLatex( 0.27, 0.124, "#nu_{e}" );
396  latex15.SetTextSize(0.02);
397  latex15.SetNDC()
398  latex15.SetTextFont(42);
399  latex15.Draw()
400  latex16 = ROOT.TLatex( 0.295, 0.124, "#bar{#nu}_{e}" );
401  latex16.SetTextSize(0.02);
402  latex16.SetNDC()
403  latex16.SetTextFont(42);
404  latex16.Draw()
405 
406  latex17 = ROOT.TLatex( 0.34, 0.124, "#nu_{#mu}" );
407  latex17.SetTextSize(0.02);
408  latex17.SetNDC()
409  latex17.SetTextFont(42);
410  latex17.Draw()
411  latex18 = ROOT.TLatex( 0.405, 0.124, "#bar{#nu}_{#mu}" );
412  latex18.SetTextSize(0.02);
413  latex18.SetNDC()
414  latex18.SetTextFont(42);
415  latex18.Draw()
416  latex19 = ROOT.TLatex( 0.45, 0.124, "#nu_{e}" );
417  latex19.SetTextSize(0.02);
418  latex19.SetNDC()
419  latex19.SetTextFont(42);
420  latex19.Draw()
421  latex20 = ROOT.TLatex( 0.475, 0.124, "#bar{#nu}_{e}" );
422  latex20.SetTextSize(0.02);
423  latex20.SetNDC()
424  latex20.SetTextFont(42);
425  latex20.Draw()
426 
427  latex21 = ROOT.TLatex( 0.52, 0.124, "#nu_{#mu}" );
428  latex21.SetTextSize(0.02);
429  latex21.SetNDC()
430  latex21.SetTextFont(42);
431  latex21.Draw()
432  latex22 = ROOT.TLatex( 0.585, 0.124, "#bar{#nu}_{#mu}" );
433  latex22.SetTextSize(0.02);
434  latex22.SetNDC()
435  latex22.SetTextFont(42);
436  latex22.Draw()
437  latex23 = ROOT.TLatex( 0.631, 0.124, "#nu_{e}" );
438  latex23.SetTextSize(0.02);
439  latex23.SetNDC()
440  latex23.SetTextFont(42);
441  latex23.Draw()
442  latex24 = ROOT.TLatex( 0.654, 0.124, "#bar{#nu}_{e}" );
443  latex24.SetTextSize(0.02);
444  latex24.SetNDC()
445  latex24.SetTextFont(42);
446  latex24.Draw()
447 
448  latex25 = ROOT.TLatex( 0.70, 0.124, "#nu_{#mu}" );
449  latex25.SetTextSize(0.02);
450  latex25.SetNDC()
451  latex25.SetTextFont(42);
452  latex25.Draw()
453  latex26 = ROOT.TLatex( 0.77, 0.124, "#bar{#nu}_{#mu}" );
454  latex26.SetTextSize(0.02);
455  latex26.SetNDC()
456  latex26.SetTextFont(42);
457  latex26.Draw()
458  latex27 = ROOT.TLatex( 0.81, 0.124, "#nu_{e}" );
459  latex27.SetTextSize(0.02);
460  latex27.SetNDC()
461  latex27.SetTextFont(42);
462  latex27.Draw()
463  latex28 = ROOT.TLatex( 0.835, 0.124, "#bar{#nu}_{e}" );
464  latex28.SetTextSize(0.02);
465  latex28.SetNDC()
466  latex28.SetTextFont(42);
467  latex28.Draw()
468 
469  latex29 = ROOT.TLatex( 0.11,0.18, "#nu_{#mu}" );
470  latex29.SetTextAngle(90)
471  latex29.SetTextSize(0.02);
472  latex29.SetNDC()
473  latex29.SetTextFont(42);
474  latex29.Draw()
475  latex30 = ROOT.TLatex( 0.11,0.25, "#bar{#nu}_{#mu}" );
476  latex30.SetTextAngle(90)
477  latex30.SetTextSize(0.02);
478  latex30.SetNDC()
479  latex30.SetTextFont(42);
480  latex30.Draw()
481  latex31 = ROOT.TLatex( 0.11,0.295, "#nu_{e}" );
482  latex31.SetTextAngle(90)
483  latex31.SetTextSize(0.02);
484  latex31.SetNDC()
485  latex31.SetTextFont(42);
486  latex31.Draw()
487  latex32 = ROOT.TLatex( 0.11,0.32, "#bar{#nu}_{e}" );
488  latex32.SetTextAngle(90)
489  latex32.SetTextSize(0.02);
490  latex32.SetNDC()
491  latex32.SetTextFont(42);
492  latex32.Draw()
493 
494  latex33 = ROOT.TLatex( 0.11,0.37, "#nu_{#mu}" );
495  latex33.SetTextAngle(90)
496  latex33.SetTextSize(0.02);
497  latex33.SetNDC()
498  latex33.SetTextFont(42);
499  latex33.Draw()
500  latex34 = ROOT.TLatex( 0.11,0.435, "#bar{#nu}_{#mu}" );
501  latex34.SetTextAngle(90)
502  latex34.SetTextSize(0.02);
503  latex34.SetNDC()
504  latex34.SetTextFont(42);
505  latex34.Draw()
506  latex35 = ROOT.TLatex( 0.11,0.48, "#nu_{e}" );
507  latex35.SetTextAngle(90)
508  latex35.SetTextSize(0.02);
509  latex35.SetNDC()
510  latex35.SetTextFont(42);
511  latex35.Draw()
512  latex36 = ROOT.TLatex( 0.11,0.505, "#bar{#nu}_{e}" );
513  latex36.SetTextAngle(90)
514  latex36.SetTextSize(0.02);
515  latex36.SetNDC()
516  latex36.SetTextFont(42);
517  latex36.Draw()
518 
519  latex37 = ROOT.TLatex( 0.11,0.56, "#nu_{#mu}" );
520  latex37.SetTextAngle(90)
521  latex37.SetTextSize(0.02);
522  latex37.SetNDC()
523  latex37.SetTextFont(42);
524  latex37.Draw()
525  latex38 = ROOT.TLatex( 0.11,0.62, "#bar{#nu}_{#mu}" );
526  latex38.SetTextAngle(90)
527  latex38.SetTextSize(0.02);
528  latex38.SetNDC()
529  latex38.SetTextFont(42);
530  latex38.Draw()
531  latex39 = ROOT.TLatex( 0.11,0.67, "#nu_{e}" );
532  latex39.SetTextAngle(90)
533  latex39.SetTextSize(0.02);
534  latex39.SetNDC()
535  latex39.SetTextFont(42);
536  latex39.Draw()
537  latex40 = ROOT.TLatex( 0.11,0.70, "#bar{#nu}_{e}" );
538  latex40.SetTextAngle(90)
539  latex40.SetTextSize(0.02);
540  latex40.SetNDC()
541  latex40.SetTextFont(42);
542  latex40.Draw()
543 
544  latex41 = ROOT.TLatex( 0.11,0.75, "#nu_{#mu}" );
545  latex41.SetTextAngle(90)
546  latex41.SetTextSize(0.02);
547  latex41.SetNDC()
548  latex41.SetTextFont(42);
549  latex41.Draw()
550  latex42 = ROOT.TLatex( 0.11,0.81, "#bar{#nu}_{#mu}" );
551  latex42.SetTextAngle(90)
552  latex42.SetTextSize(0.02);
553  latex42.SetNDC()
554  latex42.SetTextFont(42);
555  latex42.Draw()
556  latex43 = ROOT.TLatex( 0.11,0.855, "#nu_{e}" );
557  latex43.SetTextAngle(90)
558  latex43.SetTextSize(0.02);
559  latex43.SetNDC()
560  latex43.SetTextFont(42);
561  latex43.Draw()
562  latex44 = ROOT.TLatex( 0.11,0.88, "#bar{#nu}_{e}" );
563  latex44.SetTextAngle(90)
564  latex44.SetTextSize(0.02);
565  latex44.SetNDC()
566  latex44.SetTextFont(42);
567  latex44.Draw()
568 
569  latex29 = ROOT.TLatex( 0.105,0.18, "#nu_{#mu}" );
570  latex29.SetTextSize(0.02);
571  latex29.SetNDC()
572  latex29.SetTextFont(42);
573  latex29.Draw()
574  latex30 = ROOT.TLatex( 0.105,0.25, "#bar{#nu}_{#mu}" );
575  latex30.SetTextSize(0.02);
576  latex30.SetNDC()
577  latex30.SetTextFont(42);
578  latex30.Draw()
579  latex31 = ROOT.TLatex( 0.105,0.295, "#nu_{e}" );
580  latex31.SetTextSize(0.02);
581  latex31.SetNDC()
582  latex31.SetTextFont(42);
583  latex31.Draw()
584  latex32 = ROOT.TLatex( 0.105,0.32, "#bar{#nu}_{e}" );
585  latex32.SetTextSize(0.02);
586  latex32.SetNDC()
587  latex32.SetTextFont(42);
588  latex32.Draw()
589 
590  latex33 = ROOT.TLatex( 0.105,0.37, "#nu_{#mu}" );
591  latex33.SetTextSize(0.02);
592  latex33.SetNDC()
593  latex33.SetTextFont(42);
594  latex33.Draw()
595  latex34 = ROOT.TLatex( 0.105,0.435, "#bar{#nu}_{#mu}" );
596  latex34.SetTextSize(0.02);
597  latex34.SetNDC()
598  latex34.SetTextFont(42);
599  latex34.Draw()
600  latex35 = ROOT.TLatex( 0.105,0.48, "#nu_{e}" );
601  latex35.SetTextSize(0.02);
602  latex35.SetNDC()
603  latex35.SetTextFont(42);
604  latex35.Draw()
605  latex36 = ROOT.TLatex( 0.105,0.505, "#bar{#nu}_{e}" );
606  latex36.SetTextSize(0.02);
607  latex36.SetNDC()
608  latex36.SetTextFont(42);
609  latex36.Draw()
610 
611  latex37 = ROOT.TLatex( 0.105,0.56, "#nu_{#mu}" );
612  latex37.SetTextSize(0.02);
613  latex37.SetNDC()
614  latex37.SetTextFont(42);
615  latex37.Draw()
616  latex38 = ROOT.TLatex( 0.105,0.62, "#bar{#nu}_{#mu}" );
617  latex38.SetTextSize(0.02);
618  latex38.SetNDC()
619  latex38.SetTextFont(42);
620  latex38.Draw()
621  latex39 = ROOT.TLatex( 0.105,0.67, "#nu_{e}" );
622  latex39.SetTextSize(0.02);
623  latex39.SetNDC()
624  latex39.SetTextFont(42);
625  latex39.Draw()
626  latex40 = ROOT.TLatex( 0.105,0.70, "#bar{#nu}_{e}" );
627  latex40.SetTextSize(0.02);
628  latex40.SetNDC()
629  latex40.SetTextFont(42);
630  latex40.Draw()
631 
632  latex41 = ROOT.TLatex( 0.105,0.75, "#nu_{#mu}" );
633  latex41.SetTextSize(0.02);
634  latex41.SetNDC()
635  latex41.SetTextFont(42);
636  latex41.Draw()
637  latex42 = ROOT.TLatex( 0.105,0.81, "#bar{#nu}_{#mu}" );
638  latex42.SetTextSize(0.02);
639  latex42.SetNDC()
640  latex42.SetTextFont(42);
641  latex42.Draw()
642  latex43 = ROOT.TLatex( 0.105,0.855, "#nu_{e}" );
643  latex43.SetTextSize(0.02);
644  latex43.SetNDC()
645  latex43.SetTextFont(42);
646  latex43.Draw()
647  latex44 = ROOT.TLatex( 0.105,0.88, "#bar{#nu}_{e}" );
648  latex44.SetTextSize(0.02);
649  latex44.SetNDC()
650  latex44.SetTextFont(42);
651  latex44.Draw()
652 
653  line_locations = [19.5,38.5,45.5,52.5,71.5,90.5,97.5,104.5,123.5,142.5,149.5,156.5,175.5,194.5,201.5]
654  if small_bins:
655  line_locations = [19.5,38.5,57.5,76.5,95.5,114.5,133.5,152.5,171.5,190.5,209.5,228.5,247.5,266.5,285.5]
656 
657  line_1 = ROOT.TLine(19.5,0.5,19.5,208)
658  line_2 = ROOT.TLine(38.5,0.5,38.5,208)
659  line_3 = ROOT.TLine(45.5,0.5,45.5,208)
660  line_4 = ROOT.TLine(52.5,0.5,52.5,208)
661  line_5 = ROOT.TLine(71.5,0.5,71.5,208)
662  line_6 = ROOT.TLine(90.5,0.5,90.5,208)
663  line_7 = ROOT.TLine(97.5,0.5,97.5,208)
664  line_8 = ROOT.TLine(104.5,0.5,104.5,208)
665  line_9 = ROOT.TLine(123.5,0.5,123.5,208)
666  line_10 = ROOT.TLine(142.5,0.5,142.5,208)
667  line_11 = ROOT.TLine(149.5,0.5,149.5,208)
668  line_12 = ROOT.TLine(156.5,0.5,156.5,208)
669  line_13 = ROOT.TLine(175.5,0.5,175.5,208)
670  line_14 = ROOT.TLine(194.5,0.5,194.5,208)
671  line_15 = ROOT.TLine(201.5,0.5,201.5,208)
672  lineb_1 = ROOT.TLine(0.5,19.5,208,19.5)
673  lineb_2 = ROOT.TLine(0.5,38.5,208,38.5)
674  lineb_3 = ROOT.TLine(0.5,45.5,208,45.5)
675  lineb_4 = ROOT.TLine(0.5,52.5,208,52.5)
676  lineb_5 = ROOT.TLine(0.5,71.5,208,71.5)
677  lineb_6 = ROOT.TLine(0.5,90.5,208,90.5)
678  lineb_7 = ROOT.TLine(0.5,97.5,208,97.5)
679  lineb_8 = ROOT.TLine(0.5,104.5,208,104.5)
680  lineb_9 = ROOT.TLine(0.5,123.5,208,123.5)
681  lineb_10 = ROOT.TLine(0.5,142.5,208,142.5)
682  lineb_11 = ROOT.TLine(0.5,149.5,208,149.5)
683  lineb_12 = ROOT.TLine(0.5,156.5,208,156.5)
684  lineb_13 = ROOT.TLine(0.5,175.5,208,175.5)
685  lineb_14 = ROOT.TLine(0.5,194.5,208,194.5)
686  lineb_15 = ROOT.TLine(0.5,201.5,208,201.5)
687 
688  line_8.SetLineWidth(4)
689  lineb_8.SetLineWidth(4)
690 
691  line_4.SetLineWidth(2)
692  lineb_4.SetLineWidth(2)
693  line_12.SetLineWidth(2)
694  lineb_12.SetLineWidth(2)
695 
696  for line in [line_1,line_2,line_3,line_4,line_5,line_6,line_7,line_8,line_9,line_10,line_11,line_12,line_13,line_14,line_15,lineb_1,lineb_2,lineb_3,lineb_4,lineb_5,lineb_6,lineb_7,lineb_8,lineb_9,lineb_10,lineb_11,lineb_12,lineb_13,lineb_14,lineb_15]:
697  line.SetLineColor(13)
698  line.Draw()
699 
700 c1.Print("focusing_covariance_"+beam_config+".eps")
701 c1.Print("focusing_covariance_"+beam_config+".png")
702 
703 
704 
705 NRGBs = 3;
706 NCont = 999;
707 ROOT.gStyle.SetNumberContours(NCont);
708 stops = [ 0.00, 0.50, 1.00];
709 red = [ 0.00, 1.00, 1.00];
710 green = [ 0.00, 1.00, 0.00];
711 blue = [ 1.00, 1.00, 0.00];
712 ROOT.TColor.CreateGradientColorTable(NRGBs, array.array('d',stops), array.array('d',red), array.array('d',green), array.array('d',blue), NCont);
713 
714 correlation = GetCorrelation(covariance)
715 correlation.Draw("COLZ")
716 
717 for line in [line_1,line_2,line_3,line_4,line_5,line_6,line_7,line_8,line_9,line_10,line_11,line_12,line_13,line_14,line_15,lineb_1,lineb_2,lineb_3,lineb_4,lineb_5,lineb_6,lineb_7,lineb_8,lineb_9,lineb_10,lineb_11,lineb_12,lineb_13,lineb_14,lineb_15]:
718  line.SetLineColor(13)
719  line.Draw()
720 for latex in [latex1,latex2,latex3,latex4,latex5,latex6,latex7,latex8,latex9,latex10,latex11,latex12,latex13,latex14,latex15,latex16,latex17,latex18,latex19,latex20,latex21,latex22,latex23,latex24,latex25,latex26,latex27,latex28,latex29,latex30,latex31,latex32,latex33,latex34,latex35]:
721  latex.Draw();
722 
723 small_bins_string = ""
724 if small_bins: small_bins_string = "_small_bins"
725 c1.Print("focusing_correlation_"+beam_config+small_bins_string+".eps")
726 c1.Print("focusing_correlation_"+beam_config+small_bins_string+".png")
727 
728 f = ROOT.TFile("focusing_covariance_"+beam_config+small_bins_string+".root","RECREATE")
729 covariance.Write()
730 correlation.Write()
731 
732 
733 # Draw plots of individual sources of uncertainty for each flavor/mode/detector
734 plots = ["numu_neutrino_ND","numubar_neutrino_ND","nue_neutrino_ND","nuebar_neutrino_ND","numu_antineutrino_ND","numubar_antineutrino_ND","nue_antineutrino_ND","nuebar_antineutrino_ND","numu_neutrino_FD","numubar_neutrino_FD","nue_neutrino_FD","nuebar_neutrino_FD","numu_antineutrino_FD","numubar_antineutrino_FD","nue_antineutrino_FD","nuebar_antineutrino_FD"]
735 
736 start_indices = [0,19,38,45,52,71,90,97,104,123,142,149,156,175,194,201]
737 for plot in plots:
738  f = ROOT.TFile("1sigma_histograms_"+beam_config+small_bins_string+"_"+plot+".root","RECREATE")
739  #error_bins = [0,2,4,6,8,10,20,100]
740  #if plot.startswith("numu") or small_bins:
741  # error_bins = [0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,7 ,8,12,16,20,40,100]
742 
743  error_bins = [0,2,4,6,8,10]
744  if plot.startswith("numu") or small_bins:
745  error_bins = [0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,7 ,8,12]
746  error_hists = []
747  total_uncertainties = [0]*(len(error_bins)-1)
748  for i in range(len(uncertainties)):
749 
750  temp_hist = ROOT.TH1D(uncertainty_titles[i],uncertainty_titles[i],len(error_bins)-1,array.array('d',error_bins))
751  for j in range(len(error_bins)-1):
752  # Temporarily turning off abs to plot 1 sigma shifts instead of errors
753 # temp_hist.SetBinContent(j+1,abs(shifts[i][start_indices[plots.index(plot)]+j]))
754  temp_hist.SetBinContent(j+1,(shifts[i][start_indices[plots.index(plot)]+j]))
755  total_uncertainties[j] = total_uncertainties[j]+math.pow(shifts[i][start_indices[plots.index(plot)]+j],2)
756  error_hists.append(temp_hist)
757 
758 
759  temp_hist = ROOT.TH1D("Total","Total",len(error_bins)-1,array.array('d',error_bins))
760  for j in range(len(error_bins)-1):
761  temp_hist.SetBinContent(j+1,math.sqrt(total_uncertainties[j]))
762  error_hists.append(temp_hist)
763 
764 
765 # error_hists[-1].SetMinimum(0);
766  error_hists[-1].SetMinimum(-0.09);
767  error_hists[-1].SetMaximum(0.09);
768  error_hists[-1].SetLineWidth(3)
769  error_hists[-1].Draw("hist")
770  error_hists[-1].GetXaxis().SetTitle("Neutrino Energy (GeV)")
771 # error_hists[-1].GetYaxis().SetTitle("Fractional Uncertainty")
772  error_hists[-1].GetYaxis().SetTitle("1-Sigma Shift")
773  error_hists[-1].GetYaxis().SetTitleOffset(1.3)
774 
775  for i in range(0,len(error_hists)-1):
776  colors = [2,ROOT.kOrange+3,6,ROOT.kOrange+7,ROOT.kGreen+2,3,ROOT.kAzure+7,ROOT.kBlue-4,7,ROOT.kViolet+1,13,15,17,19]
777  error_hists[i].SetLineColor(colors[i])
778  error_hists[i].SetLineWidth(3)
779  error_hists[i].SetMinimum(-0.09);
780  error_hists[i].SetMaximum(0.09);
781  error_hists[i].GetXaxis().SetTitle("Neutrino Energy (GeV)")
782 # error_hists[i].GetYaxis().SetTitle("Fractional Uncertainty")
783  error_hists[i].GetYaxis().SetTitle("1-Sigma Shift")
784  error_hists[i].GetYaxis().SetTitleOffset(1.3)
785  error_hists[i].Draw("histsame")
786  error_hists[i].Write();
787 
788  leg = ROOT.TLegend(0.57,0.5,0.92,0.85)
789  leg.AddEntry(error_hists[-1],error_hists[-1].GetTitle(),"l")
790  for i in range(0,len(error_hists)-1):
791  leg.AddEntry(error_hists[i],error_hists[i].GetTitle(),"l")
792  leg.SetFillStyle(0)
793  leg.Draw()
794 
795  c1.Print("focusing_error_overlay_"+plot+"_"+beam_config+small_bins_string+".eps")
796  c1.Print("focusing_error_overlay_"+plot+"_"+beam_config+small_bins_string+".png")
797 
798 
799  f.Close()
open
int open(const char *, int)
Opens a file descriptor.
cet::split
void split(std::string const &s, char c, OutIter dest)
Definition: split.h:35
if
if(!yymsg) yymsg
Generated by   doxygen 1.8.11