Tune cspull scripts and associated scripts

(still need to remove some debug statements later)

Tune cspull scripts and associated scripts
(still need to remove some debug statements later)
Jürgen Knödlseder
1 parent 7de3c72c
Showing 5 changed files with 289 additions and 209 deletions
cscripts/cspull.py
cscripts/obsutils.py
src/ctexpcube/ctexpcube.cpp
src/ctlike/ctlike.cpp
src/ctlike/ctlike.hpp
@@ -30,22 +30,16 @@ from cscripts import obsutils
 # ============ #
 class cspull(ctools.cscript):
     """
-    Generates pull distributions for a model.
-    
-    This class implements the pull distribution generation script. It derives
-    from the ctools.cscript class which provides support for parameter
-    files, command line arguments, and logging. In that way the Python
-    script behaves just as a regular ctool.
+    Generates pull distributions for a model
     """
-
     # Constructor
     def __init__(self, *argv):
         """
         Constructor.
         """
         # Set name
-        self._name    = "cspull"
-        self._version = "1.1.0"
+        self._name    = 'cspull'
+        self._version = '1.1.0'
  
         # Initialise some members
         self._edisp       = False
@@ -54,15 +48,15 @@ class cspull(ctools.cscript):
         self._bckcube     = None
         self._edispcube   = None
         self._stackmodels = None
-        self._coordsys    = "CEL"
-        self._proj        = "TAN"
-        self._log_clients = False
+        self._coordsys    = 'CEL'
+        self._proj        = 'TAN'
+        self._log_clients = False  # Static parameter
         self._offset      = 0.0
         self._ntrials     = 0
         self._outfile     = gammalib.GFilename()
         self._npix        = 0
         self._binsz       = 0.0
-        self._pattern     = "single"
+        self._pattern     = 'single'
         self._enumbins    = 0
         self._seed        = 1
         self._chatter     = 2
@@ -81,7 +75,7 @@ class cspull(ctools.cscript):
     # Private methods
     def _get_parameters(self):
         """
-        Get parameters from parfile.
+        Get parameters from parfile
         """
         # If there are no observations in container then get some ...
         if self._obs.size() == 0:
@@ -89,8 +83,8 @@ class cspull(ctools.cscript):
  
             # Check for requested pattern and use above
             # observation parameters to set wobble pattern
-            self._pattern = self["pattern"].string()
-            if self._pattern == "four":
+            self._pattern = self['pattern'].string()
+            if self._pattern == 'four':
                 self._obs = self._set_obs()
  
         # ... otherwise add response information and energy boundaries
@@ -99,153 +93,173 @@ class cspull(ctools.cscript):
             self._setup_observations(self._obs)
  
         # Get number of energy bins
-        self._enumbins = self["enumbins"].integer()
+        self._enumbins = self['enumbins'].integer()
  
         # Read parameters for binned if requested
         if self._enumbins != 0:
-            self._npix     = self["npix"].integer()
-            self._binsz    = self["binsz"].real()
-            self._coordsys = self["coordsys"].string()
-            self._proj     = self["proj"].string()
+            self._npix     = self['npix'].integer()
+            self._binsz    = self['binsz'].real()
+            self._coordsys = self['coordsys'].string()
+            self._proj     = self['proj'].string()
  
         # Set models if we have none
         if self._obs.models().size() == 0:
-            self._obs.models(self["inmodel"].filename())
+            self._obs.models(self['inmodel'].filename())
  
         # Read other parameters    
-        self._outfile = self["outfile"].filename()
-        self._ntrials = self["ntrials"].integer()
-        self._edisp   = self["edisp"].boolean()
-        self._offset  = self["offset"].real()
-        self._seed    = self["seed"].integer()
-        self._chatter = self["chatter"].integer()
+        self._outfile = self['outfile'].filename()
+        self._ntrials = self['ntrials'].integer()
+        self._edisp   = self['edisp'].boolean()
+        self._offset  = self['offset'].real()
+        self._seed    = self['seed'].integer()
+        self._chatter = self['chatter'].integer()
  
         # Query some parameters
