Pierrick Martin / ctools

Commit cbe8c7475ba27dcda1a44fec99f72ecfb9735a11

Authored by Jürgen Knödlseder
1 parent cbe733dd

Write migration matrix instead of ERES and EBIAS into FITS file

Inline Side-by-side
Showing 1 changed file with 143 additions and 4 deletions
scripts/cta_root2caldb.py
  View file @cbe8c74
... ... @@ -3,7 +3,7 @@
3 3 # This script generates IRFs in the CALDB format of HEASARC using ROOT 2D
4 4 # performance files.
5 5 # -------------------------------------------------------------------------
6   -# Copyright (C) 2011-2014 Juergen Knoedlseder
  6 +# Copyright (C) 2011-2015 Juergen Knoedlseder
7 7 # -------------------------------------------------------------------------
8 8 #
9 9 # This program is free software: you can redistribute it and/or modify
... ... @@ -24,7 +24,7 @@
24 24 # - gammalib
25 25 # - pyroot
26 26 # ==========================================================================
27   -from ROOT import TFile, TH1F, TH2F
  27 +from ROOT import TFile, TH1F, TH2F, TH3F
28 28 from gammalib import *
29 29 from datetime import datetime
30 30 import math
... ... @@ -818,6 +818,113 @@ class caldb():
818 818  
819 819 # Return boundary information
820 820 return bounds
  821 +
  822 + def make_3D_migra(self, array, hdu, name, unit, scale=1.0):
  823 + """
  824 + Make 3D cube as function of ETRUE, MIGRA and THETA form ROOT 3D
  825 + histogram. If the HDU has already the columns, this method will
  826 + simply add another data column.
  827 + If name==None, the method will not append any data column.
  828 +
  829 + Parameters:
  830 + array - ROOT 3D histogram.
  831 + hdu - FITS HDU.
  832 + name - Data column name.
  833 + unit - Data unit.
  834 + Keywords:
  835 + scale - Scaling factor for histogram values.
  836 + """
  837 + # Extract Etrue, Eobs/Etrue and offset angle vectors
  838 + etrue = array.GetXaxis()
  839 + netrue = etrue.GetNbins()
  840 + migra = array.GetYaxis()
  841 + nmigra = migra.GetNbins()
  842 + offsets = array.GetZaxis()
  843 + noffset = offsets.GetNbins()
  844 + ewidth = [] # in MeV
  845 + for ieng in range(netrue):
  846 + ewidth.append(pow(10.0, etrue.GetBinUpEdge(ieng+1) +6.0) - \
  847 + pow(10.0, etrue.GetBinLowEdge(ieng+1)+6.0))
  848 +
  849 + # Attach columns to HDU. We do this only if they do not yet
  850 + # exist. This allows adding columns to the same HDU in successive
  851 + # calls of this method
  852 + #
  853 + # ETRUE_LO
  854 + if not hdu.contains("ETRUE_LO"):
  855 + hdu.append(GFitsTableFloatCol("ETRUE_LO", 1, netrue))
  856 + hdu["ETRUE_LO"].unit("TeV")
  857 + for ieng in range(netrue):
  858 + e_lo = pow(10.0, etrue.GetBinLowEdge(ieng+1))
  859 + hdu["ETRUE_LO"][0,ieng] = e_lo
  860 + #
  861 + # ETRUE_HI
  862 + if not hdu.contains("ETRUE_HI"):
  863 + hdu.append(GFitsTableFloatCol("ETRUE_HI", 1, netrue))
  864 + hdu["ETRUE_HI"].unit("TeV")
  865 + for ieng in range(netrue):
  866 + e_hi = pow(10.0, etrue.GetBinUpEdge(ieng+1))
  867 + hdu["ETRUE_HI"][0,ieng] = e_hi
  868 + #
  869 + # MIGRA_LO
  870 + if not hdu.contains("MIGRA_LO"):
  871 + hdu.append(GFitsTableFloatCol("MIGRA_LO", 1, nmigra))
  872 + hdu["MIGRA_LO"].unit("")
  873 + for ieng in range(nmigra):
  874 + e_lo = pow(10.0, migra.GetBinLowEdge(ieng+1))
  875 + hdu["MIGRA_LO"][0,ieng] = e_lo
  876 + #
  877 + # MIGRA_HI
  878 + if not hdu.contains("MIGRA_HI"):
  879 + hdu.append(GFitsTableFloatCol("MIGRA_HI", 1, nmigra))
  880 + hdu["MIGRA_HI"].unit("")
  881 + for ieng in range(nmigra):
  882 + e_hi = pow(10.0, migra.GetBinUpEdge(ieng+1))
  883 + hdu["MIGRA_HI"][0,ieng] = e_hi
  884 + #
  885 + # THETA_LO
  886 + if not hdu.contains("THETA_LO"):
  887 + hdu.append(GFitsTableFloatCol("THETA_LO", 1, noffset))
  888 + hdu["THETA_LO"].unit("deg")
  889 + for ioff in range(noffset):
  890 + o_lo = offsets.GetBinLowEdge(ioff+1)
  891 + hdu["THETA_LO"][0,ioff] = o_lo
  892 + #
  893 + # THETA_LO
  894 + if not hdu.contains("THETA_HI"):
  895 + hdu.append(GFitsTableFloatCol("THETA_HI", 1, noffset))
  896 + hdu["THETA_HI"].unit("deg")
  897 + for ioff in range(noffset):
  898 + o_hi = offsets.GetBinUpEdge(ioff+1)
  899 + hdu["THETA_HI"][0,ioff] = o_hi
  900 + #
  901 + # "NAME"
  902 + if name != None and not hdu.contains(name):
  903 + hdu.append(GFitsTableFloatCol(name, 1, netrue*nmigra*noffset))
  904 + hdu[name].unit(unit)
  905 + hdu[name].dim([netrue, nmigra, noffset])
  906 + for ioff in range(noffset):
  907 + for imigra in range(nmigra):
  908 + for ieng in range(netrue):
  909 + index = ieng + (imigra + ioff * nmigra) * netrue
  910 + value = array.GetBinContent(ieng+1,imigra+1,ioff+1)
  911 + hdu[name][0,index] = value * scale
  912 +
  913 + # Collect boundary information
  914 + bd_eng = "ETRUE(%.4f-%.2f)TeV" % \
  915 + (pow(10.0, etrue.GetBinLowEdge(1)), \
  916 + pow(10.0, etrue.GetBinUpEdge(netrue)))
  917 + bd_migra = "MIGRA(%.4f-%.2f)TeV" % \
  918 + (pow(10.0, migra.GetBinLowEdge(1)), \
  919 + pow(10.0, migra.GetBinUpEdge(nmigra)))
  920 + bd_off = "THETA(%.2f-%.2f)deg" % \
  921 + (offsets.GetBinLowEdge(1), \
  922 + offsets.GetBinUpEdge(noffset))
  923 + bd_phi = "PHI(0-360)deg"
  924 + bounds = [bd_eng, bd_migra, bd_off, bd_phi]
  925 +
  926 + # Return boundary information
  927 + return bounds
