diff --git a/ChangeLog b/ChangeLog index 02b87b6..1f8832d 100644 --- a/ChangeLog +++ b/ChangeLog @@ -26,6 +26,7 @@ Add "NODES" method to csspec to support non-CTA analysis (#1927) Use t-type parameters for all "tmin" and "tmax" parameters (#1864) Set OBS_ID in ctobssim (#1836) + Revise cstsdist parameters (#1567) Add ctbutterfly::butterfly() method (#1512) diff --git a/NEWS b/NEWS index 0e77805..a83fb3f 100644 --- a/NEWS +++ b/NEWS @@ -319,8 +319,14 @@ The script now derives from the csobservation base class. A "statistic" parameter was added to allow setting the optimizer statistic for the Test Statistic distribution generation (#2237). -Use t-type parameters for "tmin" and "tmax" to all specifying times in various -formats (#1864). +The parameter interface was considerably revised, in particular the following +parameters were removed: "ra", "dec", "rad", "coordsys", "proj", "emin", +"emax", "enumbins", "tmin", "tmax", "npix" and "binsz". The "inobs" parameter +is now always queried. The script now reads in an observation or an observation +container and derives the TS distribution for a source based on that observation +container. So far the script supports unbinned, binned and stacked observations. +If an observation definition XML file is specified the script assumes that all +observations are event lists (#1567). cstsmapmerge - Merge Test Statistic maps diff --git a/cscripts/cstsdist.par b/cscripts/cstsdist.par index 8a023fd..3e968a8 100644 --- a/cscripts/cstsdist.par +++ b/cscripts/cstsdist.par @@ -16,7 +16,7 @@ # # Input and output parameters #============================ -inobs, f, h, NONE,,, "Input event list, counts cube, or observation definition XML file" +inobs, f, a, events.fits,,, "Input event list, counts cube, or observation definition XML file" inmodel, f, a, $CTOOLS/share/models/crab.xml,,, "Input model definition XML file" srcname, s, a, Crab,,, "Source name" expcube, f, a, NONE,,, "Input exposure cube file" @@ -33,18 +33,6 @@ outfile, f, a, ts.dat,,, "Output Test Statistic distribution file" # Script parameters #================== ntrials, i, a, 10,,, "Number of trials" -ra, r, a, 83.63,0,360, "RA of pointing (degrees)" -dec, r, a, 22.01,-90,90, "Dec of pointing (degrees)" -rad, r, h, 5.0,0,180, "Radius of ROI (degrees)" -coordsys, s, a, CEL,CEL|GAL,, "Coordinate system (CEL - celestial, GAL - galactic)" -proj, s, a, TAN,AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method" -emin, r, a, 0.1,,, "Lower energy limit (TeV)" -emax, r, a, 100.0,,, "Upper energy limit (TeV)" -enumbins, i, a, 0,,, "Number of energy bins (0=unbinned)" -tmin, t, a, 2020-01-01T00:00:00,,, "Start time (UTC string, JD, MJD or MET in seconds)" -tmax, t, a, 2020-01-01T00:30:00,,, "Stop time (UTC string, JD, MJD or MET in seconds)" -npix, i, a, 200,,, "Number of pixels for binned" -binsz, r, a, 0.05,,, "Pixel size for binned (degrees/pixel)" statistic, s, h, DEFAULT,DEFAULT|CSTAT|WSTAT|CHI2,, "Optimization statistic" # diff --git a/cscripts/cstsdist.py b/cscripts/cstsdist.py index 41351d7..e31af55 100755 --- a/cscripts/cstsdist.py +++ b/cscripts/cstsdist.py @@ -2,7 +2,7 @@ # ========================================================================== # Generates the TS distribution for a particular model # -# Copyright (C) 2011-2017 Juergen Knoedlseder +# Copyright (C) 2011-2018 Juergen Knoedlseder # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -37,7 +37,7 @@ class cstsdist(ctools.csobservation): # Constructor def __init__(self, *argv): """ - Constructor. + Constructor """ # Initialise application by calling the appropriate class constructor self._init_csobservation(self.__class__.__name__, ctools.__version__, argv) @@ -45,6 +45,7 @@ class cstsdist(ctools.csobservation): # Initialise some members self._srcname = '' self._log_clients = False + self._model = None # Return return @@ -65,13 +66,6 @@ class cstsdist(ctools.csobservation): # Get source name self._srcname = self['srcname'].string() - # Query parameters for binned if requested - if self['enumbins'].integer() != 0: - self['npix'].integer() - self['binsz'].real() - self['coordsys'].string() - self['proj'].string() - # Set models if we have none if self.obs().models().size() == 0: self.obs().models(self['inmodel'].filename()) @@ -79,15 +73,73 @@ class cstsdist(ctools.csobservation): # Query parameters self['edisp'].boolean() self['ntrials'].integer() - self['outfile'].filename() + self['debug'].boolean() + # Read ahead output parameters + #if self._read_ahead(): + # self['outfile'].filename() + self['outfile'].filename() + # Write input parameters into logger self._log_parameters(gammalib.TERSE) # Return return + def _sim(self, seed): + """ + Return a simulated observation container + + Parameters + ---------- + seed : int + Random number generator seed + + Returns + ------- + sim : `~gammalib.GObservations` + Simulated observation container + """ + # If observation is a counts cube then simulate events from the counts + # cube model ... + if self.obs().size() == 1 and self.obs()[0].eventtype() == 'CountsCube': + + # If no counts cube model exists then compute it now + if self._model == None: + model = ctools.ctmodel(self.obs()) + model['debug'] = self['debug'].boolean() + model['chatter'] = self['chatter'].integer() + model.run() + self._model = model.cube().copy() # Save copy for persistence + + # Allocate random number generator + ran = gammalib.GRan() + + # Get copy of model map + counts = self._model.counts().copy() + + # Randomize counts + for i in range(counts.npix()): + counts[i] = ran.poisson(counts[i]) + + # Copy observations + sim = self.obs().copy() + + # Set counts map + sim[0].events().counts(counts) + + # ... otherwise simuate events from the observation container (works + # only for event lists + else: + sim = obsutils.sim(self.obs(), + seed = seed, + log = self._log_clients, + debug = self['debug'].boolean()) + + # Return simulated observation + return sim + def _trial(self, seed): """ Create the TS for a single trial @@ -105,29 +157,8 @@ class cstsdist(ctools.csobservation): # Write header self._log_header2(gammalib.EXPLICIT, 'Trial %d' % (seed+1)) - # Set default binned parameters - coordsys = 'CEL' - proj = 'TAN' - npix = 0 - binsz = 0.0 - - # If binned analysis is requested the read the binned parameters - if self['enumbins'].integer() != 0: - npix = self['npix'].integer() - binsz = self['binsz'].real() - coordsys = self['coordsys'].string() - proj = self['proj'].string() - # Simulate events - sim = obsutils.sim(self.obs(), - nbins = self['enumbins'].integer(), - seed = seed, - proj = proj, - coord = coordsys, - binsz = binsz, - npix = npix, - log = self._log_clients, - debug = self['debug'].boolean()) + sim = self._sim(seed) # Determine number of events in simulation nevents = 0.0 @@ -160,7 +191,7 @@ class cstsdist(ctools.csobservation): for model in models: self._log_value(gammalib.EXPLICIT, 'Model', model.name()) for par in model: - self._log_string(gammalib.EXPLICIT, str(par)+'\n') + self._log_string(gammalib.EXPLICIT, str(par)) elif self._logNormal(): name = 'Trial %d' % seed value = 'TS=%.3f Prefactor=%e +/- %e' % \ diff --git a/doc/source/users/reference_manual/cstsdist.rst b/doc/source/users/reference_manual/cstsdist.rst index 6e14e55..344e62f 100644 --- a/doc/source/users/reference_manual/cstsdist.rst +++ b/doc/source/users/reference_manual/cstsdist.rst @@ -15,9 +15,6 @@ sets. The Test Statistics is defined as twice the log-likelihood difference that is obtained when fitting simulated data with and without a given source model component. -This script either performs an unbinned or a binned simulation, stacked -analysis is not yet supported. - cstsdist will create an ASCII file in comma-separated value (CSV) format, containing one row per TS computation. The first row is a header row providing the column names. The following rows give the TS value, the log-likelihood @@ -26,13 +23,18 @@ fitted events, as well as the values and errors for all fitted parameters. From the output file, TS distribution plots can be generated using for example the ``show_ts_distribution.py`` script in the examples folder. The -script require matplotlib for plotting. +script requires matplotlib for plotting. + +.. warning:: + This script does not work for On/Off observations. If an observation + definition XML file is specified the script assumes that all observations + are event lists. General parameters ------------------ -``(inobs = NONE) [file]`` +``inobs [file]`` Input event list, counts cube or observation definition XML file. ``inmodel [file]`` @@ -72,42 +74,6 @@ General parameters ``ntrials [integer]`` Number of Monte Carlo samples. -``ra [real]`` - Right Ascension of CTA pointing (J2000, in degrees). - -``dec [real]`` - Declination of CTA pointing (J2000, in degrees). - -``(rad = 5.0) [real]`` - ROI radius (in degrees). - -``coordsys [string]`` - Coordinate system (CEL - celestial, GAL - galactic). - -``proj [string]`` - Projection method. - -``emin [real]`` - Lower energy limit (in TeV). - -``emax [real]`` - Upper energy limit (in TeV). - -``enumbins [integer]`` - Number of energy bins (0=unbinned). - -``tmin [time]`` - Start time (UTC string, JD, MJD or MET in seconds). - -``tmax [time]`` - Stop time (UTC string, JD, MJD or MET in seconds). - -``npix [integer]`` - Number of pixels for binned analysis. - -``binsz [real]`` - Pixel size for binned analysis. - ``(statistic = DEFAULT) [string]`` Optimization statistic. ``DEFAULT`` uses the default statistic for all observations, which is ``CSTAT`` or the statistic specified in the diff --git a/src/Makefile.am b/src/Makefile.am index 7c54235..1dbf1fe 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -100,7 +100,7 @@ help-install: echo "$(INSTALL_DATA) $$file $(DESTDIR)$(datadir)/help/$$helpname" ; \ $(INSTALL_DATA) "$$file" $(DESTDIR)$(datadir)/help/"$$helpname" ; \ sed -i.bak -e 's/``//g' -e 's/`//g' -e 's/:ref://g' -e 's/:doc://g' \ - -e 's/.. note::/Note/g' -e '1d' \ + -e 's/.. note::/Note/g' -e 's/.. warning::/Warning/g' -e '1d' \ $(DESTDIR)$(datadir)/help/"$$helpname" ; \ rm -rf $(DESTDIR)$(datadir)/help/"$$helpname".bak ; \ done diff --git a/test/test_cstsdist.py b/test/test_cstsdist.py index ca3524d..f81a036 100755 --- a/test/test_cstsdist.py +++ b/test/test_cstsdist.py @@ -2,7 +2,7 @@ # ========================================================================== # This scripts performs unit tests for the cstsdist script. # -# Copyright (C) 2016-2017 Juergen Knoedlseder +# Copyright (C) 2016-2018 Juergen Knoedlseder # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -42,6 +42,10 @@ class Test(test): # Call base class constructor test.__init__(self) + # Set members + self._stacked_inobs = self._datadir + '/obs_stacked.xml' + self._stacked_model = self._datadir + '/crab_bkgcube.xml' + # Return return @@ -55,7 +59,12 @@ class Test(test): # Append tests self.append(self._test_cmd, 'Test cstsdist on command line') - self.append(self._test_python, 'Test cstsdist from Python') + self.append(self._test_python_unbinned, + 'Test cstsdist from Python in unbinned mode') + self.append(self._test_python_binned, + 'Test cstsdist from Python in binned mode') + self.append(self._test_python_stacked, + 'Test cstsdist from Python in stacked mode') # Return return @@ -69,12 +78,10 @@ class Test(test): cstsdist = self._script('cstsdist') # Setup cstsdist command - cmd = cstsdist+' inmodel="'+self._model+'"'+ \ - ' srcname="Crab"'+ \ + cmd = cstsdist+' inobs="'+self._events+'"'+ \ + ' inmodel="'+self._model+'" srcname="Crab"'+ \ ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \ - ' ntrials=1 ra=83.63 dec=22.01 emin=0.1 emax=100.0'+ \ - ' enumbins=0 tmin=0.0 tmax=1800.0 rad=5.0 npix=200 npix=200'+ \ - ' binsz=0.05'+ \ + ' ntrials=2 '+ \ ' outfile="cstsdist_cmd1.dat"'+ \ ' logfile="cstsdist_cmd1.log" chatter=2' @@ -90,14 +97,12 @@ class Test(test): self._check_result_file('cstsdist_cmd1.dat') # Setup cstsdist command - cmd = cstsdist+' inmodel="model_that_does_not_exist.xml"'+ \ - ' srcname="Crab"'+ \ + cmd = cstsdist+' inobs="event_file_that_does_not_exist.xml"'+ \ + ' inmodel="'+self._model+'" srcname="Crab"'+ \ ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \ - ' ntrials=1 ra=83.63 dec=22.01 emin=0.1 emax=100.0'+ \ - ' enumbins=0 tmin=0.0 tmax=1800.0 rad=5.0 npix=200 npix=200'+ \ - ' binsz=0.05'+ \ + ' ntrials=1 '+ \ ' outfile="cstsdist_cmd2.dat"'+ \ - ' logfile="cstsdist_cmd2.log" debug=yes chatter=2' + ' logfile="cstsdist_cmd2.log" chatter=2' # Check if execution failed self.test_assert(self._execute(cmd) != 0, @@ -110,25 +115,18 @@ class Test(test): return # Test cstsdist from Python - def _test_python(self): + def _test_python_unbinned(self): """ Test cstsdist from Python """ - # Set-up cstsdist + # Set-up cstsdist for event list tsdist = cscripts.cstsdist() + tsdist['inobs'] = self._events tsdist['inmodel'] = self._model tsdist['srcname'] = 'Crab' tsdist['caldb'] = self._caldb tsdist['irf'] = self._irf - tsdist['ntrials'] = 1 - tsdist['ra'] = 83.63 - tsdist['dec'] = 22.01 - tsdist['emin'] = 0.1 - tsdist['emax'] = 100.0 - tsdist['enumbins'] = 0 - tsdist['tmin'] = 0.0 - tsdist['tmax'] = 1800.0 - tsdist['rad'] = 5.0 + tsdist['ntrials'] = 2 tsdist['outfile'] = 'cstsdist_py1.dat' tsdist['logfile'] = 'cstsdist_py1.log' tsdist['chatter'] = 3 @@ -144,8 +142,61 @@ class Test(test): # Return return + def _test_python_binned(self): + """ + Test cstsdist from Python + """ + # Set-up cstsdist for counts cube + tsdist = cscripts.cstsdist() + tsdist['inobs'] = self._cntcube + tsdist['inmodel'] = self._model + tsdist['srcname'] = 'Crab' + tsdist['expcube'] = 'NONE' + tsdist['psfcube'] = 'NONE' + tsdist['bkgcube'] = 'NONE' + tsdist['caldb'] = self._caldb + tsdist['irf'] = self._irf + tsdist['ntrials'] = 2 + tsdist['outfile'] = 'cstsdist_py2.dat' + tsdist['logfile'] = 'cstsdist_py2.log' + tsdist['chatter'] = 4 + + # Run cstsdist script + tsdist.logFileOpen() # Make sure we get a log file + tsdist.execute() + + # Check pull distribution file + self._check_result_file('cstsdist_py2.dat') + + # Return + return + + def _test_python_stacked(self): + """ + Test cstsdist from Python + """ + # Set-up cstsdist for stacked cube + tsdist = cscripts.cstsdist() + tsdist['inobs'] = self._stacked_inobs + tsdist['inmodel'] = self._stacked_model + tsdist['srcname'] = 'Crab' + tsdist['ntrials'] = 2 + tsdist['outfile'] = 'cstsdist_py3.dat' + tsdist['logfile'] = 'cstsdist_py3.log' + tsdist['chatter'] = 4 + + # Run cstsdist script + tsdist.logFileOpen() # Make sure we get a log file + tsdist.execute() + + # Check pull distribution file + self._check_result_file('cstsdist_py3.dat') + + # Return + return + # Check result file - def _check_result_file(self, filename, ncols=9): + def _check_result_file(self, filename, nrows=3, ncols=9): """ Check result file @@ -153,6 +204,8 @@ class Test(test): ---------- filename : str Name of result file + nrows : int, optional + Required number of rows ncols : int, optional Required number of columns """ @@ -160,7 +213,7 @@ class Test(test): results = gammalib.GCsv(filename, ',') # Check dimensions - self.test_value(results.nrows(), 2, 'Check rows in TS file') + self.test_value(results.nrows(), nrows, 'Check rows in TS file') self.test_value(results.ncols(), ncols, 'Check columns in TS file') # Return