Add On/Off pull distribution to science verification

Jürgen Knödlseder
1 parent d8468231
Showing 6 changed files with 78 additions and 10 deletions
cscripts/cspull.par
cscripts/cspull.py
doc/source/users/reference_manual/cspull.rst
test/data/sciver/crab_plaw_onoff.xml
test/data/sciver/obs_onoff.xml
test/science_verification.py
@@ -18,7 +18,7 @@
 #============================
 inobs,   f, h, NONE,,, "Input event list, counts cube, or observation definition XML file"
 inmodel, f, a, $CTOOLS/share/models/crab.xml,,, "Input model definition XML file"
-onsrc,   s, a, ,,, "Name of On source (only for On/Off analysis)"
+onsrc,   s, a, NONE,,, "Name of On source (only for On/Off analysis)"
 onrad,   r, a, 0.2,,, "Radius of On region (deg)"
 caldb,   s, a, prod2,,, "Calibration database"
 irf,     s, a, South_0.5h,,, "Instrument response function"
@@ -65,7 +65,7 @@ class cspull(ctools.csobservation):
         enumbins = self['enumbins'].integer()
  
         # Query parameters for On/Off observation
-        if len(self['onsrc'].string()) > 0:
+        if gammalib.toupper(self['onsrc'].string()) != 'NONE':
             self['onrad'].real()
  
         # Query parameters for binned if requested
@@ -112,7 +112,7 @@ class cspull(ctools.csobservation):
         # Extract information from observation
         if obs.classname() == 'GCTAOnOffObservation':
             emin = obs.on_spec().ebounds().emin().TeV()
-            emax = obs.on_spec().ebounds().emin().TeV()
+            emax = obs.on_spec().ebounds().emax().TeV()
             mode = 'On/Off'
         else:
             emin = obs.events().ebounds().emin().TeV()
@@ -164,7 +164,7 @@ class cspull(ctools.csobservation):
         coordsys = 'CEL'
  
         # If we have a On source name then set On region radius
-        if len(onsrc) > 0:
+        if gammalib.toupper(onsrc) != 'NONE':
             onrad = self['onrad'].real()
             emin  = self['emin'].real()
             emax  = self['emax'].real()
@@ -41,6 +41,12 @@ General parameters
 ``inmodel [file]``
     Input model XML file.
  
+``onsrc [string]``
+    Name of On source (only for On/Off analysis; specify ``NONE`` for other analyses).
+
+``onrad [real]``
+    Radius of On region (deg).
+
 ``caldb [string]``
     Calibration database.
  
+<?xml version="1.0" standalone="no"?>
+<source_library title="source library">
+  <source name="Crab" type="PointSource">
+    <spectrum type="PowerLaw">
+      <parameter name="Prefactor"   scale="1e-16" value="5.7"  min="1e-07" max="1000.0" free="1"/>
+      <parameter name="Index"       scale="-1"    value="2.48" min="0.0"   max="+5.0"   free="1"/>
+      <parameter name="PivotEnergy" scale="1e6"   value="0.3"  min="0.01"  max="1000.0" free="0"/>
+    </spectrum>
+    <spatialModel type="PointSource">
+      <parameter name="RA"  scale="1.0" value="83.6331" min="-360" max="360" free="0"/>
+      <parameter name="DEC" scale="1.0" value="22.0145" min="-90"  max="90"  free="0"/>
+    </spatialModel>
+  </source>
+  <source name="Background" type="CTAIrfBackground" instrument="CTA">
+    <spectrum type="PowerLaw">
+      <parameter name="Prefactor"   scale="1.0" value="1.0" min="1e-3" max="1e3"    free="1"/>
+      <parameter name="Index"       scale="1.0" value="0.0" min="-5.0" max="+5.0"   free="1"/>
+      <parameter name="PivotEnergy" scale="1e6" value="1.0" min="0.01" max="1000.0" free="0"/>
+    </spectrum>
+  </source>
+</source_library>
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<observation_list title="observation list">
+  <observation name="Crab" id="0001" instrument="CTA" statistic="cstat">
+    <parameter name="Pointing" ra="83.63" dec="22.51" />
+    <parameter name="GoodTimeIntervals" tmin="0" tmax="900" />
+    <parameter name="TimeReference" mjdrefi="51544" mjdreff="0.5" timeunit="s" timesys="TT" timeref="LOCAL" />
+    <parameter name="Deadtime" deadc="0.98" />
+    <parameter name="Calibration" database="Prod2" response="South_0.5h" />
+  </observation>
+  <observation name="Crab" id="0002" instrument="CTA" statistic="cstat">
+    <parameter name="Pointing" ra="83.63" dec="21.51" />
+    <parameter name="GoodTimeIntervals" tmin="900" tmax="1800" />
+    <parameter name="TimeReference" mjdrefi="51544" mjdreff="0.5" timeunit="s" timesys="TT" timeref="LOCAL" />
+    <parameter name="Deadtime" deadc="0.98" />
+    <parameter name="Calibration" database="Prod2" response="South_0.5h" />
+  </observation>
+</observation_list>
+
@@ -39,8 +39,8 @@ import cscripts
 # ========================== #
 # Generate pull distribution #
 # ========================== #
