3 from __future__
import print_function
5 from bisect
import bisect
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' 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]
34 fn = os.path.join(infp, basefn)
35 fgenie.append(TFile(fn))
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)))
48 for i
in range(0,100):
49 hl.append(fppfx_other.Get(
'xs_nominal/hXS1Dsec_ppfx{}'.
format(i)))
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)
62 up_shift = sorted_list[up_pos] - cv
63 dw_shift = cv - sorted_list[dw_pos]
67 fu = 0
if cv <= 0
else max(up_shift, dw_shift)/cv
73 up_shifts = [0
for k
in range(nsys)]
74 dw_shifts = [0
for k
in range(nsys)]
76 dwshift = cv - hother[sstart+k*2].GetBinContent(i)
77 upshift = hother[sstart+k*2+1].GetBinContent(i) - cv
79 if upshift > 0
and dwshift > 0:
80 up_shifts[k] = upshift
81 dw_shifts[k] = dwshift
82 elif upshift <= 0
and dwshift <= 0:
83 up_shifts[k] = -dwshift
84 dw_shifts[k] = -upshift
85 elif upshift > 0
and dwshift <= 0:
86 up_shifts[k] = upshift
88 dw_shifts[k] = dwshift
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]))
95 fu = 0
if cv <= 0
else max(up_shift, dw_shift)/cv
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. 109 upshift = hother[upidx].GetBinContent(b)-cv
110 dwshift = cv-hother[dwidx].GetBinContent(b)
113 if upshift > 0
and dwshift > 0:
116 elif upshift <= 0
and dwshift <= 0:
119 elif upshift > 0
and dwshift <= 0:
124 fu = 0
if cv <= 0
else larger/cv
130 Each bin is initialized to 0 to allow one-sided error bars. 135 shift = hother[idx].GetBinContent(b)-cv
141 fu = 0
if cv <= 0
else abs(shift)/cv
146 hIn.SetLineColor(col_fe[ihist])
148 hIn.SetLineStyle(col_st[ihist])
152 if __name__ ==
'__main__':
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()
164 fmc = TFile(os.path.join(infp,
'nominal_xsec.root'))
169 fppfx_other = TFile(os.path.join(infp,
'ppfx_others_v2.root'))
171 fppfx_other = TFile(os.path.join(infp,
'ppfx_others_data.root'))
179 hmc = fmc.Get(
'hXS1Dsec_nominal')
181 hdt = fppfx_other.Get(
'xs_nominal/hXS1Dsec_cv')
197 for syst_name
in syst_names:
198 hother.append(fppfx_other.Get(
'xs_nominal/hXS1Dsec_{}'.
format(syst_name)))
207 Nbins = hdt.GetNbinsX()
210 for b
in range(1,Nbins+1):
213 cv = hdt.GetBinContent(b)
242 fe_norm = 0.028422526
245 fu = 0
if cv <= 0
else fe_norm
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)
280 lfu = TLegend(0.2,0.5,0.45,0.88)
281 lhfu[0].
GetXaxis().SetNdivisions(510)
283 for i
in range(len(lhfu)):
293 lhfu[i].
Draw(
'histsame')
294 lfu.AddEntry(lhfu[i],final_fu[i],
'L')
297 cfu.SaveAs(
'{}/fe1D.pdf'.
format(dest_path))
299 raw_input(
'Press enter to continue.')
def do_multiverse_1bin(isys, univ_vals)
correl_xv GetYaxis() -> SetDecimals()
def get_ppfx_universes(fppfx_other)
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
def get_genie_universes(basefn)
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)
def do_focusing_1bin(isys, sstart, nsys)
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)
def do_single_shift_1bin(isys, idx)
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)