821 928  
822 929 def root2caldb(self, filename, rebin=False, psftype="Gauss", scale=1.0):
823 930 """
... ... @@ -838,7 +945,7 @@ class caldb():
838 945 self.root2psf(file, psftype)
839 946  
840 947 # Create energy dispersion
841   - self.root2edisp(file)
  948 + self.root2edisp_migra(file)
842 949  
843 950 # Create background
844 951 #self.root2bgd(file)
... ... @@ -1126,7 +1233,7 @@ class caldb():
1126 1233 # Return
1127 1234 return
1128 1235  
1129   - def root2edisp(self, file):
  1236 + def root2edisp_2D(self, file):
1130 1237 """
1131 1238 Translate ROOT to CALDB energy dispersion extension. The following ROOT
1132 1239 histograms are used:
... ... @@ -1164,6 +1271,38 @@ class caldb():
1164 1271 # Return
1165 1272 return
1166 1273  
  1274 + def root2edisp_migra(self, file):
  1275 + """
  1276 + Translate ROOT to CALDB energy dispersion extension. The following ROOT
  1277 + histograms are used:
  1278 +
  1279 + EestOverEtrue_offaxis -> MATRIX
  1280 +
  1281 + Parameters:
  1282 + file - ROOT file.
  1283 + Keywords:
  1284 + None
  1285 + """
  1286 + # Continue only if energy dispersion HDU has been opened
  1287 + if self.hdu_edisp != None:
  1288 +
  1289 + # Allocate ROOT 3D array
  1290 + matrix = TH3F()
  1291 + file.GetObject("EestOverEtrue_offaxis", matrix)
  1292 +
  1293 + # Set boundaries
  1294 + bounds = self.make_3D_migra(matrix, self.hdu_edisp, None, "")
  1295 + for b in bounds:
  1296 + self.edisp_bounds.append(b)
  1297 + self.set_cif_keywords(self.hdu_edisp, self.edisp_name, \
  1298 + self.edisp_bounds, self.edisp_desc)
  1299 +
  1300 + # MATRIX
  1301 + self.make_3D_migra(matrix, self.hdu_edisp, "MATRIX", "")
  1302 +
  1303 + # Return
  1304 + return
  1305 +
1167 1306 def root2bgd(self, file):
1168 1307 """
1169 1308 Translate ROOT to CALDB background extension. The following ROOT
... ...