-        self["profile"].boolean()
+        self['profile'].boolean()
  
         #  Write input parameters into logger
         if self._logTerse():
             self._log_parameters()
-            self._log("\n")
+            self._log('\n')
  
         # Return
         return
  
     def _set_obs(self):
         """
-        Set CTA observation(s).
+        Set CTA observation(s)
  
-        Returns:
-            An observation container.
+        Returns
+        -------
+        obs : `~gammalib.GObservations`
+            Observation container.
         """
         # Setup observation definition list
         obsdeflist = obsutils.set_obs_patterns(self._pattern,
-                                               ra=self["ra"].real(),
-                                               dec=self["dec"].real(),
-                                               offset=self["offset"].real())
+                                               ra=self['ra'].real(),
+                                               dec=self['dec'].real(),
+                                               offset=self['offset'].real())
  
         # Create list of observations
         obs = obsutils.set_obs_list(obsdeflist,
-                                    tstart=self["tmin"].real(),
-                                    duration=self["tmax"].real()-self["tmin"].real(),
-                                    deadc=self["deadc"].real(),
-                                    emin=self["emin"].real(),
-                                    emax=self["emax"].real(),
-                                    rad=self["rad"].real(),
-                                    irf=self["irf"].string(),
-                                    caldb=self["caldb"].string())
+                                    tstart=self['tmin'].real(),
+                                    duration=self['tmax'].real()-self['tmin'].real(),
+                                    deadc=self['deadc'].real(),
+                                    emin=self['emin'].real(),
+                                    emax=self['emax'].real(),
+                                    rad=self['rad'].real(),
+                                    irf=self['irf'].string(),
+                                    caldb=self['caldb'].string())
  
         # Return observation container
         return obs
  
     def _set_stacked_irf(self):
         """
-        For stacked analysis prepare stacked irfs
+        Prepare stacked IRFs for stacked analysis
         """
         # Write header into logger
         if self._logTerse():
-            self._log("\n")
-            self._log.header1("Compute stacked response")
+            self._log('\n')
+            self._log.header1('Compute stacked response')
  
         # Get stacked exposure
         expcube = ctools.ctexpcube(self._obs)
-        expcube["incube"]   = "NONE"
-        expcube["usepnt"]   = True
-        expcube["ebinalg"]  = "LOG"
-        expcube["binsz"]    = self._binsz
-        expcube["nxpix"]    = self._npix
-        expcube["nypix"]    = self._npix
-        expcube["enumbins"] = self._enumbins
-        expcube["emin"]     = self["emin"].real()
-        expcube["emax"]     = self["emax"].real()
-        expcube["coordsys"] = self._coordsys
-        expcube["proj"]     = self._proj
+        expcube['incube']   = 'NONE'
+        expcube['usepnt']   = True
+        expcube['ebinalg']  = 'LOG'
+        expcube['binsz']    = self._binsz
+        expcube['nxpix']    = self._npix
+        expcube['nypix']    = self._npix
+        expcube['enumbins'] = self._enumbins
+        expcube['emin']     = self['emin'].real()
+        expcube['emax']     = self['emax'].real()
+        expcube['coordsys'] = self._coordsys
+        expcube['proj']     = self._proj
+        expcube['chatter']  = self._chatter
+        expcube['debug']    = self._logDebug()
         expcube.run()
+        # JK: Debug
+        expcube['outcube'] = 'expcube.fits'
+        expcube.save()
  
         # Notify exposure computation
         if self._logTerse():
-            self._log("Computed exposure cube\n")
+            self._log('Computed exposure cube\n')
  
         # Get stacked Psf
         psfcube = ctools.ctpsfcube(self._obs)
