Johan Bregeon / ctools

Commit 4cc94b92383e48ceca76648fb3e19085724209a9

Authored by Jürgen Knödlseder
1 parent 169cd922

Implement background rate replacement by a power law above given energy threshold.

Inline Side-by-side
Showing 2 changed files with 265 additions and 42 deletions
cscripts/csroot2caldb.par
  View file @4cc94b9
@@ -29,6 +29,7 @@ rebin, b, h, no,,, "Rebin" @@ -29,6 +29,7 @@ rebin, b, h, no,,, "Rebin"
29 eascale, r, h, 1.0,,, "Scaling factor for effective areas" 29 eascale, r, h, 1.0,,, "Scaling factor for effective areas"
30 bgdscale, r, h, 1.0,,, "Scaling factor for background rates" 30 bgdscale, r, h, 1.0,,, "Scaling factor for background rates"
31 bgdoversample, i, h, 1,,, "Spatial oversampling factor for background template" 31 bgdoversample, i, h, 1,,, "Spatial oversampling factor for background template"
  32 +bgdethres, r, h, 1000.0,,, "Energy above which to replace background by power law (TeV)"
32 bgdinfill, b, h, no,,, "Infill of background template?" 33 bgdinfill, b, h, no,,, "Infill of background template?"
33 34
34 # 35 #
cscripts/csroot2caldb.py
  View file @4cc94b9
@@ -571,6 +571,11 @@ class csroot2caldb(ctools.cscript): @@ -571,6 +571,11 @@ class csroot2caldb(ctools.cscript):
571 # Continue only if effective area HDU has been opened 571 # Continue only if effective area HDU has been opened
572 if self._hdu_ea != None: 572 if self._hdu_ea != None:
573 573
  574 + # Write header
  575 + if self._logTerse():
  576 + self._log('\n')
  577 + self._log.header1('Generate effective area extension')
  578 +
574 # Get User parameters 579 # Get User parameters
575 norm1d = self['norm1d'].boolean() 580 norm1d = self['norm1d'].boolean()
576 581
@@ -642,8 +647,13 @@ class csroot2caldb(ctools.cscript): @@ -642,8 +647,13 @@ class csroot2caldb(ctools.cscript):
642 # Continue only if point spread function HDU has been opened 647 # Continue only if point spread function HDU has been opened
643 if self._hdu_psf != None: 648 if self._hdu_psf != None:
644 649
  650 + # Write header
  651 + if self._logTerse():
  652 + self._log('\n')
  653 + self._log.header1('Generate point spread function extension')
  654 +
645 # King profile PSF 655 # King profile PSF
646 - if psftype == "King": 656 + if psftype == 'King':
647 self._root2psf_king(file) 657 self._root2psf_king(file)
648 else: 658 else:
649 self._root2psf_gauss(file) 659 self._root2psf_gauss(file)
@@ -874,19 +884,24 @@ class csroot2caldb(ctools.cscript): @@ -874,19 +884,24 @@ class csroot2caldb(ctools.cscript):
874 # Continue only if energy dispersion HDU has been opened 884 # Continue only if energy dispersion HDU has been opened
875 if self._hdu_edisp != None: 885 if self._hdu_edisp != None:
876 886
  887 + # Write header
  888 + if self._logTerse():
  889 + self._log('\n')
  890 + self._log.header1('Generate energy dispersion extension')
  891 +
877 # Allocate ROOT 3D array 892 # Allocate ROOT 3D array
878 matrix = TH3F() 893 matrix = TH3F()
879 - file.GetObject("EestOverEtrue_offaxis", matrix) 894 + file.GetObject('EestOverEtrue_offaxis', matrix)
880 895
881 # Set boundaries 896 # Set boundaries
882 - bounds = self._make_3D_migra(matrix, self._hdu_edisp, None, "") 897 + bounds = self._make_3D_migra(matrix, self._hdu_edisp, None, '')
883 for b in bounds: 898 for b in bounds:
884 self._edisp_bounds.append(b) 899 self._edisp_bounds.append(b)
885 self._set_cif_keywords(self._hdu_edisp, self._edisp_name, \ 900 self._set_cif_keywords(self._hdu_edisp, self._edisp_name, \
886 self._edisp_bounds, self._edisp_desc) 901 self._edisp_bounds, self._edisp_desc)
887 902
888 # MATRIX 903 # MATRIX
889 - self._make_3D_migra(matrix, self._hdu_edisp, "MATRIX", "") 904 + self._make_3D_migra(matrix, self._hdu_edisp, 'MATRIX', '')
890 905
891 # Return 906 # Return
892 return 907 return
@@ -907,22 +922,27 @@ class csroot2caldb(ctools.cscript): @@ -907,22 +922,27 @@ class csroot2caldb(ctools.cscript):
907 # Continue only if background HDU has been opened 922 # Continue only if background HDU has been opened
908 if self._hdu_bgd != None: 923 if self._hdu_bgd != None:
909 924
  925 + # Write header
  926 + if self._logTerse():
  927 + self._log('\n')
  928 + self._log.header1('Generate 2D background extension')
  929 +
