Add "response" parameter to comobsbin

Jürgen Knödlseder
1 parent d78abf9b
Showing 4 changed files with 82 additions and 39 deletions
doc/source/users/reference_manual/comobsbin.rst
modules/comscripts/comobsbin.par
modules/comscripts/comobsbin.py
test/test_comobsbin.py
@@ -24,6 +24,10 @@ General parameters
 ``outfolder [string]``
     Output folder for files. Put here the location of your COMPTEL data store.
+``response <MODEL|SIM2|SIM3> [string]``
+    Response type specifying whether ``MODEL`` IAQs or simulated IAQs provided
+    in the COMPTEL calibration database should be used.
+
 ``ebinalg <FILE|LIN|LOG> [string]``
     Algorithm for defining energy bins.
@@ -23,6 +23,7 @@ outfolder, s, a, $COMDATASTORE,,, &quot;Output folder for files&quot;
 #
 # Script parameters
 #==================
+response,    s, a, MODEL,MODEL|SIM2|SIM3,, "Response type"
 ebinalg,     s, a, FILE,FILE|LIN|LOG,, "Algorithm for defining energy bins"
 emin,        r, a, 0.75,0.05,50.0, "Minimum energy (MeV)"
 emax,        r, a, 30.0,0.05,50.0, "Maximum energy (MeV)"
@@ -60,6 +60,7 @@ class comobsbin(ctools.csobservation):
             self.obs().load(self['inobs'].filename())
         # Query parameters
+        self['response'].string()
         self._get_ebounds()
         self._get_phases()
         if not self._phases.is_empty():
@@ -582,63 +583,99 @@ class comobsbin(ctools.csobservation):
         dpsi    = self['dpsi'].real()
         dphibar = self['dphibar'].real()
+        # Get calibration database and response type
+        caldb    = gammalib.GCaldb('cgro', 'comptel')
+        response = gammalib.toupper(self['response'].string())
+
         # Initialse list of IAQ names
         iaqnames = []
         # Generate one IAQ for each energy boundary
         for i in range(ebounds.size()):
-            # Set IAQ filename
-            iaqname = '%s/iaq_%6.6d-%6.6dkeV.fits' % \
-                      (self['outfolder'].string(),
-                       ebounds.emin(i).keV(), ebounds.emax(i).keV())
-            iaqfile = gammalib.GFilename(iaqname)
+            # Handle simulated response
+            if response == 'SIM2' or response == 'SIM3':
-            # Set global IAQ filename
-            iaqname_global = '%s/iaq_%6.6d-%6.6dkeV.fits' % \
-                             (self._global_datastore,
-                              ebounds.emin(i).keV(), ebounds.emax(i).keV())
-            iaqfile_global = gammalib.GFilename(iaqname_global)
+                # Write header
+                self._log_header3(gammalib.NORMAL,
+                                  'Use CALDB IAQ for %.3f - %.3f MeV' %
+                                  (ebounds.emin(i).MeV(), ebounds.emax(i).MeV()))
-            # Write header
-            self._log_header3(gammalib.NORMAL, 'Compute IAQ for %.3f - %.3f MeV' % \
-                              (ebounds.emin(i).MeV(), ebounds.emax(i).MeV()))
+                # Build response name
+                iaqname = '%s(%.2f-%.2f)MeV(%.0f)deg' % \
+                          (response, ebounds.emin(i).MeV(), ebounds.emax(i).MeV(),
+                           dphibar)
-            # If IAQ file exists in global datastore then do nothing
-            if iaqfile_global.exists():
-                iaqname = iaqname_global
-                self._log_value(gammalib.NORMAL, 'Global IAQ file exists', iaqfile_global.url())
+                # Raise an exception if response does not exist
+                filename = caldb.filename('', '', 'IAQ', '', '', iaqname)
+                if filename.is_empty():
+                    msg = ('Could not find IAQ "'+iaqname+'" in calibration database. '
+                           'Please select a different response type.')
+                    raise RuntimeError(msg)
-            # If IAQ file exists then do nothing
-            elif iaqfile.exists():
-                self._log_value(gammalib.NORMAL, 'Local IAQ file exists', iaqfile.url())
+                # Log IAQ usage
+                self._log_value(gammalib.NORMAL, 'Use CALDB IAQ', iaqname)
-            # ... otherwise compute and save it
+                # Append IAQ name
+                iaqnames.append(iaqname)
+
+            # ... otherwise handle modeled response
             else:
