Thierry Reposeur / ctools

Commit e27b5c3455652aba9fe756f83bd3c70cd210630b

Authored by Jürgen Knödlseder
1 parent 12b82d2a

Implement csroot2caldb background template generation.

Inline Side-by-side
Showing 4 changed files with 129 additions and 84 deletions
ChangeLog
  View file @e27b5c3
1   -2016-04-22 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
  1 +2016-04-25 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
2 2  
3 3 * Version 1.1.0 released
4 4 ========================
5 5  
  6 + Add "bgdoversample" and "bgdinfill" parameters to csroot2caldb
6 7 Add csobsdef.pntdef() method (#1771)
7 8 Add ctmapcube tool (#1760)
8 9 Add "seed" parameter to cspull (#1755)
... ...
1 1 New Features and Important Changes in ctools 1.1.0
2 2  
3   -22 April 2016
  3 +25 April 2016
4 4  
5 5  
6 6 Introduction
... ... @@ -198,7 +198,12 @@ support energy dispersion for stacked analysis (#1724).
198 198  
199 199 Add csworkflow script that executes an arbitrary analysis workflow (#1508).
200 200  
201   -Transformed cta_root2caldb into csroot2caldb script.
  201 +Transformed cta_root2caldb into csroot2caldb script. The parameters
  202 +"bgdoversample" and "bgdinfill" have been added to the parameter file to
  203 +provide more control over the background template generation. The background
  204 +template is now generated using linear interpolation, and the
  205 +"bgdoversample" allows to generate the background template at a finer
  206 +sampling than the spatial sampling in the ROOT file.
202 207  
203 208 Add csobsdef.pntdef() method to set directly a pointing definition from a
204 209 gammalib.GCsv object (#1771).
... ...
cscripts/csroot2caldb.par
  View file @e27b5c3
... ... @@ -17,16 +17,18 @@
17 17 #
18 18 # Input parameters
19 19 #=================
20   -infile, f, a, irf.root,,, "Input ROOT file"
21   -outdir, s, h, ,,, "Output CALDB directory"
22   -inst, s, a, prod3,,, "Instrument name"
23   -id, s, a, South_50h,,, "Instrument Response Function identifier"
24   -version, s, a, prod3,,, "Instrument Response Function version"
25   -analysis, s, a, DESY,,, "Instrument Response Function analysis"
26   -psftype, s, h, Gauss,Gauss|King,, "Point spread function type"
27   -rebin, b, h, no,,, "Rebin"
28   -eascale, r, h, 1.0,,, "Scaling factor for effective areas"
29   -bgdscale, r, h, 1.0,,, "Scaling factor for background rates"
  20 +infile, f, a, irf.root,,, "Input ROOT file"
  21 +outdir, s, h, ,,, "Output CALDB directory"
  22 +inst, s, a, prod3,,, "Instrument name"
  23 +id, s, a, South_50h,,, "Instrument Response Function identifier"
  24 +version, s, a, prod3,,, "Instrument Response Function version"
  25 +analysis, s, a, DESY,,, "Instrument Response Function analysis"
  26 +psftype, s, h, Gauss,Gauss|King,, "Point spread function type"
  27 +rebin, b, h, no,,, "Rebin"
  28 +eascale, r, h, 1.0,,, "Scaling factor for effective areas"
  29 +bgdscale, r, h, 1.0,,, "Scaling factor for background rates"
  30 +bgdoversample, i, h, 1,,, "Spatial oversampling factor for background template"
  31 +bgdinfill, b, h, no,,, "Infill of background template?"
30 32  
31 33 #
32 34 # Standard parameters
... ...
cscripts/csroot2caldb.py
  View file @e27b5c3
... ... @@ -73,16 +73,18 @@ class csroot2caldb(ctools.cscript):
73 73 Get parameters from parfile.
74 74 """
75 75 # Query parameters
76   - self["infile"].filename()
77   - self["outdir"].string()
78   - self["inst"].string()
79   - self["id"].string()
80   - self["analysis"].string()
81   - self["version"].string()
82   - self["psftype"].string()
83   - self["rebin"].boolean()
84   - self["eascale"].real()
85   - self["bgdscale"].real()
  76 + self['infile'].filename()
  77 + self['outdir'].string()
  78 + self['inst'].string()
  79 + self['id'].string()
  80 + self['analysis'].string()
  81 + self['version'].string()
  82 + self['psftype'].string()
  83 + self['rebin'].boolean()
  84 + self['eascale'].real()
  85 + self['bgdscale'].real()
  86 + self['bgdoversample'].integer()
  87 + self['bgdinfill'].boolean()
86 88  
87 89 # Write input parameters into logger
88 90 if self._logTerse():
... ... @@ -545,7 +547,7 @@ class csroot2caldb(ctools.cscript):
545 547  
546 548 # Create background
547 549 #self._root2bgd(file)
548   - self._root2bgd3D(file, bgdscale=self["bgdscale"].real())
  550 + self._root2bgd3D(file)
549 551  
550 552 # Return
551 553 return
... ... @@ -912,34 +914,34 @@ class csroot2caldb(ctools.cscript):
912 914 # Return
913 915 return
914 916  
915   - def _root2bgd3D(self, file, bgdscale=1.0):
  917 + def _root2bgd3D(self, file):
916 918 """
917   - Translate ROOT to CALDB background extension. The following ROOT
918   - histograms are used:
919   -
  919 + Translate ROOT to CALDB background extension.
  920 +
  921 + The following ROOT histograms are used:
920 922 BGRatePerSqDeg_offaxis -> BGD
921 923  
922   - Parameters:
923   - file - ROOT file.
924   - Keywords:
925   - scale - Background rate scaling factor
  924 + Parameters
  925 + ----------
  926 + file : `~ROOT.TFile`
  927 + ROOT file.
926 928 """
927 929 # Continue only if background HDU has been opened
928 930 if self._hdu_bgd != None:
929 931  
930 932 # Allocate ROOT 2D array
931 933 array = TH2F()
932   - file.GetObject("BGRatePerSqDeg_offaxis", array)
  934 + file.GetObject('BGRatePerSqDeg_offaxis', array)
933 935  
934 936 # Set boundaries
935   - bounds = self._make_3D(array, self._hdu_bgd, None, "deg")
  937 + bounds = self._make_3D(array, self._hdu_bgd, None, 'deg')
936 938 for b in bounds:
937 939 self._bgd_bounds.append(b)
938   - self._set_cif_keywords(self._hdu_bgd, self._bgd_name, \
  940 + self._set_cif_keywords(self._hdu_bgd, self._bgd_name,
939 941 self._bgd_bounds, self._bgd_desc)
940 942  
941 943 # BGD (reconstructed energy)
942   - self._make_3D(array, self._hdu_bgd, "BGD", "1/s/MeV/sr", scale=bgdscale)
  944 + self._make_3D(array, self._hdu_bgd, 'BGD', '1/s/MeV/sr')
943 945  
944 946 # Return
945 947 return
... ... @@ -1038,40 +1040,76 @@ class csroot2caldb(ctools.cscript):
1038 1040 # Return boundary information
1039 1041 return bounds
1040 1042  
1041   - def _make_3D(self, array, hdu, name, unit, scale=1.0):
  1043 + def _make_3D(self, array, hdu, name, unit):
1042 1044 """
1043   - Make 3D cube as function of DETX, DETY and energy from a
  1045 + Generate background 3D cube
  1046 +
  1047 + Generates 3D cube as function of DETX, DETY and energy from a
1044 1048 ROOT 2D histogram. If the HDU has already the energy and
1045 1049 offset angle columns, this method will simply add another data
1046   - column.
1047   - If name==None, the method will not append any data column.
1048   -
1049   - Parameters:
1050   - array - ROOT 2D histogram.
1051   - hdu - FITS HDU.
1052   - name - Data column name.
1053   - unit - Data unit.
1054   - Keywords:
1055   - scale - Scaling factor for histogram values.
  1050 + column. If name==None, the method will not append any data column.
  1051 +
  1052 + Parameters
  1053 + ----------
  1054 + array : `~ROOT.TH2F`
  1055 + ROOT 2D histogram.
  1056 + hdu : `~gammalib.GFitsHDU`
  1057 + FITS HDU.
  1058 + name : str
  1059 + Data column name.
  1060 + unit : str
  1061 + Data unit.
1056 1062 """
  1063 + # Get User parameters
  1064 + scale = self['bgdscale'].real()
  1065 + oversample = self['bgdoversample'].integer()
  1066 + infill = self['bgdinfill'].boolean()
  1067 +
1057 1068 # Set constants
1058 1069 deg2sr = 0.01745329*0.01745329
1059 1070  
1060 1071 # Extract energy and offset angle vectors
1061   - energies = array.GetXaxis()
1062   - neng = energies.GetNbins()
1063   - offsets = array.GetYaxis()
1064   - noffset = offsets.GetNbins()
1065   - ewidth = [] # in MeV
  1072 + energies = array.GetXaxis()
  1073 + neng = energies.GetNbins()
  1074 + offsets = array.GetYaxis()
  1075 + noffset = offsets.GetNbins()
  1076 + theta_max = offsets.GetBinUpEdge(noffset) # Maximum theta
  1077 + ewidth = [] # in MeV
1066 1078 for ieng in range(neng):
1067   - ewidth.append(pow(10.0, energies.GetBinUpEdge(ieng+1) +6.0) - \
  1079 + ewidth.append(pow(10.0, energies.GetBinUpEdge(ieng+1) +6.0) -
1068 1080 pow(10.0, energies.GetBinLowEdge(ieng+1)+6.0))
1069 1081  
  1082 + # Optionally fill all empty bins by the preceding value that occur
  1083 + # for all offset angles between the first and the last non zero
  1084 + # bin
  1085 + if infill:
  1086 + for ioff in range(noffset):
  1087 + i_start = -1
  1088 + i_stop = -1
  1089 + value = 0.0
  1090 + for ieng in range(neng):
  1091 + if array.GetBinContent(ieng+1,ioff+1) > 0.0:
  1092 + i_start = ieng
  1093 + value = array.GetBinContent(ieng+1,ioff+1)
  1094 + break
  1095 + for ieng in range(neng-1, -1, -1):
  1096 + if array.GetBinContent(ieng+1,ioff+1) > 0.0:
  1097 + i_stop = ieng
  1098 + break
  1099 + if i_start > -1 and i_stop > -1:
  1100 + for ieng in range(i_start, i_stop+1):
  1101 + if array.GetBinContent(ieng+1,ioff+1) == 0.0:
  1102 + array.SetBinContent(ieng+1,ioff+1,value)
  1103 + print('Background infill ieng=%d ioff=%d value=%f' %
  1104 + (ieng,ioff,value))
  1105 + else:
  1106 + value = array.GetBinContent(ieng+1,ioff+1)
  1107 +
1070 1108 # Build DETX and DETY axes
1071 1109 ndet = array.GetYaxis().GetNbins()
1072   - ndets = 2*ndet
  1110 + ndets = 2*ndet*oversample
1073 1111 det_max = array.GetYaxis().GetBinUpEdge(ndet)
1074   - det_bin = det_max/float(ndet)
  1112 + det_bin = det_max/float(ndet*oversample)
1075 1113 dets_lo = []
1076 1114 dets_hi = []
1077 1115 dets2 = []
... ... @@ -1084,36 +1122,36 @@ class csroot2caldb(ctools.cscript):
1084 1122 dets2.append(det_val*det_val)
1085 1123  
1086 1124 # Append DETX_LO column
1087   - if not hdu.contains("DETX_LO"):
1088   - hdu.append(gammalib.GFitsTableFloatCol("DETX_LO", 1, ndets))
1089   - hdu["DETX_LO"].unit("deg")
  1125 + if not hdu.contains('DETX_LO'):
  1126 + hdu.append(gammalib.GFitsTableFloatCol('DETX_LO', 1, ndets))
  1127 + hdu['DETX_LO'].unit('deg')
1090 1128 for i in range(ndets):
1091   - hdu["DETX_LO"][0,i] = dets_lo[i]
  1129 + hdu['DETX_LO'][0,i] = dets_lo[i]
1092 1130  
1093 1131 # Append DETX_HI column
1094   - if not hdu.contains("DETX_HI"):
1095   - hdu.append(gammalib.GFitsTableFloatCol("DETX_HI", 1, ndets))
1096   - hdu["DETX_HI"].unit("deg")
  1132 + if not hdu.contains('DETX_HI'):
  1133 + hdu.append(gammalib.GFitsTableFloatCol('DETX_HI', 1, ndets))
  1134 + hdu['DETX_HI'].unit('deg')
1097 1135 for i in range(ndets):
1098   - hdu["DETX_HI"][0,i] = dets_hi[i]
  1136 + hdu['DETX_HI'][0,i] = dets_hi[i]
1099 1137  
1100 1138 # Append DETY_LO column
1101   - if not hdu.contains("DETY_LO"):
1102   - hdu.append(gammalib.GFitsTableFloatCol("DETY_LO", 1, ndets))
1103   - hdu["DETY_LO"].unit("deg")
  1139 + if not hdu.contains('DETY_LO'):
  1140 + hdu.append(gammalib.GFitsTableFloatCol('DETY_LO', 1, ndets))
  1141 + hdu['DETY_LO'].unit('deg')
1104 1142 for i in range(ndets):
1105   - hdu["DETY_LO"][0,i] = dets_lo[i]
  1143 + hdu['DETY_LO'][0,i] = dets_lo[i]
1106 1144  
1107 1145 # Append DETY_HI column
1108   - if not hdu.contains("DETY_HI"):
1109   - hdu.append(gammalib.GFitsTableFloatCol("DETY_HI", 1, ndets))
1110   - hdu["DETY_HI"].unit("deg")
  1146 + if not hdu.contains('DETY_HI'):
  1147 + hdu.append(gammalib.GFitsTableFloatCol('DETY_HI', 1, ndets))
  1148 + hdu['DETY_HI'].unit('deg')
1111 1149 for i in range(ndets):
1112   - hdu["DETY_HI"][0,i] = dets_hi[i]
  1150 + hdu['DETY_HI'][0,i] = dets_hi[i]
1113 1151  
1114 1152 # Append ENERG_LO and ENERG_HI columns
1115   - self._append_column(hdu, "ENERG_LO", "TeV", energies, log=True)
1116   - self._append_column(hdu, "ENERG_HI", "TeV", energies, log=True)
  1153 + self._append_column(hdu, 'ENERG_LO', 'TeV', energies, log=True)
  1154 + self._append_column(hdu, 'ENERG_HI', 'TeV', energies, log=True)
1117 1155  
1118 1156 # Append array column
1119 1157 if name != None and not hdu.contains(name):
... ... @@ -1125,12 +1163,11 @@ class csroot2caldb(ctools.cscript):
1125 1163 for ieng in range(neng):
1126 1164 index = ix + (iy + ieng * ndets) * ndets
1127 1165 theta = math.sqrt(dets2[ix] + dets2[iy])
1128   - ioff = offsets.FindBin(theta)
1129   - if ioff > 0 and ioff <= noffset:
1130   - binsq = array.GetBinContent(ieng+1,ioff)
1131   - if binsq >= 0:
1132   - value = binsq / deg2sr / (ewidth[ieng])
1133   - hdu[name][0,index] = value * scale
  1166 + if theta < theta_max:
  1167 + binsq = array.Interpolate(energies.GetBinCenter(ieng+1),
  1168 + theta)
  1169 + value = binsq / deg2sr / (ewidth[ieng])
  1170 + hdu[name][0,index] = value * scale
1134 1171  
1135 1172 # Debugging
1136 1173 if False:
... ... @@ -1144,13 +1181,13 @@ class csroot2caldb(ctools.cscript):
1144 1181 if ioff == 0:
1145 1182 peak = value
1146 1183 peaksq = binsq
1147   - print(pow(10.0, energies.GetBinLowEdge(ieng+1))+ \
  1184 + print(pow(10.0, energies.GetBinLowEdge(ieng+1))+
1148 1185 pow(10.0, energies.GetBinUpEdge(ieng+1)), sum, peak, peaksq)
1149 1186  
1150 1187 # Collect boundary information
1151   - bd_detx = "DETX(%.2f-%.2f)deg" % (-det_max, det_max)
1152   - bd_dety = "DETY(%.2f-%.2f)deg" % (-det_max, det_max)
1153   - bd_eng = "ENERG(%.4f-%.2f)TeV" % \
  1188 + bd_detx = 'DETX(%.2f-%.2f)deg' % (-det_max, det_max)
  1189 + bd_dety = 'DETY(%.2f-%.2f)deg' % (-det_max, det_max)
  1190 + bd_eng = 'ENERG(%.4f-%.2f)TeV' % \
1154 1191 (pow(10.0, energies.GetBinLowEdge(1)),
1155 1192 pow(10.0, energies.GetBinUpEdge(neng)))
1156 1193 bounds = [bd_detx, bd_dety, bd_eng]
... ...