implemented stacking

Luigi Tibaldo
1 parent 25495305
Showing 2 changed files with 35 additions and 31 deletions
cscripts/csphagen.py
test/test_csphagen.py
@@ -184,6 +184,7 @@ class csphagen(ctools.cscript):
         self._log_header1(gammalib.NORMAL,
                           'Generation of source and background spectra')
         # Loop through observations and generate pha, arf, rmf files
+        # Save ON/OFF regions in ds9 reg files
         outobs = gammalib.GObservations()
         etrue = self._ebounds
         ereco = self._ebounds
@@ -195,9 +196,14 @@ class csphagen(ctools.cscript):
                 bkg_reg.append(region)
             if bkg_reg.size() >= 1 or (
                             self._bkgmethod == "REFLECTED" and bkg_reg.size() >= self._bkgregmin):
-                onoff = gammalib.GCTAOnOffObservation(obs, self._src_dir, etrue, ereco,self._src_reg, bkg_reg)
+                onoff = gammalib.GCTAOnOffObservation(obs, self._src_dir, etrue,
+                                                      ereco, self._src_reg,
+                                                      bkg_reg)
                 onoff.id(obs.id())
                 outobs.append(onoff)
+                self._src_reg.save(
+                    self._outroot + '_{}_on.reg'.format(obs.id()))
+                bkg_reg.save(self._outroot + '_{}_off.reg'.format(obs.id()))
             else:
                 if self._logNormal():
                     msg = 'Observation {} not included in spectra generation\n'.format(
@@ -205,8 +211,8 @@ class csphagen(ctools.cscript):
                     self._log(msg)
         # Save PHA, ARF and RMFs
-        for obs in outobs:
-            if outobs.size() < 2 or self._stack == False:
+        if outobs.size() < 2 or self._stack == False:
+            for obs in outobs:
                 obs.on_spec().save(
                     self._outroot + '_{}_pha_on.fits'.format(obs.id()),
                     True)
@@ -217,26 +223,24 @@ class csphagen(ctools.cscript):
                                True)
                 obs.rmf().save(self._outroot + '_{}_rmf.fits'.format(obs.id()),
                                True)
-            obs.on_regions().save(self._outroot + '_{}_on.reg'.format(obs.id()))
-            obs.off_regions().save(
-                self._outroot + '_{}_off.reg'.format(obs.id()))
-
-        if outobs.size() > 1 and self._stack == True:
+        else:
             if self._logVerbose():
                 msg = 'Stacking {} observations\n'.format(
                     outobs.size())
                 self._log(msg)
