Alexander Ziegler / ctools

Commit 7104b3c2e868f22742160fd8cb805bb2ea744a22

Authored by Jürgen Knödlseder
1 parent 5b44aa68

Normalize DESY effective area to take into account the 68% containment radius.

Inline Side-by-side
Showing 3 changed files with 131 additions and 160 deletions
doc/source/user_manual/getting_started/response.rst
  View file @7104b3c
... ... @@ -88,13 +88,13 @@ the database is installed using
88 88  
89 89 .. code-block:: bash
90 90  
91   - $ [sudo] tar -C $CALDB -zxvf cta-caldb-20140220.tar.gz
  91 + $ [sudo] tar -C $CALDB -zxvf cta-caldb-20150220.tar.gz
92 92  
93 93 or
94 94  
95 95 .. code-block:: bash
96 96  
97   - $ [sudo] tar -C $CALDB -xvf cta-caldb-20140220.tar
  97 + $ [sudo] tar -C $CALDB -xvf cta-caldb-20150220.tar
98 98  
99 99 depending on whether the database has been retrieved as a gzipped file or
100 100 not. Depending on your access rights to the ``$CALDB`` directory,
... ...
examples/show_response.py
  View file @7104b3c
... ... @@ -6,7 +6,7 @@
6 6 # - matplotlib
7 7 # - numpy
8 8 #
9   -# Copyright (C) 2014 Juergen Knoedlseder
  9 +# Copyright (C) 2014-2015 Juergen Knoedlseder
