jean perrot / ctools

Commit 4dce62140a1d78f62281a0fd11971fd4d6fb847e

Authored by Jürgen Knödlseder
2 parents 1fa95096 65346e20

Merge branch '00-07-02' into integration

Inline Side-by-side
Showing 3 changed files with 414 additions and 32 deletions
scripts/Makefile.am
  View file @4dce621
... ... @@ -4,10 +4,13 @@
4 4 binscripts = cta_kb_bkgrates.py cta_root2irf.py cta_make_nodes.py
5 5  
6 6 # Scripts to be distributed in pythondir
7   -pythonscripts = cspull.py cssens.py cstsdist.py obsutils.py
  7 +pythonscripts = cspull.py cssens.py cstsdist.py \
  8 + obsutils.py \
  9 + lsmodel.py
8 10  
9 11 # Scripts for which symbolic links will be created in bindir
10   -linkscripts = cspull cssens cstsdist
  12 +linkscripts = cspull cssens cstsdist \
  13 + lsmodel
11 14  
12 15 # Distribute scripts
13 16 EXTRA_DIST = $(binscripts) $(pythonscripts)
... ...
scripts/cta_root2caldb.py
  View file @4dce621
... ... @@ -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)
... ...
scripts/lsmodel.py 0 → 100755
  View file @4dce621
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script generates a model map from a counts map and a source model.
  4 +#
  5 +# Copyright (C) 2014 Juergen Knoedlseder
  6 +#
  7 +# This program is free software: you can redistribute it and/or modify
  8 +# it under the terms of the GNU General Public License as published by
  9 +# the Free Software Foundation, either version 3 of the License, or
  10 +# (at your option) any later version.
  11 +#
  12 +# This program is distributed in the hope that it will be useful,
  13 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  15 +# GNU General Public License for more details.
  16 +#
  17 +# You should have received a copy of the GNU General Public License
  18 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  19 +#
  20 +# ==========================================================================
  21 +import gammalib
  22 +import sys
  23 +
  24 +
  25 +# ============= #
  26 +# lsmodel class #
  27 +# ============= #
  28 +class lsmodel(gammalib.GApplication):
  29 + """
  30 + This class implements the model map generation script. It derives from
  31 + the GammaLib::GApplication class which provides support for parameter
  32 + files, command line arguments, and logging. In that way the Python
  33 + script behaves just as a regular ctool.
  34 + """
  35 + def __init__(self, *argv):
  36 + """
  37 + Constructor.
  38 + """
  39 + # Set name
  40 + self.name = "lsmodel"
  41 + self.version = "0.1.0"
  42 +
  43 + # Make sure that parfile exists
  44 + file = self.parfile()
  45 +
  46 + # Initialise application
  47 + if len(argv) == 0:
  48 + gammalib.GApplication.__init__(self, self.name, self.version)
  49 + elif len(argv) ==1:
  50 + gammalib.GApplication.__init__(self, self.name, self.version, *argv)
  51 + else:
  52 + raise TypeError("Invalid number of arguments given.")
  53 +
  54 + # Set logger properties
  55 + self.log_header()
  56 + self.log.date(True)
  57 +
  58 + # Return
  59 + return
  60 +
  61 + def __del__(self):
  62 + """
  63 + Destructor.
  64 + """
  65 + # Write separator into logger
  66 + if self.logTerse():
  67 + self.log("\n")
  68 +
  69 + # Return
  70 + return
  71 +
  72 + def parfile(self):
  73 + """
  74 + Check if parfile exists. If parfile does not exist then create a
  75 + default parfile. This kluge avoids shipping the script with a parfile.
  76 + """
  77 + # Set parfile name
  78 + parfile = self.name+".par"
  79 +
  80 + try:
  81 + pars = gammalib.GApplicationPars(parfile)
  82 + except:
  83 + # Signal if parfile was not found
  84 + sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
  85 +
  86 + # Create default parfile
  87 + pars = gammalib.GApplicationPars()
  88 + pars.append(gammalib.GApplicationPar("cntmap","f","a","cntmap.fits","","","Counts map"))
  89 + pars.append(gammalib.GApplicationPar("expmap","f","a","binned_expmap.fits","","","Exposure map"))
  90 + pars.append(gammalib.GApplicationPar("ltcube","f","a","ltcube.fits","","","Livetime cube"))
  91 + pars.append(gammalib.GApplicationPar("srcmdl","f","a","crab_model.xml","","","Source model"))
  92 + pars.append(gammalib.GApplicationPar("caldb","s","h","","","","Calibration database"))
  93 + pars.append(gammalib.GApplicationPar("irf","s","a","P7REP_SOURCE_V15","","","Instrument response function"))
  94 + pars.append(gammalib.GApplicationPar("outfile","f","a","model_cube.fits","","","Output file name"))
  95 + pars.append_standard()
  96 + pars.save(parfile)
  97 +
  98 + # Return
  99 + return
  100 +
  101 + def get_parameters(self):
  102 + """
  103 + Get parameters from parfile and setup the observation.
  104 + """
  105 + # Get parameters
  106 + self.m_cntmap = self["cntmap"].filename()
  107 + self.m_expmap = self["expmap"].filename()
  108 + self.m_ltcube = self["ltcube"].filename()
  109 + self.m_srcmdl = self["srcmdl"].filename()
  110 + self.m_caldb = self["caldb"].string()
  111 + self.m_irf = self["irf"].string()
  112 + self.m_outfile = self["outfile"].filename()
  113 +
  114 + # Set some fixed parameters
  115 + self.m_log = False # Logging in client tools
  116 + self.m_debug = False # Debugging in client tools
  117 +
  118 + # Return
  119 + return
  120 +
  121 + def execute(self):
  122 + """
  123 + Execute the script.
  124 + """
  125 + # Run the script
  126 + self.run()
  127 +
  128 + # Return
  129 + return
  130 +
  131 + def run(self):
  132 + """
  133 + Run the script.
  134 + """
  135 + # Switch screen logging on in debug mode
  136 + if self.logDebug():
  137 + self.log.cout(True)
  138 +
  139 + # Get parameters
  140 + self.get_parameters()
  141 +
  142 + # Write input parameters into logger
  143 + if self.logTerse():
  144 + self.log_parameters()
  145 + self.log("\n")
  146 +
  147 + # Write header
  148 + if self.logTerse():
  149 + self.log("\n")
  150 + self.log.header1("Generate model cube")
  151 +
  152 + # Set LAT observation
  153 + obs = gammalib.GLATObservation()
  154 + obs.load_binned(self.m_cntmap, self.m_expmap, self.m_ltcube)
  155 + obs.response(self.m_irf, self.m_caldb)
  156 +
  157 + # Set model (use first model for the moment)
  158 + models = gammalib.GModels(self.m_srcmdl)
  159 + model = models[0]
  160 +
  161 + # Get deep copy of GLATEventCube from observation
  162 + cube = gammalib.GLATEventCube(obs.events())
  163 +
  164 + # Fill cube with model values
  165 + for i in range(cube.size()):
  166 + event = cube[i]
  167 + value = model.eval(event, obs) * event.size()
  168 + event.counts(value)
  169 + cube[i] = event
  170 +
  171 + # Log cube
  172 + if self.logTerse():
  173 + self.log(str(cube))
  174 + self.log("\n")
  175 +
  176 + # Write header
  177 + if self.logTerse():
  178 + self.log("\n")
  179 + self.log.header1("Save model cube")
  180 +
  181 + # Save model cube
  182 + cube.save(self.m_outfile, self.clobber())
  183 +
  184 + # Return
  185 + return
  186 +
  187 +
  188 +# ======================== #
  189 +# Main routine entry point #
  190 +# ======================== #
  191 +if __name__ == '__main__':
  192 + """
  193 + Generates a model map from a counts map and a source model.
  194 + """
  195 + # Create instance of application
  196 + app = lsmodel(sys.argv)
  197 +
  198 + # Open logfile
  199 + app.logFileOpen()
  200 +
  201 + # Execute application
  202 + app.execute()
  203 +
0 204 \ No newline at end of file
... ...