-        psfcube["incube"]   = "NONE"
-        psfcube["usepnt"]   = True
-        psfcube["ebinalg"]  = "LOG"
-        psfcube["binsz"]    = self._binsz*10.0
-        psfcube["nxpix"]    = self._npix/10
-        psfcube["nypix"]    = self._npix/10
-        psfcube["enumbins"] = self._enumbins
-        psfcube["emin"]     = self["emin"].real()
-        psfcube["emax"]     = self["emax"].real()
-        psfcube["coordsys"] = self._coordsys
-        psfcube["proj"]     = self._proj
+        psfcube['incube']   = 'NONE'
+        psfcube['usepnt']   = True
+        psfcube['ebinalg']  = 'LOG'
+        psfcube['binsz']    = self._binsz*10.0
+        psfcube['nxpix']    = self._npix/10
+        psfcube['nypix']    = self._npix/10
+        psfcube['enumbins'] = self._enumbins
+        psfcube['emin']     = self['emin'].real()
+        psfcube['emax']     = self['emax'].real()
+        psfcube['coordsys'] = self._coordsys
+        psfcube['proj']     = self._proj
+        psfcube['chatter']  = self._chatter
+        psfcube['debug']    = self._logDebug()
         psfcube.run()
+        # JK: Debug
+        psfcube['outcube'] = 'psfcube.fits'
+        psfcube.save()
  
         # Notify Psf computation
         if self._logTerse():
-            self._log("Computed point spread function cube\n")
+            self._log('Computed point spread function cube\n')
  
         # Optionally get stacked Edisp
         if self._edisp:
             edispcube = ctools.ctedispcube(self._obs)
-            edispcube["incube"]   = "NONE"
-            edispcube["usepnt"]   = True
-            edispcube["ebinalg"]  = "LOG"
-            edispcube["binsz"]    = self._binsz*10.0
-            edispcube["nxpix"]    = self._npix/10
-            edispcube["nypix"]    = self._npix/10
-            edispcube["enumbins"] = self._enumbins
-            edispcube["emin"]     = self["emin"].real()
-            edispcube["emax"]     = self["emax"].real()
-            edispcube["coordsys"] = self._coordsys
-            edispcube["proj"]     = self._proj
-            edispcube.run()   
+            edispcube['incube']   = 'NONE'
+            edispcube['usepnt']   = True
+            edispcube['ebinalg']  = 'LOG'
+            edispcube['binsz']    = self._binsz*10.0
+            edispcube['nxpix']    = self._npix/10
+            edispcube['nypix']    = self._npix/10
+            edispcube['enumbins'] = self._enumbins
+            edispcube['emin']     = self['emin'].real()
+            edispcube['emax']     = self['emax'].real()
+            edispcube['coordsys'] = self._coordsys
+            edispcube['proj']     = self._proj
+            edispcube['chatter']  = self._chatter
+            edispcube['debug']    = self._logDebug()
+            edispcube.run()
  
             # Store result
             self._edispcube = edispcube.edispcube().copy()
  
             # Logging
             if self._logTerse():
-                self._log("Computed energy dispersion cube\n")
+                self._log('Computed energy dispersion cube\n')
  
         # Get stacked background
         bkgcube = ctools.ctbkgcube(self._obs)
-        bkgcube["incube"]   = "NONE"
-        bkgcube["usepnt"]   = True
-        bkgcube["ebinalg"]  = "LOG"
-        bkgcube["binsz"]    = self._binsz
-        bkgcube["nxpix"]    = self._npix
-        bkgcube["nypix"]    = self._npix
-        bkgcube["enumbins"] = self._enumbins
-        bkgcube["emin"]     = self["emin"].real()
-        bkgcube["emax"]     = self["emax"].real()
-        bkgcube["coordsys"] = self._coordsys
-        bkgcube["proj"]     = self._proj
+        bkgcube['incube']   = 'NONE'
+        bkgcube['usepnt']   = True
+        bkgcube['ebinalg']  = 'LOG'
+        bkgcube['binsz']    = self._binsz
+        bkgcube['nxpix']    = self._npix
+        bkgcube['nypix']    = self._npix
+        bkgcube['enumbins'] = self._enumbins
+        bkgcube['emin']     = self['emin'].real()
+        bkgcube['emax']     = self['emax'].real()
+        bkgcube['coordsys'] = self._coordsys
+        bkgcube['proj']     = self._proj
+        bkgcube['chatter']  = self._chatter
+        bkgcube['debug']    = self._logDebug()
         bkgcube.run()
