jean perrot / ctools

Commit 65346e200c814f540d919d77b73f9134ea0d9850

Authored by Jürgen Knödlseder
1 parent 698513d6

Now creates a 3D background file.

Inline Side-by-side
Showing 1 changed file with 206 additions and 30 deletions
scripts/cta_root2caldb.py
  View file @65346e2
... ... @@ -650,7 +650,6 @@ class caldb():
650 650 for ioff in range(noffset):
651 651 o_lo = offsets.GetBinLowEdge(ioff+1)
652 652 hdu["THETA_LO"][0,ioff] = o_lo
653   -
654 653 #
655 654 # THETA_LO
656 655 if not hdu.contains("THETA_HI"):
... ... @@ -683,8 +682,144 @@ class caldb():
683 682  
684 683 # Return boundary information
685 684 return bounds
  685 +
  686 + def make_3D(self, array, hdu, name, unit, scale=1.0):
  687 + """
  688 + Make 3D cube as function of DETX, DETY and energy from a
  689 + ROOT 2D histogram. If the HDU has already the energy and
  690 + offset angle columns, this method will simply add another data
  691 + column.
  692 + If name==None, the method will not append any data column.
  693 +
  694 + Parameters:
  695 + array - ROOT 2D histogram.
  696 + hdu - FITS HDU.
  697 + name - Data column name.
  698 + unit - Data unit.
  699 + Keywords:
  700 + scale - Scaling factor for histogram values.
  701 + """
  702 + # Set constants
  703 + deg2sr = 0.01745329*0.01745329
  704 +
  705 + # Extract energy and offset angle vectors
  706 + energies = array.GetXaxis()
  707 + neng = energies.GetNbins()
  708 + offsets = array.GetYaxis()
  709 + noffset = offsets.GetNbins()
  710 + ewidth = [] # in MeV
  711 + for ieng in range(neng):
  712 + ewidth.append(pow(10.0, energies.GetBinUpEdge(ieng+1) +6.0) - \
  713 + pow(10.0, energies.GetBinLowEdge(ieng+1)+6.0))
  714 +
  715 + # Build DETX and DETY axes
  716 + ndet = array.GetYaxis().GetNbins()
  717 + ndets = 2*ndet
  718 + det_max = array.GetYaxis().GetBinUpEdge(ndet)
  719 + det_bin = det_max/float(ndet)
  720 + dets_lo = []
  721 + dets_hi = []
  722 + dets2 = []
  723 + for i in range(ndets):
  724 + det_lo = -det_max + i*det_bin
  725 + det_hi = det_lo + det_bin
  726 + det_val = 0.5*(det_lo + det_hi)
  727 + dets_lo.append(det_lo)
  728 + dets_hi.append(det_hi)
  729 + dets2.append(det_val*det_val)
  730 +
  731 + # Attach columns to HDU. We do this only if they do not yet
  732 + # exist. This allows adding columns to the same HDU in successive
  733 + # calls of this method
  734 + #
  735 + # DETX_LO
  736 + if not hdu.contains("DETX_LO"):
  737 + hdu.append(GFitsTableFloatCol("DETX_LO", 1, ndets))
  738 + hdu["DETX_LO"].unit("deg")
  739 + for i in range(ndets):
  740 + hdu["DETX_LO"][0,i] = dets_lo[i]
  741 + #
  742 + # DETX_HI
  743 + if not hdu.contains("DETX_HI"):
  744 + hdu.append(GFitsTableFloatCol("DETX_HI", 1, ndets))
  745 + hdu["DETX_HI"].unit("deg")
  746 + for i in range(ndets):
  747 + hdu["DETX_HI"][0,i] = dets_hi[i]
  748 + #
  749 + # DETY_LO
  750 + if not hdu.contains("DETY_LO"):
  751 + hdu.append(GFitsTableFloatCol("DETY_LO", 1, ndets))
  752 + hdu["DETY_LO"].unit("deg")
  753 + for i in range(ndets):
  754 + hdu["DETY_LO"][0,i] = dets_lo[i]
  755 + #
  756 + # DETY_HI
  757 + if not hdu.contains("DETY_HI"):
  758 + hdu.append(GFitsTableFloatCol("DETY_HI", 1, ndets))
  759 + hdu["DETY_HI"].unit("deg")
  760 + for i in range(ndets):
  761 + hdu["DETY_HI"][0,i] = dets_hi[i]
  762 + #
  763 + # ENERG_LO
  764 + if not hdu.contains("ENERG_LO"):
  765 + hdu.append(GFitsTableFloatCol("ENERG_LO", 1, neng))
  766 + hdu["ENERG_LO"].unit("TeV")
  767 + for ieng in range(neng):
  768 + e_lo = pow(10.0, energies.GetBinLowEdge(ieng+1))
  769 + hdu["ENERG_LO"][0,ieng] = e_lo
  770 + #
  771 + # ENERG_HI
  772 + if not hdu.contains("ENERG_HI"):
  773 + hdu.append(GFitsTableFloatCol("ENERG_HI", 1, neng))
  774 + hdu["ENERG_HI"].unit("TeV")
  775 + for ieng in range(neng):
  776 + e_hi = pow(10.0, energies.GetBinUpEdge(ieng+1))
  777 + hdu["ENERG_HI"][0,ieng] = e_hi
  778 + #
  779 + # "NAME"
  780 + if name != None and not hdu.contains(name):
  781 + hdu.append(GFitsTableFloatCol(name, 1, ndets*ndets*neng))
  782 + hdu[name].unit(unit)
  783 + hdu[name].dim([ndets,ndets,neng])
  784 + for ix in range(ndets):
  785 + for iy in range(ndets):
  786 + for ieng in range(neng):
  787 + index = ix + (iy + ieng * ndets) * ndets
  788 + theta = math.sqrt(dets2[ix] + dets2[iy])
  789 + ioff = offsets.FindBin(theta)
  790 + if ioff > 0 and ioff <= noffset:
  791 + binsq = array.GetBinContent(ieng+1,ioff)
  792 + if binsq >= 0:
  793 + value = binsq / deg2sr / (ewidth[ieng])
  794 + hdu[name][0,index] = value * scale
  795 + #
  796 + # Debugging
  797 + if False:
  798 + for ieng in range(neng):
  799 + sum = 0.0
  800 + for ioff in range(noffset):
  801 + index = ieng + ioff * neng
  802 + binsq = array.GetBinContent(ieng+1,ioff+1)
  803 + value = binsq / deg2sr / ewidth[ieng]
  804 + sum += value
  805 + if ioff == 0:
  806 + peak = value
  807 + peaksq = binsq
  808 + print(pow(10.0, energies.GetBinLowEdge(ieng+1))+ \
  809 + pow(10.0, energies.GetBinUpEdge(ieng+1)), sum, peak, peaksq)
  810 +
  811 + # Collect boundary information
  812 + bd_detx = "DETX(%.2f-%.2f)deg" % (-det_max, det_max)
  813 + bd_dety = "DETY(%.2f-%.2f)deg" % (-det_max, det_max)
  814 + bd_eng = "ENERG(%.4f-%.2f)TeV" % \
  815 + (pow(10.0, energies.GetBinLowEdge(1)), \
  816 + pow(10.0, energies.GetBinUpEdge(neng)))
  817 + bounds = [bd_detx, bd_dety, bd_eng]
  818 +
  819 + # Return boundary information
  820 + return bounds
