diff --git a/cscripts/cspull.py b/cscripts/cspull.py index f1bcf65..a49e478 100755 --- a/cscripts/cspull.py +++ b/cscripts/cspull.py @@ -36,25 +36,19 @@ class cspull(ctools.cscript): # Constructor def __init__(self, *argv): """ - Constructor. + Constructor """ # Set name self._name = 'cspull' - self._version = '1.1.0' + self._version = '1.2.0' # Initialise some members self._edisp = False - self._exposure = None - self._psfcube = None - self._bckcube = None - self._edispcube = None - self._stackmodels = None 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' @@ -62,11 +56,10 @@ class cspull(ctools.cscript): self._seed = 1 self._chatter = 2 - # Initialise observation container from constructor arguments. + # Initialise observation container from constructor arguments self._obs, argv = self._set_input_obs(argv) - # Initialise application by calling the appropriate class - # constructor. + # Initialise application by calling the appropriate class constructor self._init_cscript(argv) # Return @@ -108,7 +101,6 @@ class cspull(ctools.cscript): 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() @@ -116,6 +108,7 @@ class cspull(ctools.cscript): self._chatter = self['chatter'].integer() # Query some parameters + self['outfile'].filename() self['profile'].boolean() # Write input parameters into logger @@ -135,6 +128,9 @@ class cspull(ctools.cscript): obs : `~gammalib.GObservations` Observation container. """ + # Set duration + duration = self['tmax'].real() - self['tmin'].real() + # Setup observation definition list obsdeflist = obsutils.set_obs_patterns(self._pattern, ra=self['ra'].real(), @@ -144,7 +140,7 @@ class cspull(ctools.cscript): # Create list of observations obs = obsutils.set_obs_list(obsdeflist, tstart=self['tmin'].real(), - duration=self['tmax'].real()-self['tmin'].real(), + duration=duration, deadc=self['deadc'].real(), emin=self['emin'].real(), emax=self['emax'].real(), @@ -155,128 +151,6 @@ class cspull(ctools.cscript): # Return observation container return obs - def _set_stacked_irf(self): - """ - Prepare stacked IRFs for stacked analysis - """ - # Write header into logger - if self._logTerse(): - self._log('\n') - self._log.header1('Compute stacked response') - - # Set number of energy bins to at least 30 per energy decade - enumbins = int((math.log10(self['emax'].real()) - - math.log10(self['emin'].real())) * 30.0) - if self._enumbins > enumbins: - enumbins = self._enumbins - - # Compute spatial binning for point spread function and - # energy dispersion cubes - binsz = 10.0 * self['binsz'].real() - nxpix = self['npix'].integer() // 10 # Make sure result is int - nypix = self['npix'].integer() // 10 # Make sure result is int - if nxpix < 2: - nxpix = 2 - if nypix < 2: - nypix = 2 - - # 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'] = 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() - - # Notify exposure computation - if self._logTerse(): - 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'] = binsz - psfcube['nxpix'] = nxpix - psfcube['nypix'] = nypix - psfcube['enumbins'] = 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() - - # Notify Psf computation - if self._logTerse(): - 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'] = binsz - edispcube['nxpix'] = nxpix - edispcube['nypix'] = nypix - edispcube['enumbins'] = 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') - - # 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'] = 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() - - # Notify background cube computation - if self._logTerse(): - self._log('Computed background cube\n') - - # Store results - self._exposure = expcube.expcube().copy() - self._psfcube = psfcube.psfcube().copy() - self._bckcube = bkgcube.bkgcube().copy() - self._stackmodels = bkgcube.models().copy() - - # Return - return - def _obs_string(self, obs): """ Generate summary string for observation @@ -310,7 +184,7 @@ class cspull(ctools.cscript): # Return summary string return text - def _trial(self, seed, stacked): + def _trial(self, seed): """ Compute the pull for a single trial @@ -318,8 +192,6 @@ class cspull(ctools.cscript): ---------- seed : int Random number generator seed - stacked : bool - Use stacked analysis Returns ------- @@ -354,16 +226,6 @@ class cspull(ctools.cscript): debug=self._logDebug(), chatter=self._chatter) - # Set response for a stacked observation - if stacked: - if self._edisp: - obs[0].response(self._exposure, self._psfcube, - self._edispcube, self._bckcube) - else: - obs[0].response(self._exposure, self._psfcube, - self._bckcube) - obs.models(self._stackmodels) - # Determine number of events in simulation nevents = 0.0 for run in obs: @@ -427,10 +289,10 @@ class cspull(ctools.cscript): # Set parameter name name = model_name+'_'+par.name() - # Append parameter, Pull_parameter and Unc_parameter + # Append parameter, Pull_parameter and e_parameter colnames.append(name) colnames.append('Pull_'+name) - colnames.append('Unc_'+name) + colnames.append('e_'+name) # Compute pull fitted_value = par.value() @@ -444,7 +306,7 @@ class cspull(ctools.cscript): # Store results values[name] = fitted_value values['Pull_'+name] = pull - values['Unc_'+name] = error + values['e_'+name] = error # Write result if self._logNormal(): @@ -481,13 +343,10 @@ class cspull(ctools.cscript): self._log.header1(gammalib.number('Observation',len(self._obs))) self._log(str(self._obs)) self._log('\n') - - # Signal if stacked analysis is requested - stacked = self._obs.size() > 1 and self._enumbins > 0 - - # If several observations and binned: prepare stacked irfs - if stacked: - self._set_stacked_irf() + if self._logExplicit(): + for obs in self._obs: + self._log(str(obs)) + self._log('\n') # Write header if self._logTerse(): @@ -498,18 +357,21 @@ class cspull(ctools.cscript): for seed in range(self._ntrials): # Make a trial and add initial seed - result = self._trial(seed + self._seed, stacked) + result = self._trial(seed + self._seed) + + # Set output filename + outfile = self['outfile'].filename().url() # Write out result immediately if seed == 0: - f = open(self._outfile.url(), 'w') + f = open(outfile, 'w') writer = csv.DictWriter(f, result['colnames']) headers = {} for n in result['colnames']: headers[n] = n writer.writerow(headers) else: - f = open(self._outfile.url(), 'a') + f = open(outfile, 'a') writer = csv.DictWriter(f, result['colnames']) writer.writerow(result['values']) f.close() diff --git a/cscripts/cstsdist.par b/cscripts/cstsdist.par index c8922c9..9656ce4 100644 --- a/cscripts/cstsdist.par +++ b/cscripts/cstsdist.par @@ -5,7 +5,7 @@ # ------------------------------------------------------------------------- # # # # File: cstsdist.par # -# Version: 1.1.0 # +# Version: 1.2.0 # # # # Author: Juergen Knoedlseder # # Jurgen.Knodlseder@irap.omp.eu # @@ -36,6 +36,8 @@ outfile, f, a, ts.dat,,, "Output Test Statistic distribution file" ntrials, i, a, 10,,, "Number of trials" ra, r, a, 83.63,0,360, "RA of pointing (deg)" dec, r, a, 22.01,-90,90, "Dec of pointing (deg)" +coordsys, s, a, CEL,CEL|GAL,, "Coordinate system (CEL - celestial, GAL - galactic)" +proj, s, a, TAN,AIT|AZP|CAR|MER|MOL|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)" diff --git a/cscripts/cstsdist.py b/cscripts/cstsdist.py index ae1c6cd..7f3e2f7 100755 --- a/cscripts/cstsdist.py +++ b/cscripts/cstsdist.py @@ -30,13 +30,7 @@ from cscripts import obsutils # ============== # class cstsdist(ctools.cscript): """ - Generates Test Statistic distribution for a model. - - - This class implements the TS 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 Test Statistic distribution for a model """ # Constructor @@ -45,24 +39,17 @@ class cstsdist(ctools.cscript): Constructor. """ # Set name - self._name = "cstsdist" - self._version = "1.1.0" + self._name = 'cstsdist' + self._version = '1.2.0' # Initialise some members - self._obs = gammalib.GObservations() - self._pattern = "single" - self._srcname = "" - self._enumbins = 0 - self._npix = 0 - self._binsz = 0.0 - self._outfile = gammalib.GFilename("ts.dat") - self._ntrials = 10 - self._edisp = False - self._debug = False - self._log_clients = False - - # Initialise application by calling the appropriate class - # constructor. + self._srcname = '' + self._log_clients = False + + # Initialise observation container from constructor arguments + self._obs, argv = self._set_input_obs(argv) + + # Initialise application by calling the appropriate class constructor self._init_cscript(argv) # Return @@ -74,165 +61,116 @@ class cstsdist(ctools.cscript): """ Get parameters from parfile and setup the observation. """ - # Set observation if not done before if self._obs == None or self._obs.size() == 0: self._obs = self._get_observations() - # Check for requested pattern and use above - # observation parameters to set wobble pattern - self._pattern = self["pattern"].string() - if self._pattern == "four": + # Check for requested pattern and use above observation parameters + # to set wobble pattern + if self['pattern'].string() == 'four': self._obs = self._set_obs() # Get source name - self._srcname = self["srcname"].string() - - # Get number of energy bins - self._enumbins = self["enumbins"].integer() + self._srcname = self['srcname'].string() - # Read parameters for binned if requested - if not self._enumbins == 0: - self._npix = self["npix"].integer() - self._binsz = self["binsz"].real() + # 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()) + self._obs.models(self['inmodel'].filename()) - # Get other parameters - self._outfile = self["outfile"].filename() - self._ntrials = self["ntrials"].integer() - self._edisp = self["edisp"].boolean() - self._debug = self["debug"].boolean() + # Query parameters + self['edisp'].boolean() + self['ntrials'].integer() + self['outfile'].filename() + self['debug'].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 an observation container. + Set an observation container - Returns: - Observation container. + Returns + ------- + obs : `~gammalib.GObservations` + Observation container """ + # Set duration + duration = self['tmax'].real() - self['tmin'].real() # Setup observation definition list - obsdeflist = obsutils.set_obs_patterns(self._pattern, - ra=self["ra"].real(), - dec=self["dec"].real(), - offset=self["offset"].real()) + obsdeflist = obsutils.set_obs_patterns(self['pattern'].string(), + 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 = duration, + 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_models(self, fitspat=False, fitspec=False): + def _trial(self, seed): """ - Set full model and background model. + Create the TS for a single trial - Kwargs: - fitspat: Fit spatial parameter (default: False). - fitspec: Fit spectral parameters (default: False). - - Returns: - Tuple containing full model and background model. - """ - # Retrieve full model from observation container - full_model = self._obs.models().copy() - - # Get source model - model = full_model[self._srcname] - - # Check that model has a Prefactor - if not model.has_par("Prefactor"): - msg = "Model \""+self._srcname+"\" has no parameter "+\ - "\"Prefactor\". Only spectral models with a "+\ - "\"Prefactor\" parameter are supported." - raise RuntimeError(msg) - - # Fit or fix spatial parameters - if fitspat: - if model.has_par("RA"): - model["RA"].free() - if model.has_par("DEC"): - model["DEC"].free() - if model.has_par("Sigma"): - model["Sigma"].free() - if model.has_par("Radius"): - model["Radius"].free() - if model.has_par("Width"): - model["Width"].free() - else: - if model.has_par("RA"): - model["RA"].fix() - if model.has_par("DEC"): - model["DEC"].fix() - if model.has_par("Sigma"): - model["Sigma"].fix() - if model.has_par("Radius"): - model["Radius"].fix() - if model.has_par("Width"): - model["Width"].fix() - - # Fit or fix spectral parameters - if fitspec: - if model.has_par("Index"): - model["Index"].free() - if model.has_par("Cutoff"): - model["Cutoff"].free() - else: - if model.has_par("Index"): - model["Index"].fix() - if model.has_par("Cutoff"): - model["Cutoff"].fix() - - # Create background model - bkg_model = full_model.copy() - bkg_model.remove(self._srcname) - - # Return models - return full_model, bkg_model - - def _trial(self, seed, full_model, bkg_model): - """ - Create the TS for a single trial. + Parameters + ---------- + seed : int + Random number generator seed - Args: - seed: Random number generator seed - full_model: Full model - bkg_model: Background model - - Returns: + Returns + ------- + result : dict Result dictionary """ # Write header if self._logExplicit(): - self._log.header2("Trial "+str(seed+1)) + self._log.header2('Trial '+str(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, - seed=seed, - binsz=self._binsz, - npix=self._npix, - log=self._log_clients, - debug=self._debug) + nbins = self['enumbins'].integer(), + seed = seed, + proj = proj, + coord = coordsys, + binsz = binsz, + npix = npix, + log = self._log_clients, + debug = self['debug'].boolean()) # Determine number of events in simulation nevents = 0.0 @@ -241,110 +179,70 @@ class cstsdist(ctools.cscript): # Write simulation results if self._logExplicit(): - self._log.header3("Simulation") - self._log.parformat("Number of simulated events") + self._log.header3('Simulation') + self._log.parformat('Number of simulated events') self._log(nevents) - self._log("\n") - - # Fit background only - sim.models(bkg_model) - like_bgm = obsutils.fit(sim, - log=self._log_clients, - debug=self._debug) - result_bgm = like_bgm.obs().models().copy() - LogL_bgm = like_bgm.opt().value() - npred_bgm = like_bgm.obs().npred() - - # Write background fit results - if self._logExplicit(): - self._log.header3("Background model fit") - self._log.parformat("log likelihood") - self._log(LogL_bgm) - self._log("\n") - self._log.parformat("Number of predicted events") - self._log(npred_bgm) - self._log("\n") - for model in result_bgm: - self._log.parformat("Model") - self._log(model.name()) - self._log("\n") - for par in model: - self._log(str(par)+"\n") - - # Fit background and test source - sim.models(full_model) - like_all = obsutils.fit(sim, - log=self._log_clients, - debug=self._debug) - result_all = like_all.obs().models().copy() - LogL_all = like_all.opt().value() - npred_all = like_all.obs().npred() - ts = 2.0*(LogL_bgm-LogL_all) - - # Write background and test source fit results + self._log('\n') + + # Fit model + fit = obsutils.fit(sim, log = self._log_clients, + debug = self['debug'].boolean()) + + # Get model fitting results + logL = fit.opt().value() + npred = fit.obs().npred() + models = fit.obs().models() + model = models[self._srcname] + ts = model.ts() + + # Write fit results, either explicit or normal if self._logExplicit(): - self._log.header3("Background and test source model fit") - self._log.parformat("Test statistics") + self._log.header3('Test source model fit') + self._log.parformat('Test statistics') self._log(ts) - self._log("\n") - self._log.parformat("log likelihood") - self._log(LogL_all) - self._log("\n") - self._log.parformat("Number of predicted events") - self._log(npred_all) - self._log("\n") - for model in result_all: - self._log.parformat("Model") + self._log('\n') + self._log.parformat('log likelihood') + self._log(logL) + self._log('\n') + self._log.parformat('Number of predicted events') + self._log(npred) + self._log('\n') + for model in models: + self._log.parformat('Model') self._log(model.name()) - self._log("\n") + self._log('\n') for par in model: - self._log(str(par)+"\n") - - # Write result - elif self._logTerse(): - self._log.parformat("Trial "+str(seed)) - self._log("TS=") + self._log(str(par)+'\n') + elif self._logNormal(): + self._log.parformat('Trial '+str(seed)) + self._log('TS=') self._log(ts) - self._log(" Prefactor=") - self._log(result_all[self._srcname]["Prefactor"].value()) - self._log("+/-") - self._log(result_all[self._srcname]["Prefactor"].error()) - self._log("\n") + self._log(' Prefactor=') + self._log(model['Prefactor'].value()) + self._log('+/-') + self._log(model['Prefactor'].error()) + self._log('\n') # Initialise results colnames = [] values = {} # Set TS value - colnames.append("TS") - values["TS"] = ts - - # Set logL for background fit - colnames.append("LogL_bgm") - values["LogL_bgm"] = LogL_bgm - - # Set logL for full fit - colnames.append("LogL_all") - values["LogL_all"] = LogL_all + colnames.append('TS') + values['TS'] = ts # Set Nevents - colnames.append("Nevents") - values["Nevents"] = nevents - - # Set Npred for background fit - colnames.append("Npred_bkg") - values["Npred_bkg"] = npred_bgm + colnames.append('Nevents') + values['Nevents'] = nevents - # Set Npred for full fit - colnames.append("Npred_all") - values["Npred_all"] = npred_all + # Set Npred + colnames.append('Npred') + values['Npred'] = npred # Gather free full fit parameters - for i in range(result_all.size()): - model = result_all[i] + for model in models: model_name = model.name() - for k in range(model.size()): - par = model[k] + for par in model: if par.is_free(): # Set parameter name @@ -355,7 +253,7 @@ class cstsdist(ctools.cscript): values[name] = par.value() # Append error - name = "Unc_"+name + name = 'e_'+name colnames.append(name) values[name] = par.error() @@ -369,7 +267,7 @@ class cstsdist(ctools.cscript): # Public methods def run(self): """ - Run the script. + Run the script """ # Switch screen logging on in debug mode if self._logDebug(): @@ -378,55 +276,51 @@ class cstsdist(ctools.cscript): # Get parameters self._get_parameters() - # Initialise models - full_model, bkg_model = self._set_models() + # Set models for this observation + self.models(obsutils.set_ts_model(self._obs.models(), self._srcname)) # Write models into logger if self._logTerse(): - self._log("\n") - self._log.header1("Models") - self._log.header2("Background model") - self._log(str(bkg_model)) - self._log("\n\n") - self._log.header2("Full model") - self._log(str(full_model)) - self._log("\n") + self._log('\n') + self._log.header1('Models') + self._log(str(self._obs.models())) + self._log('\n') # Write observation into logger if self._logTerse(): - self._log("\n") - self._log.header1("Observation") + self._log('\n') + self._log.header1(gammalib.number('Observation',len(self._obs))) self._log(str(self._obs)) - self._log("\n") - - # Write models into logger - if self._logTerse(): - self._log("\n") - self._log.header1("Test model") - self._log(str(full_model)) - self._log("\n") + self._log('\n') + if self._logExplicit(): + for obs in self._obs: + self._log(str(obs)) + self._log('\n') # Write header if self._logTerse(): - self._log("\n") - self._log.header1("Generate TS distribution") + self._log('\n') + self._log.header1('Generate TS distribution') # Loop over trials - for seed in range(self._ntrials): + for seed in range(self['ntrials'].integer()): # Make a trial - result = self._trial(seed, full_model, bkg_model) + result = self._trial(seed) + + # Set output filename + outfile = self['outfile'].filename().url() # Write out result immediately if seed == 0: - f = open(self._outfile.url(), 'w') + f = open(outfile, 'w') writer = csv.DictWriter(f, result['colnames']) headers = {} for n in result['colnames']: headers[n] = n writer.writerow(headers) else: - f = open(self._outfile.url(), 'a') + f = open(outfile, 'a') writer = csv.DictWriter(f, result['colnames']) writer.writerow(result['values']) f.close() @@ -436,17 +330,22 @@ class cstsdist(ctools.cscript): def models(self, models): """ - Set model. + Set model + + Parameters + ---------- + models : `~gammalib.GModels` + Model container """ # Copy models - self._obs.models(models.copy()) + self._obs.models(models) # Return return def execute(self): """ - Execute the script. + Execute the script """ # Open logfile self.logFileOpen() diff --git a/cscripts/obsutils.py b/cscripts/obsutils.py index 0904b26..d87cfb5 100755 --- a/cscripts/obsutils.py +++ b/cscripts/obsutils.py @@ -1,5 +1,5 @@ # ========================================================================== -# CTA observation handling support functions. +# CTA observation handling support functions # # Copyright (C) 2011-2016 Juergen Knoedlseder # @@ -30,6 +30,14 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, """ Simulate events for all observations in the container + Simulate events for all observations using ctobssim. If the number of + energy bins is positive, the events are then binned in a counts cube + using ctbin. If multiple observations are simulated, the counts cube is + a stacked cube and the corresponding response cubes are computed using + ctexpcube, ctpsfcube, ctbkgcube and optionally ctedispcube. The response + cubes are attached to the first observation in the container, which + normally is the observation with the counts cube. + Parameters ---------- obs : `~gammalib.GObservations` @@ -42,7 +50,7 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, Chatter level edisp : bool, optional Apply energy dispersion? - seed : int, integer + seed : int, optional Seed value for simulations emin : float, optional Minimum energy of counts cube for binned (TeV) @@ -64,7 +72,6 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, obs : `~gammalib.GObservations` Observation container filled with simulated events """ - # Allocate ctobssim application and set parameters sim = ctools.ctobssim(obs) sim['seed'] = seed @@ -125,6 +132,48 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, # the container and bin the events in counts maps bin.run() + # If we have multiple input observations then create stacked response + # cubes and append them to the observation + if len(sim.obs()) > 1: + + # First get the cube. We need this because ctbin builds a cube on + # the fly. We should avoid this and make the cube part of ctbin + cube = bin.cube() + + # Get xref, yref, coordsys and proj from counts cube + _xref = cube.counts().projection().crval(0) + _yref = cube.counts().projection().crval(1) + _coordsys = cube.counts().projection().coordsys() + _proj = cube.counts().projection().code() + + # Kluge: map "EQU" to "CEL" coordinate system (not sure why the + # GSkyProjection::coordsys returns "EQU" instead of "CEL" + if _coordsys == 'EQU': + _coordsys = 'CEL' + + # Get stacked response + response = get_stacked_response(sim.obs(), _xref, _yref, + binsz=binsz, nxpix=npix, nypix=npix, + emin=emin, emax=emax, enumbins=nbins, + edisp=edisp, + coordsys=_coordsys, proj=_proj, + log=log, debug=debug, + chatter=chatter) + + # Set stacked response + if edisp: + bin.obs()[0].response(response['expcube'], + response['psfcube'], + response['edispcube'], + response['bkgcube']) + else: + bin.obs()[0].response(response['expcube'], + response['psfcube'], + response['bkgcube']) + + # Set new models + bin.obs().models(response['models']) + # Make a deep copy of the observation that will be returned # (the ctbin object will go out of scope one the function is # left) @@ -149,33 +198,36 @@ def sim(obs, log=False, debug=False, chatter=2, edisp=False, seed=0, # ================ # def fit(obs, log=False, debug=False, chatter=2, edisp=False): """ - Perform maximum likelihood fitting of observations in the container. - - Parameters: - obs - Observation container - Keywords: - log - Create log file(s) - debug - Create screen dump - chatter - Chatter level - edisp - Apply energy dispersion? + Perform maximum likelihood fitting of observations in the container + + Parameters + ---------- + obs : `~gammalib.GObservations` + Observation container + log : bool, optional + Create log file(s) + debug : bool, optional + Create console dump? + chatter : int, optional + Chatter level + edisp : bool, optional + Apply energy dispersion? + + Returns + ------- + like : `~ctools.ctlike` + ctlike application """ - # Allocate ctlike application + # Allocate ctlike application and set parameters like = ctools.ctlike(obs) + like['debug'] = debug + like['chatter'] = chatter + like['edisp'] = edisp # Optionally open the log file if log: like.logFileOpen() - # Optionally switch-on debugging model - if debug: - like["debug"] = True - - # Set chatter level - like["chatter"] = chatter - - # Optionally apply energy dispersion - like["edisp"] = edisp - # Run ctlike application. like.run() @@ -188,34 +240,36 @@ def fit(obs, log=False, debug=False, chatter=2, edisp=False): # ============================================================== # def cterror(obs, srcname, log=False, debug=False, chatter=2): """ - Perform maximum likelihood fitting of observations in the container. - - Parameters: - obs - Observation container - srcname - Source name - Keywords: - log - Create log file(s) - debug - Create screen dump - chatter - Chatter level - edisp - Apply energy dispersion? + Fit observations and determine errors using likelihood profile + + Parameters + ---------- + obs : `~gammalib.GObservations` + Observation container + srcname : str + Source name + log : bool, optional + Create log file(s) + debug : bool, optional + Create console dump? + chatter : int, optional + Chatter level + + Returns + ------- + error : `~ctools.cterror` + cterror application """ - # Allocate cterror application + # Allocate cterror application and set parameters error = ctools.cterror(obs) - - # Set cterror parameters - error["srcname"] = srcname + error['srcname'] = srcname + error['debug'] = debug + error['chatter'] = chatter # Optionally open the log file if log: error.logFileOpen() - # Optionally switch-on debugging model - if debug: - error["debug"] = True - - # Set chatter level - error["chatter"] = chatter - # Run cterror application. error.run() @@ -226,24 +280,39 @@ def cterror(obs, srcname, log=False, debug=False, chatter=2): # ================= # # Create counts map # # ================= # -def cntmap(obs, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \ +def cntmap(obs, proj='TAN', coord='GAL', xval=0.0, yval=0.0, \ binsz=0.05, nxpix=200, nypix=200, \ - outname="cntmap.fits"): + outname='cntmap.fits'): """ - Creates a counts map by combining the events of all observations. - The counts map will be a summed map over all energies. - - Parameters: - obs - Observation container - Keywords: - proj - Projection type (e.g. TAN, CAR, STG, ...) (default: TAN) - coord - Coordinate type (GAL, CEL) (default: GAL) - xval - Reference longitude value [deg] (default: 0.0) - yval - Reference latitude value [deg] (default: 0.0) - binsz - Pixel size [deg/pixel] (default: 0.05) - nxpix - Number of pixels in X direction (default: 200) - nypix - Number of pixels in Y direction (default: 200) - outname - Counts map FITS filename (default: cntmap.fits) + Creates a counts map by combining the events of all observations + + The counts map will be a summed map over all energies + + Parameters + ---------- + obs : `~gammalib.GObservations` + Observation container without events + proj : str, optional + Projection for binned simulation + coord : str, optional + Coordinate system for binned simulation + xval : float, optional + Counts map centre of X direction (deg) + yval : float, optional + Counts map centre of Y direction (deg) + binsz : float, optional + Pixel size for binned simulation (deg/pixel) + nxpix : int, optional + Number of pixels in X direction + nypix : int, optional + Number of pixels in Y direction + outname : str, optional + Counts map file name + + Returns + ------- + map : `~gammalib.GSkyMap` + Counts map """ # Allocate counts map map = gammalib.GSkyMap(proj, coord, xval, yval, -binsz, binsz, nxpix, nypix, 1) @@ -276,26 +345,42 @@ def cntmap(obs, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \ # ================ # # Create model map # # ================ # -def modmap(obs, eref=0.1, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \ +def modmap(obs, eref=0.1, proj='TAN', coord='GAL', xval=0.0, yval=0.0, \ binsz=0.05, nxpix=200, nypix=200, \ - outname="modmap.fits"): + outname='modmap.fits'): """ - Make model map for a given reference energy by combining all observations. + Make model map for a given reference energy by combining all observations + The model map will be evaluated for a given reference energy 'eref' and will be given in units of [counts/(sr MeV s)]. - Parameters: - obs - Observation container - Keywords: - eref - Reference energy for which model is created [TeV] (default: 0.1) - proj - Projection type (e.g. TAN, CAR, STG, ...) (default: TAN) - coord - Coordinate type (GAL, CEL) (default: GAL) - xval - Reference longitude value [deg] (default: 0.0) - yval - Reference latitude value [deg] (default: 0.0) - binsz - Pixel size [deg/pixel] (default: 0.05) - nxpix - Number of pixels in X direction (default: 200) - nypix - Number of pixels in Y direction (default: 200) - outname - Model map FITS filename (default: modmap.fits) + Parameters + ---------- + obs : `~gammalib.GObservations` + Observation container without events + eref : float, optional + Reference energy for which model is created (TeV) + proj : str, optional + Projection for binned simulation + coord : str, optional + Coordinate system for binned simulation + xval : float, optional + Counts map centre of X direction (deg) + yval : float, optional + Counts map centre of Y direction (deg) + binsz : float, optional + Pixel size for binned simulation (deg/pixel) + nxpix : int, optional + Number of pixels in X direction + nypix : int, optional + Number of pixels in Y direction + outname : str, optional + Model map file name + + Returns + ------- + map : `~gammalib.GSkyMap` + Model map """ # Allocate model map map = gammalib.GSkyMap(proj, coord, xval, yval, -binsz, binsz, nxpix, nypix, 1) @@ -338,53 +423,46 @@ def modmap(obs, eref=0.1, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \ return map -# ===================================== # -# Set one CTA observation (old version) # -# ===================================== # -def set(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \ - emin=0.1, emax=100.0, rad=5.0, \ - irf="South_50h", caldb="prod2"): - """ - Obsolete function, use set_obs instead. - """ - # Print warning - print("Warning: obsutils.set is obsolete, use obsutils.set_obs instead.") - - # Call new function - obs_cta = set_obs(pntdir, tstart=tstart, duration=duration, deadc=deadc, \ - emin=emin, emax=emax, rad=rad, \ - irf=irf, caldb=caldb) - - # Return CTA observation - return obs_cta - - # ======================= # # Set one CTA observation # # ======================= # def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \ emin=0.1, emax=100.0, rad=5.0, \ - irf="South_50h", caldb="prod2", id="000000"): + irf='South_50h', caldb='prod2', id='000000'): """ - Set a single CTA observation. + Set a single CTA observation The function sets a single CTA observation containing an empty CTA event list. By looping over this function you can add CTA observations to the observation container. - Args: - pntdir: Pointing direction [GSkyDir] - - Kwargs: - tstart: Start time (seconds) (default: 0.0) - duration: Duration of observation (seconds) (default: 1800.0) - deadc: Deadtime correction factor (default: 0.95) - emin: Minimum event energy (TeV) (default: 0.1) - emax: Maximum event energy (TeV) (default: 100.0) - rad: ROI radius used for analysis (deg) (default: 5.0) - irf: Instrument response function (default: "South_50h") - caldb: Calibration database path (default: "prod2") - id: Run identifier (default: "000000") + Parameters + ---------- + pntdir : `~gammalib.GSkyDir` + Pointing direction + tstart : float, optional + Start time (seconds) + duration : float, optional + Duration of observation (seconds) + deadc : float, optional + Deadtime correction factor + emin : float, optional + Minimum event energy (TeV) + emax : float, optional + Maximum event energy (TeV) + rad : float, optional + ROI radius used for analysis (deg) + irf : str, optional + Instrument response function + caldb : str, optional + Calibration database path + id : str, optional + Observation identifier + + Returns + ------- + obs : `~gammalib.GCTAObservation` + CTA observation """ # Allocate CTA observation obs_cta = gammalib.GCTAObservation() @@ -394,7 +472,7 @@ def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \ if (gammalib.dir_exists(caldb)): db.rootdir(caldb) else: - db.open("cta", caldb) + db.open('cta', caldb) # Set pointing direction pnt = gammalib.GCTAPointing() @@ -413,8 +491,8 @@ def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \ gti.append(gammalib.GTime(tstart), gammalib.GTime(tstart+duration)) # Set energy boundaries - ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, "TeV"), - gammalib.GEnergy(emax, "TeV")) + ebounds = gammalib.GEbounds(gammalib.GEnergy(emin, 'TeV'), + gammalib.GEnergy(emax, 'TeV')) # Allocate event list events = gammalib.GCTAEventList() @@ -440,10 +518,11 @@ def set_obs(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \ # Set list of CTA observations # # ============================ # def set_obs_list(obsdeflist, tstart=0.0, duration=1800.0, deadc=0.95, \ - emin=0.1, emax=100.0, rad=5.0, \ - irf="South_50h", caldb="prod2"): + emin=0.1, emax=100.0, rad=5.0, \ + irf='South_50h', caldb='prod2'): """ - Returns an observation container filled with a list of CTA observations. + Returns an observation container filled with a list of CTA observations + The list is defined by the obsdeflist parameter which is a dictionnary containing the mandatory keywords 'ra' and 'dec' that specify the pointing direction for a given observation. Optional keyword give control @@ -451,27 +530,31 @@ def set_obs_list(obsdeflist, tstart=0.0, duration=1800.0, deadc=0.95, \ energy range, etc. If an optional keyword is not specified, the function keyword is used instead. - Parameters: - obsdeflist - Observation definition list [{'ra': x.xx, 'dec': x.xx}] - The directory can take the following optional keywords: - - duration: Duration of CTA observation [seconds] - - deadc: Deadtime correction factor - - emin: Minimum event energy [TeV] - - emax: Maximum event energy [TeV] - - rad: ROI radius used for analysis [deg] - - irf: Instrument response function - - caldb: Calibration database path - Optional keywords overwrite keywords specified in the - function call. - Keywords: - tstart - Start time [seconds] (default: 0.0) - duration - Duration of one CTA observation [seconds] (default: 1800.0) - deadc - Deadtime correction factor (default: 0.95) - emin - Minimum event energy [TeV] (default: 0.1) - emax - Maximum event energy [TeV] (default: 100.0) - rad - ROI radius used for analysis [deg] (default: 5.0) - irf - Instrument response function (default: cta_dummy_irf) - caldb - Calibration database path (default: "dummy") + Parameters + ---------- + obsdeflist : list of dict + Observation definition list + tstart : float, optional + Start time (seconds) + duration : float, optional + Duration of observation (seconds) + deadc : float, optional + Deadtime correction factor + emin : float, optional + Minimum event energy (TeV) + emax : float, optional + Maximum event energy (TeV) + rad : float, optional + ROI radius used for analysis (deg) + irf : str, optional + Instrument response function + caldb : str, optional + Calibration database path + + Returns + ------- + obs : `~gammalib.GObservations` + Observation container filled with CTA observation """ # Initialise empty observation container obs = gammalib.GObservations() @@ -530,7 +613,7 @@ def set_obs_list(obsdeflist, tstart=0.0, duration=1800.0, deadc=0.95, \ obs_irf = irf # Generate identifier string - id = "%6.6d" % obs_id + id = '%6.6d' % obs_id # Set CTA observation obs_cta = set_obs(pntdir, tstart=obs_start, duration=obs_duration, \ @@ -553,26 +636,38 @@ def set_obs_list(obsdeflist, tstart=0.0, duration=1800.0, deadc=0.95, \ # ============================ # def set_obs_patterns(pattern, ra=83.6331, dec=22.0145, offset=1.5): """ - Sets a number of standard patterns. - - Parameters: - pattern - Observation pattern. Possible options are: - - "single": single pointing - - "four": four pointings 'offset' around pattern centre + Sets a number of standard patterns + + Parameters + ---------- + pattern : str + Observation pattern ("single", "four") + ra : float, optional + Right Ascension of pattern centre (deg) + dec : float, optional + Declination of pattern centre (deg) + offset : float, optional + Offset from pattern centre (deg) - Keywords: - ra - Right Ascension of pattern centre [deg] (default: 83.6331) - dec - Declination of pattern centre [deg] (default: 22.0145) - offset - Offset from pattern centre [deg] (default: 1.5) + Returns + ------- + obsdeflist : list + Observation definition list """ # Initialise observation definition list obsdeflist = [] - # Add patterns - if pattern == "single": + # If the pattern is a single observation then append the Right Ascension + # and Declination to the observation definition list + if pattern == 'single': obsdef = {'ra': ra, 'dec': dec} obsdeflist.append(obsdef) - elif pattern == "four": + + # ... otherwise, if the pattern is four observations then append four + # observations offset by a certain amount from the pattern centre to the + # observation definition list + elif pattern == 'four': + # Set pattern centre centre = gammalib.GSkyDir() centre.radec_deg(ra, dec) @@ -582,8 +677,232 @@ def set_obs_patterns(pattern, ra=83.6331, dec=22.0145, offset=1.5): pntdir = centre.copy() pntdir.rotate_deg(phi, offset) obsdeflist.append({'ra': pntdir.ra_deg(), 'dec': pntdir.dec_deg()}) + + # ... otherwise we have an unknown pattern else: - print("Warning: Observation pattern '"+str(pattern)+"' not recognized.") + print('Warning: Observation pattern "'+str(pattern)+'" not recognized.') # Return observation definition list return obsdeflist + + +# ======================== # +# Set model for TS fitting # +# ======================== # +def set_ts_model(models, srcname, ra=None, dec=None, fitspat=False, fitspec=False): + """ + Set model for TS fitting + + Parameters + ---------- + models : `~gammalib.GModels` + Input model container + srcname : str + Test source name + ra : float, optional + Right Ascension of test source (deg) + dec : float, optional + Declination of test source (deg) + fitspat : bool, optional + Fit spatial parameter? + fitspec : bool, optional + Fit spectral parameters? + + Returns + ------- + model : `~gammalib.GModels` + Model container for TS fitting + """ + # Create a clone of the input model + outmodels = models.copy() + + # Disable TS computation for all model components (will enable the + # test source later) + for model in outmodels: + model.tscalc(False) + + # Get source model and enable TS computation + model = outmodels[srcname] + model.tscalc(True) + + # If source model has no "Prefactor" parameter then raise an exception + if not model.has_par('Prefactor'): + msg = ('Model "%s" has no parameter "Prefactor". Only spectral ' + 'models with a "Prefactor" parameter are supported.' % + srcname) + raise RuntimeError(msg) + + # Set position of test source + if ra != None and dec != None: + if model.has_par('RA') and model.has_par('DEC'): + model['RA'].value(ra) + model['DEC'].value(dec) + + # Set possible spatial and spectral parameters + spatial = ['RA', 'DEC', 'Sigma', 'Radius', 'Width', 'PA', + 'MinorRadius', 'MajorRadius'] + spectral = ['Index', 'Index1', 'Index2', 'BreakEnergy', 'CutoffEnergy', + 'InverseCutoffEnergy'] + + # Fit or fix spatial parameters + for par in spatial: + if model.has_par(par): + if fitspat: + model[par].free() + else: + model[par].fix() + + # Fit or fix spectral parameters + for par in spectral: + if model.has_par(par): + if fitspec: + model[par].free() + else: + model[par].fix() + + # Return model container + return outmodels + + +# ==================== # +# Get stacked response # +# ==================== # +def get_stacked_response(obs, xref, yref, binsz=0.05, nxpix=200, nypix=200, + emin=0.1, emax=100.0, enumbins=20, edisp=False, + coordsys='GAL', proj='TAN', + log=False, debug=False, chatter=2): + """ + Get stacked response cubes + + Parameters + ---------- + obs : `~gammalib.GObservations` + Observation container + xref : float + Right Ascension or Galactic longitude of response centre (deg) + yref : float + Declination or Galactic latitude of response centre (deg) + binsz : float, optional + Pixel size (deg/pixel) + nxpix : int, optional + Number of pixels in X direction + nypix : int, optional + Number of pixels in Y direction + emin : float, optional + Minimum energy (TeV) + emax : float, optional + Maximum energy (TeV) + enumbins : int, optional + Number of energy bins + edisp : bool, optional + Apply energy dispersion? + coordsys : str, optional + Coordinate system + proj : str, optional + Projection + log : bool, optional + Create log file(s) + debug : bool, optional + Create console dump? + chatter : int, optional + Chatter level + + Returns + ------- + result : dict + Dictionary of response cubes + """ + # Compute spatial binning for point spread function and + # energy dispersion cubes + psf_binsz = 10.0 * binsz + psf_nxpix = min(nxpix // 10, 2) # Make sure result is int + psf_nypix = min(nypix // 10, 2) # Make sure result is int + + # Create exposure cube + expcube = ctools.ctexpcube(obs) + expcube['incube'] = 'NONE' + expcube['usepnt'] = False + expcube['ebinalg'] = 'LOG' + expcube['xref'] = xref + expcube['yref'] = yref + expcube['binsz'] = binsz + expcube['nxpix'] = nxpix + expcube['nypix'] = nypix + expcube['enumbins'] = enumbins + expcube['emin'] = emin + expcube['emax'] = emax + expcube['coordsys'] = coordsys + expcube['proj'] = proj + if log: + expcube.logFileOpen() + expcube.run() + + # Create point spread function cube + psfcube = ctools.ctpsfcube(obs) + psfcube['incube'] = 'NONE' + psfcube['usepnt'] = False + psfcube['ebinalg'] = 'LOG' + psfcube['xref'] = xref + psfcube['yref'] = yref + psfcube['binsz'] = psf_binsz + psfcube['nxpix'] = psf_nxpix + psfcube['nypix'] = psf_nypix + psfcube['enumbins'] = enumbins + psfcube['emin'] = emin + psfcube['emax'] = emax + psfcube['coordsys'] = coordsys + psfcube['proj'] = proj + if log: + psfcube.logFileOpen() + psfcube.run() + + # Create background cube + bkgcube = ctools.ctbkgcube(obs) + bkgcube['incube'] = 'NONE' + bkgcube['usepnt'] = False + bkgcube['ebinalg'] = 'LOG' + bkgcube['xref'] = xref + bkgcube['yref'] = yref + bkgcube['binsz'] = binsz + bkgcube['nxpix'] = nxpix + bkgcube['nypix'] = nypix + bkgcube['enumbins'] = enumbins + bkgcube['emin'] = emin + bkgcube['emax'] = emax + bkgcube['coordsys'] = coordsys + bkgcube['proj'] = proj + if log: + bkgcube.logFileOpen() + bkgcube.run() + + # If energy dispersion is requested then create energy dispersion cube + if edisp: + edispcube = ctools.ctedispcube(obs) + edispcube['incube'] = 'NONE' + edispcube['usepnt'] = False + edispcube['ebinalg'] = 'LOG' + edispcube['xref'] = xref + edispcube['yref'] = yref + edispcube['binsz'] = psf_binsz + edispcube['nxpix'] = psf_nxpix + edispcube['nypix'] = psf_nypix + edispcube['enumbins'] = enumbins + edispcube['emin'] = emin + edispcube['emax'] = emax + edispcube['coordsys'] = coordsys + edispcube['proj'] = proj + if log: + edispcube.logFileOpen() + edispcube.run() + + # Build response dictionary + response = {} + response['expcube'] = expcube.expcube().copy() + response['psfcube'] = psfcube.psfcube().copy() + response['bkgcube'] = bkgcube.bkgcube().copy() + response['models'] = bkgcube.models().copy() + if edisp: + response['edispcube'] = edispcube.expcube().copy() + + # Return response cubes + return response diff --git a/doc/source/users/reference_manual/cstsdist.rst b/doc/source/users/reference_manual/cstsdist.rst index aa705dd..9884871 100644 --- a/doc/source/users/reference_manual/cstsdist.rst +++ b/doc/source/users/reference_manual/cstsdist.rst @@ -75,6 +75,12 @@ General parameters ``dec [real]`` Declination of CTA pointing (J2000, in degrees). +``coordsys [string]`` + Coordinate system (CEL - celestial, GAL - galactic). + +``proj [string]`` + Projection method. + ``emin [real]`` Lower energy limit (in TeV). diff --git a/test/test_cstsdist.py b/test/test_cstsdist.py index a0e9a81..a8a47cd 100755 --- a/test/test_cstsdist.py +++ b/test/test_cstsdist.py @@ -76,7 +76,7 @@ class Test(test): ' enumbins=0 tmax=1800.0 rad=5.0 npix=200 npix=200'+ \ ' binsz=0.05'+ \ ' outfile="cstsdist_cmd1.dat"'+ \ - ' logfile="cstsdist_cmd1.log" chatter=1' + ' logfile="cstsdist_cmd1.log" chatter=2' # Check if execution of wrong command fails self.test_assert(self._execute('command_that_does_not_exist') != 0, @@ -117,7 +117,6 @@ class Test(test): tsdist['srcname'] = 'Crab' tsdist['caldb'] = self._caldb tsdist['irf'] = self._irf - tsdist['outfile'] = 'cstsdist_py1.dat' tsdist['ntrials'] = 1 tsdist['ra'] = 83.63 tsdist['dec'] = 22.01 @@ -126,10 +125,9 @@ class Test(test): tsdist['enumbins'] = 0 tsdist['tmax'] = 1800.0 tsdist['rad'] = 5.0 - tsdist['npix'] = 200 - tsdist['binsz'] = 0.05 + tsdist['outfile'] = 'cstsdist_py1.dat' tsdist['logfile'] = 'cstsdist_py1.log' - tsdist['chatter'] = 2 + tsdist['chatter'] = 3 # Run cstsdist script tsdist.logFileOpen() # Make sure we get a log file @@ -139,20 +137,56 @@ class Test(test): # Check pull distribution file self._check_result_file('cstsdist_py1.dat') + # Set-up binned cstsdist in four observation wobble mode + tsdist = cscripts.cstsdist() + 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['coordsys'] = 'CEL' + tsdist['proj'] = 'TAN' + tsdist['emin'] = 0.1 + tsdist['emax'] = 100.0 + tsdist['enumbins'] = 10 + tsdist['tmax'] = 1800.0 + tsdist['rad'] = 5.0 + tsdist['npix'] = 50 + tsdist['binsz'] = 0.1 + tsdist['pattern'] = 'four' + tsdist['outfile'] = 'cstsdist_py2.dat' + tsdist['logfile'] = 'cstsdist_py2.log' + tsdist['chatter'] = 4 + + # Execute cstsdist script + tsdist.execute() + + # Check pull distribution file + self._check_result_file('cstsdist_py2.dat', ncols=9) + # Return return # Check result file - def _check_result_file(self, filename): + def _check_result_file(self, filename, ncols=11): """ Check result file + + Parameters + ---------- + filename : str + Name of result file + ncols : int, optional + Required number of columns """ # Open result file as CSV file results = gammalib.GCsv(filename, ',') # Check dimensions - self.test_value(results.nrows(), 2, 'Check for 2 rows in TS file') - self.test_value(results.ncols(), 14, 'Check for 14 columns in TS file') + self.test_value(results.nrows(), 2, 'Check rows in TS file') + self.test_value(results.ncols(), ncols, 'Check columns in TS file') # Return return