+        # JK: Debug
+        bkgcube['outcube']  = 'bkgcube.fits'
+        bkgcube['outmodel'] = 'cube_model.xml'
+        bkgcube.save()
  
         # Notify background cube computation
         if self._logTerse():
-            self._log("Computed background cube\n")
+            self._log('Computed background cube\n')
  
         # Store results
         self._exposure    = expcube.expcube().copy()
@@ -256,19 +270,63 @@ class cspull(ctools.cscript):
         # Return
         return
  
-    def _trial(self, seed):
+    def _obs_string(self, obs):
         """
-        Compute the pull for a single trial.
+        Generate summary string for observation
+
+        Parameters
+        ----------
+        obs : `~gammalib.GCTAObservation`
+            Observation
  
-        Args:
-            seed: Random number generator seed.
+        Returns
+        -------
+        text : str
+            Summary string
+        """
+        # Extract information from observation
+        emin   = obs.events().ebounds().emin().TeV()
+        emax   = obs.events().ebounds().emax().TeV()
+        events = obs.events().number()
+        binned = (obs.events().classname() == 'GCTAEventCube')
+        if binned:
+            mode = 'binned'
+        else:
+            mode = 'unbinned'
+
+        # Compose summary string
+        if events > 0:
+            text = '%d events, %.3f-%.3f TeV, %s' % (events, emin, emax, mode)
+        else:
+            text = '%.3f-%.3f TeV, %s' % (emin, emax, mode)
+
+        # Return summary string
+        return text
+
+    def _trial(self, seed, stacked):
+        """
+        Compute the pull for a single trial
+
+        Parameters
+        ----------
+        seed : int
+            Random number generator seed
+        stacked : bool
+            Use stacked analysis
+
+        Returns
+        -------
+        result : dict
+            Dictionary of results
         """
         # Write header
         if self._logNormal():
-            self._log.header2("Trial "+str(seed+1))
+            self._log.header2('Trial '+str(seed-self._seed+1))
  
         # Simulate events
         obs = obsutils.sim(self._obs,
+                           emin=self['emin'].real(),
+                           emax=self['emax'].real(),
                            nbins=self._enumbins,
                            seed=seed,
                            binsz=self._binsz,
@@ -281,7 +339,7 @@ class cspull(ctools.cscript):
                            chatter=self._chatter)
  
         # Set response for a stacked observation
-        if self._obs.size() > 1 and self._enumbins > 0:
+        if stacked:
             if self._edisp:
                 obs[0].response(self._exposure,  self._psfcube,
                                 self._edispcube, self._bckcube)
@@ -289,6 +347,8 @@ class cspull(ctools.cscript):
                 obs[0].response(self._exposure, self._psfcube,
                                 self._bckcube)
             obs.models(self._stackmodels)
+        # JK: Test
+        obs[0].save('stacked_cube.fits', True)
  
         # Determine number of events in simulation
         nevents = 0.0
@@ -297,13 +357,21 @@ class cspull(ctools.cscript):
  
         # Write simulation results
         if self._logNormal():
-            self._log.header3("Simulation")
-            self._log.parformat("Number of simulated events")
+            self._log.header3('Simulation')
+            for run in self._obs:
+                self._log.parformat('Input observation %s' % (run.id()))
+                self._log(self._obs_string(run))
+                self._log('\n')
+            for run in obs:
+                self._log.parformat('Output observation %s' % (run.id()))
+                self._log(self._obs_string(run))
+                self._log('\n')
+            self._log.parformat('Number of simulated events')
             self._log(nevents)
-            self._log("\n")
+            self._log('\n')
  
         # Fit model
-        if self["profile"].boolean():
+        if self['profile'].boolean():
             models = self._obs.models()
             for i in range(models.size()):
                 model_name = models[i].name()
@@ -324,17 +392,17 @@ class cspull(ctools.cscript):
  
         # Write result header
         if self._logNormal():
-            self._log.header3("Pulls")
+            self._log.header3('Pulls')
  
         # Gather results
         colnames = []
         values   = {}
