plot_xsec_2d.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 
3 from __future__ import print_function
4 
5 from ROOT import * # import pyroot
6 from bisect import bisect
7 import argparse
8 import math
9 import numpy as np
10 import os
11 
12 # define constants used in this script
13 infp = './'
14 syst_names = [
15  'lightdown','lightup','ckv','calibneg','calibpos','calibshape','MuES_onlyMuKEDw','MuES_onlyMuKEUp',
16  '2kA_Dw','2kA_Up','02mmBeamSpotSize_Dw','02mmBeamSpotSize_Up','1mmBeamShiftX_Dw','1mmBeamShiftX_Up',
17  '1mmBeamShiftY_Dw','1mmBeamShiftY_Up','3mmHorn1X_Dw','3mmHorn1X_Up','3mmHorn1Y_Dw','3mmHorn1Y_Up',
18  '3mmHorn2X_Dw','3mmHorn2X_Up','3mmHorn2Y_Dw','3mmHorn2Y_Up','7mmTargetZ_Dw','7mmTargetZ_Up',
19  'MagneticFieldinDecayPipe_Dw','MagneticFieldinDecayPipe_Up','1mmHornWater_Dw','1mmHornWater_Up'
20 ]
21 Nfinal_fe = 9
22 final_fe = ["Flux (HP)","GENIE","Light","Calibration","Cherenkov","Calib. Shape","Mu Energy Scale","Normalization","Focusing"]
23 final_fe_name = ["Flux_HP","GENIE","Light","Calibration","Cherenkov","CalibShape","MuEScale","Normalization","Focusing"]
24 
25 Nabins2X = 13 # pre calculated N bins X
26 leg_angle = [
27  "0.50<Cos#theta_{#mu}<0.56",
28  "0.56<Cos#theta_{#mu}<0.62",
29  "0.62<Cos#theta_{#mu}<0.68",
30  "0.68<Cos#theta_{#mu}<0.74",
31  "0.74<Cos#theta_{#mu}<0.80",
32  "0.80<Cos#theta_{#mu}<0.85",
33  "0.85<Cos#theta_{#mu}<0.88",
34  "0.88<Cos#theta_{#mu}<0.91",
35  "0.91<Cos#theta_{#mu}<0.94",
36  "0.94<Cos#theta_{#mu}<0.96",
37  "0.96<Cos#theta_{#mu}<0.98",
38  "0.98<Cos#theta_{#mu}<0.99",
39  "0.99<Cos#theta_{#mu}<1.00"
40 ]
41 col_fe = [2,4,2,4,8,8,kOrange-3,kOrange-3,1]
42 col_st = [1,1,2,2,1,2,1,2,2]
43 
44 tit_ke = 'Muon Kinetic Energy (GeV)'
45 tit_nu = 'Neutrino Energy (GeV)'
46 tit_fe = 'Fractional Uncertainties'
47 tit_xs = '#sigma(cm^{2}/nucleon)'
48 
49 
50 def get_genie_universes(basefn):
51 
52  hl = []
53  fn = os.path.join(infp, basefn)
54  fgenie.append(TFile(fn))
55  # find the genie universe IDs from filename
56  sIdx = int(basefn.split('_')[1])
57  eIdx = int(basefn.split('_')[2])+1
58  for i in range(sIdx,eIdx):
59  hl.append(fgenie[len(fgenie)-1].Get('xs_nominal/hXS2Dsec_genie{}'.format(i)))
60 
61  return hl
62 
63 def get_ppfx_universes(fppfx_other):
64 
65  hl = []
66  # by default we have 100 ppfx universes
67  for i in range(0,100):
68  hl.append(fppfx_other.Get('xs_nominal/hXS2Dsec_ppfx{}'.format(i)))
69 
70  return hl
71 
72 
73 def up_down_counts_1sigma(count_list, cv):
74  sorted_list = sorted(count_list)
75  cv_pos = bisect(sorted_list, cv)
76  up_pos = int((len(sorted_list) - cv_pos)*.68 + cv_pos)
77  down_pos = int(cv_pos - cv_pos*.68)
78  return sorted_list[up_pos]-cv, cv-sorted_list[down_pos]
79 
80 
81 def process_multiuniverse(isys, huniv):
82 
83  for i in range(1, NbinsX+1):
84  for j in range(1, NbinsY+1):
85  dt_cv = hdt.GetBinContent(i,j)
86 
87  if dt_cv <= 0:
88  lhup[isys].SetBinContent(i,j,0)
89  lhdw[isys].SetBinContent(i,j,0)
90  lhfe[isys].SetBinContent(i,j,0)
91  else:
92  # accumulate bin content from all ppfx universes
93  binc_all_uni = [h.GetBinContent(i,j) for h in huniv]
94  cup, cdown = up_down_counts_1sigma(binc_all_uni, dt_cv)
95  lhup[isys].SetBinContent(i,j,cup)
96  lhdw[isys].SetBinContent(i,j,cdown)
97  # calculate conservative fractional uncertainty
98  conserv_error = max(cup,cdown)
99  lhfe[isys].SetBinContent(i,j,conserv_error/dt_cv)
100 
101 
102 def process_light_calibration_mue(isys, hidxdown, hidxup):
103 
104  for i in range(1, NbinsX+1):
105  for j in range(1, NbinsY+1):
106  dt_cv = hdt.GetBinContent(i,j)
107 
108  lhup[isys].SetBinContent(i,j,0)
109  lhdw[isys].SetBinContent(i,j,0)
110  lhfe[isys].SetBinContent(i,j,0)
111 
112  sigma_dw = dt_cv - hother[hidxdown].GetBinContent(i,j)
113  sigma_up = dt_cv - hother[hidxup].GetBinContent(i,j)
114  sigma = max(abs(sigma_dw), abs(sigma_up))
115 
116  if sigma_dw < 0 and sigma_up < 0:
117  lhup[isys].SetBinContent(i,j,sigma)
118  elif sigma_dw > 0 and sigma_up > 0:
119  lhdw[isys].SetBinContent(i,j,sigma)
120  else:
121  lhup[isys].SetBinContent(i,j,abs(sigma_up))
122  lhdw[isys].SetBinContent(i,j,abs(sigma_dw))
123 
124  if dt_cv > 0:
125  lhfe[isys].SetBinContent(i,j,sigma/dt_cv)
126  else:
127  lhfe[isys].SetBinContent(i,j,0)
128 
129 
130 def process_ckv_calibshape(isys, hidx):
131 
132  for i in range(1, NbinsX+1):
133  for j in range(1, NbinsY+1):
134  dt_cv = hdt.GetBinContent(i,j)
135 
136  lhup[isys].SetBinContent(i,j,0)
137  lhdw[isys].SetBinContent(i,j,0)
138  lhfe[isys].SetBinContent(i,j,0)
139 
140  sigma = dt_cv - hother[hidx].GetBinContent(i,j)
141 
142  # if sigma > 0:
143  # lhdw[isys].SetBinContent(i,j,abs(sigma))
144  # if sigma <= 0:
145  # lhup[isys].SetBinContent(i,j,abs(sigma))
146  lhdw[isys].SetBinContent(i,j,abs(sigma))
147  lhup[isys].SetBinContent(i,j,abs(sigma))
148 
149  if dt_cv > 0:
150  lhfe[isys].SetBinContent(i,j,sigma/dt_cv)
151  else:
152  lhfe[isys].SetBinContent(i,j,0)
153 
154 
156 
157  for i in range(1, NbinsX+1):
158  for j in range(1, NbinsY+1):
159  dt_cv = hdt.GetBinContent(i,j)
160 
161  sigma = dt_cv * fe_norm
162  # lhup[isys].SetBinContent(i,j,dt_cv+sigma)
163  # lhdw[isys].SetBinContent(i,j,dt_cv-sigma)
164  lhup[isys].SetBinContent(i,j,sigma)
165  lhdw[isys].SetBinContent(i,j,sigma)
166  if dt_cv > 0:
167  lhfe[isys].SetBinContent(i,j,fe_norm)
168  else:
169  lhfe[isys].SetBinContent(i,j,0)
170 
171 
172 def process_focusing(isys, sstart, nsys):
173 
174  # debug information
175  print('Focusing systematics:')
176  for k in range(nsys):
177  print(syst_names[sstart+k*2], syst_names[sstart+k*2+1])
178 
179  for i in range(1, NbinsX+1):
180  for j in range(1, NbinsY+1):
181  dt_cv = hdt.GetBinContent(i,j)
182 
183  lhup[isys].SetBinContent(i,j,0)
184  lhdw[isys].SetBinContent(i,j,0)
185  lhfe[isys].SetBinContent(i,j,0)
186 
187  larger_sigma = []
188  for k in range(nsys):
189  sigma_dw = dt_cv - hother[sstart+k*2].GetBinContent(i,j)
190  sigma_up = dt_cv - hother[sstart+k*2+1].GetBinContent(i,j)
191  larger_sigma.append(max(abs(sigma_dw), abs(sigma_up)))
192 
193  qsum = math.sqrt(sum([x**2 for x in larger_sigma]))
194  lhup[isys].SetBinContent(i,j,qsum)
195  lhdw[isys].SetBinContent(i,j,qsum)
196 
197  if dt_cv > 0:
198  lhfe[isys].SetBinContent(i,j,qsum/dt_cv)
199  else:
200  lhfe[isys].SetBinContent(i,j,0)
201 
202 
203 def databin(xbin, ybin):
204  answer = False
205  if xbin==1 and ybin<=4: answer = True
206  if (xbin==2 or xbin==3) and ybin<=5: answer = True
207  if xbin==4 and ybin<=7: answer = True
208  if xbin==5 and ybin<=9: answer = True
209  if xbin==6 and ybin<=10: answer = True
210  if xbin==7 and ybin<=12: answer = True
211  if xbin==8 and ybin<=15: answer = True
212  if xbin==9 and ybin<=19: answer = True
213  if (xbin==10 or xbin==11 or xbin==12 or xbin==13) and ybin<=20: answer = True
214 
215  return answer
216 
217 
219 
220  for i in range(1, NbinsX+1):
221  for j in range(1, NbinsY+1):
222  kin = hdt.GetYaxis().GetBinCenter(j)
223  binw = hdt.GetYaxis().GetBinWidth(j)
224  cont = hdt.GetBinContent(i,j)
225  errUp = 0
226  errDw = 0
227 
228  for p in range(Nfinal_fe):
229  up = lhup[p].GetBinContent(i,j)
230  dw = lhdw[p].GetBinContent(i,j)
231  errUp += up**2
232  errDw += dw**2
233 
234  errUp = math.sqrt(errUp)
235  errDw = math.sqrt(errDw)
236 
237  if databin(i,j):
238  count[i-1] += 1
239  gdt[i-1].SetPoint(count[i-1], kin, cont)
240  gdt[i-1].SetPointError(count[i-1],binw/2, binw/2, errDw, errUp)
241 
242 
244 
245  for i in range(1, NbinsX+1):
246  for j in range(1, NbinsY+1):
247  if databin(i,j):
248  binw = hdt.GetYaxis().GetBinWidth(j)
249  cont = hmc.GetBinContent(i,j)
250  h1MC[i-1].SetBinContent(j,cont)
251  for p in range(Nfinal_fe):
252  ferr = lhfe[p].GetBinContent(i,j)
253  h1fe[i-1][p].SetBinContent(j,ferr)
254 
255 def setMC1D(hIn):
256  hIn.SetLineColor(2)
257  hIn.SetLineWidth(2)
258  hIn.SetStats(0)
259 
260 def setData1D(gIn):
261  gIn.SetMarkerColor(kBlack)
262  gIn.SetMarkerStyle(20)
263  gIn.SetLineColor(kBlack)
264  gIn.SetLineWidth(3)
265 
266 def setUncer1D(hIn, ihist):
267  hIn.SetLineColor(col_fe[ihist])
268  hIn.SetLineWidth(3)
269  hIn.SetLineStyle(col_st[ihist])
270  hIn.SetStats(0)
271 
272 
273 if __name__ == '__main__':
274  '''
275  Remember:
276  The unit of the output plots from this script has to be
277  rescaled by bin size when showing the final results!
278  '''
279 
280  # fake data / real data switch
281  parser = argparse.ArgumentParser()
282  parser.add_argument('-d', '--data', action='store_true', help='Switch for using data instead of fake data.')
283  args = parser.parse_args()
284  use_data = args.data
285 
286  # read histograms into containers
287 
288  # file containers
289  fgenie = []
290  fmc = TFile(os.path.join(infp, 'nominal_xsec.root'))
291  if not use_data:
292  fppfx_other = TFile(os.path.join(infp, 'ppfx_others_v2.root'))
293  else:
294  fppfx_other = TFile(os.path.join(infp, 'ppfx_others_data.root'))
295 
296  # historam containers
297  hgenie = []
298  hppfx = []
299  hother = []
300 
301  # get nominal cross section
302  hmc = fmc.Get('hXS2Dsec_nominal')
303  # get (fake) data cross section
304  hdt = fppfx_other.Get('xs_nominal/hXS2Dsec_cv')
305 
306  # store genie histograms
307  if not use_data:
308  hgenie += get_genie_universes('genie_0_99_v2.root')
309  hgenie += get_genie_universes('genie_100_299_v2.root')
310  hgenie += get_genie_universes('genie_300_499_v2.root')
311  else:
312  hgenie += get_genie_universes('genie_0_99_data.root')
313  hgenie += get_genie_universes('genie_100_299_data.root')
314  hgenie += get_genie_universes('genie_300_499_data.root')
315 
316  # store ppfx histograms
317  hppfx += get_ppfx_universes(fppfx_other)
318 
319  # load all other systematics
320  for syst_name in syst_names:
321  hother.append(fppfx_other.Get('xs_nominal/hXS2Dsec_{}'.format(syst_name)))
322 
323  # below starts to calculate up and down histograms
324  lhup = []
325  lhdw = []
326  lhfe = []
327  for i in range(Nfinal_fe):
328  lhup.append(hmc.Clone('hUp_{}'.format(i)))
329  lhdw.append(hmc.Clone('hDw_{}'.format(i)))
330  lhfe.append(hmc.Clone('hFE_{}'.format(i)))
331 
332  # some global variables common to all systematics
333  # number of bins
334  NbinsX = hdt.GetNbinsX()
335  NbinsY = hdt.GetNbinsY()
336 
337  # ppfx
338  process_multiuniverse(0, hppfx)
339 
340  # genie
341  process_multiuniverse(1, hgenie)
342 
343  # light
345 
346  # calibration
348 
349  # cherenkov
350  # this does not look right
352 
353  # calibshape
354  # this does not look right
356 
357  # muon energy scale
359 
360  # normalization
361  # Normalization (TN, doc-):
362  # POT: 0.28% (Section 9.1)
363  # Intensity: 2% (Section 9.2)
364  # Particle propagation: 1-2%, so we take 2% (Section 9.5).
365  # fe_norm = sqrt(pow(0.0028,2)+pow(0.02,2)+pow(0.02,2));
366  fe_norm = 0.028422526
368 
369  # focusing
370  # process_light_calibration_mue_foc(8, 8, 9)
371  process_focusing(8, 8, 11)
372 
373  # data
374  gdt = [TGraphAsymmErrors() for i in range(Nabins2X)]
375  count = [0 for i in range(Nabins2X)]
376  process_data()
377 
378  # MC
379  h1MC = [TH1D('h1MC_{}'.format(i),'',20,0.5,2.5) for i in range(Nabins2X)]
380  h1fe = [[TH1D('h1fe_{}_{}'.format(i,j),'',20,0.5,2.5) for j in range(Nfinal_fe)] for i in range(Nabins2X)]
381  process_mc()
382 
383  # make sure the output directory exists
384  out_subdir = 'real_data' if use_data else 'fake_data'
385  dest_path = os.path.join('output', out_subdir)
386  if not os.path.exists(dest_path):
387  os.makedirs(dest_path)
388 
389  # make cross section plots
390  cxs = []
391  leg = [TLegend(0.67,0.68,0.87,0.88) for i in range(Nabins2X)]
392  lat = [TLatex() for i in range(Nabins2X)]
393  maxY = -1.
394  for i in range(NbinsX):
395  maxY = h1MC[i].GetMaximum()
396  h1MC[i].SetMaximum(maxY*2)
397  cxs.append(TCanvas())
398  lat[i].SetNDC()
399  leg[i].SetBorderSize(0)
400  setMC1D(h1MC[i])
401  setData1D(gdt[i])
402  h1MC[i].GetXaxis().SetTitle('Muon Kinetic Energy (GeV)')
403  h1MC[i].GetXaxis().CenterTitle()
404  h1MC[i].GetYaxis().SetTitle('d^{2}#sigma/dT_{#mu}dcos#theta_{#mu} (cm^{2}/GeV/nucleon)')
405  h1MC[i].GetYaxis().CenterTitle()
406  h1MC[i].Draw("hist")
407  gdt[i].Draw("p1 same")
408  leg[i].SetTextSize(0.05)
409  leg[i].AddEntry(h1MC[i], "MC", "fl")
410  leg[i].AddEntry(gdt[i],"Data (Syst.)", "lp")
411  leg[i].Draw()
412  lat[i].DrawLatex(0.2,0.83,'{}'.format(leg_angle[i]))
413  cxs[i].SaveAs('{}/xs2D_anglebin{}.png'.format(dest_path,i))
414  cxs[i].SaveAs('{}/xs2D_anglebin{}.pdf'.format(dest_path,i))
415 
416  # make fractional uncertainty plots
417  cfe = []
418  latfe = [TLatex() for i in range(Nabins2X)]
419  lfe = []
420  for i in range(NbinsX):
421  if i > 6:
422  lfe = TLegend(0.2,0.5,0.45,0.88)
423  else:
424  lfe = TLegend(0.6,0.4,0.85,0.78)
425  lfe.SetBorderSize(0)
426  cfe.append(TCanvas())
427  latfe[i].SetNDC()
428  for j in range(Nfinal_fe):
429  setUncer1D(h1fe[i][j],j)
430  h1fe[i][j].GetXaxis().SetTitle('Muon Kinetic Energy (GeV)')
431  h1fe[i][j].GetXaxis().CenterTitle()
432  h1fe[i][j].GetYaxis().SetTitle('Fractional Uncertainty')
433  h1fe[i][j].GetYaxis().CenterTitle()
434  if j == 0:
435  h1fe[i][j].Draw('hist')
436  else:
437  h1fe[i][j].Draw('histsame')
438  h1fe[i][j].SetMaximum(.4)
439  h1fe[i][j].SetMinimum(0)
440  lfe.AddEntry(h1fe[i][j],final_fe[j],'l')
441  latfe[i].DrawLatex(0.63,0.83,'{}'.format(leg_angle[i]))
442  lfe.Draw()
443  cfe[i].SaveAs('{}/fe2D_anglebin{}.png'.format(dest_path,i))
444  cfe[i].SaveAs('{}/fe2D_anglebin{}.pdf'.format(dest_path,i))
445 
446  # storing the final histograms
447  fileOut = TFile('{}/numucc_inc_{}.root'.format(dest_path,out_subdir), 'recreate')
448  fileOut.cd()
449  fileOut.mkdir('CrossSections')
450  for i in range(NbinsX):
451  fileOut.cd('CrossSections')
452  gdt[i].GetXaxis().SetTitle(tit_ke)
453  gdt[i].GetYaxis().SetTitle(tit_xs)
454  gdt[i].Write('mock_data_bin{}'.format(i))
455  h1MC[i].Write('mc_bin{}'.format(i))
456  fileOut.mkdir('FractUncer_bin{}'.format(i))
457  fileOut.cd('FractUncer_bin{}'.format(i))
458  for j in range(Nfinal_fe):
459  h1fe[i][j].Write('FractUncer_bin{}_{}'.format(i,final_fe_name[j]))
460 
461  raw_input('Press enter to continue.')
def process_light_calibration_mue(isys, hidxdown, hidxup)
tree Draw("slc.nhit")
h7 SetMaximum(maxY)
correl_xv GetYaxis() -> SetDecimals()
def process_focusing(isys, sstart, nsys)
prelim SetTextSize(2/30.)
def process_data()
void abs(TH1 *hist)
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
def setMC1D(hIn)
def process_ckv_calibshape(isys, hidx)
gargamelle SetTitle("Gargamelle #nu_{e} CC data")
def process_normalization(isys)
legend SetBorderSize(0)
correl_xv GetXaxis() -> SetDecimals()
bool print
leg AddEntry(GRdata,"data","p")
tex DrawLatex(28.4, 12.8,"x = 0.015(5*F2)")
def databin(xbin, ybin)
std::string format(const int32_t &value, const int &ndigits=8)
Definition: HexUtils.cpp:14
cout<< "--"<< endl;for(Int_t iP=1;iP<=hyz->GetNbinsX();iP++){for(Int_t iC=1;iC<=hyz->GetNbinsY();iC++){if(hyv->GetBinContent(iP, iC)>-999){goal_hyv-> SetBinContent(iP, iC,-(dy[iP-1][iC-1]))
def get_ppfx_universes(fppfx_other)
Definition: plot_xsec_2d.py:63
def process_multiuniverse(isys, huniv)
Definition: plot_xsec_2d.py:81
def get_genie_universes(basefn)
Definition: plot_xsec_2d.py:50
histo SetMinimum(1.0e-9)
Double_t sum
Definition: plot.C:31
cosmicTree SaveAs("cosmicTree.root")
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
def up_down_counts_1sigma(count_list, cv)
Definition: plot_xsec_2d.py:73
def setUncer1D(hIn, ihist)
prelim SetNDC()
def setData1D(gIn)
def process_mc()
gm Write()