-def generate_pull_distribution(model, obs='NONE', trials=100, \
-                               caldb='prod2', irf='South_50h', \
+def generate_pull_distribution(model, obs='NONE', onsrc='NONE', onrad=0.2, \
+                               trials=100, caldb='prod2', irf='South_50h', \
                                deadc=0.98, edisp=False, \
                                ra=83.63, dec=22.01, rad=5.0, \
                                emin=0.1, emax=100.0, enumbins=0, \
@@ -57,6 +57,10 @@ def generate_pull_distribution(model, obs=&#39;NONE&#39;, trials=100, \
         Model XML filename (without .xml extension)
     obs : str, optional
         Input observation definition XML filename
+    onsrc : str, optional
+        Name of On source for On/Off analysis
+    onrad : float, optional
+        Radius of On region
     trials : int, optional
         Number of trials
     caldb : str, optional
@@ -108,6 +112,8 @@ def generate_pull_distribution(model, obs=&#39;NONE&#39;, trials=100, \
     pull = cscripts.cspull()
     pull['inobs']    = obs
     pull['inmodel']  = inmodel
+    pull['onsrc']    = onsrc
+    pull['onrad']    = onrad
     pull['outfile']  = outfile
     pull['caldb']    = caldb
     pull['irf']      = irf
@@ -164,7 +170,6 @@ def analyse_pull_distribution(filename):
  
     # Read rows
     first   = True
-    #index   = -1
     samples = 0.0
     for row in reader:
  
@@ -237,6 +242,7 @@ class sciver(gammalib.GPythonTestSuite):
         self.append(self.spec_plaw, 'Test power law model')
         self.append(self.spec_plaw_edisp, 'Test power law model with energy dispersion')
         self.append(self.spec_plaw_stacked, 'Test power law model with stacked analysis')
+        self.append(self.spec_plaw_onoff, 'Test power law model with On/Off analysis')
         self.append(self.spec_plaw2, 'Test power law 2 model')
         self.append(self.spec_smoothbplaw, 'Test smoothly broken power law model')
         self.append(self.spec_eplaw, 'Test exponentially cut off power law model')
@@ -263,7 +269,8 @@ class sciver(gammalib.GPythonTestSuite):
         return
  
     # Generate and analyse pull distributions
-    def pull(self, model, obs='NONE', trials=100, duration=1800.0,
+    def pull(self, model, obs='NONE', onsrc='NONE',
+             trials=100, duration=1800.0,
              ra=83.63, dec=22.01, rad=5.0,
              emin=0.1, emax=100.0, enumbins=0,
              edisp=False, debug=False):
@@ -276,6 +283,8 @@ class sciver(gammalib.GPythonTestSuite):
             Model XML filename (without .xml extension)
         obs : str, optional
             Input observation definition XML filename
+        onsrc : str, optional
+            Name of On source for On/Off analysis
         trials : int, optional
             Number of trials
         duration : float, optional
@@ -298,8 +307,8 @@ class sciver(gammalib.GPythonTestSuite):
             Enable debugging?
         """
         # Generate pull distribution
-        outfile = generate_pull_distribution(model, obs=obs, trials=trials,
-                                             duration=duration,
+        outfile = generate_pull_distribution(model, obs=obs, onsrc=onsrc,
+                                             trials=trials, duration=duration,
                                              ra=ra, dec=dec, rad=rad,
                                              emin=emin, emax=emax,
                                              enumbins=enumbins,
@@ -402,6 +411,20 @@ class sciver(gammalib.GPythonTestSuite):
         self.test('Pull_BackgroundModel_Index')
         return
  
+    # Test power law model with On/Off analysis
+    def spec_plaw_onoff(self):
+        """
+        Test power law model with On/Off analysis
+        """
+        self.pull('data/sciver/crab_plaw_onoff',
+                  obs='data/sciver/obs_onoff.xml', onsrc='Crab',
+                  emin=0.1, emax=100.0, enumbins=20)
+        self.test('Pull_Crab_Prefactor')
+        self.test('Pull_Crab_Index')
+        self.test('Pull_Background_Prefactor')
+        self.test('Pull_Background_Index')
+        return
+
     # Test power law 2 model
     def spec_plaw2(self):
         """