-        colnames.append("LogL")
-        colnames.append("Sim_Events")
-        colnames.append("Npred_Events")
-        values["LogL"]         = logL
-        values["Sim_Events"]   = nevents
-        values["Npred_Events"] = npred
+        colnames.append('LogL')
+        colnames.append('Sim_Events')
+        colnames.append('Npred_Events')
+        values['LogL']         = logL
+        values['Sim_Events']   = nevents
+        values['Npred_Events'] = npred
         for i in range(models.size()):
             model      = models[i]
             model_name = model.name()
@@ -343,12 +411,12 @@ class cspull(ctools.cscript):
                 if par.is_free():
  
                     # Set parameter name
-                    name = model_name+"_"+par.name()
+                    name = model_name+'_'+par.name()
  
                     # Append parameter, Pull_parameter and Unc_parameter
                     colnames.append(name)
-                    colnames.append("Pull_"+name)
-                    colnames.append("Unc_"+name)
+                    colnames.append('Pull_'+name)
+                    colnames.append('Unc_'+name)
  
                     # Compute pull
                     fitted_value = par.value()
@@ -361,18 +429,18 @@ class cspull(ctools.cscript):
  
                     # Store results
                     values[name] = fitted_value
-                    values["Pull_"+name] = pull
-                    values["Unc_"+name] = error
+                    values['Pull_'+name] = pull
+                    values['Unc_'+name]  = error
  
                     # Write result
                     if self._logNormal():
                         self._log.parformat(name)
                         self._log(pull)
-                        self._log(" (")
+                        self._log(' (')
                         self._log(fitted_value)
-                        self._log(" +/- ")
+                        self._log(' +/- ')
                         self._log(error)
-                        self._log(")\n")
+                        self._log(')\n')
  
         # Bundle together results
         result = {'colnames': colnames, 'values': values}
@@ -384,7 +452,7 @@ class cspull(ctools.cscript):
     # Public methods
     def run(self):
         """
-        Run the script.
+        Run the script
         """
         # Switch screen logging on in debug mode
         if self._logDebug():
@@ -395,25 +463,29 @@ class cspull(ctools.cscript):
  
         # Write observation into logger
         if self._logTerse():
-            self._log("\n")
-            self._log.header1("Observation(s)")
+            self._log('\n')
+            self._log.header1(gammalib.number('Observation',len(self._obs)))
             self._log(str(self._obs))
-            self._log("\n")
+            self._log('\n')
+
+        # Signal if stacked analysis is requested
+        #stacked = self._obs.size() > 1 and self._enumbins > 0
+        stacked = self._enumbins > 0
  
         # If several observations and binned: prepare stacked irfs
-        if self._obs.size() > 1 and self._enumbins > 0:
+        if stacked:
             self._set_stacked_irf()
  
         # Write header
         if self._logTerse():
-            self._log("\n")
-            self._log.header1("Generate pull distribution")
+            self._log('\n')
+            self._log.header1('Generate pull distribution')
  
         # Loop over trials
         for seed in range(self._ntrials):
  
             # Make a trial and add initial seed
-            result = self._trial(seed + self._seed)
+            result = self._trial(seed + self._seed, stacked)
  
             # Write out result immediately
             if seed == 0:
@@ -434,7 +506,7 @@ class cspull(ctools.cscript):
  
     def execute(self):
         """
-        Execute the script.
+        Execute the script
         """
         # Open logfile
         self.logFileOpen()
@@ -447,10 +519,15 @@ class cspull(ctools.cscript):
  
     def models(self, models):
         """
-        Set model.
+        Set model
+
+        Parameters
+        ----------
+        models : `~gammalib.GModels`
+            Set model container
         """
-        # Copy models
-        self.obs.models(models.copy())
+        # Set model container
+        self._obs.models(models)
  
         # Return
         return
@@ -24,46 +24,58 @@ import ctools
 # ===================== #
 # Simulate observations #
 # ===================== #
