working on stacking, broken by gammalib update

Luigi Tibaldo
1 parent 40d2e152
Showing 2 changed files with 64 additions and 19 deletions
cscripts/csphagen.py
test/test_csphagen.py
@@ -195,8 +195,7 @@ 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, 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)
             else:
@@ -207,20 +206,38 @@ class csphagen(ctools.cscript):
  
         # Save PHA, ARF and RMFs
         for obs in outobs:
-            obs.on_spec().save(
-                self._outroot + '_{}_pha_on.fits'.format(obs.id()),
-                True)
-            obs.off_spec().save(
-                self._outroot + '_{}_pha_off.fits'.format(obs.id()),
-                True)
-            obs.arf().save(self._outroot + '_{}_arf.fits'.format(obs.id()),
-                           True)
-            obs.rmf().save(self._outroot + '_{}_rmf.fits'.format(obs.id()),
-                           True)
+            if outobs.size() < 2 or self._stack == False:
+                obs.on_spec().save(
+                    self._outroot + '_{}_pha_on.fits'.format(obs.id()),
+                    True)
+                obs.off_spec().save(
+                    self._outroot + '_{}_pha_off.fits'.format(obs.id()),
+                    True)
+                obs.arf().save(self._outroot + '_{}_arf.fits'.format(obs.id()),
+                               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:
+            if self._logVerbose():
+                msg = 'Stacking {} observations\n'.format(
+                    outobs.size())
+                self._log(msg)
+            outobs = gammalib.GCTAOnOffObservation(outobs)
+            outobs.on_spec().save(
+                self._outroot + '_stacked_pha_on.fits',
+                True)
+            outobs.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)
+
         # Save On/Off observations
         outname = self._outroot + '.xml'
         outobs.save(outname)
@@ -204,6 +204,33 @@ class Test(test):
                                self._nreg_mul[s])
         self._check_outobs('phagen_py3', 2)
  
+        # Test with multiple input observations and stacking
+
+        # Setup csphagen
+        phagen = cscripts.csphagen(obs)
+        phagen['ebinalg'] = 'LOG'
+        phagen['emin'] = 0.1
+        phagen['emax'] = 100.
+        phagen['enumbins'] = nbins
+        phagen['coordsys'] = 'CEL'
+        phagen['ra'] = 83.633
+        phagen['dec'] = 22.0145
+        phagen['rad'] = 0.2
+        phagen['stack'] = True
+        phagen['exclusion'] = self._exclusion
+        phagen['outroot'] = 'phagen_py4'
+        phagen['logfile'] = 'csphagen_py4.log'
+        phagen['chatter'] = 1
+
+        # Run script
+        phagen.execute()
+
+        # Check outout
+        for s in range(2):
+            self._check_output('phagen_py4_stacked', nbins,
+                               0,check_regions=False)
+        self._check_outobs('phagen_py4', 1)
+
         return
  
     def _check_ebounds(self, table, bins):
@@ -341,7 +368,7 @@ class Test(test):
         # Return
         return
  
-    def _check_output(self, filenameroot, bins, nreg):
+    def _check_output(self, filenameroot, bins, nreg, check_regions = True):
         """
         Check the output from a csphagen run
         """
@@ -356,12 +383,13 @@ class Test(test):
  
  
         # Regions
-        reg = gammalib.GSkyRegions(filenameroot + "_on.reg")
-        self.test_value(reg.size(), 1, 'Check for ' + str(
-            1) + ' region in source region file')
-        reg = gammalib.GSkyRegions(filenameroot + "_off.reg")
-        self.test_value(reg.size(), nreg, 'Check for ' + str(
-            nreg) + ' region in background region file')
+        if check_regions:
+            reg = gammalib.GSkyRegions(filenameroot + "_on.reg")
+            self.test_value(reg.size(), 1, 'Check for ' + str(
+                1) + ' region in source region file')
+            reg = gammalib.GSkyRegions(filenameroot + "_off.reg")
+            self.test_value(reg.size(), nreg, 'Check for ' + str(
+                nreg) + ' region in background region file')
  
         return