10 10 #
11 11 # This program is free software: you can redistribute it and/or modify
12 12 # it under the terms of the GNU General Public License as published by
... ... @@ -190,7 +190,7 @@ def show_one_sensitivity(rsp, name, color="r", duration=180000.0, alpha=0.2, sig
190 190 Non = Non_lima(sigma, Noff)
191 191 src_counts = Non - bgd_counts
192 192 flux[i] = src_counts / (0.68*aeff * duration) * E[i]*1.0e6 * gammalib.MeV2erg
193   - print logE[i], flux[i]
  193 + #print logE[i], flux[i]
194 194  
195 195 # Plot data
196 196 plt.loglog(E, flux, color+'-', label=name)
... ... @@ -210,12 +210,14 @@ def show_one_sensitivity(rsp, name, color="r", duration=180000.0, alpha=0.2, sig
210 210 # ==================== #
211 211 # Show one sensitivity #
212 212 # ==================== #
213   -def show_one_response(rspname, dbname, name, color="r"):
  213 +def show_one_response(rspname, dbname, name, rootdir=None, color="r"):
214 214 """
215 215 Show one response.
216 216 """
217 217 # Set-up calibration database
218 218 caldb = gammalib.GCaldb()
  219 + if rootdir != None:
  220 + caldb.rootdir(rootdir)
219 221 if gammalib.dir_exists(dbname):
220 222 caldb.rootdir(dbname)
221 223 else:
... ... @@ -242,23 +244,15 @@ def show_one_response(rspname, dbname, name, color="r"):
242 244 # ======================== #
243 245 if __name__ == '__main__':
244 246 """
245   - Display spectrum in PHA format.
  247 + Display response information.
246 248 """
247   - # Print usage information
248   - #usage = "Usage: show_pha filename" \
249   - # " [-n bins] [-c column] [-t title] [-p plot]"
250   - #if len(sys.argv) < 2:
251   - # print(usage)
252   - # sys.exit()
253   -
254   - # Extract parameters
255   - #filename = sys.argv[1]
256   -
257 249 # Create figures
258 250 plt.figure(1)
259 251 plt.title("Effective area")
260 252 plt.figure(2)
261 253 plt.title("Background rate")
  254 + plt.figure(3)
  255 + plt.title("Sensitivity")
262 256  
263 257 # Set response dictionary
264 258 rsps = [
... ... @@ -284,9 +278,13 @@ if __name__ == &#39;__main__&#39;:
284 278 rspname = rsp['rspname']
285 279 color = rsp['color']
286 280 name = rsp['name']
  281 + if rsp.has_key('rootdir'):
  282 + rootdir = rsp['rootdir']
  283 + else:
  284 + rootdir = None
287 285  
288 286 # Show response
289   - show_one_response(rspname, dbname, name, color=color)
  287 + show_one_response(rspname, dbname, name, rootdir=rootdir, color=color)
290 288  
291 289 # Show plot
292 290 plt.show()
... ...
scripts/cta_root2caldb.py
  View file @7104b3c
1 1 #! /usr/bin/env python
2 2 # ==========================================================================
3   -# This script generates IRFs in the CALDB format of HEASARC using ROOT 2D
4   -# performance files.
  3 +# This script generates IRFs in the CALDB format of HEASARC using ROOT
  4 +# offaxis performance files.
5 5 # -------------------------------------------------------------------------
6 6 # Copyright (C) 2011-2015 Juergen Knoedlseder
7 7 # -------------------------------------------------------------------------
... ... @@ -25,7 +25,8 @@
25 25 # - pyroot
26 26 # ==========================================================================
27 27 from ROOT import TFile, TH1F, TH2F, TH3F
28   -from gammalib import *
  28 +#from gammalib import *
  29 +import gammalib
29 30 from datetime import datetime
30 31 import math
31 32 import sys
... ... @@ -205,7 +206,7 @@ class caldb():
205 206 self.bgd_file = "irf_"+version+".fits"
206 207  
207 208 # Open calibration database index
208   - self.cif = GFits(self.base_path+"/caldb.indx", True)
  209 + self.cif = gammalib.GFits(self.base_path+"/caldb.indx", True)
209 210  
210 211 # If file has no CIF extension than create it now
211 212 try:
... ... @@ -222,12 +223,12 @@ class caldb():
222 223  
223 224 # Open files
224 225 if split:
225   - self.ea = GFits(ea_filename, True)
226   - self.psf = GFits(psf_filename, True)
227   - self.edisp = GFits(edisp_filename, True)
228   - self.bgd = GFits(bgd_filename, True)
  226 + self.ea = gammalib.GFits(ea_filename, True)
  227 + self.psf = gammalib.GFits(psf_filename, True)
  228 + self.edisp = gammalib.GFits(edisp_filename, True)
  229 + self.bgd = gammalib.GFits(bgd_filename, True)
229 230 else:
230   - self.irf = GFits(ea_filename, True)
  231 + self.irf = gammalib.GFits(ea_filename, True)
231 232 self.ea = self.irf
232 233 self.psf = self.irf
233 234 self.edisp = self.irf
... ... @@ -314,29 +315,29 @@ class caldb():
314 315 None
315 316 """
316 317 # Create binary table
317   - table = GFitsBinTable()
  318 + table = gammalib.GFitsBinTable()
318 319  
319 320 # Set boundary
320 321  
321 322 # Attach columns. Reference: CAL/GEN/92-008
322   - table.append(GFitsTableStringCol("TELESCOP", 0, 10))
323   - table.append(GFitsTableStringCol("INSTRUME", 0, 10))
324   - table.append(GFitsTableStringCol("DETNAM", 0, 20))
325   - table.append(GFitsTableStringCol("FILTER", 0, 10))
326   - table.append(GFitsTableStringCol("CAL_DEV", 0, 20))
327   - table.append(GFitsTableStringCol("CAL_DIR", 0, 70))
328   - table.append(GFitsTableStringCol("CAL_FILE", 0, 40))
329   - table.append(GFitsTableStringCol("CAL_CLAS", 0, 3))
330   - table.append(GFitsTableStringCol("CAL_DTYP", 0, 4))
331   - table.append(GFitsTableStringCol("CAL_CNAM", 0, 20))
332   - table.append(GFitsTableStringCol("CAL_CBD", 0, 70, 9))
333   - table.append(GFitsTableShortCol("CAL_XNO", 0))
334   - table.append(GFitsTableStringCol("CAL_VSD", 0, 10))
335   - table.append(GFitsTableStringCol("CAL_VST", 0, 8))
336   - table.append(GFitsTableDoubleCol("REF_TIME", 0))
337   - table.append(GFitsTableShortCol("CAL_QUAL", 0))
338   - table.append(GFitsTableStringCol("CAL_DATE", 0, 8))
339   - table.append(GFitsTableStringCol("CAL_DESC", 0, 70))
  323 + table.append(gammalib.GFitsTableStringCol("TELESCOP", 0, 10))
  324 + table.append(gammalib.GFitsTableStringCol("INSTRUME", 0, 10))
  325 + table.append(gammalib.GFitsTableStringCol("DETNAM", 0, 20))
  326 + table.append(gammalib.GFitsTableStringCol("FILTER", 0, 10))
  327 + table.append(gammalib.GFitsTableStringCol("CAL_DEV", 0, 20))
  328 + table.append(gammalib.GFitsTableStringCol("CAL_DIR", 0, 70))
  329 + table.append(gammalib.GFitsTableStringCol("CAL_FILE", 0, 40))
  330 + table.append(gammalib.GFitsTableStringCol("CAL_CLAS", 0, 3))
  331 + table.append(gammalib.GFitsTableStringCol("CAL_DTYP", 0, 4))
  332 + table.append(gammalib.GFitsTableStringCol("CAL_CNAM", 0, 20))
  333 + table.append(gammalib.GFitsTableStringCol("CAL_CBD", 0, 70, 9))
  334 + table.append(gammalib.GFitsTableShortCol("CAL_XNO", 0))
  335 + table.append(gammalib.GFitsTableStringCol("CAL_VSD", 0, 10))
  336 + table.append(gammalib.GFitsTableStringCol("CAL_VST", 0, 8))
  337 + table.append(gammalib.GFitsTableDoubleCol("REF_TIME", 0))
  338 + table.append(gammalib.GFitsTableShortCol("CAL_QUAL", 0))
  339 + table.append(gammalib.GFitsTableStringCol("CAL_DATE", 0, 8))
  340 + table.append(gammalib.GFitsTableStringCol("CAL_DESC", 0, 70))
340 341  
341 342 # Set keywords. Reference: CAL/GEN/92-008
342 343 table.extname("CIF")
... ... @@ -454,7 +455,7 @@ class caldb():
454 455 self.hdu_ea = self.ea.table("EFFECTIVE AREA")
455 456 except:
456 457 # Create binary table
457   - table = GFitsBinTable()
  458 + table = gammalib.GFitsBinTable()
458 459  
459 460 # Set extension name
460 461 table.extname("EFFECTIVE AREA")
... ... @@ -488,7 +489,7 @@ class caldb():
488 489 self.hdu_psf = self.psf.table("POINT SPREAD FUNCTION")
489 490 except:
490 491 # Create binary table
491   - table = GFitsBinTable()
  492 + table = gammalib.GFitsBinTable()
492 493  
493 494 # Set extension name
494 495 table.extname("POINT SPREAD FUNCTION")
... ... @@ -522,7 +523,7 @@ class caldb():
522 523 self.hdu_edisp = self.edisp.table("ENERGY DISPERSION")
523 524 except:
524 525 # Create binary table
525   - table = GFitsBinTable()
  526 + table = gammalib.GFitsBinTable()
526 527  
527 528 # Set extension name
528 529 table.extname("ENERGY DISPERSION")
... ... @@ -556,7 +557,7 @@ class caldb():
556 557 self.hdu_bgd = self.bgd.table("BACKGROUND")
557 558 except:
558 559 # Create binary table
559   - table = GFitsBinTable()
  560 + table = gammalib.GFitsBinTable()
560 561  
561 562 # Set extension name
562 563 table.extname("BACKGROUND")
... ... @@ -629,7 +630,7 @@ class caldb():
629 630 #
630 631 # ENERG_LO
631 632 if not hdu.contains("ENERG_LO"):
632   - hdu.append(GFitsTableFloatCol("ENERG_LO", 1, neng))
  633 + hdu.append(gammalib.GFitsTableFloatCol("ENERG_LO", 1, neng))
633 634 hdu["ENERG_LO"].unit("TeV")
634 635 for ieng in range(neng):
635 636 e_lo = pow(10.0, energies.GetBinLowEdge(ieng+1))
... ... @@ -637,7 +638,7 @@ class caldb():
637 638 #
638 639 # ENERG_HI
639 640 if not hdu.contains("ENERG_HI"):
640   - hdu.append(GFitsTableFloatCol("ENERG_HI", 1, neng))
  641 + hdu.append(gammalib.GFitsTableFloatCol("ENERG_HI", 1, neng))
641 642 hdu["ENERG_HI"].unit("TeV")
642 643 for ieng in range(neng):
643 644 e_hi = pow(10.0, energies.GetBinUpEdge(ieng+1))
... ... @@ -645,7 +646,7 @@ class caldb():
645 646 #
646 647 # THETA_LO
647 648 if not hdu.contains("THETA_LO"):
648   - hdu.append(GFitsTableFloatCol("THETA_LO", 1, noffset))
  649 + hdu.append(gammalib.GFitsTableFloatCol("THETA_LO", 1, noffset))
649 650 hdu["THETA_LO"].unit("deg")
650 651 for ioff in range(noffset):
651 652 o_lo = offsets.GetBinLowEdge(ioff+1)
... ... @@ -653,7 +654,7 @@ class caldb():
653 654 #
654 655 # THETA_LO
655 656 if not hdu.contains("THETA_HI"):
656   - hdu.append(GFitsTableFloatCol("THETA_HI", 1, noffset))
  657 + hdu.append(gammalib.GFitsTableFloatCol("THETA_HI", 1, noffset))
657 658 hdu["THETA_HI"].unit("deg")
658 659 for ioff in range(noffset):
659 660 o_hi = offsets.GetBinUpEdge(ioff+1)
... ... @@ -661,7 +662,7 @@ class caldb():
661 662 #
662 663 # "NAME"
663 664 if name != None and not hdu.contains(name):
664   - hdu.append(GFitsTableFloatCol(name, 1, neng*noffset))
  665 + hdu.append(gammalib.GFitsTableFloatCol(name, 1, neng*noffset))
665 666 hdu[name].unit(unit)
666 667 hdu[name].dim([neng, noffset])
667 668 for ioff in range(noffset):
... ... @@ -734,35 +735,35 @@ class caldb():
734 735 #
735 736 # DETX_LO
736 737 if not hdu.contains("DETX_LO"):
737   - hdu.append(GFitsTableFloatCol("DETX_LO", 1, ndets))
  738 + hdu.append(gammalib.GFitsTableFloatCol("DETX_LO", 1, ndets))
738 739 hdu["DETX_LO"].unit("deg")
739 740 for i in range(ndets):
740 741 hdu["DETX_LO"][0,i] = dets_lo[i]
741 742 #
742 743 # DETX_HI
743 744 if not hdu.contains("DETX_HI"):
744   - hdu.append(GFitsTableFloatCol("DETX_HI", 1, ndets))
  745 + hdu.append(gammalib.GFitsTableFloatCol("DETX_HI", 1, ndets))
745 746 hdu["DETX_HI"].unit("deg")
746 747 for i in range(ndets):
747 748 hdu["DETX_HI"][0,i] = dets_hi[i]
748 749 #
749 750 # DETY_LO
750 751 if not hdu.contains("DETY_LO"):
751   - hdu.append(GFitsTableFloatCol("DETY_LO", 1, ndets))
  752 + hdu.append(gammalib.GFitsTableFloatCol("DETY_LO", 1, ndets))
752 753 hdu["DETY_LO"].unit("deg")
753 754 for i in range(ndets):
754 755 hdu["DETY_LO"][0,i] = dets_lo[i]
755 756 #
756 757 # DETY_HI
757 758 if not hdu.contains("DETY_HI"):
758   - hdu.append(GFitsTableFloatCol("DETY_HI", 1, ndets))
  759 + hdu.append(gammalib.GFitsTableFloatCol("DETY_HI", 1, ndets))
759 760 hdu["DETY_HI"].unit("deg")
760 761 for i in range(ndets):
761 762 hdu["DETY_HI"][0,i] = dets_hi[i]
762 763 #
763 764 # ENERG_LO
764 765 if not hdu.contains("ENERG_LO"):
765   - hdu.append(GFitsTableFloatCol("ENERG_LO", 1, neng))
  766 + hdu.append(gammalib.GFitsTableFloatCol("ENERG_LO", 1, neng))
766 767 hdu["ENERG_LO"].unit("TeV")
767 768 for ieng in range(neng):
768 769 e_lo = pow(10.0, energies.GetBinLowEdge(ieng+1))
... ... @@ -770,7 +771,7 @@ class caldb():
770 771 #
771 772 # ENERG_HI
772 773 if not hdu.contains("ENERG_HI"):
773   - hdu.append(GFitsTableFloatCol("ENERG_HI", 1, neng))
  774 + hdu.append(gammalib.GFitsTableFloatCol("ENERG_HI", 1, neng))
774 775 hdu["ENERG_HI"].unit("TeV")
775 776 for ieng in range(neng):
776 777 e_hi = pow(10.0, energies.GetBinUpEdge(ieng+1))
... ... @@ -778,7 +779,7 @@ class caldb():
778 779 #
779 780 # "NAME"
780 781 if name != None and not hdu.contains(name):
781   - hdu.append(GFitsTableFloatCol(name, 1, ndets*ndets*neng))
  782 + hdu.append(gammalib.GFitsTableFloatCol(name, 1, ndets*ndets*neng))
782 783 hdu[name].unit(unit)
783 784 hdu[name].dim([ndets,ndets,neng])
784 785 for ix in range(ndets):
... ... @@ -852,7 +853,7 @@ class caldb():
852 853 #
853 854 # ETRUE_LO
854 855 if not hdu.contains("ETRUE_LO"):
855   - hdu.append(GFitsTableFloatCol("ETRUE_LO", 1, netrue))
  856 + hdu.append(gammalib.GFitsTableFloatCol("ETRUE_LO", 1, netrue))
856 857 hdu["ETRUE_LO"].unit("TeV")
857 858 for ieng in range(netrue):
858 859 e_lo = pow(10.0, etrue.GetBinLowEdge(ieng+1))
... ... @@ -860,7 +861,7 @@ class caldb():
860 861 #
861 862 # ETRUE_HI
862 863 if not hdu.contains("ETRUE_HI"):
863   - hdu.append(GFitsTableFloatCol("ETRUE_HI", 1, netrue))
  864 + hdu.append(gammalib.GFitsTableFloatCol("ETRUE_HI", 1, netrue))
864 865 hdu["ETRUE_HI"].unit("TeV")
865 866 for ieng in range(netrue):
866 867 e_hi = pow(10.0, etrue.GetBinUpEdge(ieng+1))
... ... @@ -868,7 +869,7 @@ class caldb():
868 869 #
869 870 # MIGRA_LO
870 871 if not hdu.contains("MIGRA_LO"):
871   - hdu.append(GFitsTableFloatCol("MIGRA_LO", 1, nmigra))
  872 + hdu.append(gammalib.GFitsTableFloatCol("MIGRA_LO", 1, nmigra))
872 873 hdu["MIGRA_LO"].unit("")
873 874 for ieng in range(nmigra):
874 875 e_lo = pow(10.0, migra.GetBinLowEdge(ieng+1))
... ... @@ -876,7 +877,7 @@ class caldb():
876 877 #
877 878 # MIGRA_HI
878 879 if not hdu.contains("MIGRA_HI"):
879   - hdu.append(GFitsTableFloatCol("MIGRA_HI", 1, nmigra))
  880 + hdu.append(gammalib.GFitsTableFloatCol("MIGRA_HI", 1, nmigra))
880 881 hdu["MIGRA_HI"].unit("")
881 882 for ieng in range(nmigra):
882 883 e_hi = pow(10.0, migra.GetBinUpEdge(ieng+1))
... ... @@ -884,7 +885,7 @@ class caldb():
884 885 #
885 886 # THETA_LO
886 887 if not hdu.contains("THETA_LO"):
887   - hdu.append(GFitsTableFloatCol("THETA_LO", 1, noffset))
  888 + hdu.append(gammalib.GFitsTableFloatCol("THETA_LO", 1, noffset))
888 889 hdu["THETA_LO"].unit("deg")
889 890 for ioff in range(noffset):
890 891 o_lo = offsets.GetBinLowEdge(ioff+1)
... ... @@ -892,7 +893,7 @@ class caldb():
892 893 #
893 894 # THETA_LO
894 895 if not hdu.contains("THETA_HI"):
895   - hdu.append(GFitsTableFloatCol("THETA_HI", 1, noffset))
  896 + hdu.append(gammalib.GFitsTableFloatCol("THETA_HI", 1, noffset))
896 897 hdu["THETA_HI"].unit("deg")
897 898 for ioff in range(noffset):
898 899 o_hi = offsets.GetBinUpEdge(ioff+1)
... ... @@ -900,7 +901,7 @@ class caldb():
900 901 #
901 902 # "NAME"
902 903 if name != None and not hdu.contains(name):
903   - hdu.append(GFitsTableFloatCol(name, 1, netrue*nmigra*noffset))
  904 + hdu.append(gammalib.GFitsTableFloatCol(name, 1, netrue*nmigra*noffset))
904 905 hdu[name].unit(unit)
905 906 hdu[name].dim([netrue, nmigra, noffset])
906 907 for ioff in range(noffset):
... ... @@ -926,20 +927,23 @@ class caldb():
926 927 # Return boundary information
927 928 return bounds
928 929  
929   - def root2caldb(self, filename, rebin=False, psftype="Gauss", scale=1.0):
  930 + def root2caldb(self, filename, rebin=False, psftype="Gauss", eascale=1.0, bgdscale=1.0):
930 931 """
931 932 Translate ROOT to CALDB information.
932 933  
933 934 Parameters:
934 935 filename - ROOT 2D performance filename.
935 936 Keywords:
936   - rebin - Rebin Etrue histogram (useful for Prod1 IFAE runs)
  937 + rebin - Rebin Etrue histogram (useful for Prod1 IFAE runs)
  938 + psftype - PSF type
  939 + eascale - Effective area scaling factor
  940 + bgdscale - Background scaling factor
937 941 """
938 942 # Open ROOT performance file
939 943 file = TFile(filename)
940 944  
941 945 # Create effective area
942   - self.root2ea(file, rebin=rebin)
  946 + self.root2ea(file, rebin=rebin, eascale=eascale)
943 947  
944 948 # Create point spread function
945 949 self.root2psf(file, psftype)
... ... @@ -949,12 +953,12 @@ class caldb():
949 953  
950 954 # Create background
951 955 #self.root2bgd(file)
952   - self.root2bgd3D(file, scale=scale)
  956 + self.root2bgd3D(file, bgdscale=bgdscale)
953 957  
954 958 # Return
955 959 return
956 960  
957   - def root2ea(self, file, rebin=False):
  961 + def root2ea(self, file, rebin=False, eascale=1.0):
958 962 """
959 963 Translate ROOT to CALDB effective area extension. The following ROOT
960 964 histograms are used:
... ... @@ -965,7 +969,8 @@ class caldb():
965 969 Parameters:
966 970 file - ROOT file.
967 971 Keywords:
968   - rebin - Rebin Etrue histogram (useful for Prod1 IFAE runs)
  972 + rebin - Rebin Etrue histogram (useful for Prod1 IFAE runs)
  973 + eascale - Effective area scaling factor
969 974 """
970 975 # Continue only if effective area HDU has been opened
971 976 if self.hdu_ea != None:
... ... @@ -975,8 +980,6 @@ class caldb():
975 980 ereco = TH2F()
976 981 file.GetObject("EffectiveAreaEtrue_offaxis", etrue)
977 982 file.GetObject("EffectiveArea_offaxis", ereco)
978   - #print(etrue.GetXaxis().GetBinLowEdge(1), etrue.GetXaxis().GetBinUpEdge(etrue.GetXaxis().GetNbins()))
979   - #print(ereco.GetXaxis().GetBinLowEdge(1), ereco.GetXaxis().GetBinUpEdge(ereco.GetXaxis().GetNbins()))
980 983  
981 984 # Rebin etrue histogram
982 985 if rebin:
... ... @@ -988,6 +991,21 @@ class caldb():
988 991 value = etrue.GetBinContent(ieng+1,ioff+1) / 10.0
989 992 etrue.SetBinContent(ieng+1,ioff+1,value)
990 993  
  994 + # Scale histograms (if needed)
  995 + if eascale != 1.0:
  996 + neng = etrue.GetXaxis().GetNbins()
  997 + noffset = etrue.GetYaxis().GetNbins()
  998 + for ioff in range(noffset):
  999 + for ieng in range(neng):
  1000 + value = etrue.GetBinContent(ieng+1,ioff+1) * eascale
  1001 + etrue.SetBinContent(ieng+1,ioff+1,value)
  1002 + neng = ereco.GetXaxis().GetNbins()
  1003 + noffset = ereco.GetYaxis().GetNbins()
  1004 + for ioff in range(noffset):
  1005 + for ieng in range(neng):
  1006 + value = ereco.GetBinContent(ieng+1,ioff+1) * eascale
  1007 + ereco.SetBinContent(ieng+1,ioff+1,value)
  1008 +
991 1009 # Write boundaries (use Ereco boundaries)
992 1010 bounds = self.make_2D(ereco, self.hdu_ea, None, "m2")
993 1011 for b in bounds:
... ... @@ -1339,7 +1357,7 @@ class caldb():
1339 1357 # Return
1340 1358 return
1341 1359  
1342   - def root2bgd3D(self, file, scale=1.0):
  1360 + def root2bgd3D(self, file, bgdscale=1.0):
1343 1361 """
1344 1362 Translate ROOT to CALDB background extension. The following ROOT
1345 1363 histograms are used:
... ... @@ -1367,10 +1385,10 @@ class caldb():
1367 1385 self.bgd_bounds, self.bgd_desc)
1368 1386  
1369 1387 # BGD
1370   - self.make_3D(array, self.hdu_bgd, "BGD", "1/s/MeV/sr", scale=scale)
  1388 + self.make_3D(array, self.hdu_bgd, "BGD", "1/s/MeV/sr", scale=bgdscale)
1371 1389  
1372 1390 # BGD_RECO
1373   - self.make_3D(array, self.hdu_bgd, "BGD_RECO", "1/s/MeV/sr", scale=scale)
  1391 + self.make_3D(array, self.hdu_bgd, "BGD_RECO", "1/s/MeV/sr", scale=bgdscale)
1374 1392  
1375 1393 # Return
1376 1394 return
... ... @@ -1403,9 +1421,9 @@ class caldb():
1403 1421 return
1404 1422  
1405 1423  
1406   -#===================#
  1424 +# ================= #
1407 1425 # Set test database #
1408   -#===================#
  1426 +# ================= #
1409 1427 def set_test():
1410 1428 """
1411 1429 Set Prod1 IFAE database.
... ... @@ -1413,23 +1431,22 @@ def set_test():
1413 1431 # Set database attributes
1414 1432 path = "/project-data/cta/performance/prod1/IFAEOffaxisPerformanceBEI_May2012"
1415 1433 rebin = True
1416   - scale = 1.0
1417 1434  
1418 1435 # Set database content
1419   - db = [{'inst': "e", 'id': "IFAE20120510_50h",
1420   - 'path': path, 'rebin': rebin, 'psftype': "Gauss", 'scale': scale,
  1436 + db = [{'inst': "e", 'id': "IFAE20120510_50h",
  1437 + 'path': path, 'rebin': rebin, 'psftype': "Gauss",
1421 1438 'file': "SubarrayE_IFAE_50hours_20120510_offaxis.root"},
1422   - {'inst': "e", 'id': "IFAE20120510_50h_King",
1423   - 'path': path, 'rebin': rebin, 'psftype': "King", 'scale': scale,
  1439 + {'inst': "e", 'id': "IFAE20120510_50h_King",
  1440 + 'path': path, 'rebin': rebin, 'psftype': "King",
1424 1441 'file': "SubarrayE_IFAE_50hours_20120510_offaxis.root"}]
1425 1442  
1426 1443 # Return database
1427 1444 return db
1428 1445  
1429 1446  
1430   -#=========================#
  1447 +# ======================= #
1431 1448 # Set Prod1 IFAE database #
1432   -#=========================#
  1449 +# ======================= #
1433 1450 def set_prod1_ifae():
1434 1451 """
1435 1452 Set Prod1 IFAE database.
... ... @@ -1438,76 +1455,25 @@ def set_prod1_ifae():
1438 1455 path = "/project-data/cta/performance/prod1/IFAEOffaxisPerformanceBEI_May2012"
1439 1456 rebin = True
1440 1457 psftype = "Gauss"
1441   - scale = 1.0
1442 1458  
1443 1459 # Set database content
1444 1460 db = [{'inst': "b", 'id': "IFAE20120510_50h",
1445   - 'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
  1461 + 'path': path, 'rebin': rebin, 'psftype': psftype,
1446 1462 'file': "SubarrayB_IFAE_50hours_20120510_offaxis.root"},
1447 1463 {'inst': "e", 'id': "IFAE20120510_50h",
1448   - 'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
  1464 + 'path': path, 'rebin': rebin, 'psftype': psftype,
1449 1465 'file': "SubarrayE_IFAE_50hours_20120510_offaxis.root"},
1450 1466 {'inst': "i", 'id': "IFAE20120510_50h",
1451   - 'path': path, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
  1467 + 'path': path, 'rebin': rebin, 'psftype': psftype,
1452 1468 'file': "SubarrayI_IFAE_50hours_20120510_offaxis.root"}]
1453 1469  
1454 1470 # Return database
1455 1471 return db
1456 1472  
1457 1473  
1458   -#=========================#
1459   -# Set Prod2 DESY database #
1460   -#=========================#
1461   -def set_prod2_desy_ignore():
1462   - """
1463   - Set Prod2 DESY database (ignore these files).
1464   - """
1465   - # Set database attributes
1466   - path_aar = "/project-data/cta/performance/prod2/Performance_DESY_20140131_Aar"
1467   - path_ten = "/project-data/cta/performance/prod2/Performance_DESY_20140131_Tenerife_50h"
1468   - rebin = False
1469   - psftype = "Gauss"
1470   - #scale = 1.0/180000.0
1471   - scale = 1.0
1472   -
1473   - # Set database content
1474   - db = [{'inst': "aar", 'id': "DESY20140105_50h",
1475   - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1476   - 'file': "DESY.d20140105.Erec1.V3.ID0NIM2.prod2-Aar-NS.S.2a.180000s.root"},
1477   - {'inst': "aar", 'id': "DESY20140105_50h_0deg",
1478   - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1479   - 'file': "DESY.d20140105.Erec1.V3.ID0_0degNIM2.prod2-Aar-NS.S.2a.180000s.root"},
1480   - {'inst': "aar", 'id': "DESY20140105_50h_180deg",
1481   - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1482   - 'file': "DESY.d20140105.Erec1.V3.ID0_180degNIM2.prod2-Aar-NS.S.2a.180000s.root"},
1483   -
1484   - {'inst': "aar500", 'id': "DESY20140105_50h",
1485   - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1486   - 'file': "DESY.d20140105.Erec1.V3.ID0NIM2.prod2-Aar-500m-NS.S.2a.180000s.root"},
1487   - {'inst': "aar500", 'id': "DESY20140105_50h_0deg",
1488   - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1489   - 'file': "DESY.d20140105.Erec1.V3.ID0_0degNIM2.prod2-Aar-500m-NS.S.2a.180000s.root"},
1490   - {'inst': "aar500", 'id': "DESY20140105_50h_180deg",
1491   - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1492   - 'file': "DESY.d20140105.Erec1.V3.ID0_180degNIM2.prod2-Aar-500m-NS.S.2a.180000s.root"},
1493   -
1494   - {'inst': "tenerife", 'id': "DESY20140105_50h",
1495   - 'path': path_ten, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1496   - 'file': "DESY.d20140105.Erec1.V3.ID0NIM2.prod2-Tenerife-NS.N.2NN.180000s.root"},
1497   - {'inst': "tenerife", 'id': "DESY20140105_50h_0deg",
1498   - 'path': path_ten, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1499   - 'file': "DESY.d20140105.Erec1.V3.ID0_0degNIM2.prod2-Tenerife-NS.N.2NN.180000s.root"},
1500   - {'inst': "tenerife", 'id': "DESY20140105_50h_180deg",
1501   - 'path': path_ten, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
1502   - 'file': "DESY.d20140105.Erec1.V3.ID0_180degNIM2.prod2-Tenerife-NS.N.2NN.180000s.root"}]
1503   -
1504   - # Return database
1505   - return db
1506   -
1507   -
1508   -#=========================#
  1474 +# ======================= #
1509 1475 # Set Prod2 DESY database #
1510   -#=========================#
  1476 +# ======================= #
1511 1477 def set_prod2_desy():
1512 1478 """
1513 1479 Set Prod2 DESY database.
... ... @@ -1517,37 +1483,37 @@ def set_prod2_desy():
1517 1483 path_ten = "/project-data/cta/performance/prod2/Performance_DESY_20140131_Tenerife_50h"
1518 1484 rebin = False
1519 1485 psftype = "Gauss"
1520   - scale = 1.0
  1486 + eascale = 1.0 / 0.68 # Prod 2 DESY effective areas ar for 68% containment
1521 1487  
1522 1488 # Set database content
1523 1489 db = [{'inst': "aar", 'id': "DESY20140105_50h",
1524   - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
  1490 + 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale,
1525 1491 'file': "DESY.d20140105.Erec1.V3.ID0NIM2.prod2-Aar-NS.S.2a.180000s.root"},
1526 1492 {'inst': "aar", 'id': "DESY20140105_50h_0deg",
1527   - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
  1493 + 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale,
1528 1494 'file': "DESY.d20140105.Erec1.V3.ID0_0degNIM2.prod2-Aar-NS.S.2a.180000s.root"},
1529 1495 {'inst': "aar", 'id': "DESY20140105_50h_180deg",
1530   - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
  1496 + 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale,
1531 1497 'file': "DESY.d20140105.Erec1.V3.ID0_180degNIM2.prod2-Aar-NS.S.2a.180000s.root"},
1532 1498  
1533 1499 {'inst': "aar500", 'id': "DESY20140105_50h",
1534   - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
  1500 + 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale,
1535 1501 'file': "DESY.d20140105.Erec1.V3.ID0NIM2.prod2-Aar-500m-NS.S.2a.180000s.root"},
1536 1502 {'inst': "aar500", 'id': "DESY20140105_50h_0deg",
1537   - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
  1503 + 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale,
1538 1504 'file': "DESY.d20140105.Erec1.V3.ID0_0degNIM2.prod2-Aar-NS.S.2a.180000s.root"},
1539 1505 {'inst': "aar500", 'id': "DESY20140105_50h_180deg",
1540   - 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
  1506 + 'path': path_aar, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale,
1541 1507 'file': "DESY.d20140105.Erec1.V3.ID0_180degNIM2.prod2-Aar-NS.S.2a.180000s.root"},
1542 1508  
1543 1509 {'inst': "tenerife", 'id': "DESY20140105_50h",
1544   - 'path': path_ten, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
  1510 + 'path': path_ten, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale,
1545 1511 'file': "DESY.d20140105.Erec1.V3.ID0NIM2.prod2-Tenerife-NS.N.2NN.180000s.root"},
1546 1512 {'inst': "tenerife", 'id': "DESY20140105_50h_0deg",
1547   - 'path': path_ten, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
  1513 + 'path': path_ten, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale,
1548 1514 'file': "DESY.d20140105.Erec1.V3.ID0_0degNIM2.prod2-Tenerife-NS.N.2NN.180000s.root"},
1549 1515 {'inst': "tenerife", 'id': "DESY20140105_50h_180deg",
1550   - 'path': path_ten, 'rebin': rebin, 'psftype': psftype, 'scale': scale,
  1516 + 'path': path_ten, 'rebin': rebin, 'psftype': psftype, 'eascale': eascale,
1551 1517 'file': "DESY.d20140105.Erec1.V3.ID0_180degNIM2.prod2-Tenerife-NS.N.2NN.180000s.root"}
1552 1518 ]
1553 1519  
... ... @@ -1555,9 +1521,9 @@ def set_prod2_desy():
1555 1521 return db
1556 1522  
1557 1523  
1558   -#==========================#
  1524 +# ======================== #
1559 1525 # Main routine entry point #
1560   -#==========================#
  1526 +# ======================== #
1561 1527 if __name__ == '__main__':
1562 1528 """
1563 1529 This script translates a CTA ROOT 2D performance files into a CALDB
... ... @@ -1573,8 +1539,8 @@ if __name__ == &#39;__main__&#39;:
1573 1539 - BGRatePerSqDeg_offaxis (2D) - Background rate per square degree [counts/s/deg2]
1574 1540 - BGRate_offaxis (2D) - Background rate
1575 1541 - EffectiveAreaEtrue_offaxis (2D) - Effective area in true energy (finer bins)
1576   - - MigMatrix_offaxis (3D)
1577   - - EestOverEtrue_offaxis (3D)
  1542 + - MigMatrix_offaxis (3D) - Migration matrix
  1543 + - EestOverEtrue_offaxis (3D) - Migration matrix
1578 1544 - ERes_offaxis (2D) - Energy Resolution
1579 1545 - Ebias_offaxis (2D) - Energy bias
1580 1546 """
... ... @@ -1592,7 +1558,14 @@ if __name__ == &#39;__main__&#39;:
1592 1558 id = entry['id']
1593 1559 rebin = entry['rebin']
1594 1560 psftype = entry['psftype']
1595   - scale = entry['scale']
  1561 + if entry.has_key('eascale'):
  1562 + eascale = entry['eascale']
  1563 + else:
  1564 + eascale = 1.0
  1565 + if entry.has_key('bgdscale'):
  1566 + bgdscale = entry['bgdscale']
  1567 + else:
  1568 + bgdscale = 1.0
1596 1569 filename = path+"/"+entry['file']
1597 1570  
1598 1571 # Allocate caldb entry
... ... @@ -1602,4 +1575,4 @@ if __name__ == &#39;__main__&#39;:
1602 1575 irf.open("file")
1603 1576  
1604 1577 # Translate ROOT to CALDB information
1605   - irf.root2caldb(filename, rebin=rebin, psftype=psftype, scale=scale)
  1578 + irf.root2caldb(filename, rebin=rebin, psftype=psftype, eascale=eascale, bgdscale=bgdscale)
... ...