-def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, nbins=0,
-        binsz=0.05, npix=200, proj="TAN", coord="GAL"):
+def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0,
+        emin=None, emax=None, nbins=0,
+        binsz=0.05, npix=200, proj='TAN', coord='GAL'):
     """
-    Simulate events for all observations in the container.
-
-    Args:
-        obs: Observation container without events
-
-    Kwargs:
-        log:     Create log file(s)
-        debug:   Create console dump?
-        chatter: Chatter level
-        edisp:   Apply energy dispersion?
-        seed:    Seed value for simulations (default: 0)
-        nbins:   Number of energy bins (default: 0=unbinned)
-        binsz:   Pixel size for binned simulation (deg/pixel)
-        npix:    Number of pixels in X and Y for binned simulation
-        proj:    Projection for binned simulation
-        coord:   Coordinate system for binned simulation
-
-    Returns:
-        Observation container, filled with simulated events.
+    Simulate events for all observations in the container
+
+    Parameters
+    ----------
+    obs : `~gammalib.GObservations`
+        Observation container without events
+    log : bool, optional
+        Create log file(s)
+    debug : bool, optional
+        Create console dump?
+    chatter : int, optional
+        Chatter level
+    edisp : bool, optional
+        Apply energy dispersion?
+    seed : int, integer
+        Seed value for simulations
+    emin : float, optional
+        Minimum energy of counts cube for binned (TeV)
+    emax : float, optional
+        Maximum energy of counts cube for binned (TeV)
+    nbins : int, optional
+        Number of energy bins (0=unbinned)
+    binsz : float, optional
+        Pixel size for binned simulation (deg/pixel)
+    npix : int, optional
+        Number of pixels in X and Y for binned simulation
+    proj : str, optional
+        Projection for binned simulation
+    coord : str, optional
+        Coordinate system for binned simulation
+
+    Returns
+    -------
+    obs : `~gammalib.GObservations`
+        Observation container filled with simulated events
     """
  
     # Allocate ctobssim application and set parameters
     sim = ctools.ctobssim(obs)
-    sim["seed"]  = seed
-    sim["edisp"] = edisp
+    sim['seed']    = seed
+    sim['edisp']   = edisp
+    sim['chatter'] = chatter
+    sim['debug']   = debug
  
     # Optionally open the log file
     if log:
         sim.logFileOpen()
  
-    # Optionally switch-on debugging model
-    if debug:
-        sim["debug"] = True
-
-    # Set chatter level
-    sim["chatter"] = chatter
-
     # Run ctobssim application. This will loop over all observations in the
     # container and simulation the events for each observation. Note that
     # events are not added together, they still apply to each observation
@@ -73,45 +85,42 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, nbins=0,
     # Binned option?
     if nbins > 0:
  
-        # Determine common energy boundaries for observations
-        emin = None
-        emax = None
-        for run in sim.obs():
-            run_emin = run.events().ebounds().emin().TeV()
-            run_emax = run.events().ebounds().emax().TeV()
-            if emin == None:
-                emin = run_emin
-            elif run_emin > emin:
-                emin = run_emin
-            if emax == None:
-                emax = run_emax
-            elif run_emax > emax:
-                emax = run_emax
+        # If energy boundaries are not given then determine common energy
+        # boundaries for all observations
+        if emin == None or emax == None:
+            emin = None
+            emax = None
+            for run in sim.obs():
+                run_emin = run.events().ebounds().emin().TeV()
+                run_emax = run.events().ebounds().emax().TeV()
+                if emin == None:
+                    emin = run_emin
+                elif run_emin > emin:
+                    emin = run_emin
+                if emax == None:
+                    emax = run_emax
+                elif run_emax > emax:
+                    emax = run_emax
  
         # Allocate ctbin application and set parameters
         bin = ctools.ctbin(sim.obs())