-                # Initialise IAQ
-                iaq = gammalib.GCOMIaq(55.0, dchi, 50.0, dphibar)
+                # Set IAQ filename
+                iaqname = '%s/iaq_%6.6d-%6.6dkeV.fits' % \
+                          (self['outfolder'].string(),
+                           ebounds.emin(i).keV(), ebounds.emax(i).keV())
+                iaqfile = gammalib.GFilename(iaqname)
+
+                # Set global IAQ filename
+                iaqname_global = '%s/iaq_%6.6d-%6.6dkeV.fits' % \
+                                 (self._global_datastore,
+                                  ebounds.emin(i).keV(), ebounds.emax(i).keV())
+                iaqfile_global = gammalib.GFilename(iaqname_global)
+
+                # Write header
+                self._log_header3(gammalib.NORMAL,
+                                  'Compute IAQ for %.3f - %.3f MeV' %
+                                  (ebounds.emin(i).MeV(), ebounds.emax(i).MeV()))
+
+                # If IAQ file exists in global datastore then do nothing
+                if iaqfile_global.exists():
+                    iaqname = iaqname_global
+                    self._log_value(gammalib.NORMAL,
+                                    'Global IAQ file exists', iaqfile_global.url())
+
+                # If IAQ file exists then do nothing
+                elif iaqfile.exists():
+                    self._log_value(gammalib.NORMAL,
+                                    'Local IAQ file exists', iaqfile.url())
+
+                # ... otherwise compute and save it
+                else:
-                # Set IAQ energy range
-                ebds = gammalib.GEbounds(ebounds.emin(i), ebounds.emax(i))
+                    # Initialise IAQ
+                    iaq = gammalib.GCOMIaq(55.0, dchi, 50.0, dphibar)
-                # Generate continuum IAQ
-                spectrum = gammalib.GModelSpectralPlaw(1.0,-2.0,gammalib.GEnergy(1.0,'MeV'))
+                    # Set IAQ energy range
+                    ebds = gammalib.GEbounds(ebounds.emin(i), ebounds.emax(i))
-                # Compute IAQ
-                iaq.set(spectrum, ebds)
+                    # Generate continuum IAQ
+                    spectrum = gammalib.GModelSpectralPlaw(1.0,-2.0,gammalib.GEnergy(1.0,'MeV'))
-                # Save IAQ
-                iaq.save(iaqfile, True)
+                    # Compute IAQ
+                    iaq.set(spectrum, ebds)
-                # Log creation
-                self._log_value(gammalib.NORMAL, 'Local IAQ file created', iaqfile.url())
+                    # Save IAQ
+                    iaq.save(iaqfile, True)
+
+                    # Log creation
+                    self._log_value(gammalib.NORMAL, 'Local IAQ file created', iaqfile.url())
-                # Log IAQ
-                self._log_string(gammalib.NORMAL, str(iaq))
+                    # Log IAQ
+                    self._log_string(gammalib.NORMAL, str(iaq))
-            # Append IAQ filename
-            iaqnames.append(iaqname)
+                # Append IAQ filename
+                iaqnames.append(iaqname)
         # Return
         return iaqnames
@@ -79,7 +79,7 @@ class Test(test):
         comobsbin = self._comscript('comobsbin')
         # Setup command
-        cmd = comobsbin+' inobs="'+self._obs+'" ebinalg="FILE" '+ \
+        cmd = comobsbin+' inobs="'+self._obs+'" response="MODEL" ebinalg="FILE" '+ \
                         'ebinfile="'+self._ebounds+'" phase="NONE" '+ \
                         'coordsys="GAL" proj="TAN" '+ \
                         'outfolder="comobsbin_cmd1" outobs="comobsbin_cmd1.xml" '+ \
@@ -93,7 +93,7 @@ class Test(test):
         self._check_result('comobsbin_cmd1.xml', 'comobsbin_cmd1')
         # Setup command
-        cmd = comobsbin+' inobs="file_that_does_not_exist.xml" ebinalg="FILE" '+ \
+        cmd = comobsbin+' inobs="file_that_does_not_exist.xml" response="MODEL" ebinalg="FILE" '+ \
                         'ebinfile="'+self._ebounds+'" phase="NONE" '+ \
                         'coordsys="GAL" proj="TAN" '+ \
                         'outfolder="comobsbin_cmd2" outobs="comobsbin_cmd2.xml" '+ \
@@ -117,6 +117,7 @@ class Test(test):
         # Set-up script
         bin = comscripts.comobsbin()
         bin['inobs']     = self._obs
+        bin['response']  = 'MODEL'
         bin['ebinalg']   = 'FILE'
         bin['ebinfile']  = self._ebounds
         bin['phase']     = 'NONE'