Merge branch '1828-enhance-code-quality' into integration

Jürgen Knödlseder
2 parents a552305f 77f1091b
Showing 10 changed files with 92 additions and 65 deletions
cscripts/cssens.py
examples/pipeline_binned_disk.py
examples/pipeline_binned_mem.py
examples/pipeline_stacked_disk.py
examples/pipeline_stacked_mem.py
test/test_cspull.py
test/test_cssens.py
test/test_ctmapcube.py
test/test_cttsmap.py
test/test_examples.py.in
@@ -373,7 +373,17 @@ class cssens(ctools.cscript):
             logL   = fit.opt().value()
             npred  = fit.obs().npred()
             models = fit.obs().models()
-            ts     = models[self._srcname].ts()
+            source = models[self._srcname]
+            ts     = source.ts()
+
+            # Get fitted Crab, photon and energy fluxes
+            crab_flux   = source['Prefactor'].value() / crab_prefactor
+            photon_flux = source.spectral().flux(emin, emax)
+            energy_flux = source.spectral().eflux(emin, emax)
+
+            # Compute differential sensitivity in unit erg/cm2/s
+            energy      = gammalib.GEnergy(e_mean, 'TeV')
+            sensitivity = source.spectral().eval(energy) * e_mean*erg_mean * 1.0e6
             # Assess quality based on a comparison between Npred and Nevents
             quality = npred - nevents
@@ -400,27 +410,30 @@ class cssens(ctools.cscript):
                     for par in model:
                         self._log(str(par)+'\n')
-            # Start over if TS was non-positive
+            # If TS was non-positive then increase the test flux and start
+            # over
             if ts <= 0.0:
-                if self._logExplicit():
-                    self._log('Non positive TS. Start over.\n')
-                continue
-            # Get fitted Crab, photon and energy fluxes
-            crab_flux   = models[self._srcname]['Prefactor'].value() / \
-                          crab_prefactor
-            photon_flux = models[self._srcname].spectral().flux(emin, emax)
-            energy_flux = models[self._srcname].spectral().eflux(emin, emax)
+                # If the number of iterations was exceeded then stop
+                if (iterations >= self._max_iter):
+                    if self._logTerse():
+                        self._log(' Test ended after %d iterations.\n' %
+                                  self._max_iter)
+                    break
-            # Compute differential sensitivity in unit erg/cm2/s
-            energy      = gammalib.GEnergy(e_mean, 'TeV')
-            sensitivity = models[self._srcname].spectral().eval(energy) * \
-                          e_mean*erg_mean * 1.0e6
+                # Increase test flux by a factor of 2
+                test_crab_flux = test_crab_flux * 2.0
+
+                # Signal start we start over
+                if self._logExplicit():
+                    self._log('Non positive TS, increase test flux and start '
+                              'over.\n')
+
+                # ... and start over
+                continue
             # Compute flux correction factor based on average TS
-            correct = 1.0
-            if ts > 0:
-                correct = math.sqrt(self._ts_thres/ts)
+            correct = math.sqrt(self._ts_thres/ts)
             # Compute extrapolated fluxes
             crab_flux   = correct * crab_flux
@@ -432,7 +445,7 @@ class cssens(ctools.cscript):
             energy_flux_value.append(energy_flux)
             sensitivity_value.append(sensitivity)
-            # Write background and test source fit results
+            # Write fit results
             if self._logExplicit():
                 self._log.parformat('Photon flux')
                 self._log(photon_flux)
@@ -453,7 +466,7 @@ class cssens(ctools.cscript):
                     for par in model:
                         self._log(str(par)+'\n')
             elif self._logTerse():
-                self._log.parformat('Iteration '+str(iterations))
+                self._log.parformat('Iteration %d' % iterations)
                 self._log('TS=')
                 self._log(ts)
                 self._log(' ')
@@ -589,7 +602,19 @@ class cssens(ctools.cscript):
         results  = []
         # Set test source model for this observation
-        models = modutils.test_source(self._obs.models(), self._srcname)
+        models = modutils.test_source(self._obs.models(), self._srcname,
+                                      ra=self._ra, dec=self._dec)
+
+        # Write observations into logger
+        if self._logTerse():
+            self._log('\n')
+            self._log.header1(gammalib.number('Observation',len(self._obs)))
+            self._log(str(self._obs))
+            self._log('\n')
+        if self._logExplicit():
+            for obs in self._obs:
+                self._log(str(obs))
+                self._log('\n')
         # Write models into logger
         if self._logTerse():