686 821  
687   - def root2caldb(self, filename, rebin=False, psftype="Gauss"):
  822 + def root2caldb(self, filename, rebin=False, psftype="Gauss", scale=1.0):
688 823 """
689 824 Translate ROOT to CALDB information.
690 825  
... ... @@ -706,7 +841,8 @@ class caldb():
706 841 self.root2edisp(file)
707 842  
708 843 # Create background
709   - self.root2bgd(file)
  844 + #self.root2bgd(file)
  845 + self.root2bgd3D(file, scale=scale)
710 846  
711 847 # Return
712 848 return
... ... @@ -1057,6 +1193,42 @@ class caldb():
1057 1193 # Return
1058 1194 return
1059 1195  
  1196 + def root2bgd3D(self, file, scale=1.0):
  1197 + """
  1198 + Translate ROOT to CALDB background extension. The following ROOT
  1199 + histograms are used:
  1200 +
  1201 + BGRatePerSqDeg_offaxis -> BGD
  1202 + BGRatePerSqDeg_offaxis -> BGD_RECO
  1203 +
  1204 + Parameters:
  1205 + file - ROOT file.
  1206 + Keywords:
  1207 + scale - Background rate scaling factor
  1208 + """
  1209 + # Continue only if background HDU has been opened
  1210 + if self.hdu_bgd != None:
  1211 +
  1212 + # Allocate ROOT 2D array
  1213 + array = TH2F()
  1214 + file.GetObject("BGRatePerSqDeg_offaxis", array)
  1215 +
  1216 + # Set boundaries
  1217 + bounds = self.make_3D(array, self.hdu_bgd, None, "deg")
  1218 + for b in bounds:
  1219 + self.bgd_bounds.append(b)
  1220 + self.set_cif_keywords(self.hdu_bgd, self.bgd_name, \
  1221 + self.bgd_bounds, self.bgd_desc)
  1222 +
  1223 + # BGD
  1224 + self.make_3D(array, self.hdu_bgd, "BGD", "1/s/MeV/sr", scale=scale)
  1225 +
  1226 + # BGD_RECO
  1227 + self.make_3D(array, self.hdu_bgd, "BGD_RECO", "1/s/MeV/sr", scale=scale)
  1228 +
  1229 + # Return
  1230 + return
  1231 +
1060 1232 def set_cif_keywords(self, hdu, name, bounds, desc):
1061 1233 """
1062 1234 Set standard CIF keywords for extension.
... ... @@ -1095,13 +1267,14 @@ def set_test():
1095 1267 # Set database attributes
1096 1268 path = "/project-data/cta/performance/prod1/IFAEOffaxisPerformanceBEI_May2012"
1097 1269 rebin = True
  1270 + scale = 1.0
