diff --git a/cscripts/Makefile.am b/cscripts/Makefile.am index 39fb8cb..a4ebb5f 100644 --- a/cscripts/Makefile.am +++ b/cscripts/Makefile.am @@ -164,6 +164,7 @@ testdata = $(top_srcdir)/test/data/crab.xml \ $(top_srcdir)/test/data/crab_psfcube.fits.gz \ $(top_srcdir)/test/data/crab_edispcube.fits.gz \ $(top_srcdir)/test/data/crab_bkgcube.fits.gz \ + $(top_srcdir)/test/data/crab_onoff.xml \ $(top_srcdir)/test/data/crab_offaxis1.fits.gz \ $(top_srcdir)/test/data/crab_offaxis2.fits.gz \ $(top_srcdir)/test/data/crab_exclusion.fits.gz \ diff --git a/cscripts/cslightcrv.par b/cscripts/cslightcrv.par index 3cf4918..ec4f7bf 100644 --- a/cscripts/cslightcrv.par +++ b/cscripts/cslightcrv.par @@ -16,12 +16,13 @@ # # Input and output parameters #============================ -inobs, f, a, events.fits,,, "Input event list or observation definition XML file" -inmodel, f, a, $CTOOLS/share/models/crab.xml,,, "Input model definition XML file" -srcname, s, a, Crab,,, "Source name" -caldb, s, a, prod2,,, "Calibration database" -irf, s, a, South_0.5h,,, "Instrument response function" -edisp, b, h, no,,, "Apply energy dispersion?" +inobs, f, a, events.fits,,, "Input event list or observation definition XML file" +inmodel, f, a, $CTOOLS/share/models/crab.xml,,, "Input model definition XML file" +srcname, s, a, Crab,,, "Source name" +caldb, s, a, prod2,,, "Calibration database" +irf, s, a, South_0.5h,,, "Instrument response function" +inexclusion, f, h, NONE,,, "Sky map of exclusion regions for On/Off analysis" +edisp, b, h, no,,, "Apply energy dispersion in Cube analysis?" outfile, f, a, lightcurve.fits,,, "Output light curve file" # @@ -32,16 +33,25 @@ tmin, t, a, MJD 51544.5,,,"Lightcurve start time (UTC string, JD, MJD or ME tmax, t, a, MJD 51544.6,,,"Lightcurve stop time (UTC string, JD, MJD or MET in seconds)" tbins, i, a, 5,1,10000, "Number of time bins" tbinfile, f, a, NONE,,, "File defining the time binning" +method, s, a, 3D, 3D|ONOFF,, "Analysis method" emin, r, a, 0.1,,, "Lower energy limit of events (TeV)" emax, r, a, 100.0,,, "Upper energy limit of events (TeV)" -enumbins, i, a, 0,0,100, "Number of energy bins for binned (0=unbinned)" +enumbins, i, a, 0,0,100, "Number of energy bins for binned (0=unbinned for 3D analysis only)" coordsys, s, a, CEL,CEL|GAL,, "Coordinate System" proj, s, a, CAR,AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method" -xref, r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)" -yref, r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)" +xref, r, a, 83.63,0,360, "First coordinate of image/source region center in degrees (RA or galactic l)" +yref, r, a, 22.01,-90,90, "Second coordinate of image/source region center in degrees (DEC or galactic b)" nxpix, i, a, 200,,, "Size of the X axis in pixels" nypix, i, a, 200,,, "Size of the Y axis in pixels" binsz, r, a, 0.02,,, "Pixel size (deg/pixel)" +srcshape, s, h, CIRCLE,CIRCLE,, "Shape of the source region" +rad, r, a, 0.2,0,180, "Radius of source region circle (deg)" +bkgmethod, s, h, REFLECTED,REFLECTED,, "Method for background estimation" +bkgregmin, i, h, 2,1,1000, "Minimum number of background regions" +maxoffset, r, h, 4.0,,, "Maximum offset of source from camera center to accept observation" +etruemin, r, h, 0.01,,, "Minimum true energy (TeV)" +etruemax, r, h, 300.0,,, "Maximum true energy (TeV)" +etruebins, i, h, 30,,, "Number of bins per decade for true energy bins" statistic, s, h, DEFAULT,DEFAULT|CSTAT|WSTAT|CHI2,, "Optimization statistic" calc_ts, b, h, yes,,, "Compute TS value for each bin?" calc_ulim, b, h, yes,,, "Compute upper limit for each bin?" diff --git a/cscripts/cslightcrv.py b/cscripts/cslightcrv.py index 880454b..bf99296 100755 --- a/cscripts/cslightcrv.py +++ b/cscripts/cslightcrv.py @@ -74,6 +74,7 @@ class cslightcrv(ctools.csobservation): # Initialise some members self._srcname = '' self._tbins = gammalib.GGti() + self._onoff = False self._stacked = False self._fits = gammalib.GFits() @@ -103,8 +104,13 @@ class cslightcrv(ctools.csobservation): # Get time boundaries self._tbins = self._create_tbounds() + # Set On/Off analysis flag and query relevant user parameters + self._onoff = self._is_onoff() + + # If cube analysis is selected # Set stacked analysis flag and query relevant user parameters - self._stacked = self._is_stacked() + if not self._onoff: + self._stacked = self._is_stacked() # Query the hidden parameters, just in case self['edisp'].boolean() @@ -357,40 +363,6 @@ class cslightcrv(ctools.csobservation): # Return upper limit tuple return ul_diff, ul_flux, ul_eflux - def _bin_observation(self, obs): - """ - Bin an observation if a binned analysis was requested - - Parameters - ---------- - obs : `~gammalib.GObservations` - Observation container - - Returns - ------- - obs : `~gammalib.GObservations` - Observation container where the first observation is a binned observation - """ - # Get new observation container with binned observation - new_obs = obsutils.get_stacked_obs(self, obs) - - # Extract models - models = new_obs.models() - - # Fix background models if required - if self['fix_bkg'].boolean(): - for model in models: - if model.classname() != 'GModelSky': - for par in model: - par.fix() - - # Put back models - new_obs.models(models) - - # Return new oberservation container - return new_obs - - # Public methods def run(self): """ @@ -483,31 +455,116 @@ class cslightcrv(ctools.csobservation): # Retrieve observation obs = select.obs() - # If a stacked analysis is requested then bin the events - # and compute the stacked response functions and setup - # an observation container with a single stacked observation. - if self._stacked: - obs = self._bin_observation(obs) + # Deal with stacked and On/Off Observations + if self._stacked or self._onoff: + + # If a stacked analysis is requested bin the events + # and compute the stacked response functions and setup + # an observation container with a single stacked observation. + if self._stacked: + new_obs = obsutils.get_stacked_obs(self, obs) + + # ... otherwise if On/Off analysis is requested generate + # the On/Off observations and response + elif self._onoff: + new_obs = obsutils.get_onoff_obs(self, obs) + + # Extract models + models = new_obs.models() + + # Fix background models if required + if self['fix_bkg'].boolean(): + for model in models: + if model.classname() != 'GModelSky': + for par in model: + par.fix() + + # Put back models + new_obs.models(models) + + # Continue with new oberservation container + obs = new_obs # Header self._log_header3(gammalib.EXPLICIT, 'Fitting the data') # Do maximum likelihood model fitting - like = ctools.ctlike(obs) - like['edisp'] = self['edisp'].boolean() - like.run() + if obs.size() > 0: + like = ctools.ctlike(obs) + like['edisp'] = self['edisp'].boolean() + like.run() + + # Skip bin if no event was present + if like.obs().logL() == 0.0: - # Skip bin if no event was present - if like.obs().logL() == 0.0: + # Signal skipping of bin + self._log_value(gammalib.TERSE, 'Warning', + 'No event in this time bin, skip bin.') - # Signal skipping of bin + # Set all results to 0 + for par in pars: + result['values'][par] = 0.0 + result['values']['e_'+par] = 0.0 + + # Append result + results.append(result) + + # Continue with next time bin + continue + + # Retrieve model fitting results for source of interest + source = like.obs().models()[self._srcname] + + # Extract parameter values + for par in pars: + result['values'][par] = source[par].value() + result['values']['e_'+par] = source[par].error() + + # Calculate upper limit (-1 if not computed) + #ul_diff, ul_flux, ul_eflux = self._compute_ulimit(like.obs()) + ul_diff, ul_flux, ul_eflux = self._compute_ulimit(obs) + if ul_diff > 0.0: + result['ul_diff'] = ul_diff + result['ul_flux'] = ul_flux + result['ul_eflux'] = ul_eflux + + # Extract Test Statistic value + if self['calc_ts'].boolean(): + result['ts'] = source.ts() + + # Append result to list of dictionaries + results.append(result) + + # Log results for this time bin + self._log.header3('Results') + pars = self._get_free_par_names() + for par in pars: + value = source[par].value() + error = source[par].error() + unit = source[par].unit() + self._log_value(gammalib.NORMAL, par, + str(value)+' +/- '+str(error)+' '+unit) + if ul_diff > 0.0: + self._log_value(gammalib.NORMAL, 'Upper flux limit', + str(result['ul_diff'])+' ph/cm2/s/MeV') + self._log_value(gammalib.NORMAL, 'Upper flux limit', + str(result['ul_flux'])+' ph/cm2/s') + self._log_value(gammalib.NORMAL, 'Upper flux limit', + str(result['ul_eflux'])+' erg/cm2/s') + if self['calc_ts'].boolean(): + self._log_value(gammalib.NORMAL, 'Test Statistic', result['ts']) + + # Otherwise, if observations size is 0, signal bin is skipped and + # fill results table with zeros + else: self._log_value(gammalib.TERSE, 'Warning', - 'No event in this time bin, skip bin.') + 'No observations available in this time bin', + 'Skip bin.') # Set all results to 0 for par in pars: - result['values'][par] = 0.0 - result['values']['e_'+par] = 0.0 + result['values'][par] = 0.0 + result['values']['e_' + par] = 0.0 # Append result results.append(result) @@ -515,48 +572,6 @@ class cslightcrv(ctools.csobservation): # Continue with next time bin continue - # Retrieve model fitting results for source of interest - source = like.obs().models()[self._srcname] - - # Extract parameter values - for par in pars: - result['values'][par] = source[par].value() - result['values']['e_'+par] = source[par].error() - - # Calculate upper limit (-1 if not computed) - #ul_diff, ul_flux, ul_eflux = self._compute_ulimit(like.obs()) - ul_diff, ul_flux, ul_eflux = self._compute_ulimit(obs) - if ul_diff > 0.0: - result['ul_diff'] = ul_diff - result['ul_flux'] = ul_flux - result['ul_eflux'] = ul_eflux - - # Extract Test Statistic value - if self['calc_ts'].boolean(): - result['ts'] = source.ts() - - # Append result to list of dictionaries - results.append(result) - - # Log results for this time bin - self._log.header3('Results') - pars = self._get_free_par_names() - for par in pars: - value = source[par].value() - error = source[par].error() - unit = source[par].unit() - self._log_value(gammalib.NORMAL, par, - str(value)+' +/- '+str(error)+' '+unit) - if ul_diff > 0.0: - self._log_value(gammalib.NORMAL, 'Upper flux limit', - str(result['ul_diff'])+' ph/cm2/s/MeV') - self._log_value(gammalib.NORMAL, 'Upper flux limit', - str(result['ul_flux'])+' ph/cm2/s') - self._log_value(gammalib.NORMAL, 'Upper flux limit', - str(result['ul_eflux'])+' erg/cm2/s') - if self['calc_ts'].boolean(): - self._log_value(gammalib.NORMAL, 'Test Statistic', result['ts']) - # Create FITS table from results table = self._create_fits_table(results) diff --git a/cscripts/csphasecrv.par b/cscripts/csphasecrv.par index aac533a..49720a3 100644 --- a/cscripts/csphasecrv.par +++ b/cscripts/csphasecrv.par @@ -21,7 +21,8 @@ inmodel, f, a, $CTOOLS/share/models/crab.xml,,, "Input model definition XML fi srcname, s, a, Crab,,, "Source name" caldb, s, a, prod2,,, "Calibration database" irf, s, a, South_0.5h,,, "Instrument response function" -edisp, b, h, no,,, "Apply energy dispersion?" +inexclusion, f, h, NONE,,, "Sky map of exclusion regions for On/Off analysis" +edisp, b, h, no,,, "Apply energy dispersion in Cube analysis?" outfile, f, a, phasecrv.fits,,, "Output fits file name." # @@ -30,16 +31,25 @@ outfile, f, a, phasecrv.fits,,, "Output fits file name." phbinalg, s, a, LIN,LIN|FILE,, "Algorithm for defining phase bins" phbins, i, a, 5,1,10000, "Number of phase bins" phbinfile, f, a, NONE,,, "File defining the phase binning" +method, s, a, 3D, 3D|ONOFF,, "Analysis method" emin, r, a, 0.1,,, "Lower energy limit of events (TeV)" emax, r, a, 100.0,,, "Upper energy limit of events (TeV)" -enumbins, i, a, 0,0,100, "Number of energy bins for binned (0=unbinned)" +enumbins, i, a, 0,0,100, "Number of energy bins for binned (0=unbinned for 3D analysis only))" coordsys, s, a, CEL,CEL|GAL,, "Coordinate System" proj, s, a, CAR,AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method" -xref, r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)" -yref, r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)" +xref, r, a, 83.63,0,360, "First coordinate of image/source region center in degrees (RA or galactic l)" +yref, r, a, 22.01,-90,90, "Second coordinate of image/source region center in degrees (DEC or galactic b)" nxpix, i, a, 200,,, "Size of the X axis in pixels" nypix, i, a, 200,,, "Size of the Y axis in pixels" binsz, r, a, 0.02,,, "Pixel size (deg/pixel)" +srcshape, s, h, CIRCLE,CIRCLE,, "Shape of the source region" +rad, r, a, 0.2,0,180, "Radius of source region circle (deg)" +bkgmethod, s, h, REFLECTED,REFLECTED,, "Method for background estimation" +bkgregmin, i, h, 2,1,1000, "Minimum number of background regions" +maxoffset, r, h, 4.0,,, "Maximum offset of source from camera center to accept observation" +etruemin, r, h, 0.01,,, "Minimum true energy (TeV)" +etruemax, r, h, 300.0,,, "Maximum true energy (TeV)" +etruebins, i, h, 30,,, "Number of bins per decade for true energy bins" statistic, s, h, DEFAULT,DEFAULT|CSTAT|WSTAT|CHI2,, "Optimization statistic" # diff --git a/cscripts/csphasecrv.py b/cscripts/csphasecrv.py index 194bddb..ff3804d 100755 --- a/cscripts/csphasecrv.py +++ b/cscripts/csphasecrv.py @@ -63,6 +63,7 @@ class csphasecrv(ctools.csobservation): # [[ph1min,ph1max], [ph2min,ph2max],..] self._srcname = '' self._phbins = [[0.0,1.0]] + self._onoff = False self._stacked = False self._fits = gammalib.GFits() self._fitmodels = {} @@ -92,8 +93,13 @@ class csphasecrv(ctools.csobservation): # Get phase boundaries self._phbins = self._create_tbounds() + # Set On/Off analysis flag and query relevant user parameters + self._onoff = self._is_onoff() + + # If cube analysis is selected # Set stacked analysis flag and query relevant user parameters - self._stacked = self._is_stacked() + if not self._onoff: + self._stacked = self._is_stacked() # Query the hidden parameters, just in case self['edisp'].boolean() @@ -236,11 +242,20 @@ class csphasecrv(ctools.csobservation): # Store fit results phname = str(phmin)+'-'+str(phmax) - source = self._fitmodels[phname][self._srcname] - for par in pars: - result['values'][par] = source[par].value() - result['values']['e_'+par] = source[par].error() - + + # If the model contains the source of interest fill results + try: + source = self._fitmodels[phname][self._srcname] + for par in pars: + result['values'][par] = source[par].value() + result['values']['e_'+par] = source[par].error() + + # ... otherwise fills with zeros + except: + for par in pars: + result['values'][par] = 0. + result['values']['e_'+par] = 0. + # Append result to list of dictionaries results.append(result) @@ -348,36 +363,57 @@ class csphasecrv(ctools.csobservation): select['dec'] = 'UNDEFINED' select['expr'] = 'PHASE>'+str(phbin[0])+' && PHASE<'+str(phbin[1]) select.run() - + + # Set phase string + phstr = str(phbin[0]) + '-' + str(phbin[1]) + # Add phase to observation id for i in range(0,select.obs().size()): oldid = select.obs()[i].id() - select.obs()[i].id(oldid+'_'+str(phbin[0])+'-'+str(phbin[1])) + select.obs()[i].id(oldid+'_'+phstr) obs = select.obs() - # If a stacked analysis is requested then bin the events - # and compute the stacked response functions and setup + # If an On/Off analysis is requested generate the On/Off observations + if self._onoff: + obs = obsutils.get_onoff_obs(self, select.obs()) + + # ... otherwise, if stacked analysis is requested then bin the + # events and compute the stacked response functions and setup # an observation container with a single stacked observation. - if self._stacked: + elif self._stacked: obs = obsutils.get_stacked_obs(self, select.obs()) - + # Header self._log_header3(gammalib.EXPLICIT, 'Fitting the data') - - # Do maximum likelihood model fitting - like = ctools.ctlike(obs) - like['edisp'] = self['edisp'].boolean() - like.run() - - # Renormalize models to phase selection - # TODO move the scaling from the temporal to the spectral component - for model in like.obs().models(): - scaled_norm = model['Normalization'].value()/(phbin[1]-phbin[0]) - model['Normalization'].value(scaled_norm) - - # Store fit model - self._fitmodels[str(phbin[0])+'-'+str(phbin[1])] = \ - like.obs().models().copy() + + # The On/Off analysis can produce empty observation containers, + # e.g., when on-axis observations are used. To avoid ctlike asking + # for a new observation container (or hang, if in interactive mode) + # we'll run ctlike only if the size is >0 + if obs.size() > 0: + + # Do maximum likelihood model fitting + like = ctools.ctlike(obs) + like['edisp'] = self['edisp'].boolean() + like.run() + + # Renormalize models to phase selection + # TODO move the scaling from the temporal to the spectral component + for model in like.obs().models(): + scaled_norm = model['Normalization'].value()/(phbin[1]-phbin[0]) + model['Normalization'].value(scaled_norm) + + # Store fit model + self._fitmodels[phstr] = like.obs().models().copy() + + # ... otherwise we set an empty model container + else: + self._log_string(gammalib.TERSE, + 'PHASE %f - %f: no observations available' + ' for fitting' %(phbin[0], phbin[1])) + + # Set empty models container + self._fitmodels[phstr] = gammalib.GModels() # Create FITS file self._create_fits() diff --git a/cscripts/obsutils.py b/cscripts/obsutils.py index 864d62d..10af62c 100755 --- a/cscripts/obsutils.py +++ b/cscripts/obsutils.py @@ -760,3 +760,71 @@ def get_stacked_obs(cls, obs): # Return new oberservation container return new_obs + + +# ================================ # +# Get On/Off observation container # +# ================================ # +def get_onoff_obs(cls, obs): + """ + Create On/Off observations container from given observations + + Parameters + ---------- + cls : `~ctools.cscript` + cscript class + obs : `~gammalib.GObservations` + Observation container + + Returns + ------- + onoff_obs : `~gammalib.GObservations` + Observation container with On/Off observations + """ + # Write header + if cls._logExplicit(): + cls._log.header3('Creating On/Off observations') + + # Create On/Off observations + phagen = cscripts.csphagen(obs) + phagen['inexclusion'] = cls['inexclusion'].value() + phagen['emin'] = cls['emin'].real() + phagen['emax'] = cls['emax'].real() + phagen['enumbins'] = cls['enumbins'].integer() + phagen['ebinalg'] = 'LOG' + phagen['srcshape'] = cls['srcshape'].string() + phagen['coordsys'] = cls['coordsys'].string() + if cls['coordsys'].string() == 'CEL': + phagen['ra'] = cls['xref'].real() + phagen['dec'] = cls['yref'].real() + elif cls['coordsys'].string() == 'GAL': + phagen['glon'] = cls['xref'].real() + phagen['glat'] = cls['yref'].real() + if cls['srcshape'].string() == 'CIRCLE': + phagen['rad'] = cls['rad'].real() + phagen['bkgmethod'] = cls['bkgmethod'].string() + if cls['bkgmethod'].string() == 'REFLECTED': + phagen['bkgregmin'] = cls['bkgregmin'].integer() + phagen['maxoffset'] = cls['maxoffset'].real() + phagen['stack'] = True + phagen['etruemin'] = cls['etruemin'].real() + phagen['etruemax'] = cls['etruemax'].real() + phagen['etruebins'] = cls['etruebins'].integer() + phagen['chatter'] = cls['chatter'].integer() + phagen['clobber'] = cls['clobber'].boolean() + phagen['debug'] = cls['debug'].boolean() + phagen.run() + + # Clone resulting observation container + onoff_obs = phagen.obs().copy() + + # On/Off observations are created with CSTAT as default statistic + # We will change this to the user choice + if cls['statistic'].string() == 'DEFAULT': + pass + else: + for observation in onoff_obs: + observation.statistic(cls['statistic'].string()) + + # Return On/Off oberservation container + return onoff_obs diff --git a/doc/source/users/reference_manual/cslightcrv.rst b/doc/source/users/reference_manual/cslightcrv.rst index 932da9c..3c6500e 100644 --- a/doc/source/users/reference_manual/cslightcrv.rst +++ b/doc/source/users/reference_manual/cslightcrv.rst @@ -42,13 +42,19 @@ General parameters ``irf [string]`` Instrumental response function. +``(inexclusion = NONE) [file]`` + Optional FITS file containing a WCS map in the first hdu that defines sky + regions not to be used for background estimation in ``ONOFF`` analysis (where + map value != 0). + ``(edisp = no) [boolean]`` - Applies energy dispersion to response computation. + Applies energy dispersion to response computation (for ``3D`` analysis only, + energy dispersion is always taken into account in ``ONOFF`` analysis). ``outfile [file]`` Name of the light curve output file. -``(tbinalg = GTI) [string]`` +``tbinalg [string]`` Algorithm for defining time bins. ``tmin [time]`` @@ -63,6 +69,10 @@ General parameters ``tbinfile [file]`` File defining the time binning. +``method <3D|ONOFF> [string]`` + Selects between ``3D`` analysis (3D spatial/energy likelihood) and ``ONOFF`` + analysis (1D likelihood with background from Off regions). + ``emin [real]`` Lower energy limit of events (in TeV). @@ -70,7 +80,7 @@ General parameters Upper energy limit of events (in TeV). ``enumbins [integer]`` - Number of energy bins per light curve bin (0=unbinned). + Number of energy bins per light curve bin (0=unbinned for ``3D`` analysis only). ``coordsys [string]`` Coordinate system (CEL - celestial, GAL - galactic). @@ -79,19 +89,52 @@ General parameters Projection method. ``xref [real]`` - Right Ascension / Galactic longitude of image centre (J2000, in degrees). + Right Ascension / Galactic longitude of image centre for ``3D`` analysis or + source region centre for ``ONOFF`` analysis (J2000, in degrees). ``yref [real]`` - Declination / Galactic latitude of image centre (J2000, in degrees). + Declination / Galactic latitude of image centre for ``3D`` analysis or + source region centre for ``ONOFF`` analysis (J2000, in degrees). ``nxpix [integer]`` - Size of the Right Ascension / Galactic longitude axis (in pixels). + Size of the Right Ascension / Galactic longitude axis for ``3D`` analysis (in pixels). ``nypix [integer]`` - Size of the Declination / Galactic latitude axis (in pixels). + Size of the Declination / Galactic latitude axis for ``3D`` analysis (in pixels). ``binsz [real]`` - Pixel size (in degrees/pixel). + Pixel size for ``3D`` analysis (in degrees/pixel). + +``(srcshape = CIRCLE) [string]`` + Shape of the source region for ``ONOFF`` analysis. + ``CIRCLE``: circular region around given position. + +``rad [real]`` + Radius of source region circle for ``ONOFF`` analysis (deg) + +``(bkgmethod = REFLECTED) [string]`` + Method for background estimation in ``ONOFF`` analysis. + ``REFLECTED:`` background evaluated in regions with the same shape as + source region reflected w.r.t. pointing direction for each observation. + +``(bkgregmin = 2) [integer]`` + Minimum number of background regions that are required for an observation in + ``ONOFF`` analysis. If this number of background regions is not available the + observation is skipped. + +``(maxoffset = 4.0) [real]`` + Maximum offset in degrees of source from camera center to accept the + observation for ``ONOFF`` analysis. + +``(etruemin = 0.01) [real]`` + Minimum true energy to evaluate instrumental response in ``ONOFF`` analysis (TeV). + +``(etruemax = 0.01) [real]`` + Maximum true energy to evaluate instrumental response in ``ONOFF`` analysis (TeV). + +``(etruebins = 30) [integer]`` + Number of bins per decade for true energy bins to evaluate instrumental + response in ``ONOFF`` analysis. ``(statistic = DEFAULT) [string]`` Optimization statistic. ``DEFAULT`` uses the default statistic for all diff --git a/doc/source/users/reference_manual/csphasecrv.rst b/doc/source/users/reference_manual/csphasecrv.rst index ad3d124..a7d5412 100644 --- a/doc/source/users/reference_manual/csphasecrv.rst +++ b/doc/source/users/reference_manual/csphasecrv.rst @@ -10,14 +10,16 @@ Synopsis -------- This script computes spectra by performing a maximum likelihood fit using -:doc:`ctlike` in a series of phase bins for pulsars. The phase bins can be +:doc:`ctlike` in a series of phase bins for periodic sources. The input event +list must have a ``PHASE`` column containing the phase values, that can be +assigned using :doc:`ctphase`. The phase bins can be either specified in an ASCII file or as an interval divided into equally sized phase bins. The format of the ASCII file is one row per phase bin, each specifying the start of stop value of the phase bin, separated by a whitespace. The phase goes from 0.0 to 1.0. -On output the script writes the fitting results into a FITS file. The script -also produces one XML file per phase bin that contains the fitting results for +The script writes the fit results into a FITS file. The script +also produces one XML file per phase bin that contains the best-fit model for that bin. @@ -39,8 +41,14 @@ General parameters ``irf [string]`` Instrumental response function. +``(inexclusion = NONE) [file]`` + Optional FITS file containing a WCS map in the first hdu that defines sky + regions not to be used for background estimation in ``ONOFF`` analysis (where + map value != 0). + ``(edisp = no) [boolean]`` - Applies energy dispersion to response computation. + Applies energy dispersion to response computation (for ``3D`` analysis only, + energy dispersion is always taken into account in ``ONOFF`` analysis). ``outfile [file]`` Name of the XML output file. The phase interval will be automatically @@ -55,6 +63,10 @@ General parameters ``phbinfile [file]`` File defining the phase binning. +``method <3D|ONOFF> [string]`` + Selects between 3D analysis (3D spatial/energy likelihood) and ONOFF + analysis (1D likelihood with background from Off regions). + ``emin [real]`` Lower energy limit of events (in TeV). @@ -71,19 +83,52 @@ General parameters Projection method. ``xref [real]`` - Right Ascension / Galactic longitude of image centre (J2000, in degrees). + Right Ascension / Galactic longitude of image centre for ``3D`` analysis or + source region centre for ``ONOFF`` analysis (J2000, in degrees). ``yref [real]`` - Declination / Galactic latitude of image centre (J2000, in degrees). + Declination / Galactic latitude of image centre for ``3D`` analysis or + source region centre for ``ONOFF`` analysis (J2000, in degrees). ``nxpix [integer]`` - Size of the Right Ascension / Galactic longitude axis (in pixels). + Size of the Right Ascension / Galactic longitude axis for ``3D`` analysis (in pixels). ``nypix [integer]`` - Size of the Declination / Galactic latitude axis (in pixels). + Size of the Declination / Galactic latitude axis for ``3D`` analysis (in pixels). ``binsz [real]`` - Pixel size (in degrees/pixel). + Pixel size for ``3D`` analysis (in degrees/pixel). + +``(srcshape = CIRCLE) [string]`` + Shape of the source region for ``ONOFF`` analysis. + ``CIRCLE``: circular region around given position. + +``rad [real]`` + Radius of source region circle for ``ONOFF`` analysis (deg) + +``(bkgmethod = REFLECTED) [string]`` + Method for background estimation in ``ONOFF`` analysis. + ``REFLECTED:`` background evaluated in regions with the same shape as + source region reflected w.r.t. pointing direction for each observation. + +``(bkgregmin = 2) [integer]`` + Minimum number of background regions that are required for an observation in + ``ONOFF`` analysis. If this number of background regions is not available the + observation is skipped. + +``(maxoffset = 4.0) [real]`` + Maximum offset in degrees of source from camera center to accept the + observation for ``ONOFF`` analysis. + +``(etruemin = 0.01) [real]`` + Minimum true energy to evaluate instrumental response in ``ONOFF`` analysis (TeV). + +``(etruemax = 0.01) [real]`` + Maximum true energy to evaluate instrumental response in ``ONOFF`` analysis (TeV). + +``(etruebins = 30) [integer]`` + Number of bins per decade for true energy bins to evaluate instrumental + response in ``ONOFF`` analysis. ``(statistic = DEFAULT) [string]`` Optimization statistic. ``DEFAULT`` uses the default statistic for all @@ -129,4 +174,5 @@ Standard parameters Related tools or scripts ------------------------ +:doc:`ctphase` :doc:`ctlike` diff --git a/doc/source/users/tutorials/howto/howto_lightcurve.rst b/doc/source/users/tutorials/howto/howto_lightcurve.rst index c1ae9ec..d45f456 100644 --- a/doc/source/users/tutorials/howto/howto_lightcurve.rst +++ b/doc/source/users/tutorials/howto/howto_lightcurve.rst @@ -34,12 +34,13 @@ The time interval is linearly divided into 20 time bins. Input model definition XML file [$CTOOLS/share/models/crab.xml] models.xml Source name [Crab] Src001 Algorithm for defining time bins (FILE|LIN|GTI) [GTI] LIN - Lightcurve start time (MJD) [51544.5] 59235.5 - Lightcurve stop time (MJD) [51544.6] 59236.0 + Lightcurve start time (UTC string, JD, MJD or MET in seconds) [MJD 51544.5] MJD 59235.5 + Lightcurve stop time (UTC string, JD, MJD or MET in seconds) [MJD 51544.6] MJD 59236.0 Number of time bins (1-10000) [5] 20 - Number of energy bins for binned (0=unbinned) (0-100) [0] + Analysis method (3D|ONOFF) [3D] Lower energy limit of events (TeV) [0.1] Upper energy limit of events (TeV) [100.0] + Number of energy bins for binned (0=unbinned for Cube analysis only) (0-100) [0] Output light curve file [lightcurve.fits] The resulting light curve is shown in the figure below. diff --git a/pyext/ctool.i b/pyext/ctool.i index 377ba81..0e53f5b 100644 --- a/pyext/ctool.i +++ b/pyext/ctool.i @@ -100,6 +100,7 @@ public: %rename(_get_skydir) get_skydir; %rename(_set_outfile_name) set_outfile_name; %rename(_is_stacked) is_stacked; + %rename(_is_onoff) is_onoff; %rename(_set_response) set_response; %rename(_set_edisp) set_edisp; %rename(_restore_edisp) restore_edisp; @@ -146,6 +147,7 @@ public: // Protected methods that query user parameters bool is_stacked(void); + bool is_onoff(void); // Protected support methods void set_response(GObservations& obs); diff --git a/src/support/ctool.cpp b/src/support/ctool.cpp index e4fd7e5..bd6922a 100644 --- a/src/support/ctool.cpp +++ b/src/support/ctool.cpp @@ -1195,6 +1195,54 @@ bool ctool::is_stacked(void) /***********************************************************************//** + * @brief Query user parameters for On/Off analysis + * + * @return True if an On/Off analysis should be done. + * + * Queries the parameter "method", and if this is "ONOFF" signal that an On/ + * Off analysis is requested. In that case, all necessary input parameters of + * csphagen are queried: "inexclusion", "emin", "emax", "enumbins", + * "coordsys", "xref", "yref", "srcshape", "rad", "bkgmethod", "bkgregmin", + * "maxoffset", "etruemin", "etruemax", and "etruebins". + * + * Using this method assures that the parameters are always queried in the + * same order. + ***************************************************************************/ +bool ctool::is_onoff(void) +{ + // Set onoff flag + bool onoff = ((*this)["method"].string() == "ONOFF"); + + // Query analysis method chosen + if (onoff) { + + // Query csphagen parameters + (*this)["inexclusion"].filename(); + (*this)["emin"].real(); + (*this)["emax"].real(); + (*this)["enumbins"].integer(); + (*this)["coordsys"].string(); + (*this)["xref"].real(); + (*this)["yref"].real(); + if ((*this)["srcshape"].string() == "CIRCLE"){ + (*this)["rad"].real(); + } + if ((*this)["bkgmethod"].string() == "REFLECTED"){ + (*this)["bkgregmin"].integer(); + } + (*this)["maxoffset"].real(); + (*this)["etruemin"].real(); + (*this)["etruemax"].real(); + (*this)["etruebins"].integer(); + + } // endif: method was ONOFF + + // Return onoff flag + return onoff; +} + + +/***********************************************************************//** * @brief Set response for all CTA observations in container * * @param[in,out] obs Observation container diff --git a/src/support/ctool.hpp b/src/support/ctool.hpp index 26a6aea..97c66e9 100644 --- a/src/support/ctool.hpp +++ b/src/support/ctool.hpp @@ -104,6 +104,7 @@ protected: // Protected methods that query user parameters bool is_stacked(void); + bool is_onoff(void); // Protected methods for logging void log_observations(const GChatter& chatter, diff --git a/test/data/crab_offaxis1.fits.gz b/test/data/crab_offaxis1.fits.gz index 67a8082..7ded87e 100644 --- a/test/data/crab_offaxis1.fits.gz +++ b/test/data/crab_offaxis1.fits.gz diff --git a/test/data/crab_onoff.xml b/test/data/crab_onoff.xml new file mode 100644 index 0000000..cbe3e13 --- /dev/null +++ b/test/data/crab_onoff.xml @@ -0,0 +1,24 @@ + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/test/test_cscript.py b/test/test_cscript.py index 77b6f0f..3b32b78 100755 --- a/test/test_cscript.py +++ b/test/test_cscript.py @@ -455,6 +455,9 @@ class Test(test): # Test is_stacked() method # TODO: implement + # Test is_onoff() method + # TODO: implement + # Test set_response() method # TODO: implement diff --git a/test/test_cslightcrv.py b/test/test_cslightcrv.py index 0b73668..c69b271 100755 --- a/test/test_cslightcrv.py +++ b/test/test_cslightcrv.py @@ -42,6 +42,10 @@ class Test(test): # Call base class constructor test.__init__(self) + # Add off-axis events for classical analysis + self._offaxis_events = self._datadir + '/crab_offaxis1.fits' + self._model_onoff = self._datadir + '/crab_onoff.xml' + # Return return @@ -73,7 +77,7 @@ class Test(test): ' inmodel="'+self._model+'" srcname="Crab"'+ \ ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \ ' tbinalg="LIN" tmin="MJD 51544.50" tmax="MJD 51544.53"'+ \ - ' tbins=3 enumbins=0 emin=0.1 emax=100.0'+ \ + ' tbins=3 method="3D" enumbins=0 emin=0.1 emax=100.0'+ \ ' outfile="lightcurve_cmd1.fits"'+ \ ' logfile="cslightcrv_cmd1.log" chatter=1' @@ -86,14 +90,14 @@ class Test(test): 'Check successful execution from command line') # Check light curve - self._check_light_curve('lightcurve_cmd1.fits', 3) + self._check_light_curve('lightcurve_cmd1.fits', 3, 1.6e-16) # Setup cslightcrv command cmd = cslightcrv+' inobs="events_that_do_not_exist.fits"'+ \ ' inmodel="'+self._model+'" srcname="Crab"'+ \ ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \ ' tbinalg="LIN" tmin="MJD 51544.50" tmax="MJD 51544.53"'+ \ - ' tbins=3 enumbins=0 emin=0.1 emax=100.0'+ \ + ' tbins=3 method="3D" enumbins=0 emin=0.1 emax=100.0'+ \ ' outfile="lightcurve_cmd1.fits"'+ \ ' logfile="cslightcrv_cmd2.log" debug=yes' @@ -123,6 +127,7 @@ class Test(test): lcrv['tmin'] = 'MJD 51544.50' lcrv['tmax'] = 'MJD 51544.53' lcrv['tbins'] = 3 + lcrv['method'] = '3D' lcrv['enumbins'] = 0 lcrv['emin'] = 0.1 lcrv['emax'] = 100.0 @@ -137,7 +142,7 @@ class Test(test): lcrv.save() # Check light curve - self._check_light_curve('cslightcrv_py1.fits', 3) + self._check_light_curve('cslightcrv_py1.fits', 3, 1.6e-16) # Now use FILE as time bin algorithm. For this we need first to # create an ASCII file. We use now 6 time bins. The ASCII file @@ -159,6 +164,7 @@ class Test(test): lcrv['irf'] = self._irf lcrv['tbinalg'] = 'FILE' lcrv['tbinfile'] = 'cslightcrv_py2.dat' + lcrv['method'] = '3D' lcrv['enumbins'] = 0 lcrv['emin'] = 0.1 lcrv['emax'] = 100.0 @@ -171,7 +177,7 @@ class Test(test): lcrv.execute() # Check light curve - self._check_light_curve('cslightcrv_py2.fits', 2) + self._check_light_curve('cslightcrv_py2.fits', 2, 1.6e-16) # Now we setup an observation container on input. We attached the # model to the observation container so that cslightcrv should @@ -188,6 +194,7 @@ class Test(test): lcrv['caldb'] = self._caldb lcrv['irf'] = self._irf lcrv['tbinalg'] = 'GTI' + lcrv['method'] = '3D' lcrv['enumbins'] = 0 lcrv['emin'] = 0.1 lcrv['emax'] = 100.0 @@ -199,9 +206,9 @@ class Test(test): lcrv.execute() # Check light curve - self._check_light_curve('cslightcrv_py3.fits', 1) + self._check_light_curve('cslightcrv_py3.fits', 1, 1.6e-16) - # Finally we set-up a binned cslightcrv + # Binned cslightcrv lcrv = cscripts.cslightcrv() lcrv['inobs'] = self._events lcrv['inmodel'] = self._model @@ -212,6 +219,7 @@ class Test(test): lcrv['tmin'] = 'MJD 51544.50' lcrv['tmax'] = 'MJD 51544.53' lcrv['tbins'] = 2 + lcrv['method'] = '3D' lcrv['emin'] = 0.1 lcrv['emax'] = 100.0 lcrv['enumbins'] = 10 @@ -230,13 +238,43 @@ class Test(test): lcrv.execute() # Check light curve - self._check_light_curve('cslightcrv_py4.fits', 2) + self._check_light_curve('cslightcrv_py4.fits', 2, 1.6e-16) + + # cslightcrv with classical analysis + lcrv = cscripts.cslightcrv() + lcrv['inobs'] = self._offaxis_events + lcrv['inmodel'] = self._model_onoff + lcrv['srcname'] = 'Crab' + lcrv['caldb'] = self._caldb + lcrv['irf'] = self._irf + lcrv['tbinalg'] = 'LIN' + lcrv['tmin'] = 'MJD 58849.0008' + lcrv['tmax'] = 'MJD 58849.0112' + lcrv['tbins'] = 2 + lcrv['method'] = 'ONOFF' + lcrv['emin'] = 0.1 + lcrv['emax'] = 100.0 + lcrv['enumbins'] = 10 + lcrv['coordsys'] = 'CEL' + lcrv['xref'] = 83.63 + lcrv['yref'] = 22.01 + lcrv['rad'] = 0.2 + lcrv['statistic'] = 'WSTAT' + lcrv['outfile'] = 'cslightcrv_py5.fits' + lcrv['logfile'] = 'cslightcrv_py5.log' + lcrv['chatter'] = 4 + + # Execute cslightcrv script + lcrv.execute() + + # Check light curve + self._check_light_curve('cslightcrv_py5.fits', 2, 6.e-16) # Return return # Check light curve result file - def _check_light_curve(self, filename, bins): + def _check_light_curve(self, filename, bins,pref): """ Check light curve file """ @@ -266,6 +304,13 @@ class Test(test): self.test_assert(table.contains(col), 'FITS file contains "'+col+'" column') + # Check that table has been filled + # Prefactor has right order of magnitude + for s in range(table.nrows()): + self.test_value(table['Prefactor'][s], pref, 0.2*pref, + 'Check prefactor value') + + # Close FITS file fits.close() diff --git a/test/test_csphasecrv.py b/test/test_csphasecrv.py index d696641..bb318c1 100755 --- a/test/test_csphasecrv.py +++ b/test/test_csphasecrv.py @@ -44,6 +44,8 @@ class Test(test): # Set test data self._phased_events = self._datadir + '/phased_events.fits' + self._offaxis_events = self._datadir + '/crab_offaxis1.fits' + self._model_onoff = self._datadir + '/crab_onoff.xml' # Return return @@ -75,7 +77,7 @@ class Test(test): cmd = csphasecrv+' inobs="'+self._phased_events+'"'+ \ ' inmodel="'+self._model+'" srcname="Crab"'+ \ ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \ - ' phbinalg="LIN" phbins=5'+ \ + ' phbinalg="LIN" phbins=5 method="3D"'+ \ ' enumbins=0 emin=0.1 emax=100.0'+ \ ' outfile="csphasecrv_cmd1.fits"'+ \ ' logfile="csphasecrv_cmd1.log" chatter=1' @@ -95,7 +97,7 @@ class Test(test): cmd = csphasecrv+' inobs="events_that_do_not_exist.fits"'+ \ ' inmodel="'+self._model+'" srcname="Crab"'+ \ ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \ - ' phbinalg="LIN" phbins=5'+ \ + ' phbinalg="LIN" phbins=5 method="3D"'+ \ ' enumbins=0 emin=0.1 emax=100.0'+ \ ' outfile="csphasecrv_cmd2.fits"'+ \ ' logfile="csphasecrv_cmd2.log" debug=yes chatter=1' @@ -124,6 +126,7 @@ class Test(test): pcrv['irf'] = self._irf pcrv['phbinalg'] = 'LIN' pcrv['phbins'] = 5 + pcrv['method'] = '3D' pcrv['enumbins'] = 0 pcrv['emin'] = 0.1 pcrv['emax'] = 100.0 @@ -160,6 +163,7 @@ class Test(test): pcrv['irf'] = self._irf pcrv['phbinalg'] = 'FILE' pcrv['phbinfile'] = 'csphasecrv_py2.dat' + pcrv['method'] = '3D' pcrv['enumbins'] = 0 pcrv['emin'] = 0.1 pcrv['emax'] = 100.0 @@ -190,6 +194,7 @@ class Test(test): pcrv['irf'] = self._irf pcrv['phbinalg'] = 'LIN' pcrv['phbins'] = 5 + pcrv['method'] = '3D' pcrv['enumbins'] = 0 pcrv['emin'] = 0.1 pcrv['emax'] = 100.0 @@ -204,7 +209,7 @@ class Test(test): # Check phase curve self._check_phase_curve('csphasecrv_py3.fits', 5) - # Finally we set-up a binned csphasecrv + # Binned csphasecrv pcrv = cscripts.csphasecrv() pcrv['inobs'] = self._phased_events pcrv['inmodel'] = self._model @@ -213,6 +218,7 @@ class Test(test): pcrv['irf'] = self._irf pcrv['phbinalg'] = 'LIN' pcrv['phbins'] = 2 + pcrv['method'] = '3D' pcrv['emin'] = 0.1 pcrv['emax'] = 100.0 pcrv['enumbins'] = 10 @@ -234,6 +240,35 @@ class Test(test): # Check phase curve self._check_phase_curve('csphasecrv_py4.fits', 2) + # On/Off csphasecrv + pcrv = cscripts.csphasecrv() + pcrv['inobs'] = self._offaxis_events + pcrv['inmodel'] = self._model_onoff + pcrv['srcname'] = 'Crab' + pcrv['caldb'] = self._caldb + pcrv['irf'] = self._irf + pcrv['phbinalg'] = 'LIN' + pcrv['phbins'] = 2 + pcrv['method'] = 'ONOFF' + pcrv['statistic'] = 'WSTAT' + pcrv['emin'] = 0.1 + pcrv['emax'] = 100.0 + pcrv['enumbins'] = 10 + pcrv['coordsys'] = 'CEL' + pcrv['xref'] = 83.63 + pcrv['yref'] = 22.01 + pcrv['rad'] = 0.2 + pcrv['outfile'] = 'csphasecrv_py5.fits' + pcrv['logfile'] = 'csphasecrv_py5.log' + pcrv['chatter'] = 4 + + # Execute csphasecrv script + pcrv.logFileOpen() # Make sure we get a log file + pcrv.execute() + + # Check phase curve + self._check_phase_curve('csphasecrv_py5.fits', 2) + # Return return @@ -267,6 +302,12 @@ class Test(test): self.test_assert(table.contains(col), 'FITS file contains "'+col+'" column') + # Check that table has been filled + # Prefactor has right order of magnitude + for s in range(table.nrows()): + self.test_value(table['Prefactor'][s], 4.e-16, 3.e-16, + 'Check prefactor value') + # Close FITS file fits.close()