@@ -155,4 +155,5 @@ if __name__ == &#39;__main__&#39;:
     obs.models(gammalib.GModels(datadir+'/crab.xml'))
     # Run analysis pipeline
-    run_pipeline(obs, model=datadir+'/crab.xml')
+    run_pipeline(obs, model=datadir+'/crab.xml', enumbins=10, nxpix=40,
+                 nypix=40, binsz=0.1)
@@ -126,4 +126,4 @@ if __name__ == &#39;__main__&#39;:
     obs.models(gammalib.GModels(datadir+'/crab.xml'))
     # Run analysis pipeline
-    run_pipeline(obs)
+    run_pipeline(obs, enumbins=10, nxpix=40, nypix=40, binsz=0.1)
@@ -205,4 +205,5 @@ if __name__ == &#39;__main__&#39;:
     obs.models(gammalib.GModels(datadir+'/crab.xml'))
     # Run analysis pipeline
-    run_pipeline(obs, model=datadir+'/crab.xml')
+    run_pipeline(obs, model=datadir+'/crab.xml', enumbins=10, nxpix=40,
+                 nypix=40, binsz=0.1)
@@ -178,4 +178,4 @@ if __name__ == &#39;__main__&#39;:
     obs.models(gammalib.GModels(datadir+'/crab.xml'))
     # Run analysis pipeline
-    run_pipeline(obs)
+    run_pipeline(obs, enumbins=10, nxpix=40, nypix=40, binsz=0.1)
@@ -146,7 +146,7 @@ class Test(test):
         pull = cscripts.cspull()
         pull['inmodel']  = self._model
         pull['outfile']  = 'cspull_py2.dat'
-        pull['ntrials']  = 2
+        pull['ntrials']  = 1
         pull['caldb']    = self._caldb
         pull['irf']      = self._irf
         pull['ra']       = 83.6331
@@ -157,8 +157,8 @@ class Test(test):
         pull['tmax']     = 100.0
         pull['deadc']    = 0.95
         pull['rad']      = 5.0
-        pull['npix']     = 100
-        pull['binsz']    = 0.02
+        pull['npix']     = 20
+        pull['binsz']    = 0.2
         pull['coordsys'] = 'CEL'
         pull['proj']     = 'TAN'
         pull['logfile']  = 'cspull_py2.log'
@@ -168,14 +168,14 @@ class Test(test):
         pull.execute()
         # Check pull distribution file
-        self._check_pull_file('cspull_py2.dat')
+        self._check_pull_file('cspull_py2.dat', rows=2)
         # Set-up cspull from event list
         pull = cscripts.cspull()
         pull['inobs']    = self._events
         pull['inmodel']  = self._model
         pull['outfile']  = 'cspull_py3.dat'
-        pull['ntrials']  = 2
+        pull['ntrials']  = 1
         pull['caldb']    = self._caldb
         pull['irf']      = self._irf
         pull['enumbins'] = 0
@@ -186,7 +186,7 @@ class Test(test):
         pull.execute()
         # Check pull distribution file
-        self._check_pull_file('cspull_py3.dat')
+        self._check_pull_file('cspull_py3.dat', rows=2)
         # Build observation container with unbinned observation
         cta = gammalib.GCTAObservation(self._events)
@@ -197,7 +197,7 @@ class Test(test):
         pull = cscripts.cspull(obs)
         pull['inmodel']  = self._model
         pull['outfile']  = 'cspull_py4.dat'
-        pull['ntrials']  = 2
+        pull['ntrials']  = 1
         pull['caldb']    = self._caldb
         pull['irf']      = self._irf
         pull['enumbins'] = 0
@@ -208,7 +208,7 @@ class Test(test):
         pull.execute()
         # Check pull distribution file
-        self._check_pull_file('cspull_py4.dat')
+        self._check_pull_file('cspull_py4.dat', rows=2)
         # Set-up stacked cspull with one observation where response cubes
         # are specified in the input observation
@@ -216,7 +216,7 @@ class Test(test):
         pull['inobs']    = self._inobs
         pull['inmodel']  = self._stacked_model
         pull['outfile']  = 'cspull_py5.dat'
-        pull['ntrials']  = 2
+        pull['ntrials']  = 1
         pull['enumbins'] = 0
         pull['logfile']  = 'cspull_py5.log'
         pull['chatter']  = 4
@@ -225,7 +225,7 @@ class Test(test):
         pull.execute()
         # Check pull distribution file
-        self._check_pull_file('cspull_py5.dat')
+        self._check_pull_file('cspull_py5.dat', rows=2)
         # Set-up stacked cspull with two observations from which response
         # cubes will be computed internally. The IRFs are also specified
