Belle II Software light-2406-ragdoll
analysis-lookuptable-creation-photon-efficiency-datamc-ratio.py
1#!/usr/bin/env python3
2
3
10
11import glob
12import numpy as np
13import argparse
14
15import basf2 as b2
16
17parser = argparse.ArgumentParser()
18parser.add_argument(
19 '-i',
20 '--input_path',
21 default="/afs/desy.de/user/h/hsvidras/repositories/photoneff/analysis/plots_lima/lima/",
22 help='Path to folder with efficiency ratio files generated by the photon efficiency workflow.')
23args = parser.parse_args()
24
25# Define bin constructors
26
27
28def make_1D_bin(name, min_val, max_val):
29 return {name: [min_val, max_val]}
30
31
32def make_2D_bin(bin_x, bin_y):
33 bin_2d = bin_x.copy()
34 bin_2d.update(bin_y)
35 return bin_2d
36
37
38def make_3D_bin(bin_x, bin_y, bin_z):
39 bin_3d = bin_x.copy()
40 bin_3d.update(bin_y)
41 bin_3d.update(bin_z)
42 return bin_3d
43
44# To make these one needs tables generated with the photon efficiency code
45# REPOSITORY: https://gitlab.desy.de/belle2/performance/neutrals/photon-detection-efficiency
46# They provide you with 2D theta phi tables, for different photon energy bins.
47# The energy bins are defined in the file names, but the theta, phi bins are not, so they are given in separate files.
48# The naming scheme has to be exact for this to work, but luckily the code in the repository takes care of all that.
49# So you only need to input the location where you want the code to look for these tables.
50# Essentially files are create for edges named
51# BINS: table_pRecoilPhipRecoilTheta_pRecoilfrom0p2andpRecoilto0p4_pRecoil[Theta,Phi]bins.npy
52# RATIOS: table_pRecoilPhipRecoilTheta*[values, statistical_up, statistical_down, systematic_up,systematic_down].npy
53# , where * stands for bins in the format of: from0p2andpRecoilto0p4 (meaning 0.2<pRecoil<0.4)
54# All this will be read automatically, but the naming scheme has to be exact.
55
56
57table_location = args.input_path
58
59# Define bin ranges. Bins may be of different size
60yedges = np.load(f"{table_location}/table_pRecoilPhipRecoilTheta_pRecoilfrom0p2andpRecoilto0p4_pRecoilPhibins.npy")
61xedges = np.load(f"{table_location}/table_pRecoilPhipRecoilTheta_pRecoilfrom0p2andpRecoilto0p4_pRecoilThetabins.npy")
62
63# Define values to put into the table, that correspond to previously defined bins.
64values = glob.glob(f"{table_location}/table_pRecoilPhipRecoilTheta*values.npy")
65stat_up = glob.glob(f"{table_location}/table_pRecoilPhipRecoilTheta*statistical_up.npy")
66stat_down = glob.glob(f"{table_location}/table_pRecoilPhipRecoilTheta*statistical_down.npy")
67sys_up = glob.glob(f"{table_location}/table_pRecoilPhipRecoilTheta*systematic_up.npy")
68sys_down = glob.glob(f"{table_location}/table_pRecoilPhipRecoilTheta*systematic_down.npy")
69
70bins_p = []
71bins_phi = [make_1D_bin("phi", lowbin, highbin) for lowbin, highbin in zip(yedges[:-1], yedges[1:])]
72bins_theta = [make_1D_bin("theta", lowbin, highbin) for lowbin, highbin in zip(xedges[:-1], xedges[1:])]
73
74table = []
75
76for n, valuenp in enumerate(values):
77
78 bins_with_text = valuenp.split('pRecoilfrom')[1].split('pRecoilto')
79 lowbin_str = bins_with_text[0][:-3].replace('p', '.')
80 highbin_str = bins_with_text[1][:-11].replace('p', '.')
81 pbin = make_1D_bin("E", float(lowbin_str), float(highbin_str))
82
83 value = np.load(valuenp)
84 stat_err_up_table = np.load(stat_up[n])
85 stat_err_down_table = np.load(stat_down[n])
86 sys_err_up_table = np.load(sys_up[n])
87 sys_err_down_table = np.load(sys_down[n])
88 for j, ybin in enumerate(bins_phi):
89 for i, xbin in enumerate(bins_theta):
90 weightInfo = {}
91 weightInfo["Weight"] = value[i, j]
92 weightInfo["StatErrUp"] = stat_err_up_table[i, j]
93 weightInfo["StatErrDown"] = stat_err_down_table[i, j]
94 weightInfo["SystErrUp"] = sys_err_up_table[i, j]
95 weightInfo["SystErrDown"] = sys_err_down_table[i, j]
96 weightInfo["TotalErrUp"] = (sys_err_up_table[i, j]**2 + stat_err_up_table[i, j]**2)**0.5
97 weightInfo["TotalErrDown"] = (sys_err_down_table[i, j]**2 + stat_err_down_table[i, j]**2)**0.5
98 table.append([weightInfo, make_3D_bin(xbin, ybin, pbin)])
99
100# And of course let's define out-of-range bin info
101outOfRangeWeightInfo = {}
102outOfRangeWeightInfo["Weight"] = np.nan
103outOfRangeWeightInfo["StatErrUp"] = np.nan
104outOfRangeWeightInfo["SystErrUp"] = np.nan
105outOfRangeWeightInfo["StatErrDown"] = np.nan
106outOfRangeWeightInfo["SystErrDown"] = np.nan
107outOfRangeWeightInfo["TotalErrDown"] = np.nan
108outOfRangeWeightInfo["TotalErrDown"] = np.nan
109
110# Now, let's configure table creator
111addtable = b2.register_module('ParticleWeightingLookUpCreator')
112addtable.param('tableIDNotSpec', table)
113addtable.param('outOfRangeWeight', outOfRangeWeightInfo)
114addtable.param('experimentHigh', -1)
115addtable.param('experimentLow', 0)
116addtable.param('runHigh', -1)
117addtable.param('runLow', 0)
118addtable.param('tableName', "PhotonEfficiencyDataMCRatio_June2021")
119
120eventinfosetter = b2.register_module('EventInfoSetter')
121eventinfosetter.param('evtNumList', [10])
122eventinfosetter.param('runList', [0])
123eventinfosetter.param('expList', [0])
124
125my_path = b2.create_path()
126my_path.add_module(addtable)
127my_path.add_module(eventinfosetter)
128
129b2.process(my_path)