Set default deadc value to 0.98

Jürgen Knödlseder
1 parent a1f63642
Showing 17 changed files with 54 additions and 53 deletions
ChangeLog
cscripts/csobsdef.par
cscripts/cspull.par
cscripts/cssens.par
cscripts/cstsdist.par
cscripts/obsutils.py
src/ctmodel/ctmodel.par
src/ctobssim/ctobssim.par
test/data/obs_stacked_two.xml
test/data/pntdef_maximal.dat
test/data/sciver/obs_stacked.xml
test/data/workflow.xml
test/science_verification.py
test/test_csobsdef.py
test/test_cspull.py
test/test_csworkflow.py
test/test_ctobssim.py
@@ -5,6 +5,7 @@
         Set local language to english (#1370)
         Enhance logfile output of csroot2caldb script
+        Set default deadc value to 0.98 (CTA has <2% deadtime)
 2017-02-09
@@ -29,7 +29,7 @@ irf,      s, h, ,,, &quot;Instrument response function&quot;
 emin,     r, h, UNDEF,,, "Lower energy limit (TeV)"
 emax,     r, h, UNDEF,,, "Upper energy limit (TeV)"
 rad,      r, h, UNDEF,,, "Radius of ROI (degrees)"
-deadc,    r, h, 0.95,,, "Deadtime correction factor"
+deadc,    r, h, 0.98,,, "Deadtime correction factor"
 #
 # Standard parameters
@@ -21,7 +21,7 @@ inmodel, f, a, $CTOOLS/share/models/crab.xml,,, &quot;Input model definition XML file
 caldb,   s, a, prod2,,, "Calibration database"
 irf,     s, a, South_0.5h,,, "Instrument response function"
 edisp,   b, h, no,,, "Apply energy dispersion?"
-deadc,   r, h, 0.95,0,1, "Deadtime correction factor"
+deadc,   r, h, 0.98,0,1, "Deadtime correction factor"
 outfile, f, a, pull.dat,,, "Output pull distribution file"
 #
@@ -22,7 +22,7 @@ srcname, s, a, Crab,,, &quot;Source name&quot;
 caldb,   s, a, prod2,,, "Calibration database"
 irf,     s, a, South_0.5h,,, "Instrument response function"
 edisp,   b, h, no,,, "Apply energy dispersion?"
-deadc,   r, h, 0.95,0,1, "Deadtime correction factor"
+deadc,   r, h, 0.98,0,1, "Deadtime correction factor"
 outfile, f, h, sensitivity.dat,,, "Output sensitivity ASCII file"
 #
@@ -26,7 +26,7 @@ bkgcube,   f, a, NONE,,, &quot;Input background cube file (only needed for stacked an
 caldb,     s, a, prod2,,, "Calibration database"
 irf,       s, a, South_0.5h,,, "Instrument response function"
 edisp,     b, h, no,,, "Apply energy dispersion?"
-deadc,     r, h, 0.95,0,1, "Deadtime correction factor"
+deadc,     r, h, 0.98,0,1, "Deadtime correction factor"
 outfile,   f, a, ts.dat,,, "Output Test Statistic distribution file"
 #
 # ==========================================================================
 # Utility functions for observation handling
 #
-# Copyright (C) 2011-2016 Juergen Knoedlseder
+# Copyright (C) 2011-2017 Juergen Knoedlseder
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -178,7 +178,7 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0,
 # ======================= #
 # Set one CTA observation #
 # ======================= #
-def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \
+def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.98, \
             emin=0.1, emax=100.0, rad=5.0, \
             irf='South_50h', caldb='prod2', obsid='000000'):
     """
@@ -273,7 +273,7 @@ def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \
 # ============================ #
 # Set list of CTA observations #
 # ============================ #
-def set_obs_list(obsdeflist, tstart=0.0, duration=1800.0, deadc=0.95, \
+def set_obs_list(obsdeflist, tstart=0.0, duration=1800.0, deadc=0.98, \
                  emin=0.1, emax=100.0, rad=5.0, \
                  irf='South_50h', caldb='prod2'):
     """
@@ -411,7 +411,7 @@ def set_obs_patterns(pattern, ra=83.6331, dec=22.0145, offset=1.5):
 # Set observation container filled with CTA observations #
 # ====================================================== #
 def set_observations(ra, dec, rad, tstart, duration, emin, emax, irf, caldb,
-                     deadc=0.95, pattern='single', offset=1.5):
+                     deadc=0.98, pattern='single', offset=1.5):
     """
     Set an observation container filled with CTA observations
@@ -36,7 +36,7 @@ dec,      r, a, 22.01,-90,90, &quot;Dec of pointing (degrees)&quot;
 rad,      r, a, 5.0,0,180, "Radius of FOV (degrees)"
 tmin,     r, a, 0.0,,, "Start time (MET in s)"
 tmax,     r, a, 1800.0,,, "End time (MET in s)"
-deadc,    r, h, 0.95,0,1, "Average deadtime correction factor"
+deadc,    r, h, 0.98,0,1, "Average deadtime correction factor"
 #
 # Energy binning parameters
@@ -36,7 +36,7 @@ tmin,    r, a, 0.0,,, &quot;Start time (MET in s)&quot;
 tmax,    r, a, 1800.0,,, "End time (MET in s)"
 emin,    r, a, 0.1,,, "Lower energy limit (TeV)"
 emax,    r, a, 100.0,,, "Upper energy limit (TeV)"
-deadc,   r, h, 0.95,0,1, "Average deadtime correction factor"
+deadc,   r, h, 0.98,0,1, "Average deadtime correction factor"
 maxrate, r, h, 1.0e6,,, "Maximum photon rate"
 eslices, i, h, 10,1,100, "Number of energy slices"
@@ -6,7 +6,7 @@
     <parameter name="GoodTimeIntervals" tmin="0" tmax="1800" />
     <parameter name="TimeReference" mjdrefi="51544" mjdreff="0.5" timeunit="s" timesys="TT" timeref="LOCAL" />
     <parameter name="RegionOfInterest" ra="83.6331" dec="21.0145" rad="1.5" />
-    <parameter name="Deadtime" deadc="0.95" />
+    <parameter name="Deadtime" deadc="0.98" />
     <parameter name="Calibration" database="prod2" response="South_0.5h" />
   </observation>
   <observation name="Crab" id="000002" instrument="CTA">
@@ -15,7 +15,7 @@
     <parameter name="GoodTimeIntervals" tmin="1800" tmax="3600" />
     <parameter name="TimeReference" mjdrefi="51544" mjdreff="0.5" timeunit="s" timesys="TT" timeref="LOCAL" />
     <parameter name="RegionOfInterest" ra="83.6331" dec="23.0145" rad="1.5" />
-    <parameter name="Deadtime" deadc="0.95" />
+    <parameter name="Deadtime" deadc="0.98" />
     <parameter name="Calibration" database="prod2" response="South_0.5h" />
   </observation>
 </observation_list>
 name,id,lon,lat,duration,emin,emax,rad,deadc,caldb,irf
-Crab,0001,184.58,-5.78,1800.0,0.1,100.0,5.0,0.95,Prod2,South_50h
-Crab,0002,184.58,-4.78,1800.0,0.1,100.0,5.0,0.95,Prod2,South_50h
-Crab,0003,184.58,-6.78,1800.0,0.1,100.0,5.0,0.95,Prod2,South_50h
+Crab,0001,184.58,-5.78,1800.0,0.1,100.0,5.0,0.98,Prod2,South_50h
+Crab,0002,184.58,-4.78,1800.0,0.1,100.0,5.0,0.98,Prod2,South_50h
+Crab,0003,184.58,-6.78,1800.0,0.1,100.0,5.0,0.98,Prod2,South_50h
@@ -6,7 +6,7 @@
     <parameter name="GoodTimeIntervals" tmin="0" tmax="1800" />
     <parameter name="TimeReference" mjdrefi="51544" mjdreff="0.5" timeunit="s" timesys="TT" timeref="LOCAL" />
     <parameter name="RegionOfInterest" ra="83.6331" dec="21.0145" rad="1.5" />
-    <parameter name="Deadtime" deadc="0.95" />
+    <parameter name="Deadtime" deadc="0.98" />
     <parameter name="Calibration" database="prod2" response="South_0.5h" />
   </observation>
   <observation name="Crab" id="000002" instrument="CTA">
@@ -15,7 +15,7 @@
     <parameter name="GoodTimeIntervals" tmin="1800" tmax="3600" />
     <parameter name="TimeReference" mjdrefi="51544" mjdreff="0.5" timeunit="s" timesys="TT" timeref="LOCAL" />
     <parameter name="RegionOfInterest" ra="83.6331" dec="23.0145" rad="1.5" />
-    <parameter name="Deadtime" deadc="0.95" />
+    <parameter name="Deadtime" deadc="0.98" />
     <parameter name="Calibration" database="prod2" response="South_0.5h" />
   </observation>
 </observation_list>
@@ -6,7 +6,7 @@
             <parameter name="caldb"   value="prod2" />
             <parameter name="irf"     value="South_0.5h" />
             <parameter name="edisp"   value="no" />
-            <parameter name="deadc"   value="0.95" />
+            <parameter name="deadc"   value="0.98" />
             <parameter name="ra"      value="83.63" />
             <parameter name="dec"     value="22.01" />
             <parameter name="target"  value="Crab" />
@@ -11,7 +11,7 @@
 #
 # --------------------------------------------------------------------------
 #
-# Copyright (C) 2015-2016 Juergen Knoedlseder
+# Copyright (C) 2015-2017 Juergen Knoedlseder
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -41,7 +41,7 @@ import cscripts
 # ========================== #
 def generate_pull_distribution(model, obs='NONE', trials=100, \
                                caldb='prod2', irf='South_50h', \
-                               deadc=0.95, edisp=False, \
+                               deadc=0.98, edisp=False, \
                                ra=83.63, dec=22.01, rad=5.0, \
                                emin=0.1, emax=100.0, enumbins=0, \
                                duration=1800.0, \
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This scripts performs unit tests for the csobsdef script.
 #
-# Copyright (C) 2016 Juergen Knoedlseder
+# Copyright (C) 2016-2017 Juergen Knoedlseder
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -186,7 +186,7 @@ class Test(test):
                             'Check for "CTA" instrument')
             self.test_value(o.ontime(), 1800.0, 1.0e-6,
                             'Check for ontime')
-            self.test_value(o.deadc(), 0.95, 1.0e-6,
+            self.test_value(o.deadc(), 0.98, 1.0e-6,
                             'Check for deadtime correction')
             self.test_value(o.events().ebounds().emin().TeV(), 0.1, 1.0e-6,
                             'Check for minimum energy')
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This scripts performs unit tests for the cspull script.
 #
-# Copyright (C) 2016 Juergen Knoedlseder
+# Copyright (C) 2016-2017 Juergen Knoedlseder
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -79,7 +79,7 @@ class Test(test):
                      ' ntrials=2'+ \
                      ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \
                      ' ra=83.6331 dec=22.0145 emin=0.1 emax=100.0'+ \
-                     ' enumbins=0 tmax=100.0 deadc=0.95 rad=5.0'+ \
+                     ' enumbins=0 tmax=100.0 deadc=0.98 rad=5.0'+ \
                      ' npix=200 binsz=0.05'+ \
                      ' logfile="cspull_cmd1.log" chatter=1'
@@ -100,7 +100,7 @@ class Test(test):
                      ' ntrials=2'+ \
                      ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \
                      ' ra=83.6331 dec=22.0145 emin=0.1 emax=100.0'+ \
-                     ' enumbins=0 tmax=100.0 deadc=0.95 rad=5.0'+ \
+                     ' enumbins=0 tmax=100.0 deadc=0.98 rad=5.0'+ \
                      ' npix=200 binsz=0.05'+ \
                      ' logfile="cspull_cmd2.log"'
@@ -129,7 +129,7 @@ class Test(test):
         pull['emax']     = 100.0
         pull['enumbins'] = 0
         pull['tmax']     = 100.0
-        pull['deadc']    = 0.95
+        pull['deadc']    = 0.98
         pull['rad']      = 5.0
         pull['logfile']  = 'cspull_py1.log'
         pull['chatter']  = 2
@@ -155,7 +155,7 @@ class Test(test):
         pull['emax']     = 100.0
         pull['enumbins'] = 10
         pull['tmax']     = 100.0
-        pull['deadc']    = 0.95
+        pull['deadc']    = 0.98
         pull['rad']      = 5.0
         pull['npix']     = 20
         pull['binsz']    = 0.2
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This scripts performs unit tests for the csworkflow script.
 #
-# Copyright (C) 2016 Juergen Knoedlseder
+# Copyright (C) 2016-2017 Juergen Knoedlseder
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -143,13 +143,13 @@ class Test(test):
         self.test_value(models['Crab'][2].value(),
                         5.71055e-16, 1.0e-5, 'Crab prefactor')
         self.test_value(models['Crab'][3].value(),
-                        -2.47559, 1.0e-5, 'Crab index')
+                        -2.463392, 1.0e-5, 'Crab index')
         self.test_value(models['Background'][0].value(),
-                        3.00491, 1.0e-5, 'Background model sigma')
+                        3.008828, 1.0e-5, 'Background model sigma')
         self.test_value(models['Background'][1].value(),
                         6.27131e-05, 1.0e-5, 'Background model prefactor')
         self.test_value(models['Background'][2].value(),
-                        -1.83164, 1.0e-5, 'Background model index')
+                        -1.833184, 1.0e-5, 'Background model index')
         # Return
         return
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This scripts performs unit tests for the ctobssim tool
 #
-# Copyright (C) 2014-2016 Juergen Knoedlseder
+# Copyright (C) 2014-2017 Juergen Knoedlseder
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -88,7 +88,7 @@ class Test(test):
         # Load counts cube and check content.
         evt = gammalib.GCTAEventList('ctobssim_cmd1.fits')
-        self._test_list(evt, 6881)
+        self._test_list(evt, 7105)
         # Setup ctobssim command
         cmd = ctobssim+' inmodel="model_that_does_not_exist.xml"'+ \
@@ -161,14 +161,14 @@ class Test(test):
         # Check content of observation
         self._test_observation(sim)
-        self._test_list(sim.obs()[0].events(), 6881)
+        self._test_list(sim.obs()[0].events(), 7105)
         # Save events
         sim.save()
         # Load counts cube and check content.
         evt = gammalib.GCTAEventList('ctobssim_py1.fits')
-        self._test_list(evt, 6881)
+        self._test_list(evt, 7105)
         # Set-up observation container
         pnts = [{'ra': 83.63, 'dec': 21.01},
@@ -196,10 +196,10 @@ class Test(test):
         # Retrieve observation and check content
         self._test_observation(sim, nobs=4, pnts=pnts)
-        self._test_list(sim.obs()[0].events(), 6003)
-        self._test_list(sim.obs()[1].events(), 6084)
-        self._test_list(sim.obs()[2].events(), 5955)
-        self._test_list(sim.obs()[3].events(), 6030)
+        self._test_list(sim.obs()[0].events(), 6199)
+        self._test_list(sim.obs()[1].events(), 6305)
+        self._test_list(sim.obs()[2].events(), 6092)
+        self._test_list(sim.obs()[3].events(), 6192)
         # Save events
         sim.save()
@@ -208,20 +208,20 @@ class Test(test):
         obs = gammalib.GObservations('ctobssim_py2.xml')
         # Retrieve observation and check content
-        self._test_list(obs[0].events(), 6003)
-        self._test_list(obs[1].events(), 6084)
-        self._test_list(obs[2].events(), 5955)
-        self._test_list(obs[3].events(), 6030)
+        self._test_list(obs[0].events(), 6199)
+        self._test_list(obs[1].events(), 6305)
+        self._test_list(obs[2].events(), 6092)
+        self._test_list(obs[3].events(), 6192)
         # Copy ctobssim tool
         cpy_sim = sim.copy()
         # Retrieve observation and check content of copy
         self._test_observation(cpy_sim, nobs=4, pnts=pnts)
-        self._test_list(cpy_sim.obs()[0].events(), 6003)
-        self._test_list(cpy_sim.obs()[1].events(), 6084)
-        self._test_list(cpy_sim.obs()[2].events(), 5955)
-        self._test_list(cpy_sim.obs()[3].events(), 6030)
+        self._test_list(cpy_sim.obs()[0].events(), 6199)
+        self._test_list(cpy_sim.obs()[1].events(), 6305)
+        self._test_list(cpy_sim.obs()[2].events(), 6092)
+        self._test_list(cpy_sim.obs()[3].events(), 6192)
         # Execute copy of ctobssim tool again, now with a higher chatter
         # level than before
@@ -235,10 +235,10 @@ class Test(test):
         # Check result file
         obs = gammalib.GObservations('ctobssim_py3.xml')
         self.test_value(obs.size(), 4, 'Check for number of observations')
-        self._test_list(obs[0].events(), 6003)
-        self._test_list(obs[1].events(), 6084)
-        self._test_list(obs[2].events(), 5955)
-        self._test_list(obs[3].events(), 6030)
+        self._test_list(obs[0].events(), 6199)
+        self._test_list(obs[1].events(), 6305)
+        self._test_list(obs[2].events(), 6092)
+        self._test_list(obs[3].events(), 6192)
         # Return
         return
@@ -268,8 +268,8 @@ class Test(test):
                  'Observation is CTA observation')
             self.test_value(obs.ontime(), 1800.0, 1.0e-6,
                  'Ontime is 1800 sec')
-            self.test_value(obs.livetime(), 1710.0, 1.0e-6,
-                 'Livetime is 1710 sec')
+            self.test_value(obs.livetime(), 1764.0, 1.0e-6,
+                 'Livetime is 1764 sec')
             self.test_value(pnt.dir().ra_deg(), pnts[i]['ra'], 1.0e-6,
                  'Pointing Right Ascension is '+str(pnts[i]['ra'])+' deg')
             self.test_value(pnt.dir().dec_deg(), pnts[i]['dec'], 1.0e-6,