@@ -235,7 +235,7 @@ class Test(test):
         pull['inobs']    = self._inobs_two
         pull['inmodel']  = self._model
         pull['outfile']  = 'cspull_py6.dat'
-        pull['ntrials']  = 2
+        pull['ntrials']  = 1
         pull['emin']     = 0.02
         pull['emax']     = 100.0
         pull['enumbins'] = 10
@@ -250,13 +250,13 @@ class Test(test):
         pull.execute()
         # Check pull distribution file
-        self._check_pull_file('cspull_py6.dat')
+        self._check_pull_file('cspull_py6.dat', rows=2)
         # Return
         return
     # Check pull file
-    def _check_pull_file(self, filename):
+    def _check_pull_file(self, filename, rows=3):
         """
         Check pull file
         """
@@ -264,8 +264,8 @@ class Test(test):
         pulls = gammalib.GCsv(filename)
         # Check dimensions
-        self.test_value(pulls.nrows(), 3,
-                        'Check for 3 rows in pull file')
+        self.test_value(pulls.nrows(), rows,
+                        'Check for number of rows in pull file')
         # Return
         return
@@ -72,7 +72,7 @@ class Test(test):
         cmd = cssens+' inmodel="'+self._model+'"'+ \
                      ' srcname="Crab"'+ \
                      ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \
-                     ' duration=500.0 rad=3.0 emin=1.0 emax=10.0 bins=1'+ \
+                     ' duration=1800.0 rad=3.0 emin=1.0 emax=10.0 bins=1'+ \
                      ' outfile="cssens_cmd1.dat"'+ \
                      ' logfile="cssens_cmd1.log" chatter=1'
@@ -91,7 +91,7 @@ class Test(test):
         cmd = cssens+' inmodel="model_that_does_not_exist.xml"'+ \
                      ' srcname="Crab"'+ \
                      ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \
-                     ' duration=500.0 rad=3.0 emin=1.0 emax=10.0 bins=1'+ \
+                     ' duration=1800.0 rad=3.0 emin=1.0 emax=10.0 bins=1'+ \
                      ' outfile="cssens_cmd2.dat"'+ \
                      ' logfile="cssens_cmd2.log" chatter=2'
@@ -115,7 +115,7 @@ class Test(test):
         sens['caldb']    = self._caldb
         sens['irf']      = self._irf
         sens['outfile']  = 'cssens_py1.dat'
-        sens['duration'] = 500.0
+        sens['duration'] = 1800.0
         sens['rad']      = 3.0
         sens['emin']     = 1.0
         sens['emax']     = 10.0
@@ -71,8 +71,8 @@ class Test(test):
         # Setup ctmapcube command
         cmd = ctmapcube+' inmodel="'+self._model+'"'+ \
                         ' outcube="ctmapcube_cmd1.fits"'+\
-                        ' emin=0.1 emax=100.0 enumbins=20 ebinalg=LOG'+ \
-                        ' nxpix=200 nypix=200 binsz=0.02 coordsys=CEL'+ \
+                        ' emin=0.1 emax=100.0 enumbins=10 ebinalg=LOG'+ \
+                        ' nxpix=100 nypix=100 binsz=0.04 coordsys=CEL'+ \
                         ' xref=83.63 yref=22.01 proj=CAR'+ \
                         ' logfile="ctmapcube_cmd1.log" chatter=1'
@@ -90,8 +90,8 @@ class Test(test):
         # Setup ctmapcube command
         cmd = ctmapcube+' inmodel="events_that_do_not_exist.fits"'+ \
                         ' outcube="ctmapcube_cmd2.fits"'+\
-                        ' emin=0.1 emax=100.0 enumbins=20 ebinalg=LOG'+ \
-                        ' nxpix=200 nypix=200 binsz=0.02 coordsys=CEL'+ \
+                        ' emin=0.1 emax=100.0 enumbins=10 ebinalg=LOG'+ \
+                        ' nxpix=100 nypix=100 binsz=0.04 coordsys=CEL'+ \
                         ' xref=83.63 yref=22.01 proj=CAR'+ \
                         ' logfile="ctmapcube_cmd2.log"'
@@ -114,10 +114,10 @@ class Test(test):
         mapcube['ebinalg']  = 'LOG'
         mapcube['emin']     = 0.1
         mapcube['emax']     = 100.0