1098 1271  
1099 1272 # Set database content
1100 1273 db = [{'inst': "e", 'id': "IFAE20120510_50h",
1101   - 'path': path, 'rebin': rebin, 'psftype': "Gauss",
  1274 + 'path': path, 'rebin': rebin, 'psftype': "Gauss", 'scale': scale,
1102 1275 'file': "SubarrayE_IFAE_50hours_20120510_offaxis.root"},
1103 1276 {'inst': "e", 'id': "IFAE20120510_50h_King",
1104   - 'path': path, 'rebin': rebin, 'psftype': "King",
  1277 + 'path': path, 'rebin': rebin, 'psftype': "King", 'scale': scale,
1105 1278 'file': "SubarrayE_IFAE_50hours_20120510_offaxis.root"}]
1106 1279  
1107 1280 # Return database
... ... @@ -1119,16 +1292,17 @@ def set_prod1_ifae():
1119 1292 path = "/project-data/cta/performance/prod1/IFAEOffaxisPerformanceBEI_May2012"
1120 1293 rebin = True
1121 1294 psftype = "Gauss"
  1295 + scale = 1.0
1122 1296  
1123 1297 # Set database content
1124 1298 db = [{'inst': "b", 'id': "IFAE20120510_50h",
1125   - 'path': path, 'rebin': rebin, 'psftype': psftype,
  1299 + 'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1126 1300 'file': "SubarrayB_IFAE_50hours_20120510_offaxis.root"},
1127 1301 {'inst': "e", 'id': "IFAE20120510_50h",
1128   - 'path': path, 'rebin': rebin, 'psftype': psftype,
  1302 + 'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1129 1303 'file': "SubarrayE_IFAE_50hours_20120510_offaxis.root"},
1130 1304 {'inst': "i", 'id': "IFAE20120510_50h",
1131   - 'path': path, 'rebin': rebin, 'psftype': psftype,
  1305 + 'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1132 1306 'file': "SubarrayI_IFAE_50hours_20120510_offaxis.root"}]
1133 1307  
1134 1308 # Return database
... ... @@ -1146,76 +1320,77 @@ def set_prod2_desy():
1146 1320 path = "/project-data/cta/performance/prod2/Performance_DESY_20140128"
1147 1321 rebin = False
1148 1322 psftype = "Gauss"
  1323 + scale = 1.0/180000.0
1149 1324  
1150 1325 # Set database content
1151 1326 db = [{'inst': "aar", 'id': "DESY20140105_50h",
1152   - 'path': path, 'rebin': rebin, 'psftype': psftype,
  1327 + 'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1153 1328 'file': "DESY.d20140105.Erec1.V2.ID0NIM2.prod2-Aar-NS.S.2a.180000s.root"},
1154 1329 {'inst': "aar", 'id': "DESY20140105_50h_0deg",
1155   - 'path': path, 'rebin': rebin, 'psftype': psftype,
  1330 + 'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1156 1331 'file': "DESY.d20140105.Erec1.V2.ID0_0degNIM2.prod2-Aar-NS.S.2a.180000s.root"},
1157 1332 {'inst': "aar", 'id': "DESY20140105_50h_180deg",
1158   - 'path': path, 'rebin': rebin, 'psftype': psftype,
  1333 + 'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1159 1334 'file': "DESY.d20140105.Erec1.V2.ID0_180degNIM2.prod2-Aar-NS.S.2a.180000s.root"},
1160 1335  
1161 1336 {'inst': "aar500", 'id': "DESY20140105_50h",
1162   - 'path': path, 'rebin': rebin, 'psftype': psftype,
  1337 + 'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1163 1338 'file': "DESY.d20140105.Erec1.V2.ID0NIM2.prod2-Aar-500m-NS.S.2a.180000s.root"},
1164 1339 {'inst': "aar500", 'id': "DESY20140105_50h_0deg",
1165   - 'path': path, 'rebin': rebin, 'psftype': psftype,
  1340 + 'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1166 1341 'file': "DESY.d20140105.Erec1.V2.ID0_0degNIM2.prod2-Aar-500m-NS.S.2a.180000s.root"},
1167 1342 {'inst': "aar500", 'id': "DESY20140105_50h_180deg",
1168   - 'path': path, 'rebin': rebin, 'psftype': psftype,
  1343 + 'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1169 1344 'file': "DESY.d20140105.Erec1.V2.ID0_180degNIM2.prod2-Aar-500m-NS.S.2a.180000s.root"},
1170 1345  
1171 1346 {'inst': "leoncito", 'id': "DESY20140105_50h",
1172   - 'path': path, 'rebin': rebin, 'psftype': psftype,
  1347 + 'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1173 1348 'file': "DESY.d20140105.Erec1.V2.ID0NIM2.prod2-LeoncitoPP-NS.S.2a.180000s.root"},
1174 1349 {'inst': "leoncito", 'id': "DESY20140105_50h_0deg",
1175   - 'path': path, 'rebin': rebin, 'psftype': psftype,
  1350 + 'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1176 1351 'file': "DESY.d20140105.Erec1.V2.ID0_0degNIM2.prod2-LeoncitoPP-NS.S.2a.180000s.root"},
1177 1352 {'inst': "leoncito", 'id': "DESY20140105_50h_180deg",
1178   - 'path': path, 'rebin': rebin, 'psftype': psftype,
  1353 + 'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1179 1354 'file': "DESY.d20140105.Erec1.V2.ID0_180degNIM2.prod2-LeoncitoPP-NS.S.2a.180000s.root"},
1180 1355  
1181 1356 {'inst': "sac", 'id': "DESY20140105_50h",
1182   - 'path': path, 'rebin': rebin, 'psftype': psftype,
  1357 + 'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1183 1358 'file': "DESY.d20140105.Erec1.V2.ID0NIM2.prod2-SAC100-NS.S.2a.180000s.root"},
1184 1359 {'inst': "sac", 'id': "DESY20140105_50h_0deg",
1185   - 'path': path, 'rebin': rebin, 'psftype': psftype,
  1360 + 'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1186 1361 'file': "DESY.d20140105.Erec1.V2.ID0_0degNIM2.prod2-SAC100-NS.S.2a.180000s.root"},
1187 1362 {'inst': "sac", 'id': "DESY20140105_50h_180deg",
1188   - 'path': path, 'rebin': rebin, 'psftype': psftype,
  1363 + 'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1189 1364 'file': "DESY.d20140105.Erec1.V2.ID0_180degNIM2.prod2-SAC100-NS.S.2a.180000s.root"},
1190 1365  
1191 1366 {'inst': "spm", 'id': "DESY20140105_50h",
1192   - 'path': path, 'rebin': rebin, 'psftype': psftype,
  1367 + 'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1193 1368 'file': "DESY.d20140105.Erec1.V2.ID0NIM2.prod2-SPM-NS.N.2NN.180000s.root"},
1194 1369 {'inst': "spm", 'id': "DESY20140105_50h_0deg",
1195   - 'path': path, 'rebin': rebin, 'psftype': psftype,
  1370 + 'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1196 1371 'file': "DESY.d20140105.Erec1.V2.ID0_0degNIM2.prod2-SPM-NS.N.2NN.180000s.root"},
1197 1372 {'inst': "spm", 'id': "DESY20140105_50h_180deg",
1198   - 'path': path, 'rebin': rebin, 'psftype': psftype,
  1373 + 'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1199 1374 'file': "DESY.d20140105.Erec1.V2.ID0_180degNIM2.prod2-SPM-NS.N.2NN.180000s.root"},
1200 1375  
1201 1376 {'inst': "tenerife", 'id': "DESY20140105_50h",
1202   - 'path': path, 'rebin': rebin, 'psftype': psftype,
  1377 + 'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1203 1378 'file': "DESY.d20140105.Erec1.V2.ID0NIM2.prod2-Tenerife-NS.N.2NN.180000s.root"},
1204 1379 {'inst': "tenerife", 'id': "DESY20140105_50h_0deg",
1205   - 'path': path, 'rebin': rebin, 'psftype': psftype,
  1380 + 'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1206 1381 'file': "DESY.d20140105.Erec1.V2.ID0_0degNIM2.prod2-Tenerife-NS.N.2NN.180000s.root"},
1207 1382 {'inst': "tenerife", 'id': "DESY20140105_50h_180deg",
1208   - 'path': path, 'rebin': rebin, 'psftype': psftype,
  1383 + 'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1209 1384 'file': "DESY.d20140105.Erec1.V2.ID0_180degNIM2.prod2-Tenerife-NS.N.2NN.180000s.root"},
1210 1385  
1211 1386 {'inst': "us", 'id': "DESY20140105_50h",
1212   - 'path': path, 'rebin': rebin, 'psftype': psftype,
  1387 + 'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1213 1388 'file': "DESY.d20140105.Erec1.V2.ID0NIM2.prod2-US-NS.N.2NN.180000s.root"},
1214 1389 {'inst': "us", 'id': "DESY20140105_50h_0deg",
1215   - 'path': path, 'rebin': rebin, 'psftype': psftype,
  1390 + 'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1216 1391 'file': "DESY.d20140105.Erec1.V2.ID0_0degNIM2.prod2-US-NS.N.2NN.180000s.root"},
1217 1392 {'inst': "us", 'id': "DESY20140105_50h_180deg",
1218   - 'path': path, 'rebin': rebin, 'psftype': psftype,
  1393 + 'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1219 1394 'file': "DESY.d20140105.Erec1.V2.ID0_180degNIM2.prod2-US-NS.N.2NN.180000s.root"}]
1220 1395  
1221 1396 # Return database
... ... @@ -1259,6 +1434,7 @@ if __name__ == &#39;__main__&#39;:
1259 1434 id = entry['id']
1260 1435 rebin = entry['rebin']
1261 1436 psftype = entry['psftype']
  1437 + scale = entry['scale']
1262 1438 filename = path+"/"+entry['file']
1263 1439  
1264 1440 # Allocate caldb entry
... ... @@ -1268,4 +1444,4 @@ if __name__ == &#39;__main__&#39;:
1268 1444 irf.open("file")
1269 1445  
1270 1446 # Translate ROOT to CALDB information
1271   - irf.root2caldb(filename, rebin=rebin, psftype=psftype)
  1447 + irf.root2caldb(filename, rebin=rebin, psftype=psftype, scale=scale)
... ...