910 # Allocate ROOT 2D array 930 # Allocate ROOT 2D array
911 array = TH2F() 931 array = TH2F()
912 - file.GetObject("BGRatePerSqDeg_offaxis", array) 932 + file.GetObject('BGRatePerSqDeg_offaxis', array)
913 933
914 # Set boundaries 934 # Set boundaries
915 - bounds = self._make_2D(array, self._hdu_bgd, None, "deg") 935 + bounds = self._make_2D(array, self._hdu_bgd, None, 'deg')
916 for b in bounds: 936 for b in bounds:
917 self._bgd_bounds.append(b) 937 self._bgd_bounds.append(b)
918 self._set_cif_keywords(self._hdu_bgd, self._bgd_name, \ 938 self._set_cif_keywords(self._hdu_bgd, self._bgd_name, \
919 self._bgd_bounds, self._bgd_desc) 939 self._bgd_bounds, self._bgd_desc)
920 940
921 # BGD 941 # BGD
922 - self._make_2D(array, self._hdu_bgd, "BGD", "") 942 + self._make_2D(array, self._hdu_bgd, 'BGD', '')
923 943
924 # BGD_RECO 944 # BGD_RECO
925 - self._make_2D(array, self._hdu_bgd, "BGD_RECO", "") 945 + self._make_2D(array, self._hdu_bgd, 'BGD_RECO', '')
926 946
927 # Return 947 # Return
928 return 948 return
@@ -942,14 +962,28 @@ class csroot2caldb(ctools.cscript): @@ -942,14 +962,28 @@ class csroot2caldb(ctools.cscript):
942 # Continue only if background HDU has been opened 962 # Continue only if background HDU has been opened
943 if self._hdu_bgd != None: 963 if self._hdu_bgd != None:
944 964
  965 + # Write header
  966 + if self._logTerse():
  967 + self._log('\n')
  968 + self._log.header1('Generate 3D background extension')
  969 +
945 # Get User parameters 970 # Get User parameters
946 - norm1d = self['norm1d'].boolean() 971 + norm1d = self['norm1d'].boolean()
  972 + bgdethres = self['bgdethres'].real()
  973 + bgdinfill = self['bgdinfill'].boolean()
947 974
948 # Allocate ROOT 2D array 975 # Allocate ROOT 2D array
949 array = TH2F() 976 array = TH2F()
950 file.GetObject('BGRatePerSqDeg_offaxis', array) 977 file.GetObject('BGRatePerSqDeg_offaxis', array)
951 978
952 - # Normalize 2D histogram on on-axis 1D histogram 979 + # Replace 2D histogram values by power law extrapolation
  980 + self._plaw_replace(array, bgdethres)
  981 +
  982 + # Optionally fill in empty values in 2D histogram
  983 + if bgdinfill:
  984 + self._infill_bgd(array)
  985 +
  986 + # Optionally normalize 2D histogram on on-axis 1D histogram
953 if norm1d: 987 if norm1d:
954 array_1D = TH1F() 988 array_1D = TH1F()
955 file.GetObject('BGRatePerSqDeg', array_1D) 989 file.GetObject('BGRatePerSqDeg', array_1D)
@@ -1033,14 +1067,17 @@ class csroot2caldb(ctools.cscript): @@ -1033,14 +1067,17 @@ class csroot2caldb(ctools.cscript):
1033 if neng != hist1D.GetXaxis().GetNbins(): 1067 if neng != hist1D.GetXaxis().GetNbins():
1034 # Write header 1068 # Write header
1035 if self._logTerse(): 1069 if self._logTerse():
1036 - self._log('Difference in energy binning between 2D histogram'  
1037 - ' "%s" and 1D histogram "%s".\n' %  
1038 - (hist2D.GetName(), hist1D.GetName())) 1070 + self._log(gammalib.parformat('On-axis normalisation'))
  1071 + self._log('Impossible since energy binning of 1D histogram '
  1072 + '"%s" does not match that of 2D histogram "%s".\n' %
  1073 + (hist1D.GetName(), hist2D.GetName()))