-        mapcube['enumbins'] = 20
-        mapcube['nxpix']    = 200
-        mapcube['nypix']    = 200
-        mapcube['binsz']    = 0.02
+        mapcube['enumbins'] = 10
+        mapcube['nxpix']    = 100
+        mapcube['nypix']    = 100
+        mapcube['binsz']    = 0.04
         mapcube['coordsys'] = 'CEL'
         mapcube['proj']     = 'CAR'
         mapcube['xref']     = 83.63
@@ -152,12 +152,12 @@ class Test(test):
         cube = gammalib.GModelSpatialDiffuseCube(filename)
         # Test map cube
-        #self.test_value(cube.maps(), 20, "20 maps")
-        #self.test_value(cube.pixels(), 40000, "40000 map pixels")
+        #self.test_value(cube.maps(), 10, "10 maps")
+        #self.test_value(cube.pixels(), 10000, "10000 map pixels")
         # The map cube is not loaded by default !!!! We should add the method so that
         # the attributes are known !!!!
-        self.test_value(cube.maps(), 0, '20 maps')
-        self.test_value(cube.pixels(), 0, '40000 map pixels')
+        self.test_value(cube.maps(), 0, '10 maps')
+        self.test_value(cube.pixels(), 0, '10000 map pixels')
         # Return
         return
@@ -69,7 +69,7 @@ class Test(test):
                       ' inmodel="'+self._model+'" srcname="Crab"'+ \
                       ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \
                       ' outmap="cttsmap_cmd1.fits"'+ \
-                      ' nxpix=5 nypix=5 binsz=0.02'+ \
+                      ' nxpix=3 nypix=3 binsz=0.02'+ \
                       ' coordsys="CEL" proj="CAR" xref=83.63 yref=22.01'+ \
                       ' logfile="cttsmap_cmd1.log" chatter=1'
@@ -89,7 +89,7 @@ class Test(test):
                       ' inmodel="'+self._model+'" srcname="Crab"'+ \
                       ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \
                       ' outmap="cttsmap_cmd2.fits"'+ \
-                      ' nxpix=5 nypix=5 binsz=0.02'+ \
+                      ' nxpix=3 nypix=3 binsz=0.02'+ \
                       ' coordsys="CEL" proj="CAR" xref=83.63 yref=22.01'+ \
                       ' logfile="cttsmap_cmd2.log" chatter=2'
@@ -113,8 +113,8 @@ class Test(test):
         tsmap['caldb']    = self._caldb
         tsmap['irf']      = self._irf
         tsmap['outmap']   = 'cttsmap_py1.fits'
-        tsmap['nxpix']    = 5
-        tsmap['nypix']    = 5
+        tsmap['nxpix']    = 3
+        tsmap['nypix']    = 3
         tsmap['binsz']    = 0.02
         tsmap['coordsys'] = 'CEL'
         tsmap['proj']     = 'CAR'
@@ -149,14 +149,14 @@ class Test(test):
         # Check dimensions
         self.test_value(ts.naxis(), 2, 'Check for 2 dimensions')
-        self.test_value(ts.naxes(0), 5, 'Check for 5 pixels in X')
-        self.test_value(ts.naxes(1), 5, 'Check for 5 pixels in Y')
+        self.test_value(ts.naxes(0), 3, 'Check for 3 pixels in X')
+        self.test_value(ts.naxes(1), 3, 'Check for 3 pixels in Y')
         self.test_value(prefactor.naxis(), 2, 'Check for 2 dimensions')
-        self.test_value(prefactor.naxes(0), 5, 'Check for 5 pixels in X')
-        self.test_value(prefactor.naxes(1), 5, 'Check for 5 pixels in Y')
+        self.test_value(prefactor.naxes(0), 3, 'Check for 3 pixels in X')
+        self.test_value(prefactor.naxes(1), 3, 'Check for 3 pixels in Y')
         self.test_value(index.naxis(), 2, 'Check for 2 dimensions')
-        self.test_value(index.naxes(0), 5, 'Check for 5 pixels in X')
-        self.test_value(index.naxes(1), 5, 'Check for 5 pixels in Y')
+        self.test_value(index.naxes(0), 3, 'Check for 3 pixels in X')
+        self.test_value(index.naxes(1), 3, 'Check for 3 pixels in Y')
         # Return
         return
@@ -186,7 +186,7 @@ class Test(gammalib.GPythonTestSuite):
         Test make_ts_distributions
         """
         # Set arguments
-        args = '-n 2 -e 0 -t 1800 -d %s' % self._datadir
+        args = '-n 1 -e 0 -t 100 -d %s' % self._datadir
         # Execute script
         self._execute_python('make_ts_distributions', args=args)