3 from __future__
import print_function
6 from bisect
import bisect
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' 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"]
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" 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]
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)' 53 fn = os.path.join(infp, basefn)
54 fgenie.append(TFile(fn))
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)))
67 for i
in range(0,100):
68 hl.append(fppfx_other.Get(
'xs_nominal/hXS2Dsec_ppfx{}'.
format(i)))
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]
83 for i
in range(1, NbinsX+1):
84 for j
in range(1, NbinsY+1):
85 dt_cv = hdt.GetBinContent(i,j)
93 binc_all_uni = [h.GetBinContent(i,j)
for h
in huniv]
98 conserv_error =
max(cup,cdown)
104 for i
in range(1, NbinsX+1):
105 for j
in range(1, NbinsY+1):
106 dt_cv = hdt.GetBinContent(i,j)
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))
116 if sigma_dw < 0
and sigma_up < 0:
118 elif sigma_dw > 0
and sigma_up > 0:
132 for i
in range(1, NbinsX+1):
133 for j
in range(1, NbinsY+1):
134 dt_cv = hdt.GetBinContent(i,j)
140 sigma = dt_cv - hother[hidx].GetBinContent(i,j)
157 for i
in range(1, NbinsX+1):
158 for j
in range(1, NbinsY+1):
159 dt_cv = hdt.GetBinContent(i,j)
161 sigma = dt_cv * fe_norm
175 print(
'Focusing systematics:')
176 for k
in range(nsys):
177 print(syst_names[sstart+k*2], syst_names[sstart+k*2+1])
179 for i
in range(1, NbinsX+1):
180 for j
in range(1, NbinsY+1):
181 dt_cv = hdt.GetBinContent(i,j)
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)))
193 qsum = math.sqrt(
sum([x**2
for x
in larger_sigma]))
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 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)
228 for p
in range(Nfinal_fe):
229 up = lhup[p].GetBinContent(i,j)
230 dw = lhdw[p].GetBinContent(i,j)
234 errUp = math.sqrt(errUp)
235 errDw = math.sqrt(errDw)
239 gdt[i-1].SetPoint(count[i-1], kin, cont)
240 gdt[i-1].SetPointError(count[i-1],binw/2, binw/2, errDw, errUp)
245 for i
in range(1, NbinsX+1):
246 for j
in range(1, NbinsY+1):
248 binw = hdt.GetYaxis().GetBinWidth(j)
249 cont = hmc.GetBinContent(i,j)
251 for p
in range(Nfinal_fe):
252 ferr = lhfe[p].GetBinContent(i,j)
261 gIn.SetMarkerColor(kBlack)
262 gIn.SetMarkerStyle(20)
263 gIn.SetLineColor(kBlack)
267 hIn.SetLineColor(col_fe[ihist])
269 hIn.SetLineStyle(col_st[ihist])
273 if __name__ ==
'__main__':
276 The unit of the output plots from this script has to be 277 rescaled by bin size when showing the final results! 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()
290 fmc = TFile(os.path.join(infp,
'nominal_xsec.root'))
292 fppfx_other = TFile(os.path.join(infp,
'ppfx_others_v2.root'))
294 fppfx_other = TFile(os.path.join(infp,
'ppfx_others_data.root'))
302 hmc = fmc.Get(
'hXS2Dsec_nominal')
304 hdt = fppfx_other.Get(
'xs_nominal/hXS2Dsec_cv')
320 for syst_name
in syst_names:
321 hother.append(fppfx_other.Get(
'xs_nominal/hXS2Dsec_{}'.
format(syst_name)))
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)))
334 NbinsX = hdt.GetNbinsX()
335 NbinsY = hdt.GetNbinsY()
366 fe_norm = 0.028422526
374 gdt = [TGraphAsymmErrors()
for i
in range(Nabins2X)]
375 count = [0
for i
in range(Nabins2X)]
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)]
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)
391 leg = [TLegend(0.67,0.68,0.87,0.88)
for i
in range(Nabins2X)]
392 lat = [TLatex()
for i
in range(Nabins2X)]
394 for i
in range(NbinsX):
397 cxs.append(TCanvas())
404 h1MC[i].
GetYaxis().
SetTitle(
'd^{2}#sigma/dT_{#mu}dcos#theta_{#mu} (cm^{2}/GeV/nucleon)')
407 gdt[i].
Draw(
"p1 same")
409 leg[i].
AddEntry(h1MC[i],
"MC",
"fl")
410 leg[i].
AddEntry(gdt[i],
"Data (Syst.)",
"lp")
413 cxs[i].
SaveAs(
'{}/xs2D_anglebin{}.png'.
format(dest_path,i))
414 cxs[i].
SaveAs(
'{}/xs2D_anglebin{}.pdf'.
format(dest_path,i))
418 latfe = [TLatex()
for i
in range(Nabins2X)]
420 for i
in range(NbinsX):
422 lfe = TLegend(0.2,0.5,0.45,0.88)
424 lfe = TLegend(0.6,0.4,0.85,0.78)
426 cfe.append(TCanvas())
428 for j
in range(Nfinal_fe):
435 h1fe[i][j].
Draw(
'hist')
437 h1fe[i][j].
Draw(
'histsame')
440 lfe.AddEntry(h1fe[i][j],final_fe[j],
'l')
443 cfe[i].
SaveAs(
'{}/fe2D_anglebin{}.png'.
format(dest_path,i))
444 cfe[i].
SaveAs(
'{}/fe2D_anglebin{}.pdf'.
format(dest_path,i))
447 fileOut = TFile(
'{}/numucc_inc_{}.root'.
format(dest_path,out_subdir),
'recreate')
449 fileOut.mkdir(
'CrossSections')
450 for i
in range(NbinsX):
451 fileOut.cd(
'CrossSections')
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]))
461 raw_input(
'Press enter to continue.')
def process_light_calibration_mue(isys, hidxdown, hidxup)
correl_xv GetYaxis() -> SetDecimals()
def process_focusing(isys, sstart, nsys)
prelim SetTextSize(2/30.)
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
def process_ckv_calibshape(isys, hidx)
gargamelle SetTitle("Gargamelle #nu_{e} CC data")
def process_normalization(isys)
correl_xv GetXaxis() -> SetDecimals()
leg AddEntry(GRdata,"data","p")
tex DrawLatex(28.4, 12.8,"x = 0.015(5*F2)")
std::string format(const int32_t &value, const int &ndigits=8)
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)
def process_multiuniverse(isys, huniv)
def get_genie_universes(basefn)
cosmicTree SaveAs("cosmicTree.root")
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
def up_down_counts_1sigma(count_list, cv)
def setUncer1D(hIn, ihist)