-        bin["ebinalg"]  = "LOG"
-        bin["emin"]     = emin
-        bin["emax"]     = emax
-        bin["enumbins"] = nbins
-        bin["usepnt"]   = True # Use pointing for map centre
-        bin["nxpix"]    = npix
-        bin["nypix"]    = npix
-        bin["binsz"]    = binsz
-        bin["coordsys"] = coord
-        bin["proj"]     = proj
+        bin['ebinalg']  = 'LOG'
+        bin['emin']     = emin
+        bin['emax']     = emax
+        bin['enumbins'] = nbins
+        bin['usepnt']   = True # Use pointing for map centre
+        bin['nxpix']    = npix
+        bin['nypix']    = npix
+        bin['binsz']    = binsz
+        bin['coordsys'] = coord
+        bin['proj']     = proj
+        bin['chatter']  = chatter
+        bin['debug']    = debug
  
         # Optionally open the log file
         if log:
             bin.logFileOpen()
  
-        # Optionally switch-on debugging model
-        if debug:
-            bin["debug"] = True
-
-        # Set chatter level
-        bin["chatter"] = chatter
-
         # Run ctbin application. This will loop over all observations in
         # the container and bin the events in counts maps
         bin.run()
@@ -220,12 +220,7 @@ void ctexpcube::run(void)
     // Write observation(s) into logger
     if (logTerse()) {
         log << std::endl;
-        if (m_obs.size() > 1) {
-            log.header1("Observations");
-        }
-        else {
-            log.header1("Observation");
-        }
+        log.header1(gammalib::number("Observation",m_obs.size()));
         log << m_obs << std::endl;
     }
  
@@ -230,13 +230,8 @@ void ctlike::run(void)
     // Write observation(s) into logger
     if (logTerse()) {
         log << std::endl;
-        if (m_obs.size() > 1) {
-            log.header1("Observations");
-        }
-        else {
-            log.header1("Observation");
-        }
-        log << m_obs << std::endl;
+        log.header1(gammalib::number("Observation",m_obs.size()));
+        log << m_obs.print(m_chatter) << std::endl;
     }
  
     // Optimize model parameters using LM optimizer
@@ -311,7 +306,7 @@ void ctlike::run(void)
     if (logTerse()) {
         log << std::endl;
         log.header1("Maximum likelihood optimisation results");
-        log << *m_opt << std::endl;
+        log << m_opt->print(m_chatter) << std::endl;
         log << gammalib::parformat("Maximum log likelihood");
         log << gammalib::str(m_logL,3) << std::endl;
         log << gammalib::parformat("Observed events  (Nobs)");
@@ -321,7 +316,7 @@ void ctlike::run(void)
         log << " (Nobs - Npred = ";
         log << gammalib::str(num_events-npred);
         log << ")" << std::endl;
-        log << m_obs.models() << std::endl;
+        log << m_obs.models().print(m_chatter) << std::endl;
     }
  
     // Restore energy dispersion flag for all CTA observations
@@ -430,6 +425,7 @@ void ctlike::get_parameters(void)
     m_refit           = (*this)["refit"].boolean();
     m_apply_edisp     = (*this)["edisp"].boolean();
     m_fix_spat_for_ts = (*this)["fix_spat_for_ts"].boolean();
+    m_chatter         = static_cast<GChatter>((*this)["chatter"].integer());
  
     // Optionally read ahead parameters so that they get correctly
     // dumped into the log file
@@ -592,6 +588,7 @@ void ctlike::init_members(void)
     m_opt             = NULL;
     m_apply_edisp     = false;
     m_fix_spat_for_ts = false;
+    m_chatter         = static_cast<GChatter>(2);
  
     // Set logger properties
     log.date(true);
@@ -628,6 +625,7 @@ void ctlike::copy_members(const ctlike&amp; app)
     m_opt             = app.m_opt->clone();
     m_apply_edisp     = app.m_apply_edisp;
     m_fix_spat_for_ts = app.m_fix_spat_for_ts;
+    m_chatter         = app.m_chatter;
  
     // Return
     return;
@@ -77,6 +77,7 @@ protected:
     GFilename     m_outmodel;        //!< Source model output XML file name
     bool          m_apply_edisp;     //!< Apply energy dispersion?
     bool          m_fix_spat_for_ts; //!< Fix spatial parameters for TS computation?
+    GChatter      m_chatter;         //!< Chattiness
  
     // Members
     GObservations m_obs;        //!< Observations