1039 else: 1074 else:
1040 # Write header 1075 # Write header
1041 if self._logTerse(): 1076 if self._logTerse():
1042 - self._log('Renormalize 2D histogram "%s" on 1D histogram "%s".\n' %  
1043 - (hist2D.GetName(), hist1D.GetName())) 1077 + self._log(gammalib.parformat('On-axis normalisation'))
  1078 + self._log('On-axis values of 1D histogram "%s" imposed on '
  1079 + '2D histogram "%s".\n' %
  1080 + (hist1D.GetName(), hist2D.GetName()))
1044 1081
1045 # Get a copy of 2D histogram 1082 # Get a copy of 2D histogram
1046 hist2D_copy = hist2D.Clone() 1083 hist2D_copy = hist2D.Clone()
@@ -1087,6 +1124,218 @@ class csroot2caldb(ctools.cscript): @@ -1087,6 +1124,218 @@ class csroot2caldb(ctools.cscript):
1087 # Return 1124 # Return
1088 return 1125 return
1089 1126
  1127 + def _plaw_replace(self, hist2D, ethres):
  1128 + """
  1129 + Replace background rate values by power law values
  1130 +
  1131 + Parameters
  1132 + ----------
  1133 + hist2D : `~ROOT.TH2F`
  1134 + ROOT 2D histogram.
  1135 + ethres : float
  1136 + Energy threshold.
  1137 + """
  1138 + # Extract energy and offset angle vectors
  1139 + energies = hist2D.GetXaxis()
  1140 + neng = energies.GetNbins()
  1141 + offsets = hist2D.GetYaxis()
  1142 + noffset = offsets.GetNbins()
  1143 +
  1144 + # Continue only if threshold is low enough
  1145 + if ethres < math.pow(10.0, energies.GetBinUpEdge(neng)):
  1146 +
  1147 + # Determine mean energy values
  1148 + engs = []
  1149 + for ieng in range(neng):
  1150 + energy = 0.5 * (math.pow(10.0, energies.GetBinLowEdge(ieng+1)) +
  1151 + math.pow(10.0, energies.GetBinUpEdge(ieng+1)))
  1152 + engs.append(energy)
  1153 +
  1154 + # Loop over all offaxis angles
  1155 + for ioff in range(noffset):
  1156 +
  1157 + # Determine offset angle (for logging)
  1158 + offset = 0.5 * (offsets.GetBinLowEdge(ioff+1) +
  1159 + offsets.GetBinUpEdge(ioff+1))
  1160 +
  1161 + # Initialise vectors
  1162 + axis = []
  1163 + rates = []
  1164 +
  1165 + # Extract values
  1166 + for ieng in range(neng):
  1167 + energy = engs[ieng]
  1168 + rate = hist2D.GetBinContent(ieng+1,ioff+1)
  1169 + axis.append(energy)
  1170 + rates.append(rate)
  1171 +
  1172 + # Determine energy range
  1173 + emin, emax = self._plaw_energy_range(axis, rates)
  1174 +
  1175 + # Get power law coefficients
  1176 + coeff = self._plaw_coefficients(axis, rates, emin, ethres)
  1177 +
  1178 + # Replace histogram values by power law values
  1179 + for ieng in range(neng):
  1180 + energy = engs[ieng]
  1181 + if energy > ethres:
  1182 + value = self._plaw_value(coeff, energy)
  1183 + hist2D.SetBinContent(ieng+1,ioff+1,value)
  1184 + if self._logTerse():
  1185 + self._log(gammalib.parformat('Power law replacement'))
  1186 + self._log('%e cts/s/deg2 (E=%8.4f TeV (%d) '
  1187 + 'Off=%.2f deg (%d)\n' %
  1188 + (value,engs[ieng],ieng,offset,ioff))
  1189 +
  1190 + # Return
  1191 + return
  1192 +
  1193 + def _plaw_energy_range(self, energies, rates):
  1194 + """
  1195 + Determine the energy range over which the power law is valid
  1196 +
  1197 + Parameters
  1198 + ----------
  1199 + energies : list of float
  1200 + Energy values (TeV)
  1201 + rates : list of float
  1202 + Background rate values
  1203 + """
  1204 + # Find energy at which background rate takes a maximum
  1205 + min_energy = 0.0
  1206 + max_rate = 0.0
  1207 + for ieng, energy in enumerate(energies):
  1208 + if rates[ieng] > max_rate:
  1209 + max_rate = rates[ieng]
  1210 + min_energy = energy
  1211 +
  1212 + # Find last data point
  1213 + max_energy = 0.0
  1214 + for ieng in range(len(energies)-1,-1,-1):
  1215 + if rates[ieng] > 0.0:
  1216 + max_energy = energies[ieng]
  1217 + break
  1218 +
  1219 + # Return
  1220 + return (min_energy, max_energy)
  1221 +
  1222 + def _plaw_coefficients(self, energies, rates, emin, emax):
  1223 + """
  1224 + Determines the power law coefficients
  1225 +
  1226 + Parameters
  1227 + ----------
  1228 + energies : list of float
  1229 + Energy values (TeV)
  1230 + rates : list of float
  1231 + Background rate values
  1232 + emin : float
  1233 + Minimum energy (TeV)
  1234 + emax : float
  1235 + Maximum energy (TeV)
  1236 + """
  1237 + # Determine nearest energy indices
  1238 + i_emin = 0
  1239 + i_emax = 0
  1240 + delta_emin = 1.0e30
  1241 + delta_emax = 1.0e30
  1242 + for ieng, energy in enumerate(energies):
  1243 + if abs(energy-emin) < delta_emin:
  1244 + delta_emin = abs(energy-emin)
  1245 + i_emin = ieng
  1246 + if abs(energy-emax) < delta_emax:
  1247 + delta_emax = abs(energy-emax)
  1248 + i_emax = ieng
  1249 +
  1250 + # Determine coefficients of power law
  1251 + x1 = math.log10(energies[i_emin])
  1252 + x2 = math.log10(energies[i_emax])
  1253 + y1 = math.log10(rates[i_emin])
  1254 + y2 = math.log10(rates[i_emax])
  1255 + m = (y2-y1)/(x2-x1)
  1256 + t = y1 - m * x1
  1257 +
  1258 + # Return coefficients
  1259 + return {'m': m, 't': t}
  1260 +
  1261 + def _plaw_value(self, coeff, energy):
  1262 + """
  1263 + Determine the power law value
  1264 +
  1265 + Parameters
  1266 + ----------
  1267 + coeff : dict
  1268 + Dictionary of coefficients
  1269 + energy : float
  1270 + Energy value (TeV)
  1271 + """
  1272 + # Compute value
  1273 + value = math.pow(10.0, coeff['m'] * math.log10(energy) + coeff['t'])
  1274 +
  1275 + # Return value
  1276 + return value
  1277 +
  1278 + def _infill_bgd(self, hist2D):
  1279 + """
  1280 + Fill all empty bins in 2D histogram by preceeding values in energy
  1281 +
  1282 + Parameters
  1283 + ----------
  1284 + hist2D : `~ROOT.TH2F`
  1285 + ROOT 2D histogram.
  1286 + """
  1287 + # Extract energy and offset angle vectors
  1288 + energies = hist2D.GetXaxis()
  1289 + neng = energies.GetNbins()
  1290 + offsets = hist2D.GetYaxis()
  1291 + noffset = offsets.GetNbins()
  1292 +
  1293 + # Determine mean energy values (for logging)
  1294 + engs = []
  1295 + for ieng in range(neng):
  1296 + energy = 0.5 * (math.pow(10.0, energies.GetBinLowEdge(ieng+1)) +
  1297 + math.pow(10.0, energies.GetBinUpEdge(ieng+1)))
  1298 + engs.append(energy)
  1299 +
  1300 + # Loop over all offaxis angles
  1301 + for ioff in range(noffset):
  1302 +
  1303 + # Determine offset angle (for logging)
  1304 + offset = 0.5 * (offsets.GetBinLowEdge(ioff+1) +
  1305 + offsets.GetBinUpEdge(ioff+1))
  1306 +
  1307 + # Initialise indices and value
  1308 + i_start = -1
  1309 + i_stop = -1
  1310 + value = 0.0
  1311 +
  1312 + # First start and stop indices for infill
  1313 + for ieng in range(neng):
  1314 + if hist2D.GetBinContent(ieng+1,ioff+1) > 0.0:
  1315 + i_start = ieng
  1316 + value = hist2D.GetBinContent(ieng+1,ioff+1)
  1317 + break
  1318 + for ieng in range(neng-1, -1, -1):
  1319 + if hist2D.GetBinContent(ieng+1,ioff+1) > 0.0:
  1320 + i_stop = ieng
  1321 + break
  1322 +
  1323 + # If indices are valid then fill in background rates
  1324 + if i_start > -1 and i_stop > -1:
  1325 + for ieng in range(i_start, i_stop+1):
  1326 + if hist2D.GetBinContent(ieng+1,ioff+1) == 0.0:
  1327 + hist2D.SetBinContent(ieng+1,ioff+1,value)
  1328 + if self._logTerse():
  1329 + self._log(gammalib.parformat('Background infill'))
  1330 + self._log('%e cts/s/deg2 (E=%8.4f TeV (%d) '
  1331 + 'Off=%.2f deg (%d)\n' %
  1332 + (value,engs[ieng],ieng,offset,ioff))
  1333 + else:
  1334 + value = hist2D.GetBinContent(ieng+1,ioff+1)
  1335 +
  1336 + # Return
  1337 + return
  1338 +
