plot_xsec_1d.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 
3 from __future__ import print_function
4 from ROOT import *
5 from bisect import bisect
6 import argparse
7 import math
8 import os
9 
10 
11 # global variables
12 
13 # file directory
14 infp = './'
15 
16 
17 syst_names = [
18  'lightdown','lightup','ckv','calibneg','calibpos','calibshape','MuES_onlyMuKEDw','MuES_onlyMuKEUp',
19  '2kA_Dw','2kA_Up','02mmBeamSpotSize_Dw','02mmBeamSpotSize_Up','1mmBeamShiftX_Dw','1mmBeamShiftX_Up',
20  '1mmBeamShiftY_Dw','1mmBeamShiftY_Up','3mmHorn1X_Dw','3mmHorn1X_Up','3mmHorn1Y_Dw','3mmHorn1Y_Up',
21  '3mmHorn2X_Dw','3mmHorn2X_Up','3mmHorn2Y_Dw','3mmHorn2Y_Up','7mmTargetZ_Dw','7mmTargetZ_Up',
22  'MagneticFieldinDecayPipe_Dw','MagneticFieldinDecayPipe_Up','1mmHornWater_Dw','1mmHornWater_Up'
23 ]
24 
25 # names of the fractional uncertainties
26 fu_names = ['Flux_HP','GENIE','Light','Calibration','Cherenkov','CalibShape','MuEScale','Normalization','Focusing']
27 final_fu = ['Flux (HP)','GENIE','Light','Calibration','Cherenkov','Calib. Shape','Mu Energy Scale','Normalization','Focusing']
28 col_fe = [2,4,2,4,8,8,kOrange-3,kOrange-3,1]
29 col_st = [1,1,2,2,1,2,1,2,2]
30 
31 def get_genie_universes(basefn):
32 
33  hl = []
34  fn = os.path.join(infp, basefn)
35  fgenie.append(TFile(fn))
36  # find the genie universe IDs from filename
37  sIdx = int(basefn.split('_')[1])
38  eIdx = int(basefn.split('_')[2])+1
39  for i in range(sIdx,eIdx):
40  hl.append(fgenie[len(fgenie)-1].Get('xs_nominal/hXS1Dsec_genie{}'.format(i)))
41 
42  return hl
43 
44 def get_ppfx_universes(fppfx_other):
45 
46  hl = []
47  # by default we have 100 ppfx universes
48  for i in range(0,100):
49  hl.append(fppfx_other.Get('xs_nominal/hXS1Dsec_ppfx{}'.format(i)))
50 
51  return hl
52 
53 
54 def do_multiverse_1bin(isys, univ_vals):
55 
56  sorted_list = sorted(univ_vals)
57  cv_pos = bisect(sorted_list, cv)
58  nlist = len(univ_vals)
59  up_pos = nlist-1 if cv_pos >= nlist else int((nlist - cv_pos)*.68 + cv_pos)
60  dw_pos = 0 if cv_pos <= 0 else int(cv_pos - cv_pos*.68)
61 
62  up_shift = sorted_list[up_pos] - cv
63  dw_shift = cv - sorted_list[dw_pos]
64 
65  lhup[isys].SetBinContent(b, up_shift)
66  lhdw[isys].SetBinContent(b, dw_shift)
67  fu = 0 if cv <= 0 else max(up_shift, dw_shift)/cv
68  lhfu[isys].SetBinContent(b, fu)
69 
70 
71 def do_focusing_1bin(isys, sstart, nsys):
72 
73  up_shifts = [0 for k in range(nsys)]
74  dw_shifts = [0 for k in range(nsys)]
75  for k in range(nsys):
76  dwshift = cv - hother[sstart+k*2].GetBinContent(i)
77  upshift = hother[sstart+k*2+1].GetBinContent(i) - cv
78 
79  if upshift > 0 and dwshift > 0:
80  up_shifts[k] = upshift
81  dw_shifts[k] = dwshift
82  elif upshift <= 0 and dwshift <= 0: # up down flip
83  up_shifts[k] = -dwshift
84  dw_shifts[k] = -upshift
85  elif upshift > 0 and dwshift <= 0:
86  up_shifts[k] = upshift
87  else:
88  dw_shifts[k] = dwshift
89 
90  up_shift = math.sqrt(sum([x**2 for x in up_shifts]))
91  dw_shift = math.sqrt(sum([x**2 for x in dw_shifts]))
92  lhup[isys].SetBinContent(b, up_shift)
93  lhdw[isys].SetBinContent(b, dw_shift)
94 
95  fu = 0 if cv <= 0 else max(up_shift, dw_shift)/cv
96  lhfu[isys].SetBinContent(b, fu)
97 
98 
99 def do_up_down_1bin(isys, upidx, dwidx):
100  '''
101  It could happen that in some regions both systematics go the same direction.
102  In that case we take the larger value for that direction and asign zero to the other.
103  '''
104 
105  # initialize
106  lhup[isys].SetBinContent(b, 0)
107  lhdw[isys].SetBinContent(b, 0)
108 
109  upshift = hother[upidx].GetBinContent(b)-cv
110  dwshift = cv-hother[dwidx].GetBinContent(b)
111  larger = max(abs(upshift), abs(dwshift))
112 
113  if upshift > 0 and dwshift > 0:
114  lhup[isys].SetBinContent(b, upshift)
115  lhdw[isys].SetBinContent(b, dwshift)
116  elif upshift <= 0 and dwshift <= 0: # up down flip
117  lhup[isys].SetBinContent(b, -dwshift)
118  lhdw[isys].SetBinContent(b, -upshift)
119  elif upshift > 0 and dwshift <= 0:
120  lhup[isys].SetBinContent(b, upshift)
121  else:
122  lhdw[isys].SetBinContent(b, dwshift)
123 
124  fu = 0 if cv <= 0 else larger/cv
125  lhfu[isys].SetBinContent(b, fu)
126 
127 
128 def do_single_shift_1bin(isys, idx):
129  '''
130  Each bin is initialized to 0 to allow one-sided error bars.
131  '''
132 
133  lhup[isys].SetBinContent(b, 0)
134  lhdw[isys].SetBinContent(b, 0)
135  shift = hother[idx].GetBinContent(b)-cv
136  if shift > 0:
137  lhup[isys].SetBinContent(b, shift)
138  else:
139  lhdw[isys].SetBinContent(b, -shift)
140 
141  fu = 0 if cv <= 0 else abs(shift)/cv
142  lhfu[isys].SetBinContent(b, fu)
143 
144 
145 def setUncer1D(hIn, ihist):
146  hIn.SetLineColor(col_fe[ihist])
147  hIn.SetLineWidth(3)
148  hIn.SetLineStyle(col_st[ihist])
149  hIn.SetStats(0)
150 
151 
152 if __name__ == '__main__':
153 
154  # command line arguments
155  # fake data / real data switch
156  parser = argparse.ArgumentParser()
157  parser.add_argument('-d', '--data', action='store_true', help='Switch for using data instead of fake data.')
158  args = parser.parse_args()
159  use_data = args.data
160 
161  # load files into file containers
162  # mc is the Monte Carlo histogram in the final plot
163  fgenie = []
164  fmc = TFile(os.path.join(infp, 'nominal_xsec.root'))
165  # Here I load (fake) data and all other systematics.
166  # It is an unforunate design to put data together with other systematics.
167  # However it is too late. Get over!
168  if not use_data:
169  fppfx_other = TFile(os.path.join(infp, 'ppfx_others_v2.root'))
170  else:
171  fppfx_other = TFile(os.path.join(infp, 'ppfx_others_data.root'))
172 
173  # historam containers
174  hgenie = []
175  hppfx = []
176  hother = []
177 
178  # get nominal cross section
179  hmc = fmc.Get('hXS1Dsec_nominal')
180  # get (fake) data cross section
181  hdt = fppfx_other.Get('xs_nominal/hXS1Dsec_cv')
182 
183  # store genie histograms
184  if not use_data:
185  hgenie += get_genie_universes('genie_0_99_v2.root')
186  hgenie += get_genie_universes('genie_100_299_v2.root')
187  hgenie += get_genie_universes('genie_300_499_v2.root')
188  else:
189  hgenie += get_genie_universes('genie_0_99_data.root')
190  hgenie += get_genie_universes('genie_100_299_data.root')
191  hgenie += get_genie_universes('genie_300_499_data.root')
192 
193  # store ppfx histograms
194  hppfx += get_ppfx_universes(fppfx_other)
195 
196  # store all other systematics
197  for syst_name in syst_names:
198  hother.append(fppfx_other.Get('xs_nominal/hXS1Dsec_{}'.format(syst_name)))
199 
200  # histogram containers: one for each systematics
201  nfu = len(fu_names)
202  lhup = [hmc.Clone('hUp_{}'.format(i)) for i in range(nfu)]
203  lhdw = [hmc.Clone('hDw_{}'.format(i)) for i in range(nfu)]
204  lhfu = [hmc.Clone('hFU_{}'.format(i)) for i in range(nfu)]
205 
206  # number of bins to loop through
207  Nbins = hdt.GetNbinsX()
208 
209  # loop through each bin
210  for b in range(1,Nbins+1):
211 
212  # central value of nominal MC
213  cv = hdt.GetBinContent(b)
214 
215  # ppfx
216  do_multiverse_1bin(0, [h.GetBinContent(b) for h in hppfx])
217 
218  # genie
219  do_multiverse_1bin(1, [h.GetBinContent(b) for h in hgenie])
220 
221  # light
222  do_up_down_1bin(2,1,0)
223 
224 
225  # calibration
226  do_up_down_1bin(3,4,3)
227 
228  # Cherenkov
230 
231  # calibshape
233 
234  # muon energy scale
235  do_up_down_1bin(6,7,6)
236 
237  # Normalization (TN, doc-):
238  # POT: 0.28% (Section 9.1)
239  # Intensity: 2% (Section 9.2)
240  # Particle propagation: 1-2%, so we take 2% (Section 9.5).
241  # fe_norm = sqrt(pow(0.0028,2)+pow(0.02,2)+pow(0.02,2));
242  fe_norm = 0.028422526
243  lhup[7].SetBinContent(b, abs(cv*fe_norm))
244  lhdw[7].SetBinContent(b, abs(cv*fe_norm))
245  fu = 0 if cv <= 0 else fe_norm
246  lhfu[7].SetBinContent(b, fu)
247 
248  # focusing
249  # first argument starting index, second number of up-down groups
250  do_focusing_1bin(8,8,11)
251 
252 
253  # make sure the output directory exists
254  out_subdir = 'real_data' if use_data else 'fake_data'
255  dest_path = os.path.join('output1d', out_subdir)
256  if not os.path.exists(dest_path):
257  os.makedirs(dest_path)
258 
259 
260  # ia = 8
261  # ib = 9
262  # lhup = lhup[ia:ib]
263  # lhdw = lhdw[ia:ib]
264  # # plot data
265  # gdt = TGraphAsymmErrors()
266  # set up output graph points
267  # for i in range(1,Nbins+1):
268  # # sum up and down errors in quadrature
269  # qsum_up = math.sqrt(sum([h.GetBinContent(i)**2 for h in lhup]))
270  # qsum_dw = math.sqrt(sum([h.GetBinContent(i)**2 for h in lhdw]))
271 
272  # binw = hdt.GetXaxis().GetBinWidth(i)
273  # gdt.SetPoint(i, hdt.GetXaxis().GetBinCenter(i), hdt.GetBinContent(i))
274  # gdt.SetPointError(i, binw/2, binw/2, qsum_dw, qsum_up)
275 
276  # gdt.Draw('AP1')
277 
278  # plot fractional uncertainty
279  cfu = TCanvas()
280  lfu = TLegend(0.2,0.5,0.45,0.88)
281  lhfu[0].GetXaxis().SetNdivisions(510)
282  lfu.SetBorderSize(0)
283  for i in range(len(lhfu)):
284  setUncer1D(lhfu[i],i)
285  if i == 0:
286  lhfu[i].Draw('hist')
287  lhfu[i].SetMaximum(1)
288  lhfu[i].GetXaxis().CenterTitle()
289  lhfu[i].GetXaxis().SetTitle('Neutrino Energy (GeV)')
290  lhfu[i].GetYaxis().CenterTitle()
291  lhfu[i].GetYaxis().SetTitle('Fractional Uncertainty')
292  else:
293  lhfu[i].Draw('histsame')
294  lfu.AddEntry(lhfu[i],final_fu[i],'L')
295  lfu.Draw()
296  # save fractional uncertainty to file
297  cfu.SaveAs('{}/fe1D.pdf'.format(dest_path))
298 
299  raw_input('Press enter to continue.')
tree Draw("slc.nhit")
h7 SetMaximum(maxY)
def do_multiverse_1bin(isys, univ_vals)
Definition: plot_xsec_1d.py:54
correl_xv GetYaxis() -> SetDecimals()
def get_ppfx_universes(fppfx_other)
Definition: plot_xsec_1d.py:44
void abs(TH1 *hist)
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
def get_genie_universes(basefn)
Definition: plot_xsec_1d.py:31
gargamelle SetTitle("Gargamelle #nu_{e} CC data")
correl_xv GetXaxis() -> SetDecimals()
def setUncer1D(hIn, ihist)
std::string format(const int32_t &value, const int &ndigits=8)
Definition: HexUtils.cpp:14
def do_focusing_1bin(isys, sstart, nsys)
Definition: plot_xsec_1d.py:71
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 do_up_down_1bin(isys, upidx, dwidx)
Definition: plot_xsec_1d.py:99
def do_single_shift_1bin(isys, idx)
Double_t sum
Definition: plot.C:31
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68