-            outobs = gammalib.GCTAOnOffObservation(outobs)
-            outobs.on_spec().save(
+            stackobs = gammalib.GCTAOnOffObservation(outobs)
+            stackobs.on_spec().save(
                 self._outroot + '_stacked_pha_on.fits',
                 True)
-            outobs.off_spec().save(
+            stackobs.off_spec().save(
                 self._outroot + '_stacked_pha_off.fits',
                 True)
-            outobs.arf().save(self._outroot + '_stacked_arf.fits',
-                              True)
-            outobs.rmf().save(self._outroot + '_stacked_rmf.fits',
-                              True)
+            stackobs.arf().save(self._outroot + '_stacked_arf.fits',
+                                True)
+            stackobs.rmf().save(self._outroot + '_stacked_rmf.fits',
+                                True)
+            outobs = gammalib.GObservations()
+            outobs.append(stackobs)
         # Save On/Off observations
         outname = self._outroot + '.xml'
@@ -198,7 +198,7 @@ class Test(test):
         # Run script
         phagen.execute()
-        # Check outout
+        # Check output
         for s in range(2):
             self._check_output('phagen_py3_' + str(s + 1), nbins,
                                self._nreg_mul[s])
@@ -225,10 +225,10 @@ class Test(test):
         # Run script
         phagen.execute()
-        # Check outout
+        # Check output
         for s in range(2):
             self._check_output('phagen_py4_stacked', nbins,
-                               0,check_regions=False)
+                               0, check_regions=False)
         self._check_outobs('phagen_py4', 1)
         return
@@ -240,10 +240,10 @@ class Test(test):
         cols = ['E_MIN', 'E_MAX']
         self.test_value(table.ncols(), len(cols),
-                        'Check for %d columns in light curve FITS table' % len(
+                        'Check for %d columns in energy bounds table' % len(
                             cols))
         self.test_value(table.nrows(), bins,
-                        'Check for %d rows in light curve FITS table' % bins)
+                        'Check for %d rows in energy bounds table' % bins)
         for col in cols:
             self.test_assert(table.contains(col),
                              'FITS file contains "' + col + '" column')
@@ -263,7 +263,7 @@ class Test(test):
         # Check FITS file structure
         self.test_value(fits.size(), 3,
-                        'Check for 3 extensions in light curve FITS file')
+                        'Check for 3 extensions in PHA file')
         self.test_assert(fits.contains('SPECTRUM'),
                          'FITS file contains "SPECTRUM" extension')
         self.test_assert(fits.contains('EBOUNDS'),
@@ -274,10 +274,10 @@ class Test(test):
         # Check FITS table structure
         self.test_value(table.ncols(), len(cols),
-                        'Check for %d columns in light curve FITS table' % len(
+                        'Check for %d columns in PHA table' % len(
                             cols))
         self.test_value(table.nrows(), bins,
-                        'Check for %d rows in light curve FITS table' % bins)
+                        'Check for %d rows in PHA table' % bins)
         for col in cols:
             self.test_assert(table.contains(col),
                              'FITS file contains "' + col + '" column')
@@ -297,14 +297,14 @@ class Test(test):
         Check ARF file
         """
         # Expected column names
-        cols = ['ENERG_LO', 'ENERG_HI', 'SPECRESP', 'BACKGROUND', 'ALPHA']
+        cols = ['ENERG_LO', 'ENERG_HI', 'SPECRESP', 'BACKRESP']
         # Open FITS file
         fits = gammalib.GFits(filename)
         # Check FITS file structure
         self.test_value(fits.size(), 2,
-                        'Check for 2 extensions in light curve FITS file')
+                        'Check for 2 extensions in ARF file')
         self.test_assert(fits.contains('SPECRESP'),
                          'FITS file contains "SPECRESP" extension')
@@ -313,10 +313,10 @@ class Test(test):
         # Check FITS table structure
         self.test_value(table.ncols(), len(cols),
-                        'Check for %d columns in light curve FITS table' % len(
+                        'Check for %d columns in ARF table' % len(
                             cols))
         self.test_value(table.nrows(), bins,
-                        'Check for %d rows in light curve FITS table' % bins)
+                        'Check for %d rows in ARF table' % bins)
         for col in cols:
             self.test_assert(table.contains(col),
                              'FITS file contains "' + col + '" column')
@@ -339,7 +339,7 @@ class Test(test):
         # Check FITS file structure
         self.test_value(fits.size(), 3,
-                        'Check for 3 extensions in light curve FITS file')
+                        'Check for 3 extensions in RMF file')
         self.test_assert(fits.contains('MATRIX'),
                          'FITS file contains "MATRIX" extension')
         self.test_assert(fits.contains('EBOUNDS'),
@@ -350,10 +350,10 @@ class Test(test):
         # Check FITS table structure
         self.test_value(table.ncols(), len(cols),
-                        'Check for %d columns in light curve FITS table' % len(
+                        'Check for %d columns in RMF table' % len(
                             cols))
         self.test_value(table.nrows(), bins,
-                        'Check for %d rows in light curve FITS table' % bins)
+                        'Check for %d rows in RMF table' % bins)
         for col in cols:
             self.test_assert(table.contains(col),
                              'FITS file contains "' + col + '" column')
@@ -368,7 +368,7 @@ class Test(test):
         # Return
         return
-    def _check_output(self, filenameroot, bins, nreg, check_regions = True):
+    def _check_output(self, filenameroot, bins, nreg, check_regions=True):
         """
         Check the output from a csphagen run
         """