1090 def _make_2D(self, array, hdu, name, unit, scale=1.0): 1339 def _make_2D(self, array, hdu, name, unit, scale=1.0):
1091 """ 1340 """
1092 Make 2D matrix as function of energy and offset angle from a 1341 Make 2D matrix as function of energy and offset angle from a
@@ -1162,7 +1411,6 @@ class csroot2caldb(ctools.cscript): @@ -1162,7 +1411,6 @@ class csroot2caldb(ctools.cscript):
1162 # Get User parameters 1411 # Get User parameters
1163 scale = self['bgdscale'].real() 1412 scale = self['bgdscale'].real()
1164 oversample = self['bgdoversample'].integer() 1413 oversample = self['bgdoversample'].integer()
1165 - infill = self['bgdinfill'].boolean()  
1166 1414
1167 # Set constants 1415 # Set constants
1168 deg2sr = 0.01745329*0.01745329 1416 deg2sr = 0.01745329*0.01745329
@@ -1178,32 +1426,6 @@ class csroot2caldb(ctools.cscript): @@ -1178,32 +1426,6 @@ class csroot2caldb(ctools.cscript):
1178 ewidth.append(pow(10.0, energies.GetBinUpEdge(ieng+1) +6.0) - 1426 ewidth.append(pow(10.0, energies.GetBinUpEdge(ieng+1) +6.0) -
1179 pow(10.0, energies.GetBinLowEdge(ieng+1)+6.0)) 1427 pow(10.0, energies.GetBinLowEdge(ieng+1)+6.0))
1180 1428
1181 - # Optionally fill all empty bins by the preceding value that occur  
1182 - # for all offset angles between the first and the last non zero  
1183 - # bin  
1184 - if infill:  
1185 - for ioff in range(noffset):  
1186 - i_start = -1  
1187 - i_stop = -1  
1188 - value = 0.0  
1189 - for ieng in range(neng):  
1190 - if array.GetBinContent(ieng+1,ioff+1) > 0.0:  
1191 - i_start = ieng  
1192 - value = array.GetBinContent(ieng+1,ioff+1)  
1193 - break  
1194 - for ieng in range(neng-1, -1, -1):  
1195 - if array.GetBinContent(ieng+1,ioff+1) > 0.0:  
1196 - i_stop = ieng  
1197 - break  
1198 - if i_start > -1 and i_stop > -1:  
1199 - for ieng in range(i_start, i_stop+1):  
1200 - if array.GetBinContent(ieng+1,ioff+1) == 0.0:  
1201 - array.SetBinContent(ieng+1,ioff+1,value)  
1202 - print('Background infill ieng=%d ioff=%d value=%f' %  
1203 - (ieng,ioff,value))  
1204 - else:  
1205 - value = array.GetBinContent(ieng+1,ioff+1)  
1206 -  
1207 # Build DETX and DETY axes 1429 # Build DETX and DETY axes
1208 ndet = array.GetYaxis().GetNbins() 1430 ndet = array.GetYaxis().GetNbins()
1209 ndets = 2*ndet*oversample 1431 ndets = 2*ndet*oversample