diff --git a/.gitignore b/.gitignore index 23f8f1a..98079f5 100644 --- a/.gitignore +++ b/.gitignore @@ -71,6 +71,8 @@ src/ctcubemask/ctcubemask src/ctbutterfly/ctbutterfly src/ctulimit/ctulimit src/cterror/cterror +src/ctprob/ctprob +src/ctphase/ctphase test/pfiles/ test/caldb/ test/prod3 diff --git a/AUTHORS b/AUTHORS index 1f8f45a..4b6cdc1 100644 --- a/AUTHORS +++ b/AUTHORS @@ -2,6 +2,7 @@ Johan Bregeon (johan.bregeon@lupm.in2p3.fr) Rolf Buehler (rolf.buehler@desy.de) Jean-Baptiste Cayrou (cayrou@etud.insa-toulouse.fr) Christoph Deil (Christoph.Deil@mpi-hd.mpg.de) +Leonardo Di Venere (leonardo.divenere@ba.infn.it) Florent Forest (florent.forest@isae.fr) Lucie Gerard (lgerard@wgs18.zeuthen.desy.de) Maria Haupt (maria.haupt@desy.de) diff --git a/ChangeLog b/ChangeLog index ebe243c..7f305e8 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,13 @@ -2017-04-05 +2017-04-11 * Version 1.3.0 released ======================== + Add ctprob tool (#1991) + ctlike stores covariance matrix in either CSV or FITS format (#1863) + Use correct ROOT histograms in csroot2caldb + Implement regression scheme in cssens and remove num_avg parameter (#1989) + Implement Gaussian error propagation in ctbutterfly (#1862) Copy cscripts into bin directory instead of creating a symbolic links (#1987) If code originates from Git, always compile the Python bindings (#1985) Raise exception instead of exiting with return code in Python tests (#1983) diff --git a/NEWS b/NEWS index 60abd03..d41c52c 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in ctools 1.3.0.dev1 -5 April 2017 +11 April 2017 Introduction @@ -57,7 +57,10 @@ None ctlike - Maximum likelihood fitting ----------------------------------- -None +The tool now stores the covariance matrix in either FITS or CSV format, +depending on the file type specified in the file name. For file names +terminating with ".fits" or ".fit" a FITS file will be written, otherwise +a CSV file will be written (#1863). cttsmap - Generation of Test Statistic map @@ -113,7 +116,11 @@ None ctbutterfly - Generate butterfly diagrams ----------------------------------------- -None +The tool now uses Gaussian error propagation as the standard method for the +computation of the butterfly diagram. The envelope method that existed before +is still available for power law spectral models as an option. The "method" +parameter has been added to allow the selection of the computation method +(#1862). ctulimit - Compute upper limit @@ -126,6 +133,15 @@ cterror - Likelihood profile error computation None +ctprob - Event probability computation +-------------------------------------- +The ctprob tool has been added that computates for each event the probability +that the event originates from a given model component. The tool takes on input +either an event list or an observation definition XML file. It appends to each +event list a column for each model component that will contain the probability +that the event originates from that component (#1991). + + cscripts -------- The csebins script has been added that generates an energy boundaries @@ -133,6 +149,15 @@ extension that can be used for a stacked analysis. The energy boundaries are set as function of the variation in the effective area and background rates of the Instrument Response Function (#1976). +A regression scheme now replaces the averaging scheme in cssens as this +allows more accurately the prediction of the flux value at the threshold +Test Statistic. The num_avg parameter has been removed (#1989). + +The names of some ROOT histograms has been changed in the Monte Carlo data +and the "EffectiveAreaEtrue_offaxis" needed to be changed to +"EffectiveAreaEtrueNoTheta2cut_offaxis" and "EffectiveAreaEtrue" to +"EffectiveAreaEtrueNoTheta2cut" in csroot2caldb. + Examples -------- diff --git a/cscripts/csroot2caldb.py b/cscripts/csroot2caldb.py index 1a36e74..28c3df4 100755 --- a/cscripts/csroot2caldb.py +++ b/cscripts/csroot2caldb.py @@ -536,7 +536,8 @@ class csroot2caldb(ctools.cscript): Translate ROOT to CALDB effective area extension The following ROOT histograms are used: - - EffectiveAreaEtrue_offaxis -> EFFAREA + - EffectiveAreaEtrueNoTheta2cut_offaxis -> EFFAREA + - EffectiveAreaEtrueNoTheta2cut (for normalization) Parameters ---------- @@ -551,13 +552,13 @@ class csroot2caldb(ctools.cscript): self._log_header1(gammalib.TERSE, 'Generate effective area extension') # Get relevant ROOT histograms - etrue = tfile.Get('EffectiveAreaEtrue_offaxis') + etrue = tfile.Get('EffectiveAreaEtrueNoTheta2cut_offaxis') # If requested then normalize the 2D histogram on the on-axis 1D # histogram. This assures that the 2D histogram has the same on-axis # effective area dependence as the 1D histogram. if self['norm1d'].boolean(): - etrue_1D = tfile.Get('EffectiveAreaEtrue') + etrue_1D = tfile.Get('EffectiveAreaEtrueNoTheta2cut') self._renorm_onaxis(etrue, etrue_1D) # If requested then rebin the effective area Etrue histogram. This diff --git a/cscripts/cssens.par b/cscripts/cssens.par index e79fbb4..ec14a9f 100644 --- a/cscripts/cssens.par +++ b/cscripts/cssens.par @@ -20,7 +20,7 @@ inobs, f, h, NONE,,, "Input event list, counts cube, or observation definition 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" +irf, s, a, South_50h,,, "Instrument response function" edisp, b, h, no,,, "Apply energy dispersion?" deadc, r, h, 0.98,0,1, "Deadtime correction factor" outfile, f, h, sensitivity.dat,,, "Output sensitivity ASCII file" @@ -31,16 +31,15 @@ outfile, f, h, sensitivity.dat,,, "Output sensitivity ASCII file" offset, r, h, 0.0,,, "Source offset angle (deg)" duration, r, a, 180000.0,,, "Effective exposure time (s)" rad, r, a, 5.0,,, "Radius of ROI (deg)" -emin, r, a, 0.020,,, "Lower energy limit (TeV)" -emax, r, a, 200.0,,, "Upper energy limit (TeV)" -bins, i, a, 21,,, "Number of energy bins for differential sensitivity computation" +emin, r, a, 0.03163,,, "Lower energy limit (TeV)" +emax, r, a, 125.8925,,, "Upper energy limit (TeV)" +bins, i, a, 18,,, "Number of energy bins for differential sensitivity computation" enumbins, i, h, 0,,, "Number of energy bins for binned analysis (0=unbinned)" npix, i, h, 200,,, "Number of pixels for binned analysis" binsz, r, h, 0.05,,, "Pixel size for binned analysis (deg/pixel)" type, s, h, Differential,Differential|Integral,, "Sensitivity type" sigma, r, h, 5.0,,, "Significance threshold (Gaussian sigma)" max_iter, i, h, 50,,, "Maximum number of iterations" -num_avg, i, h, 3,,, "Number of iterations for sliding average" # # Standard parameters @@ -49,6 +48,6 @@ chatter, i, h, 2,0,4, "Chattiness of output" clobber, b, h, yes,,, "Overwrite existing output files with new output files?" debug, b, h, no,,, "Debugging mode activated" mode, s, h, ql,,, "Mode of automatic parameters" -logfile, f, h, cspull.log,,, "Log filename" +logfile, f, h, cssens.log,,, "Log filename" diff --git a/cscripts/cssens.py b/cscripts/cssens.py index 5d3a8b0..f16590c 100755 --- a/cscripts/cssens.py +++ b/cscripts/cssens.py @@ -1,8 +1,8 @@ #! /usr/bin/env python # ========================================================================== -# Computes the array sensitivity using a test source +# Computes the array sensitivity using the Test Statistic for a test source # -# Copyright (C) 2011-2016 Juergen Knoedlseder +# Copyright (C) 2011-2017 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 @@ -40,7 +40,7 @@ class cssens(ctools.cscript): and the flux level is determined that leads to a particular significance. The significance is determined using the Test Statistic, defined as twice - the likelihood difference between fitting with and without the test source. + the likelihood difference between fitting with and without the test source. """ # Constructor @@ -50,31 +50,20 @@ class cssens(ctools.cscript): """ # Set name self._name = 'cssens' - self._version = '1.2.0' + self._version = '1.3.0' # Initialise class members - self._obs = gammalib.GObservations() self._ebounds = gammalib.GEbounds() self._obs_ebounds = [] self._srcname = '' - self._outfile = gammalib.GFilename() self._ra = None self._dec = None - self._edisp = False - self._emin = 0.020 - self._emax = 200.0 - self._bins = 21 - self._enumbins = 0 - self._npix = 200 - self._binsz = 0.05 - self._type = 'Differential' - self._ts_thres = 25.0 - self._max_iter = 50 - self._num_avg = 3 self._log_clients = False - - # Initialise application by calling the appropriate base class - # constructor + + # 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 @@ -87,8 +76,8 @@ class cssens(ctools.cscript): Get user parameters from parfile """ # Set observation if not done before - if self._obs == None or self._obs.size() == 0: - self._obs = self._set_obs() + if self._obs.size() == 0: + self._obs = self._set_obs(self['emin'].real(), self['emax'].real()) # Set models if we have none if self._obs.models().size() == 0: @@ -98,31 +87,28 @@ class cssens(ctools.cscript): self._srcname = self['srcname'].string() # Read further parameters - self._outfile = self['outfile'].filename() - self._emin = self['emin'].real() - self._emax = self['emax'].real() - self._bins = self['bins'].integer() - - # Read parameters for binned if requested - self._enumbins = self['enumbins'].integer() - if not self._enumbins == 0: - self._npix = self['npix'].integer() - self._binsz = self['binsz'].real() - - # Read remaining parameters - self._edisp = self['edisp'].boolean() - self._ts_thres = self['sigma'].real() * self['sigma'].real() - self._max_iter = self['max_iter'].integer() - self._num_avg = self['num_avg'].integer() - self._type = self['type'].string() - - # Query remaining parameters + emin = self['emin'].real() + emax = self['emax'].real() + bins = self['bins'].integer() + + # Query parameters for binned if requested + enumbins = self['enumbins'].integer() + if not enumbins == 0: + self['npix'].integer() + self['binsz'].real() + + # Query input parameters + self['sigma'].real() + self['max_iter'].integer() + self['type'].string() + self['outfile'].filename() + self['edisp'].boolean() self['debug'].boolean() # Derive some parameters - self._ebounds = gammalib.GEbounds(self._bins, - gammalib.GEnergy(self._emin, 'TeV'), - gammalib.GEnergy(self._emax, 'TeV')) + self._ebounds = gammalib.GEbounds(bins, + gammalib.GEnergy(emin, 'TeV'), + gammalib.GEnergy(emax, 'TeV')) # Write input parameters into logger self._log_parameters(gammalib.TERSE) @@ -130,20 +116,20 @@ class cssens(ctools.cscript): # Return return - def _set_obs(self, lpnt=0.0, bpnt=0.0, emin=0.1, emax=100.0): + def _set_obs(self, emin, emax, lpnt=0.0, bpnt=0.0): """ Set an observation container Parameters ---------- + emin : float + Minimum energy (TeV) + emax : float + Maximum energy (TeV) lpnt : float, optional Galactic longitude of pointing (deg) bpnt : float, optional Galactic latitude of pointing (deg) - emin : float, optional - Minimum energy (TeV) - emax : float, optional - Maximum energy (TeV) Returns ------- @@ -205,23 +191,23 @@ class cssens(ctools.cscript): # Get energy boundaries of the observation obs_ebounds = self._obs_ebounds[i] - + # Get minimum and maximum energy of the observation obs_emin = obs_ebounds.emin() obs_emax = obs_ebounds.emax() - + # If [emin,emax] is fully contained in the observation energy range # the use [emin,emax] as energy boundaries if obs_emin <= emin and obs_emax >= emax: ebounds = gammalib.GEbounds(emin, emax) - + # ... otherwise, if [emin,emax] is completely outside the # observation energy range then set the energy boundaries to the # zero-width interval [0,0] elif emax < obs_emin or emin > obs_emax: e0 = gammalib.GEnergy(0.0, 'TeV') ebounds = gammalib.GEbounds(e0, e0) - + # ... otherwise, if [emin,emax] overlaps partially with the # observation energy range then set the energy boundaries to the # overlapping part @@ -229,10 +215,10 @@ class cssens(ctools.cscript): # Set overlapping energy range set_emin = max(emin, obs_emin) set_emax = min(emax, obs_emax) - + # Set energy boundaries ebounds = gammalib.GEbounds(set_emin, set_emax) - + # Set the energy boundaries as the boundaries of the observation obs.events().ebounds(ebounds) @@ -286,13 +272,27 @@ class cssens(ctools.cscript): # Set TeV->erg conversion factor tev2erg = 1.6021764 + # Set parameters + ts_thres = self['sigma'].real() * self['sigma'].real() + max_iter = self['max_iter'].integer() + enumbins = self['enumbins'].integer() + if not enumbins == 0: + npix = self['npix'].integer() + binsz = self['binsz'].real() + else: + npix = 200 + binsz = 0.05 + + # Set flux ratio precision required for convergence to 5% + ratio_precision = 0.05 + # Set energy boundaries self._set_obs_ebounds(emin, emax) - # Determine energy boundaries from first observation in the container - loge = math.log10(math.sqrt(emin.TeV()*emax.TeV())) - e_mean = math.pow(10.0, loge) - erg_mean = e_mean * tev2erg + # Determine mean energy for energy boundary + e_mean = math.sqrt(emin.TeV()*emax.TeV()) + loge = math.log10(e_mean) + erg_mean = e_mean * tev2erg # Compute Crab unit. This is the factor with which the Prefactor needs # to be multiplied to get 1 Crab. @@ -300,20 +300,27 @@ class cssens(ctools.cscript): src_flux = test_model[self._srcname].spectral().flux(emin, emax) crab_unit = crab_flux/src_flux - # Write header and some introductory information into logger + # Initialise regression coefficient + regcoeff = 0.0 + + # Write header for energy bin + self._log_string(gammalib.TERSE, '') self._log_header2(gammalib.TERSE, 'Energies: '+str(emin)+' - '+str(emax)) + + # Write initial parameters + self._log_header3(gammalib.TERSE, 'Initial parameters') self._log_value(gammalib.TERSE, 'Crab flux', str(crab_flux)+' ph/cm2/s') - self._log_value(gammalib.TERSE, 'Source model flux', str(src_flux)+ - ' ph/cm2/s') + self._log_value(gammalib.TERSE, 'Source model flux', str(src_flux)+' ph/cm2/s') self._log_value(gammalib.TERSE, 'Crab unit factor', crab_unit) # Initialise loop - crab_flux_value = [] - photon_flux_value = [] - energy_flux_value = [] - sensitivity_value = [] - iterations = 0 - test_crab_flux = 0.1 # Initial test flux in Crab units (100 mCrab) + results = [] + iterations = 0 + test_crab_flux = 0.1 # Initial test flux in Crab units (100 mCrab) + + # Write header for iterations for terse chatter level + if self._logTerse(): + self._log_header3(gammalib.TERSE, 'Iterations') # Loop until we break while True: @@ -335,11 +342,11 @@ class cssens(ctools.cscript): # Simulate events for the models. "sim" holds an observation # container with observations containing the simulated events. - sim = obsutils.sim(self._obs, nbins=self._enumbins, seed=iterations, - binsz=self._binsz, npix=self._npix, + sim = obsutils.sim(self._obs, nbins=enumbins, seed=iterations, + binsz=binsz, npix=npix, log=self._log_clients, debug=self['debug'].boolean(), - edisp=self._edisp) + edisp=self['edisp'].boolean()) # Determine number of events in simulation by summing the events # over all observations in the observation container @@ -354,7 +361,7 @@ class cssens(ctools.cscript): # Fit test source to the simulated events in the observation # container fit = ctools.ctlike(sim) - fit['edisp'] = edisp=self._edisp + fit['edisp'] = self['edisp'].boolean() fit['debug'] = self['debug'].boolean() fit['chatter'] = self['chatter'].integer() fit.run() @@ -383,29 +390,24 @@ class cssens(ctools.cscript): # Assess quality based on a comparison between Npred and Nevents quality = npred - nevents - # Write test source fit results into logger - self._log_header3(gammalib.EXPLICIT, 'Test source model fit') - self._log_value(gammalib.EXPLICIT, 'Test statistics', ts) - self._log_value(gammalib.EXPLICIT, 'log likelihood', logL) - self._log_value(gammalib.EXPLICIT, 'Number of predicted events', quality) - self._log_value(gammalib.EXPLICIT, 'Fit quality', npred) - for model in models: - self._log_value(gammalib.EXPLICIT, 'Model', model.name()) - for par in model: - self._log_string(gammalib.EXPLICIT, str(par)) - - # If TS was non-positive then increase the test flux and start - # over + # Write fit results into logger + name = 'Iteration %d' % iterations + value = ('TS=%10.4f Sim=%9.4f mCrab Fit=%9.4f mCrab ' + 'Sens=%e erg/cm2/s' % + (ts, test_crab_flux*1000.0, crab_flux*1000.0, sensitivity)) + self._log_value(gammalib.TERSE, name, value) + + # If TS was non-positive then increase the test flux and start over if ts <= 0.0: # If the number of iterations was exceeded then stop - if (iterations >= self._max_iter): + if (iterations >= max_iter): self._log_string(gammalib.TERSE, - ' Test ended after %d iterations.' % self._max_iter) + ' Test ended after %d iterations.' % max_iter) break - # Increase test flux by a factor of 2 - test_crab_flux = test_crab_flux * 2.0 + # Increase test flux + test_crab_flux *= 3.0 # Signal start we start over self._log_string(gammalib.EXPLICIT, @@ -414,95 +416,64 @@ class cssens(ctools.cscript): # ... and start over continue - # Compute flux correction factor based on average TS. We - # assume here that the significance is proportional to the source - # flux, and that the significance scales as the square root of - # the Test Statistic - correct = math.sqrt(self._ts_thres/ts) + # Append result entry to result list + result = {'ts': ts, 'crab_flux': crab_flux, + 'photon_flux': photon_flux, + 'energy_flux': energy_flux} + results.append(result) + + # Predict Crab flux at threshold TS using a linear regression of + # the log(TS) and log(crab_flux) values that have so far been + # computed. If not enough results are available than use a simple + # TS scaling relation. + if len(results) > 1: + pred_crab_flux, regcoeff = self._predict_flux(results, ts_thres) + correct = pred_crab_flux / crab_flux + else: + correct = math.sqrt(ts_thres/ts) # Compute extrapolated fluxes based on the flux correction factor crab_flux = correct * crab_flux photon_flux = correct * photon_flux energy_flux = correct * energy_flux sensitivity = correct * sensitivity - crab_flux_value.append(crab_flux) - photon_flux_value.append(photon_flux) - energy_flux_value.append(energy_flux) - sensitivity_value.append(sensitivity) - - # Write fit results - if self._logExplicit(): - self._log_value(gammalib.EXPLICIT, 'Photon flux', - str(photon_flux)+' ph/cm2/s') - self._log_value(gammalib.EXPLICIT, 'Energy flux', - str(energy_flux)+' erg/cm2/s') - self._log_value(gammalib.EXPLICIT, 'Crab flux', - str(crab_flux*1000.0)+' mCrab') - self._log_value(gammalib.EXPLICIT, 'Differential sensitivity', - str(sensitivity)+' erg/cm2/s') - for model in models: - self._log_value(gammalib.EXPLICIT, 'Model', model.name()) - for par in model: - self._log_string(gammalib.EXPLICIT, str(par)) - elif self._logTerse(): - name = 'Iteration %d' % iterations - value = ('TS=%.4f corr=%e %e ph/cm2/s = %e erg/cm2/s = %.2f ' - 'mCrab = %e erg/cm2/s' % (ts, correct, photon_flux, - energy_flux, crab_flux*1000.0, sensitivity)) - self._log_value(gammalib.TERSE, name, value) - - # Compute sliding average of extrapolated fitted prefactor, - # photon and energy flux. This damps out fluctuations and - # improves convergence. - crab_flux = 0.0 - photon_flux = 0.0 - energy_flux = 0.0 - sensitivity = 0.0 - num = 0.0 - for k in range(self._num_avg): - inx = len(crab_flux_value) - k - 1 - if inx >= 0: - crab_flux += crab_flux_value[inx] - photon_flux += photon_flux_value[inx] - energy_flux += energy_flux_value[inx] - sensitivity += sensitivity_value[inx] - num += 1.0 - crab_flux /= num - photon_flux /= num - energy_flux /= num - sensitivity /= num - - # Compare average flux to last average - if iterations > self._num_avg: + + # If we have at least 3 results then check if the flux determination + # at the TS threshold has converged + if len(results) > 3: if test_crab_flux > 0: + + # Compute fractional change in the Crab flux between two + # iterations ratio = crab_flux/test_crab_flux - # We have 2 convergence criteria: - # 1. The average flux does not change - # 2. The flux correction factor is small - if ratio >= 0.99 and ratio <= 1.01 and \ - correct >= 0.9 and correct <= 1.1: - if self._logTerse(): - self._log(' Converged ('+str(ratio)+')\n') + # If fractional change is smaller than the required position + # the iterations are stopped + if ratio > 1.0-ratio_precision and \ + ratio < 1.0+ratio_precision: + value = ('TS=%10.4f Sim=%9.4f mCrab ' + ' Sens=%e erg/cm2/s' % + (ts, crab_flux*1000.0, sensitivity)) + self._log_value(gammalib.TERSE, 'Converged result', value) + self._log_value(gammalib.TERSE, 'Converged flux ratio', ratio) + self._log_value(gammalib.TERSE, 'Regression coefficient', + regcoeff) break else: - if self._logTerse(): - self._log(' Flux is zero.\n') + self._log_value(gammalib.TERSE, 'Not converged', 'Flux is zero') break - # Use average for next iteration + # Set test flux for next iteration test_crab_flux = crab_flux # Exit loop if number of trials exhausted - if (iterations >= self._max_iter): - if self._logTerse(): - self._log(' Test ended after '+str(self._max_iter)+ - ' iterations.\n') + if (iterations >= max_iter): + self._log_string(gammalib.TERSE, + ' Test ended after %d iterations.' % max_iter) break # Write fit results into logger self._log_header3(gammalib.TERSE, 'Fit results') - self._log_value(gammalib.TERSE, 'Test statistics', ts) self._log_value(gammalib.TERSE, 'Photon flux', str(photon_flux)+' ph/cm2/s') self._log_value(gammalib.TERSE, 'Energy flux', @@ -523,16 +494,84 @@ class cssens(ctools.cscript): # Restore energy boundaries of observation container for i, obs in enumerate(self._obs): obs.events().ebounds(self._obs_ebounds[i]) - + # Store result result = {'loge': loge, 'emin': emin.TeV(), 'emax': emax.TeV(), \ 'crab_flux': crab_flux, 'photon_flux': photon_flux, \ 'energy_flux': energy_flux, \ - 'sensitivity': sensitivity} + 'sensitivity': sensitivity, 'regcoeff': regcoeff, \ + 'nevents': nevents, 'npred': npred} # Return result return result + def _predict_flux(self, results, ts): + """ + Predict Crab flux for a given TS value + + The Crab flux at a given Test Statistic value is predicted by doing a + linear regression of the log(TS) and log(crab_flux) values in a results + list. + + See https://en.wikipedia.org/wiki/Simple_linear_regression + + Parameters + ---------- + results : list of dict + List of results + ts : float + Test Statistic value + + Returns + ------- + crab_flux_prediction, regcoeff : tuple of float + Predicted Crab flux and regression coefficient + """ + # Compute means and regression coefficient + mean_x = 0.0 + mean_y = 0.0 + mean_xy = 0.0 + mean_xx = 0.0 + mean_yy = 0.0 + for result in results: + x = math.log(float(result['ts'])) + y = math.log(float(result['crab_flux'])) + mean_x += x + mean_y += y + mean_xy += x * y + mean_xx += x * x + mean_yy += y * y + norm = 1.0 / float(len(results)) + mean_x *= norm + mean_y *= norm + mean_xy *= norm + mean_xx *= norm + mean_yy *= norm + rxy = (mean_xy - mean_x * mean_y) / \ + math.sqrt((mean_xx - mean_x*mean_x) * + (mean_yy - mean_y*mean_y)) + regcoeff = rxy*rxy + + # Compute regression line slope + beta_nom = 0.0 + beta_denom = 0.0 + for result in results: + x = math.log(float(result['ts'])) + y = math.log(float(result['crab_flux'])) + beta_nom += (x - mean_x) * (y - mean_y) + beta_denom += (x - mean_x) * (x - mean_x) + beta = beta_nom / beta_denom + + # Compute regression line offset + alpha = mean_y - beta * mean_x + + # Predict Crab flux at TS threshold + log_ts_thres = math.log(ts) + crab_flux_prediction = math.exp(alpha + beta * log_ts_thres) + + # Return + return (crab_flux_prediction, regcoeff) + # Public methods def run(self): @@ -545,15 +584,16 @@ class cssens(ctools.cscript): # Get parameters self._get_parameters() - + # Loop over observations and store a deep copy of the energy # boundaries for later use for obs in self._obs: self._obs_ebounds.append(obs.events().ebounds().copy()) - + # Initialise script colnames = ['loge', 'emin', 'emax', 'crab_flux', 'photon_flux', - 'energy_flux', 'sensitivity'] + 'energy_flux', 'sensitivity', 'regcoeff', 'nevents', + 'npred'] results = [] # Set test source model for this observation @@ -566,30 +606,35 @@ class cssens(ctools.cscript): # Write models into logger self._log_models(gammalib.NORMAL, models, 'Model') + # Get sensitivity type + sensitivity_type = self['type'].string() + # Write header self._log_header1(gammalib.TERSE, 'Sensitivity determination') - self._log_value(gammalib.TERSE, 'Type', str(self._type)) + self._log_value(gammalib.TERSE, 'Type', sensitivity_type) # Loop over energy bins for ieng in range(self._ebounds.size()): # Set energies - if self._type == 'Differential': + if sensitivity_type == 'Differential': emin = self._ebounds.emin(ieng) emax = self._ebounds.emax(ieng) - elif self._type == 'Integral': + elif sensitivity_type == 'Integral': emin = self._ebounds.emin(ieng) emax = self._ebounds.emax() else: msg = ('Invalid sensitivity type "%s" encountered. Either ' - 'specify "Differential" or "Integral".' % self._type) + 'specify "Differential" or "Integral".' % + sensitivity_type) raise RuntimeError(msg) # Determine sensitivity result = self._get_sensitivity(emin, emax, models) # Write out trial result - ioutils.write_csv_row(self._outfile.url(), ieng, colnames, result) + ioutils.write_csv_row(self['outfile'].filename().url(), ieng, + colnames, result) # Append results results.append(result) diff --git a/doc/source/users/reference_manual/cssens.rst b/doc/source/users/reference_manual/cssens.rst index bf22f90..2cec9c5 100644 --- a/doc/source/users/reference_manual/cssens.rst +++ b/doc/source/users/reference_manual/cssens.rst @@ -100,9 +100,6 @@ General parameters ``(max_iter = 50) [integer]`` Maximum number of iterations. -``(num_avg = 3) [integer]`` - Number of iterations for sliding average. - Standard parameters ------------------- diff --git a/doc/source/users/reference_manual/ctbutterfly.rst b/doc/source/users/reference_manual/ctbutterfly.rst index cd8cbd2..b15314d 100644 --- a/doc/source/users/reference_manual/ctbutterfly.rst +++ b/doc/source/users/reference_manual/ctbutterfly.rst @@ -3,21 +3,25 @@ ctbutterfly =========== -Computes butterfly diagram for a power law model. +Computes butterfly diagram for a given spectral model. Synopsis -------- -This tool calculates a butterfly diagram for a specific source with power law -spectral model. The butterfly diagram is the envelope of all power law models -that are within a given confidence limit compatible with the data. By default -a confidence level of 68% is used, but the level can be adjusted using the -hidden ``confidence`` parameter. ctbutterfly computes this envelope by -evaluating for each energy the minimum and maximum intensity of all power law -models that fall within the error ellipse of the prefactor and index parameters. -The error ellipse is derived from the covariance matrix of a maximum likelihood -fit. +This tool calculates a butterfly diagram for a specific source according to +its spectral model. The butterfly diagram is the envelope of all spectral +models that are within a given confidence limit compatible with the data. +The default method used for the calculation is Gaussian error propagation +using the covariance matrix from a maximum likelihood fit. By default a +confidence level of 68% is used, but the level can be adjusted using the +hidden ``confidence`` parameter. For power law models, an alternative +calculation method can be specified by setting the hidden parameter +``method=ENVELOPE``. By using this method, the envelope is computed by +evaluating for each energy the minimum and maximum intensity of all power +law models that fall within the error ellipse of the prefactor and index +parameters. The error ellipse is derived from the covariance matrix of a +maximum likelihood fit. ctbutterfly assumes that the input model (parameter ``inmodel``) has been adjusted using :doc:`ctlike` to the data, but if this is not the case you @@ -68,11 +72,14 @@ General parameters ``outfile [file]`` Output butterfly ASCII file name. - + ``(fit = no) [boolean]`` Performs maximum likelihood fitting of input model ignoring any provided covariance matrix. - + +``(method = "GAUSSIAN") [string]`` + Computation method. + ``(confidence = 0.68) [real]`` Confidence level for error computation. diff --git a/doc/source/users/reference_manual/ctlike.rst b/doc/source/users/reference_manual/ctlike.rst index 95d76d2..dcdcd4e 100644 --- a/doc/source/users/reference_manual/ctlike.rst +++ b/doc/source/users/reference_manual/ctlike.rst @@ -48,13 +48,13 @@ the spatial model parameters can be kept fixed during that computation ctlike generates an output model XML file that contains the values of the best fitting model parameters. For all free parameters, an ``error`` attribute is added that provides the statistical uncertainty in the parameter estimate. -In addition, the entire output of the covariance matrix to a separate FITS file -can be queried (parameter ``outcovmat``). The parameter order given in the file -corresponds to the arrangement of the covariance matrix entries. If computation -of the Test Statistics has been requested for a model component, a ``ts`` -attribute providing the Test Statistics value is added. The output model can be -used as an input model for other ctools. +If computation of the Test Statistics has been requested for a model component, +a ``ts`` attribute providing the Test Statistics value is added. The output model +can be used as an input model for other ctools. +In addition, the covariance matrix can be saved into a file by specifying a +file name for the ``outcovmat`` parameter. If the file name ends with ``.fits`` +or ``.fit`` a FITS file is written, otherwise a CSV file is written. General parameters @@ -89,7 +89,7 @@ General parameters the maximum likelihood fit. ``(outcovmat = NONE) [string]`` - Output FITS file to store covariance matrix. + Output FITS or CSV file to store covariance matrix. ``(stat = POISSON) [string]`` Fitting statistics (POISSON or GAUSSIAN; only affects binned analysis). diff --git a/doc/source/users/reference_manual/ctprob.rst b/doc/source/users/reference_manual/ctprob.rst new file mode 100644 index 0000000..5da08e3 --- /dev/null +++ b/doc/source/users/reference_manual/ctprob.rst @@ -0,0 +1,85 @@ +.. _ctprob: + +ctprob +======== + +Compute event probability for a given model. + + +Synopsis +-------- + +This tool takes an event list file (or an XML observation definition file) and +a model as input and computes for each event the probability that it originates +from a specific model component. The probability is computed by by evaluating +the differential event probability for a given model and normalizing these +probabilities so that the sum of the probabilities for all model components is +unity for each event. + +The tool appends for each model component a single precision column to the event +list that contains the event probabilities. The name of the column is build from +the model name prefixed with ``PROB_``. + + +General parameters +------------------ + +``inobs [file]`` + Input event list or observation definition XML file + +``outobs [file]`` + Output event list or observation definition XML file. + +``(prefix = "prob_") [string]`` + Prefix for output event lists in observation definition file. + +``inmodel [string]`` + Input model definition XML file. + +``caldb [string]`` + Calibration database. + +``irf [string]`` + Instrument response function. + +``(edisp = no) [boolean]`` + Applies energy dispersion to response computation. + + +Standard parameters +------------------- + +``(publish = no) [boolean]`` + Specifies whether the event list(s) should be published on VO Hub. + +``(chatter = 2) [integer]`` + Verbosity of the executable: + chatter = 0: no information will be logged + + chatter = 1: only errors will be logged + + chatter = 2: errors and actions will be logged + + chatter = 3: report about the task execution + + chatter = 4: detailed report about the task execution + +``(clobber = yes) [boolean]`` + Specifies whether existing output files should be overwritten. + +``(debug = no) [boolean]`` + Enables debug mode. In debug mode the executable will dump any log file output to the console. + +``(mode = ql) [string]`` + Mode of automatic parameters (default is "ql", i.e. "query and learn"). + +``(logfile = ctprob.log) [string]`` + Name of log file. + + +Related tools or scripts +------------------------ + +:doc:`ctphase` +:doc:`csphasecrv` + diff --git a/doc/source/users/reference_manual/index.rst b/doc/source/users/reference_manual/index.rst index 332e84b..d4a2c2c 100644 --- a/doc/source/users/reference_manual/index.rst +++ b/doc/source/users/reference_manual/index.rst @@ -27,6 +27,7 @@ ctools ctmodel --- Computes model counts cube ctobssim --- Simulate observations ctpsfcube --- Generates point spread function cube + ctprob --- Computes event probability for a given model ctselect --- Selects event data ctskymap --- Generates sky map cttsmap --- Generates Test Statistic map diff --git a/pyext/Makefile.am b/pyext/Makefile.am index e38a57a..82faf59 100755 --- a/pyext/Makefile.am +++ b/pyext/Makefile.am @@ -59,6 +59,7 @@ testscripts = ctools/tests/__init__.py \ $(top_srcdir)/test/test_ctbutterfly.py \ $(top_srcdir)/test/test_ctulimit.py \ $(top_srcdir)/test/test_cterror.py \ + $(top_srcdir)/test/test_ctprob.py \ $(top_srcdir)/test/test_pipelines.py # Test data to be distributed in $(pyexecdir)/ctools/tests/data diff --git a/pyext/ctool.i b/pyext/ctool.i index 0d80e84..e99778f 100644 --- a/pyext/ctool.i +++ b/pyext/ctool.i @@ -1,7 +1,7 @@ /*************************************************************************** * ctool - ctool base class * * ----------------------------------------------------------------------- * - * copyright (C) 2014-2016 by Juergen Knoedlseder * + * copyright (C) 2014-2017 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -105,6 +105,8 @@ public: %rename(_insert_energy_boundaries) insert_energy_boundaries; %rename(_cube_layer_usage) cube_layer_usage; %rename(_is_valid_filename) is_valid_filename; + %rename(_get_gtiname) get_gtiname; + %rename(_save_event_list) save_event_list; %rename(_warn_too_few_energies) warn_too_few_energies; %rename(_warn_xml_suffix) warn_xml_suffix; @@ -147,6 +149,13 @@ public: std::vector cube_layer_usage(const GEbounds& cube_ebounds, const GEbounds& list_ebounds) const; bool is_valid_filename(const GFilename& filename) const; + std::string get_gtiname(const std::string& filename, + const std::string& evtname) const; + void save_event_list(const GCTAObservation* obs, + const std::string& infile, + const std::string& evtname, + const std::string& gtiname, + const std::string& outfile) const; // Protected warning strings std::string warn_too_few_energies(const GEnergies& energies) const; diff --git a/pyext/ctools/tools.i b/pyext/ctools/tools.i index e15ffba..60229e3 100755 --- a/pyext/ctools/tools.i +++ b/pyext/ctools/tools.i @@ -73,6 +73,7 @@ namespace std { %include "ctbutterfly.i" %include "ctulimit.i" %include "cterror.i" +%include "ctprob.i" /* __ Test function ______________________________________________________ */ %pythoncode %{ diff --git a/pyext/ctphase.i b/pyext/ctphase.i index 13e681b..b141aba 100644 --- a/pyext/ctphase.i +++ b/pyext/ctphase.i @@ -1,7 +1,7 @@ /*************************************************************************** * ctphase - Append phase information to CTA events file * * ----------------------------------------------------------------------- * - * copyright (C) 2017 by Leonardo Di Venere * + * copyright (C) 2017 by Joshua Cardenzana * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -21,7 +21,7 @@ /** * @file ctphase.i * @brief Append phase information to CTA events file - * @author Leonardo Di Venere + * @author Joshua Cardenzana */ %{ /* Put headers and other declarations here that are needed for compilation */ diff --git a/pyext/ctprob.i b/pyext/ctprob.i new file mode 100644 index 0000000..eec5d4c --- /dev/null +++ b/pyext/ctprob.i @@ -0,0 +1,61 @@ +/*************************************************************************** + * ctprob - Computes probability for a given model * + * ----------------------------------------------------------------------- * + * copyright (C) 2012-2016 by Leonardo Di Venere * + * ----------------------------------------------------------------------- * + * * + * 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 * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file ctprob.i + * @brief Computes probability for a given model + * @author Leonardo Di Venere + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "ctprob.hpp" +%} + + +/***********************************************************************//** + * @class ctprob + * + * @brief Data selection tool + ***************************************************************************/ +class ctprob : public ctobservation { +public: + // Constructors and destructors + ctprob(void); + explicit ctprob(const GObservations& obs); + ctprob(int argc, char *argv[]); + ctprob(const ctprob& app); + virtual ~ctprob(void); + + // Methods + void clear(void); + void run(void); + void save(void); + void publish(const std::string& name = ""); +}; + + +/***********************************************************************//** + * @brief Data selection tool Python extension + ***************************************************************************/ +%extend ctprob { + ctprob copy() { + return (*self); + } +} diff --git a/src/Makefile.am b/src/Makefile.am index b30eae9..138cdd1 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,7 +1,7 @@ ############################################################################# # Makefile.am for ctools module # # ------------------------------------------------------------------------- # -# Copyright (C) 2016 Juergen Knoedlseder # +# Copyright (C) 2016-2017 Juergen Knoedlseder # # ------------------------------------------------------------------------- # # # # This program is free software: you can redistribute it and/or modify # @@ -41,7 +41,8 @@ SUBDIRS = support \ ctbutterfly \ ctulimit \ cterror \ - ctphase + ctphase \ + ctprob # Scripts to be distributed scripts = ctools-setup \ @@ -70,7 +71,8 @@ libctools_la_LIBADD = support/libsupport.la \ ctbutterfly/libctbutterfly.la \ ctulimit/libctulimit.la \ cterror/libcterror.la \ - ctphase/libctphase.la + ctphase/libctphase.la \ + ctprob/libctprob.la libctools_la_LDFLAGS = -version-info $(CTOOLS_LT_VERSION) @OPENMP_LIB_LDFLAGS@ diff --git a/src/ctbin/ctbin.cpp b/src/ctbin/ctbin.cpp index 56a8975..585c614 100644 --- a/src/ctbin/ctbin.cpp +++ b/src/ctbin/ctbin.cpp @@ -55,7 +55,7 @@ * * Constructs an empty event binning tool. ***************************************************************************/ -ctbin::ctbin(void) : ctobservation(CTBIN_NAME, CTBIN_VERSION) +ctbin::ctbin(void) : ctobservation(CTBIN_NAME, VERSION) { // Initialise members init_members(); @@ -73,7 +73,7 @@ ctbin::ctbin(void) : ctobservation(CTBIN_NAME, CTBIN_VERSION) * Constructs event binning tool from an observation container. ***************************************************************************/ ctbin::ctbin(const GObservations& obs) : - ctobservation(CTBIN_NAME, CTBIN_VERSION, obs) + ctobservation(CTBIN_NAME, VERSION, obs) { // Initialise members init_members(); @@ -94,7 +94,7 @@ ctbin::ctbin(const GObservations& obs) : * parameter setting. ***************************************************************************/ ctbin::ctbin(int argc, char *argv[]) : - ctobservation(CTBIN_NAME, CTBIN_VERSION, argc, argv) + ctobservation(CTBIN_NAME, VERSION, argc, argv) { // Initialise members init_members(); diff --git a/src/ctbin/ctbin.hpp b/src/ctbin/ctbin.hpp index 5b3cf03..afb0cc3 100644 --- a/src/ctbin/ctbin.hpp +++ b/src/ctbin/ctbin.hpp @@ -31,8 +31,7 @@ #include "ctobservation.hpp" /* __Definitions _________________________________________________________ */ -#define CTBIN_NAME "ctbin" -#define CTBIN_VERSION "1.2.0" +#define CTBIN_NAME "ctbin" /***********************************************************************//** diff --git a/src/ctbkgcube/ctbkgcube.cpp b/src/ctbkgcube/ctbkgcube.cpp index 80b5b4a..ab48178 100644 --- a/src/ctbkgcube/ctbkgcube.cpp +++ b/src/ctbkgcube/ctbkgcube.cpp @@ -50,7 +50,7 @@ /***********************************************************************//** * @brief Void constructor ***************************************************************************/ -ctbkgcube::ctbkgcube(void) : ctobservation(CTBKGCUBE_NAME, CTBKGCUBE_VERSION) +ctbkgcube::ctbkgcube(void) : ctobservation(CTBKGCUBE_NAME, VERSION) { // Initialise members init_members(); @@ -69,7 +69,7 @@ ctbkgcube::ctbkgcube(void) : ctobservation(CTBKGCUBE_NAME, CTBKGCUBE_VERSION) * observations container. ***************************************************************************/ ctbkgcube::ctbkgcube(const GObservations& obs) : - ctobservation(CTBKGCUBE_NAME, CTBKGCUBE_VERSION, obs) + ctobservation(CTBKGCUBE_NAME, VERSION, obs) { // Initialise members init_members(); @@ -87,7 +87,7 @@ ctbkgcube::ctbkgcube(const GObservations& obs) : * @param[in] argv Array of command line arguments. ***************************************************************************/ ctbkgcube::ctbkgcube(int argc, char *argv[]) : - ctobservation(CTBKGCUBE_NAME, CTBKGCUBE_VERSION, argc, argv) + ctobservation(CTBKGCUBE_NAME, VERSION, argc, argv) { // Initialise members init_members(); diff --git a/src/ctbkgcube/ctbkgcube.hpp b/src/ctbkgcube/ctbkgcube.hpp index a2793f0..85d2cc6 100644 --- a/src/ctbkgcube/ctbkgcube.hpp +++ b/src/ctbkgcube/ctbkgcube.hpp @@ -33,8 +33,7 @@ #include "ctobservation.hpp" /* __Definitions _________________________________________________________ */ -#define CTBKGCUBE_NAME "ctbkgcube" -#define CTBKGCUBE_VERSION "1.1.0" +#define CTBKGCUBE_NAME "ctbkgcube" /***********************************************************************//** diff --git a/src/ctbutterfly/ctbutterfly.cpp b/src/ctbutterfly/ctbutterfly.cpp index cd462bb..e4cb59d 100644 --- a/src/ctbutterfly/ctbutterfly.cpp +++ b/src/ctbutterfly/ctbutterfly.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctbutterfly - butterfly calculation tool * * ----------------------------------------------------------------------- * - * copyright (C) 2014-2016 by Michael Mayer * + * copyright (C) 2014-2017 by Michael Mayer * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -52,7 +52,7 @@ /***********************************************************************//** * @brief Void constructor ***************************************************************************/ -ctbutterfly::ctbutterfly(void) : ctlikelihood(CTBUTTERFLY_NAME, CTBUTTERFLY_VERSION) +ctbutterfly::ctbutterfly(void) : ctlikelihood(CTBUTTERFLY_NAME, VERSION) { // Initialise members init_members(); @@ -71,7 +71,7 @@ ctbutterfly::ctbutterfly(void) : ctlikelihood(CTBUTTERFLY_NAME, CTBUTTERFLY_VERS * observations container. ***************************************************************************/ ctbutterfly::ctbutterfly(const GObservations& obs) : - ctlikelihood(CTBUTTERFLY_NAME, CTBUTTERFLY_VERSION, obs) + ctlikelihood(CTBUTTERFLY_NAME, VERSION, obs) { // Initialise members init_members(); @@ -89,7 +89,7 @@ ctbutterfly::ctbutterfly(const GObservations& obs) : * @param[in] argv Array of command line arguments. ***************************************************************************/ ctbutterfly::ctbutterfly(int argc, char *argv[]) : - ctlikelihood(CTBUTTERFLY_NAME, CTBUTTERFLY_VERSION, argc, argv) + ctlikelihood(CTBUTTERFLY_NAME, VERSION, argc, argv) { // Initialise members init_members(); @@ -243,11 +243,14 @@ void ctbutterfly::run(void) m_obs.optimize(m_opt); m_obs.errors(m_opt); - // Get covariance matrix by inverting the curvature matrix. Note that - // the curvature matrix is for the unscaled parameters, hence we have - // to scale them before we can proceed. + // Get covariance matrix. GObservations::likelihood likelihood = m_obs.function(); - m_covariance = likelihood.curvature()->invert(); + if (m_method == "GAUSSIAN") { + m_covariance = likelihood.curvature()->invert(); + } + else { + m_covariance = likelihood.covariance(); + } // Write optimizer, optimised model and covariance matrix into logger log_string(NORMAL, m_opt.print(m_chatter)); @@ -271,157 +274,26 @@ void ctbutterfly::run(void) GObservations::likelihood likelihood = m_obs.function(); likelihood.eval(pars); - // Get covariance matrix by inverting the curvature matrix. Note that - // the curvature matrix is for the unscaled parameters, hence we have - // to scale them before we can proceed. - m_covariance = likelihood.curvature()->invert(); + // Get covariance matrix. + if (m_method == "GAUSSIAN") { + m_covariance = likelihood.curvature()->invert(); + } + else { + m_covariance = likelihood.covariance(); + } // Write covariance matrix into logger log_string(EXPLICIT, m_covariance.print(m_chatter)); } // endif: covariance computation needed - // Write header into logger - log_header1(TERSE, "Prepare butterfly computation"); - - // Find parameter indices. We do this by the memory location of the - // model parameters. - GModelSky* skymodel = dynamic_cast(models[m_srcname]); - GModelPar* par_prefactor = &(*skymodel)["Prefactor"]; - GModelPar* par_index = &(*skymodel)["Index"]; - int inx_prefactor = -1; - int inx_index = -1; - for (int i = 0; i < models.npars(); ++i) { - if (models.pars()[i] == par_prefactor) { - inx_prefactor = i; - } - if (models.pars()[i] == par_index) { - inx_index = i; - } + // Compute butterfly depending on selected method + if (m_method == "GAUSSIAN") { + gaussian_error_propagation(models); + } + else { + ellipsoid_boundary(models); } - - // Store model values and scales - double prefactor_mean = (*skymodel)["Prefactor"].value(); - double index_mean = (*skymodel)["Index"].value(); - double pivot_mean = (*skymodel)["PivotEnergy"].value(); - double prefactor_scale = (*skymodel)["Prefactor"].scale(); - double index_scale = (*skymodel)["Index"].scale(); - GEnergy pivot(pivot_mean, "MeV"); - - // Write parameter indices into logger - log_value(NORMAL, "Prefactor", gammalib::str(prefactor_mean)+ - " ph/cm2/s/MeV ("+gammalib::str(inx_prefactor)+")"); - log_value(NORMAL, "Index", gammalib::str(index_mean)+ - " ("+gammalib::str(inx_index)+")"); - log_value(NORMAL, "Pivot energy", pivot.print()); - - // Extract covariance elements - double cov_pp = m_covariance(inx_prefactor,inx_prefactor); - double cov_pg = m_covariance(inx_prefactor,inx_index); - double cov_gg = m_covariance(inx_index,inx_index); - - // Compute eigenvectors and eigenvalues - double lambda1; - double lambda2; - GVector vector1(2); - GVector vector2(2); - eigenvectors(cov_pp, cov_pg, cov_pg, cov_gg, - &lambda1, &lambda2, &vector1, &vector2); - - // Write eigenvectors and eigenvalues into logger - log_value(NORMAL, "Eigenvalue 1", lambda1); - log_value(NORMAL, "Eigenvalue 2", lambda2); - log_value(NORMAL, "Eigenvector 1", vector1.print()); - log_value(NORMAL, "Eigenvector 2", vector2.print()); - - // Confidence scaling - double sigma = gammalib::erfinv(m_confidence) * gammalib::sqrt_two; - double scale = (sigma*sigma); - - // Write confidence level information into logger - log_value(NORMAL, "Confidence level", m_confidence); - log_value(NORMAL, "Corresponding scaling", scale); - - // Compute minor and major axes of error ellipse - double major = scale*std::sqrt(lambda1); - double minor = scale*std::sqrt(lambda2); - - // Write header into logger - log_header1(TERSE, "Generate butterfly"); - - // Initialise result arrays - m_energies.assign(m_ebounds.size(), 0.0); - m_intensities.assign(m_ebounds.size(), 0.0); - m_min_intensities.assign(m_ebounds.size(), 1e30); - m_max_intensities.assign(m_ebounds.size(), 0.0); - - // Walk around the error ellipse - int nt = 360; - double dt = gammalib::twopi/double(nt); - double t = 0.0; - for (int i = 0; i < nt; ++i, t += dt) { - - // Pre-compute sin and cosine - double sin_t = std::sin(t); - double cos_t = std::cos(t); - - // Compute prefactor and index - double prefactor = (major * cos_t * vector1[0] + - minor * sin_t * vector2[0]) * prefactor_scale + - prefactor_mean; - double index = (major * cos_t * vector1[1] + - minor * sin_t * vector2[1]) * index_scale + - index_mean; - - // Write prefactor and index into logger - log_value(EXPLICIT, "Angle "+gammalib::str(t*gammalib::rad2deg), - gammalib::str(prefactor)+" ph/cm2/s/MeV; "+ - gammalib::str(index)); - - // Setup power law model - GModelSpectralPlaw plaw(prefactor, index, pivot); - - // Evaluate power law at each energy and extract minimum and - // maximum intensity - for (int k = 0; k < m_ebounds.size(); ++k) { - - // Get the energy of current bin - GTime time; - GEnergy energy = m_ebounds.elogmean(k); - - // Evaluate power law - double intensity = plaw.eval(energy, time); - - // Get extremes - if (intensity < m_min_intensities[k]) { - m_min_intensities[k] = intensity; - } - if (intensity > m_max_intensities[k]) { - m_max_intensities[k] = intensity; - } - - } // endfor: looped over all energies - - } // endfor: looped over all angles - - // Setup power law model for mean prefactor and index - GModelSpectralPlaw plaw(prefactor_mean, index_mean, pivot); - - // Compute energy and intensity vectors for mean power law - for (int k = 0; k < m_ebounds.size(); ++k) { - - // Get the energy of current bin - GTime time; - GEnergy energy = m_ebounds.elogmean(k); - - // Evaluate power law - double intensity = plaw.eval(energy, time); - - // Store results - m_energies[k] = energy.MeV(); - m_intensities[k] = intensity; - - } // endfor: looped over energies // Write header into logger log_header1(TERSE, "Results"); @@ -431,11 +303,11 @@ void ctbutterfly::run(void) // Write results into logger log_value(NORMAL, "Intensity at "+gammalib::str(m_energies[k])+" MeV", - gammalib::str(m_intensities[k])+" ("+ - gammalib::str(m_min_intensities[k])+", "+ - gammalib::str(m_max_intensities[k])+") ph/cm2/s/MeV"); + gammalib::str(m_intensities[k])+" ("+ + gammalib::str(m_min_intensities[k])+", "+ + gammalib::str(m_max_intensities[k])+") ph/cm2/s/MeV"); - } // endfor: looped over energies + } //endfor: loop over energy bin // Restore energy dispersion flag for all CTA observations for (int i = 0; i < m_obs.size(); ++i) { @@ -504,6 +376,7 @@ void ctbutterfly::init_members(void) { // Initialise members m_srcname.clear(); + m_method.clear(); m_confidence = 0.68; m_max_iter = 50; m_apply_edisp = false; @@ -537,6 +410,7 @@ void ctbutterfly::copy_members(const ctbutterfly& app) { // Copy attributes m_srcname = app.m_srcname; + m_method = app.m_method; m_confidence = app.m_confidence; m_max_iter = app.m_max_iter; m_apply_edisp = app.m_apply_edisp; @@ -544,7 +418,7 @@ void ctbutterfly::copy_members(const ctbutterfly& app) m_chatter = app.m_chatter; m_ebounds = app.m_ebounds; m_outfile = app.m_outfile; - + // Copy protected members m_covariance = app.m_covariance; m_energies = app.m_energies; @@ -583,16 +457,13 @@ void ctbutterfly::get_parameters(void) { // Setup observations from "inobs" parameter setup_observations(m_obs); - + // Setup models from "inmodel" parameter setup_models(m_obs, (*this)["srcname"].string()); // Get name of test source m_srcname = (*this)["srcname"].string(); - // Check model name and type - check_model(); - // Create energy boundaries from user parameters m_ebounds = create_ebounds(); @@ -600,18 +471,22 @@ void ctbutterfly::get_parameters(void) std::string matrixfilename = (*this)["matrix"].filename(); if (matrixfilename != "NONE") { std::string msg = "Loading of matrix from file not implemented yet. " - "Use filename = \"NONE\" to induce a recomputation of " - "the matrix internally."; + "Use filename = \"NONE\" to induce a recomputation of " + "the matrix internally."; throw GException::feature_not_implemented(G_GET_PARAMETERS, msg); // m_covariance.load(matrixfilename); } // Get remaining parameters m_apply_edisp = (*this)["edisp"].boolean(); + m_method = (*this)["method"].string(); m_confidence = (*this)["confidence"].real(); m_fit = (*this)["fit"].boolean(); m_chatter = static_cast((*this)["chatter"].integer()); + // Check model name and type + check_model(); + // Optionally read ahead parameters so that they get correctly // dumped into the log file if (read_ahead()) { @@ -627,6 +502,279 @@ void ctbutterfly::get_parameters(void) /***********************************************************************//** + * @brief Compute butterfly using the ellipsoid boundary method + ***************************************************************************/ +void ctbutterfly::ellipsoid_boundary(GModels& models) +{ + // Write header into logger + log_header1(TERSE, "Prepare butterfly computation"); + + // Find parameter indices. We do this by the memory location of the + // model parameters. + GModelSky* skymodel = dynamic_cast(models[m_srcname]); + GModelPar* par_prefactor = &(*skymodel)["Prefactor"]; + GModelPar* par_index = &(*skymodel)["Index"]; + int inx_prefactor = -1; + int inx_index = -1; + for (int i = 0; i < models.npars(); ++i) { + if (models.pars()[i] == par_prefactor) { + inx_prefactor = i; + } + if (models.pars()[i] == par_index) { + inx_index = i; + } + } + + // Store model values + double prefactor_mean = (*skymodel)["Prefactor"].value(); + double index_mean = (*skymodel)["Index"].value(); + double pivot_mean = (*skymodel)["PivotEnergy"].value(); + double prefactor_scale = (*skymodel)["Prefactor"].scale(); + double index_scale = (*skymodel)["Index"].scale(); + GEnergy pivot(pivot_mean, "MeV"); + + // Write parameter indices into logger + log_value(NORMAL, "Prefactor", gammalib::str(prefactor_mean)+ + " ph/cm2/s/MeV ("+ + gammalib::str(inx_prefactor)+")"); + log_value(NORMAL, "Index", gammalib::str(index_mean)+ + " ("+gammalib::str(inx_index)+")"); + log_value(NORMAL, "Pivot energy", pivot.print()); + + // Extract covariance elements + double cov_pp = m_covariance(inx_prefactor,inx_prefactor); + double cov_pg = m_covariance(inx_prefactor,inx_index); + double cov_gg = m_covariance(inx_index,inx_index); + + // Compute eigenvectors and eigenvalues + double lambda1; + double lambda2; + GVector vector1(2); + GVector vector2(2); + eigenvectors(cov_pp, cov_pg, cov_pg, cov_gg, + &lambda1, &lambda2, &vector1, &vector2); + + // Write eigenvectors and eigenvalues into logger + log_value(NORMAL, "Eigenvalue 1", lambda1); + log_value(NORMAL, "Eigenvalue 2", lambda2); + log_value(NORMAL, "Eigenvector 1", vector1.print()); + log_value(NORMAL, "Eigenvector 2", vector2.print()); + + // Confidence scaling, calculated from a Chi-Sqr. + // distribution with 2 dof + double scale = std::sqrt( -2. * std::log(1.-m_confidence) ); + + // Write confidence level information into logger + log_value(NORMAL, "Confidence level", m_confidence); + log_value(NORMAL, "Corresponding scaling", scale); + + // Compute minor and major axes of error ellipse + double major = scale*std::sqrt(lambda1); + double minor = scale*std::sqrt(lambda2); + + // Write header into logger + log_header1(TERSE, "Generate butterfly"); + + // Initialise result arrays + m_energies.assign(m_ebounds.size(), 0.0); + m_intensities.assign(m_ebounds.size(), 0.0); + m_min_intensities.assign(m_ebounds.size(), 1e30); + m_max_intensities.assign(m_ebounds.size(), 0.0); + + // Walk around the error ellipse + int nt = 360; + double dt = gammalib::twopi/double(nt); + double t = 0.0; + for (int i = 0; i < nt; ++i, t += dt) { + + // Pre-compute sin and cosine + double sin_t = std::sin(t); + double cos_t = std::cos(t); + + // Compute prefactor and index + double prefactor = (major * cos_t * vector1[0] + + minor * sin_t * vector2[0]) * prefactor_scale + + prefactor_mean; + double index = (major * cos_t * vector1[1] + + minor * sin_t * vector2[1]) * index_scale + + index_mean; + + // Write prefactor and index into logger + log_value(EXPLICIT, "Angle "+gammalib::str(t*gammalib::rad2deg), + gammalib::str(prefactor)+" ph/cm2/s/MeV; "+ + gammalib::str(index)); + + // Setup power law model + GModelSpectralPlaw plaw(prefactor, index, pivot); + + // Evaluate power law at each energy and extract minimum and + // maximum intensity + for (int k = 0; k < m_ebounds.size(); ++k) { + + // Get the energy of current bin + GTime time; + GEnergy energy = m_ebounds.elogmean(k); + + // Evaluate power law + double intensity = plaw.eval(energy, time); + + // Get extremes + if (intensity < m_min_intensities[k]) { + m_min_intensities[k] = intensity; + } + if (intensity > m_max_intensities[k]) { + m_max_intensities[k] = intensity; + } + + } // endfor: looped over all energies + + } // endfor: looped over all angles + + // Setup power law model for mean prefactor and index + GModelSpectralPlaw plaw(prefactor_mean, index_mean, pivot); + + // Compute energy and intensity vectors for mean power law + for (int k = 0; k < m_ebounds.size(); ++k) { + + // Get the energy of current bin + GTime time; + GEnergy energy = m_ebounds.elogmean(k); + + // Evaluate power law + double intensity = plaw.eval(energy, time); + + // Store results + m_energies[k] = energy.MeV(); + m_intensities[k] = intensity; + + } // endfor: looped over energies + + // Return + return; +} + + +/***********************************************************************//** + * @brief Compute butterfly using Gaussian error propagation + * + * Computes the butterfly diagram using Gaussian error propagation. This + * works for any kind of model. + ***************************************************************************/ +void ctbutterfly::gaussian_error_propagation(GModels& models) +{ + // Write header into logger + log_header1(TERSE, "Prepare butterfly computation"); + + // Confidence scaling + double scale = gammalib::erfinv(m_confidence) * gammalib::sqrt_two; + + // Write confidence level information into logger + log_value(NORMAL, "Confidence level", m_confidence); + log_value(NORMAL, "Corresponding scaling", scale); + + // Write header into logger + log_header1(TERSE, "Generate butterfly"); + + // Initialise result arrays + m_energies.clear(); + m_intensities.clear(); + m_min_intensities.clear(); + m_max_intensities.clear(); + + // Initialise dummy time to evaluate spectral model + GTime time = GTime(); + + // Initialise vector of gradients + GVector grad = GVector(m_obs.models().npars()); + + // Loop over energy bins + for (int i = 0; i < m_ebounds.size(); ++i) { + + // Initialise number of gradients + int num_gradient = 0; + + // Get the mean logarithmic energy of current bin + GEnergy energy = m_ebounds.elogmean(i); + + // Initialise model flux value + double model_flux = 0.0; + + // Loop over models + for (int j = 0; j < models.size(); ++j) { + + // Check wether model is a skymodel + GModelSky* skymodel = dynamic_cast(models[j]); + + // Yes ... + if (skymodel != NULL) { + + // Set flux value of source of interest, evaluate gradients + if (skymodel->name() == m_srcname) { + + // Get spectral component of model + GModelSpectral *spectral = skymodel->spectral(); + + // Get model flux + model_flux = spectral->eval(energy, time, true); + + // Loop over all model parameters + for (int k = 0; k < skymodel->size(); ++k) { + + // Set gradient if we deal with spectral component + bool parIsSpectral = false; + for (int l = 0; l < spectral->size(); ++l) { + if ((&((*spectral)[l])) == (&((*skymodel)[k]))) { + grad[num_gradient] = (*spectral)[l].gradient(); + num_gradient++; + parIsSpectral = true; + break; + } + } // endfor: looped over spectral parameters + + // If parameter is spatial or temporal, count up + if (!parIsSpectral) { + num_gradient++; + } + + } // endfor: looped over all model parameters + + } // endif: model was source of interest + else { + num_gradient += skymodel->size(); + } + + } // endif: model was sky model + + // Skip other models (e.g. GModelData instances) + else { + num_gradient += models[j]->size(); + } + + } // endfor: Looped over models + + // Multiply covariance to the gradient vector + GVector vector = m_covariance * grad; + + // Get the error from the scalar product + double error = std::sqrt(grad * vector); + + // Multiply with confidence scaling + error *= scale; + + // Store flux, value and energy for saving + m_intensities.push_back(model_flux); + m_energies.push_back(energy.MeV()); + m_min_intensities.push_back(model_flux-error); + m_max_intensities.push_back(model_flux+error); + + } // endfor: Looped over energy bins + + // Return + return; +} + + +/***********************************************************************//** * @brief Check if sky model is valid * * @exception GException::invalid_value @@ -645,17 +793,22 @@ void ctbutterfly::check_model(void) GModelSky* model = dynamic_cast(models[m_srcname]); if (model == NULL) { std::string msg = "Source \""+m_srcname+"\" is not a sky model. " - "Please specify the name of a sky model for " - "butterfly computation."; + "Please specify the name of a sky model for " + "butterfly computation."; throw GException::invalid_value(G_CHECK_MODEL, msg); } - // Check that the spectral model is a power law - if (model->spectral()->type() != "PowerLaw") { - std::string msg = "\""+model->spectral()->type()+"\" cannot be " - "used as spectral model for an butterfly " - "computation. Please specify a power law model."; - throw GException::invalid_value(G_CHECK_MODEL, msg); + // Check that the spectral model is a power law for the "ENVELOPE" + // method + if (m_method == "ENVELOPE") { + if (model->spectral()->type() != "PowerLaw") { + std::string msg = "\""+model->spectral()->type()+"\" cannot be " + "used as spectral model for an butterfly " + "computation with method=ENVELOPE. Please " + "specify a power law model or switch to " + "method=GAUSSIAN."; + throw GException::invalid_value(G_CHECK_MODEL, msg); + } } } // endif: source model existed @@ -663,8 +816,8 @@ void ctbutterfly::check_model(void) // Otherwise throw an exception else { std::string msg = "Source \""+m_srcname+"\" not found in model " - "container. Please add a source with that name " - "or check for possible typos."; + "container. Please add a source with that name " + "or check for possible typos."; throw GException::invalid_value(G_CHECK_MODEL, msg); } @@ -689,13 +842,13 @@ void ctbutterfly::check_model(void) * http://www.math.harvard.edu/archive/21b_fall_04/exhibits/2dmatrices/index.html ***************************************************************************/ void ctbutterfly::eigenvectors(const double& a, - const double& b, - const double& c, - const double& d, - double* lambda1, - double* lambda2, - GVector* vector1, - GVector* vector2) + const double& b, + const double& c, + const double& d, + double* lambda1, + double* lambda2, + GVector* vector1, + GVector* vector2) { // Compute trace and determinant double trace = a + d; @@ -730,7 +883,7 @@ void ctbutterfly::eigenvectors(const double& a, // Normalize eigenvectors *vector1 /= norm(*vector1); *vector2 /= norm(*vector2); - + // Return return; } diff --git a/src/ctbutterfly/ctbutterfly.hpp b/src/ctbutterfly/ctbutterfly.hpp index 09ce0fc..060ac57 100644 --- a/src/ctbutterfly/ctbutterfly.hpp +++ b/src/ctbutterfly/ctbutterfly.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctbutterfly - butterfly calculation tool * * ----------------------------------------------------------------------- * - * copyright (C) 2014-2016 by Michael Mayer * + * copyright (C) 2014-2017 by Michael Mayer * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -33,8 +33,7 @@ #include "ctlikelihood.hpp" /* __Definitions _________________________________________________________ */ -#define CTBUTTERFLY_NAME "ctbutterfly" -#define CTBUTTERFLY_VERSION "1.1.0" +#define CTBUTTERFLY_NAME "ctbutterfly" /***********************************************************************//** @@ -77,6 +76,10 @@ protected: void copy_members(const ctbutterfly& app); void free_members(void); void get_parameters(void); + void gaussian_error_propagation(GModels& models); + void ellipsoid_boundary(GModels& models); + + void check_model(void); void eigenvectors(const double& a, const double& b, @@ -89,6 +92,7 @@ protected: // User parameters std::string m_srcname; //!< Name of source to compute butterfly + std::string m_method; //!< Computation method double m_confidence; //!< Confidence level int m_max_iter; //!< Maximum number of iterations bool m_apply_edisp; //!< Apply energy dispersion? diff --git a/src/ctbutterfly/ctbutterfly.par b/src/ctbutterfly/ctbutterfly.par index fbce06e..5c3f228 100644 --- a/src/ctbutterfly/ctbutterfly.par +++ b/src/ctbutterfly/ctbutterfly.par @@ -27,6 +27,7 @@ caldb, s, a, "prod2",,, "Calibration database" irf, s, a, "South_0.5h",,, "Instrument response function" edisp, b, h, no,,, "Apply energy dispersion?" outfile, f, a, "butterfly.txt",,, "Output ASCII file" +method, s, h, GAUSSIAN, GAUSSIAN|ENVELOPE,,"Butterfly computation method" # # Butterfly parameters diff --git a/src/ctcubemask/ctcubemask.cpp b/src/ctcubemask/ctcubemask.cpp index dc98944..900f5dd 100644 --- a/src/ctcubemask/ctcubemask.cpp +++ b/src/ctcubemask/ctcubemask.cpp @@ -51,7 +51,7 @@ /***********************************************************************//** * @brief Void constructor ***************************************************************************/ -ctcubemask::ctcubemask(void) : ctobservation(CTCUBEMASK_NAME, CTCUBEMASK_VERSION) +ctcubemask::ctcubemask(void) : ctobservation(CTCUBEMASK_NAME, VERSION) { // Initialise members init_members(); @@ -70,7 +70,7 @@ ctcubemask::ctcubemask(void) : ctobservation(CTCUBEMASK_NAME, CTCUBEMASK_VERSION * an observation container. ***************************************************************************/ ctcubemask::ctcubemask(const GObservations& obs) : - ctobservation(CTCUBEMASK_NAME, CTCUBEMASK_VERSION, obs) + ctobservation(CTCUBEMASK_NAME, VERSION, obs) { // Initialise members init_members(); @@ -89,7 +89,7 @@ ctcubemask::ctcubemask(const GObservations& obs) : * @param[in] argv Array of command line arguments. ***************************************************************************/ ctcubemask::ctcubemask(int argc, char *argv[]) : - ctobservation(CTCUBEMASK_NAME, CTCUBEMASK_VERSION, argc, argv) + ctobservation(CTCUBEMASK_NAME, VERSION, argc, argv) { // Initialise members init_members(); diff --git a/src/ctcubemask/ctcubemask.hpp b/src/ctcubemask/ctcubemask.hpp index 0124aac..bea60b8 100644 --- a/src/ctcubemask/ctcubemask.hpp +++ b/src/ctcubemask/ctcubemask.hpp @@ -35,8 +35,7 @@ #include "ctobservation.hpp" /* __Definitions _________________________________________________________ */ -#define CTCUBEMASK_NAME "ctcubemask" -#define CTCUBEMASK_VERSION "1.2.0" +#define CTCUBEMASK_NAME "ctcubemask" /***********************************************************************//** diff --git a/src/ctedispcube/ctedispcube.cpp b/src/ctedispcube/ctedispcube.cpp index 160f4ab..fc0288f 100644 --- a/src/ctedispcube/ctedispcube.cpp +++ b/src/ctedispcube/ctedispcube.cpp @@ -52,7 +52,7 @@ * * Constructs an empty energy dispersion tool. ***************************************************************************/ -ctedispcube::ctedispcube(void) : ctobservation(CTEDISPCUBE_NAME, CTEDISPCUBE_VERSION) +ctedispcube::ctedispcube(void) : ctobservation(CTEDISPCUBE_NAME, VERSION) { // Initialise members init_members(); @@ -71,7 +71,7 @@ ctedispcube::ctedispcube(void) : ctobservation(CTEDISPCUBE_NAME, CTEDISPCUBE_VER * in an observation container @p obs. ***************************************************************************/ ctedispcube::ctedispcube(const GObservations& obs) : - ctobservation(CTEDISPCUBE_NAME, CTEDISPCUBE_VERSION, obs) + ctobservation(CTEDISPCUBE_NAME, VERSION, obs) { // Initialise members init_members(); @@ -92,7 +92,7 @@ ctedispcube::ctedispcube(const GObservations& obs) : * on the command line. ***************************************************************************/ ctedispcube::ctedispcube(int argc, char *argv[]) : - ctobservation(CTEDISPCUBE_NAME, CTEDISPCUBE_VERSION, argc, argv) + ctobservation(CTEDISPCUBE_NAME, VERSION, argc, argv) { // Initialise members init_members(); diff --git a/src/ctedispcube/ctedispcube.hpp b/src/ctedispcube/ctedispcube.hpp index 5a89c25..58d6cfe 100644 --- a/src/ctedispcube/ctedispcube.hpp +++ b/src/ctedispcube/ctedispcube.hpp @@ -33,8 +33,7 @@ #include "ctobservation.hpp" /* __Definitions _________________________________________________________ */ -#define CTEDISPCUBE_NAME "ctedispcube" -#define CTEDISPCUBE_VERSION "1.1.0" +#define CTEDISPCUBE_NAME "ctedispcube" /***********************************************************************//** diff --git a/src/cterror/cterror.cpp b/src/cterror/cterror.cpp index 8afce1b..9e9ef4a 100644 --- a/src/cterror/cterror.cpp +++ b/src/cterror/cterror.cpp @@ -52,7 +52,7 @@ * * Constructs an empty cterror tool. ***************************************************************************/ -cterror::cterror(void) : ctlikelihood(CTERROR_NAME, CTERROR_VERSION) +cterror::cterror(void) : ctlikelihood(CTERROR_NAME, VERSION) { // Initialise members init_members(); @@ -70,7 +70,7 @@ cterror::cterror(void) : ctlikelihood(CTERROR_NAME, CTERROR_VERSION) * Constructs cterror tool from an observations container. ***************************************************************************/ cterror::cterror(const GObservations& obs) : - ctlikelihood(CTERROR_NAME, CTERROR_VERSION, obs) + ctlikelihood(CTERROR_NAME, VERSION, obs) { // Initialise members init_members(); @@ -90,7 +90,7 @@ cterror::cterror(const GObservations& obs) : * parameters that are provided as command line arguments. ***************************************************************************/ cterror::cterror(int argc, char *argv[]) : - ctlikelihood(CTERROR_NAME, CTERROR_VERSION, argc, argv) + ctlikelihood(CTERROR_NAME, VERSION, argc, argv) { // Initialise members init_members(); diff --git a/src/cterror/cterror.hpp b/src/cterror/cterror.hpp index 24293f0..98b101e 100644 --- a/src/cterror/cterror.hpp +++ b/src/cterror/cterror.hpp @@ -31,8 +31,7 @@ #include "ctlikelihood.hpp" /* __Definitions _________________________________________________________ */ -#define CTERROR_NAME "cterror" -#define CTERROR_VERSION "1.2.0" +#define CTERROR_NAME "cterror" /***********************************************************************//** diff --git a/src/ctexpcube/ctexpcube.cpp b/src/ctexpcube/ctexpcube.cpp index 9cfe847..6d8d354 100644 --- a/src/ctexpcube/ctexpcube.cpp +++ b/src/ctexpcube/ctexpcube.cpp @@ -50,7 +50,7 @@ /***********************************************************************//** * @brief Void constructor ***************************************************************************/ -ctexpcube::ctexpcube(void) : ctobservation(CTEXPCUBE_NAME, CTEXPCUBE_VERSION) +ctexpcube::ctexpcube(void) : ctobservation(CTEXPCUBE_NAME, VERSION) { // Initialise members init_members(); @@ -69,7 +69,7 @@ ctexpcube::ctexpcube(void) : ctobservation(CTEXPCUBE_NAME, CTEXPCUBE_VERSION) * observations container. ***************************************************************************/ ctexpcube::ctexpcube(const GObservations& obs) : - ctobservation(CTEXPCUBE_NAME, CTEXPCUBE_VERSION, obs) + ctobservation(CTEXPCUBE_NAME, VERSION, obs) { // Initialise members init_members(); @@ -87,7 +87,7 @@ ctexpcube::ctexpcube(const GObservations& obs) : * @param[in] argv Array of command line arguments. ***************************************************************************/ ctexpcube::ctexpcube(int argc, char *argv[]) : - ctobservation(CTEXPCUBE_NAME, CTEXPCUBE_VERSION, argc, argv) + ctobservation(CTEXPCUBE_NAME, VERSION, argc, argv) { // Initialise members init_members(); diff --git a/src/ctexpcube/ctexpcube.hpp b/src/ctexpcube/ctexpcube.hpp index c2e159f..19bec02 100644 --- a/src/ctexpcube/ctexpcube.hpp +++ b/src/ctexpcube/ctexpcube.hpp @@ -33,8 +33,7 @@ #include "ctobservation.hpp" /* __Definitions _________________________________________________________ */ -#define CTEXPCUBE_NAME "ctexpcube" -#define CTEXPCUBE_VERSION "1.1.0" +#define CTEXPCUBE_NAME "ctexpcube" /***********************************************************************//** diff --git a/src/ctlike/ctlike.cpp b/src/ctlike/ctlike.cpp index 8199bdc..5229a94 100644 --- a/src/ctlike/ctlike.cpp +++ b/src/ctlike/ctlike.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctlike - Maximum likelihood fitting tool * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2016 by Juergen Knoedlseder * + * copyright (C) 2010-2017 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -48,7 +48,7 @@ /***********************************************************************//** * @brief Void constructor ***************************************************************************/ -ctlike::ctlike(void) : ctlikelihood(CTLIKE_NAME, CTLIKE_VERSION) +ctlike::ctlike(void) : ctlikelihood(CTLIKE_NAME, VERSION) { // Initialise members init_members(); @@ -66,7 +66,7 @@ ctlike::ctlike(void) : ctlikelihood(CTLIKE_NAME, CTLIKE_VERSION) * Constructs ctlike tool from an observations container. ***************************************************************************/ ctlike::ctlike(const GObservations& obs) : - ctlikelihood(CTLIKE_NAME, CTLIKE_VERSION, obs) + ctlikelihood(CTLIKE_NAME, VERSION, obs) { // Initialise members init_members(); @@ -86,7 +86,7 @@ ctlike::ctlike(const GObservations& obs) : * that are provided as command line arguments. ***************************************************************************/ ctlike::ctlike(int argc, char *argv[]) : - ctlikelihood(CTLIKE_NAME, CTLIKE_VERSION, argc, argv) + ctlikelihood(CTLIKE_NAME, VERSION, argc, argv) { // Initialise members init_members(); @@ -321,8 +321,13 @@ void ctlike::run(void) /***********************************************************************//** * @brief Save results * - * This method saves the fit results into a XML file and a FITS file. If - * the filename parameters are "NONE", the files are not saved. + * This method saves the fit results and the covariance matrix. The fit + * results are written into a XML file while the covariance matrix is written + * into either a FITS or a CSV file, depending on the file type extension + * (an extension of `.fits` or `.fit` produce a FITS file, any other + * extension produces a CSV file). + * + * If the filenames are `NONE` no information is saved. ***************************************************************************/ void ctlike::save(void) { diff --git a/src/ctlike/ctlike.hpp b/src/ctlike/ctlike.hpp index 3bad996..496364a 100644 --- a/src/ctlike/ctlike.hpp +++ b/src/ctlike/ctlike.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctlike - Maximum likelihood fitting tool * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2016 by Juergen Knoedlseder * + * copyright (C) 2010-2017 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -31,8 +31,7 @@ #include "ctlikelihood.hpp" /* __Definitions _________________________________________________________ */ -#define CTLIKE_NAME "ctlike" -#define CTLIKE_VERSION "1.2.0" +#define CTLIKE_NAME "ctlike" /***********************************************************************//** @@ -68,17 +67,17 @@ protected: double reoptimize_lm(void); // User parameters - bool m_refit; //!< Refitting - GFilename m_outmodel; //!< Source model output XML file name - GFilename m_outcovmat; //!< Covariance matrix output FITS file name - bool m_apply_edisp; //!< Apply energy dispersion? - bool m_fix_spat_for_ts; //!< Fix spatial parameters for TS computation? - GChatter m_chatter; //!< Chattiness + bool m_refit; //!< Refitting + GFilename m_outmodel; //!< Source model output XML file name + GFilename m_outcovmat; //!< Covariance matrix output file name + bool m_apply_edisp; //!< Apply energy dispersion? + bool m_fix_spat_for_ts; //!< Fix spatial parameters for TS computation? + GChatter m_chatter; //!< Chattiness // Members - int m_max_iter; //!< Maximum number of iterations - int m_max_stall; //!< Maximum number of stalls - double m_logL; //!< Maximum log likelihood + int m_max_iter; //!< Maximum number of iterations + int m_max_stall; //!< Maximum number of stalls + double m_logL; //!< Maximum log likelihood }; #endif /* CTLIKE_HPP */ diff --git a/src/ctmapcube/ctmapcube.cpp b/src/ctmapcube/ctmapcube.cpp index 7167b70..6c3af23 100644 --- a/src/ctmapcube/ctmapcube.cpp +++ b/src/ctmapcube/ctmapcube.cpp @@ -50,7 +50,7 @@ /***********************************************************************//** * @brief Void constructor ***************************************************************************/ -ctmapcube::ctmapcube(void) : ctool(CTMAPCUBE_NAME, CTMAPCUBE_VERSION) +ctmapcube::ctmapcube(void) : ctool(CTMAPCUBE_NAME, VERSION) { // Initialise members init_members(); @@ -67,7 +67,7 @@ ctmapcube::ctmapcube(void) : ctool(CTMAPCUBE_NAME, CTMAPCUBE_VERSION) * @param[in] argv Array of command line arguments. ***************************************************************************/ ctmapcube::ctmapcube(int argc, char *argv[]) : - ctool(CTMAPCUBE_NAME, CTMAPCUBE_VERSION, argc, argv) + ctool(CTMAPCUBE_NAME, VERSION, argc, argv) { // Initialise members init_members(); diff --git a/src/ctmapcube/ctmapcube.hpp b/src/ctmapcube/ctmapcube.hpp index 0f06341..23619a7 100644 --- a/src/ctmapcube/ctmapcube.hpp +++ b/src/ctmapcube/ctmapcube.hpp @@ -33,8 +33,7 @@ #include "ctool.hpp" /* __Definitions _________________________________________________________ */ -#define CTMAPCUBE_NAME "ctmapcube" -#define CTMAPCUBE_VERSION "1.3.0" +#define CTMAPCUBE_NAME "ctmapcube" /***********************************************************************//** diff --git a/src/ctmodel/ctmodel.cpp b/src/ctmodel/ctmodel.cpp index 0decffb..d4185b4 100644 --- a/src/ctmodel/ctmodel.cpp +++ b/src/ctmodel/ctmodel.cpp @@ -53,7 +53,7 @@ const GEnergy g_energy_margin(1.0e-12, "TeV"); /***********************************************************************//** * @brief Void constructor ***************************************************************************/ -ctmodel::ctmodel(void) : ctobservation(CTMODEL_NAME, CTMODEL_VERSION) +ctmodel::ctmodel(void) : ctobservation(CTMODEL_NAME, VERSION) { // Initialise members init_members(); @@ -72,7 +72,7 @@ ctmodel::ctmodel(void) : ctobservation(CTMODEL_NAME, CTMODEL_VERSION) * observations container. ***************************************************************************/ ctmodel::ctmodel(const GObservations& obs) : - ctobservation(CTMODEL_NAME, CTMODEL_VERSION, obs) + ctobservation(CTMODEL_NAME, VERSION, obs) { // Initialise members init_members(); @@ -90,7 +90,7 @@ ctmodel::ctmodel(const GObservations& obs) : * @param[in] argv Array of command line arguments. ***************************************************************************/ ctmodel::ctmodel(int argc, char *argv[]) : - ctobservation(CTMODEL_NAME, CTMODEL_VERSION, argc, argv) + ctobservation(CTMODEL_NAME, VERSION, argc, argv) { // Initialise members init_members(); diff --git a/src/ctmodel/ctmodel.hpp b/src/ctmodel/ctmodel.hpp index c7824f6..2c06f54 100644 --- a/src/ctmodel/ctmodel.hpp +++ b/src/ctmodel/ctmodel.hpp @@ -33,8 +33,7 @@ #include "ctobservation.hpp" /* __Definitions _________________________________________________________ */ -#define CTMODEL_NAME "ctmodel" -#define CTMODEL_VERSION "1.1.0" +#define CTMODEL_NAME "ctmodel" /***********************************************************************//** diff --git a/src/ctobssim/ctobssim.cpp b/src/ctobssim/ctobssim.cpp index 7e232a6..d66c8c3 100644 --- a/src/ctobssim/ctobssim.cpp +++ b/src/ctobssim/ctobssim.cpp @@ -67,7 +67,7 @@ const double g_roi_margin = 0.5; //!< Simulation radius margin (degrees) * * Constructs an empty ctobssim tool. ***************************************************************************/ -ctobssim::ctobssim(void) : ctobservation(CTOBSSIM_NAME, CTOBSSIM_VERSION) +ctobssim::ctobssim(void) : ctobservation(CTOBSSIM_NAME, VERSION) { // Initialise members init_members(); @@ -85,7 +85,7 @@ ctobssim::ctobssim(void) : ctobservation(CTOBSSIM_NAME, CTOBSSIM_VERSION) * Constructs ctobssim tool from an observation container. ***************************************************************************/ ctobssim::ctobssim(const GObservations& obs) : - ctobservation(CTOBSSIM_NAME, CTOBSSIM_VERSION, obs) + ctobservation(CTOBSSIM_NAME, VERSION, obs) { // Initialise members init_members(); @@ -105,7 +105,7 @@ ctobssim::ctobssim(const GObservations& obs) : * setting. ***************************************************************************/ ctobssim::ctobssim(int argc, char *argv[]) : - ctobservation(CTOBSSIM_NAME, CTOBSSIM_VERSION, argc, argv) + ctobservation(CTOBSSIM_NAME, VERSION, argc, argv) { // Initialise members init_members(); @@ -1471,7 +1471,8 @@ void ctobssim::simulate_background(GCTAObservation* obs, // Debug code: signal that we step into the response MC method #if defined(G_BACKGROUND_DEBUG) - std::cout << "ctobssim::simulate_background: model->mc in" << std::endl; + std::cout << "ctobssim::simulate_background: model->mc in" + << std::endl; #endif // Get simulated CTA event list. Note that this method @@ -1481,7 +1482,8 @@ void ctobssim::simulate_background(GCTAObservation* obs, // Debug code: signal that we came back #if defined(G_BACKGROUND_DEBUG) - std::cout << "ctobssim::simulate_background: model->mc out" << std::endl; + std::cout << "ctobssim::simulate_background: model->mc out" + << std::endl; #endif // Continue only if we got a CTA event list diff --git a/src/ctobssim/ctobssim.hpp b/src/ctobssim/ctobssim.hpp index 2dd0bd1..9982549 100644 --- a/src/ctobssim/ctobssim.hpp +++ b/src/ctobssim/ctobssim.hpp @@ -33,8 +33,7 @@ #include "ctobservation.hpp" /* __Definitions _________________________________________________________ */ -#define CTOBSSIM_NAME "ctobssim" -#define CTOBSSIM_VERSION "1.2.0" +#define CTOBSSIM_NAME "ctobssim" /***********************************************************************//** diff --git a/src/ctphase/Makefile.am b/src/ctphase/Makefile.am index a31e611..d95ec3d 100644 --- a/src/ctphase/Makefile.am +++ b/src/ctphase/Makefile.am @@ -1,7 +1,7 @@ ############################################################################# # Makefile.am for ctphase tool # # ------------------------------------------------------------------------- # -# Copyright (C) 2017 Leonardo Di Venere # +# Copyright (C) 2017 Joshua Cardenzana # # ------------------------------------------------------------------------- # # # # This program is free software: you can redistribute it and/or modify # diff --git a/src/ctphase/ctphase.cpp b/src/ctphase/ctphase.cpp index cac164e..a6c1e2c 100644 --- a/src/ctphase/ctphase.cpp +++ b/src/ctphase/ctphase.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctphase - Append phase information to CTA events file * * ----------------------------------------------------------------------- * - * copyright (C) 2017 by Leonardo Di Venere * + * copyright (C) 2017 by Joshua Cardenzana * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -21,7 +21,7 @@ /** * @file ctphase.cpp * @brief Append phase information to CTA events file - * @author Leonardo Di Venere + * @author Joshua Cardenzana */ /* __ Includes ___________________________________________________________ */ @@ -33,15 +33,12 @@ #include "GTools.hpp" /* __ Method name definitions ____________________________________________ */ -#define G_RUN "ctphase::run()" -#define G_PHASE_EVENTS "ctphase::phase_events(GCTAObservation*, "\ - "std::string&, std::string&, std::string&)" -#define G_SET_EBOUNDS "ctphase::set_ebounds(GCTAObservation*, GEbounds&)" +#define G_GET_PARAMETERS "ctphase::get_parameters()" +#define G_PHASE_EVENTS "ctphase::phase_events(GCTAObservation*)" /* __ Debug definitions __________________________________________________ */ /* __ Coding definitions _________________________________________________ */ -//#define G_USE_MKSTEMP //!< Use mkstemp for temporary filename /*========================================================================== @@ -53,7 +50,7 @@ /***********************************************************************//** * @brief Void constructor ***************************************************************************/ -ctphase::ctphase(void) : ctobservation(CTPHASE_NAME, CTPHASE_VERSION) +ctphase::ctphase(void) : ctobservation(CTPHASE_NAME, VERSION) { // Initialise members init_members(); @@ -72,7 +69,7 @@ ctphase::ctphase(void) : ctobservation(CTPHASE_NAME, CTPHASE_VERSION) * provided in an observation container. ***************************************************************************/ ctphase::ctphase(const GObservations& obs) : - ctobservation(CTPHASE_NAME, CTPHASE_VERSION, obs) + ctobservation(CTPHASE_NAME, VERSION, obs) { // Initialise members init_members(); @@ -91,7 +88,7 @@ ctphase::ctphase(const GObservations& obs) : * @param[in] argv Array of command line arguments. ***************************************************************************/ ctphase::ctphase(int argc, char *argv[]) : - ctobservation(CTPHASE_NAME, CTPHASE_VERSION, argc, argv) + ctobservation(CTPHASE_NAME, VERSION, argc, argv) { // Initialise members init_members(); @@ -203,10 +200,11 @@ void ctphase::clear(void) /***********************************************************************//** - * @brief Append phase information to event data + * @brief Runs the ctprob tool * * This method reads in the application parameters and loops over all - * observations that were found to append phase information to the data. + * unbinned observations to compute the phase information for each event + * based on an input model. ***************************************************************************/ void ctphase::run(void) { @@ -222,118 +220,22 @@ void ctphase::run(void) log_observations(NORMAL, m_obs, "Input observation"); // Write header into logger - log_header1(TERSE, "Appending phase informaiton"); + log_header1(TERSE, "Compute event phase"); // Initialise counters int n_observations = 0; - // Loop over all observation in the container - for (int i = 0; i < m_obs.size(); ++i) { - - // Write header for the current observation - log_header3(TERSE, get_obs_header(m_obs[i])); - - // Initialise event input and output filenames and the event - // and GTI extension names - m_infiles.push_back(""); - m_evtname.push_back(gammalib::extname_cta_events); - m_gtiname.push_back(gammalib::extname_gti); - - // Get CTA observation - GCTAObservation* obs = dynamic_cast(m_obs[i]); - - // Skip observation if it's not CTA - if (obs == NULL) { - std::string msg = " Skipping "+m_obs[i]->instrument()+ - " observation"; - log_string(NORMAL, msg); - continue; - } - - // Skip observation if we have a binned observation - if (obs->eventtype() == "CountsCube") { - std::string msg = " Skipping binned "+m_obs[i]->instrument()+ - " observation"; - log_string(NORMAL, msg); - continue; - } + // Loop over all unbinned CTA observations in the container + for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL; + obs = next_unbinned_observation()) { // Increment counter n_observations++; - // Save event file name (for possible saving) - m_infiles[i] = obs->eventfile(); - - // Extract event and GTI extension names from input FITS file - GFilename fname(m_infiles[i]); - if (fname.has_extname()) { - m_evtname[i] = fname.extname(); - } - m_gtiname[i] = get_gtiname(fname.url(), m_evtname[i]); - - // Write input file information into logger - log_value(NORMAL, "Input filename", m_infiles[i]); - log_value(NORMAL, "Event extension name", m_evtname[i]); - log_value(NORMAL, "GTI extension name", m_gtiname[i]); - - // Fall through in case that the event file is empty - if (obs->events()->size() == 0) { - log_string(NORMAL, " Warning: No events in event file \""+ - m_infiles[i]+"\". Event phasing skipped."); - continue; - } + // Compute event phase + phase_events(obs); - // If we have an input file then check it - if (!m_infiles[i].empty()) { - std::string message = check_infile(m_infiles[i], m_evtname[i]); - if (!message.empty()) { - throw GException::invalid_value(G_RUN, message); - } - } - - // Get temporary file name - #if G_USE_MKSTEMP - char tpl[] = "ctphaseXXXXXX"; - int fileid = mkstemp(tpl); - std::string filename(tpl); - #else - std::string filename = std::tmpnam(NULL); - #endif - - // Save observation in temporary file. We add here the events and - // GTI extension name so that the GCTAObservation::save method can - // use this information for writing the proper extension names into - // the temporary file - obs->save(filename+"["+m_evtname[i]+";"+ m_gtiname[i]+"]", true); - - // Log saved FITS file. - if (logExplicit()) { - GFits tmpfile(filename); - log.header3("FITS file content of temporary file"); - log << tmpfile << std::endl; - tmpfile.close(); - } - - // If we have a temporary file then check it - if (!filename.empty()) { - std::string message = check_infile(filename, m_evtname[i]); - if (!message.empty()) { - throw GException::invalid_value(G_RUN, message); - } - } - - // Load observation from temporary file and compute phase - phase_events(obs, filename, m_evtname[i], m_gtiname[i]); - - // Close temporary file - #if G_USE_MKSTEMP - close(fileid); - #endif - - // Remove temporary file - std::remove(filename.c_str()); - - } // endfor: looped over all observations + } // endfor: looped over observations // If more than a single observation has been handled then make sure that // an XML file will be used for storage @@ -380,12 +282,12 @@ void ctphase::save(void) // Case A: Save event file(s) and XML metadata information if (m_use_xml) { - save_xml(); + save_events_xml(); } // Case B: Save event file as FITS file else { - save_fits(); + save_events_fits(); } // Return @@ -459,15 +361,10 @@ void ctphase::publish(const std::string& name) void ctphase::init_members(void) { // Initialise parameters - m_outobs.clear(); - m_prefix.clear(); - m_phase.clear(); m_chatter = static_cast(2); // Initialise protected members - m_infiles.clear(); - m_evtname.clear(); - m_gtiname.clear(); + m_phase.clear(); // Return return; @@ -482,15 +379,10 @@ void ctphase::init_members(void) void ctphase::copy_members(const ctphase& app) { // Copy parameters - m_outobs = app.m_outobs; - m_prefix = app.m_prefix; - m_chatter = app.m_chatter; + m_chatter = app.m_chatter; // Copy protected members - m_infiles = app.m_infiles; - m_evtname = app.m_evtname; - m_gtiname = app.m_gtiname; - m_phase = app.m_phase; + m_phase = app.m_phase; // Return return; @@ -502,8 +394,6 @@ void ctphase::copy_members(const ctphase& app) ***************************************************************************/ void ctphase::free_members(void) { - m_phase.clear(); - // Return return; } @@ -526,381 +416,81 @@ void ctphase::get_parameters(void) // information and do not accept counts cubes. setup_observations(m_obs, false, true, false); - // Check whether a model XML file and source name were supplied - if (!((*this)["inmodel"].filename().file().empty()) && - !((*this)["source"].string().empty())) { - - // Get the information from an xml file - read_phase_info_from_xml(); - } - - // Check for sanity of frequency and phase information - else if ((*this)["p0"].is_valid() && - (*this)["mjd"].is_valid() && - (*this)["f0"].is_valid()&& - (*this)["f1"].is_valid()&& - (*this)["f2"].is_valid()) { - - // The following creates an XML object responsible for initializing - // 'm_phase'. This was done as the default values for the parameters in - // 'GModelTemporalPhaseCurve' are such that the ranges are predefined. - // Initializing 'm_phase' from an XML object circumvents this restriction. - - GXmlElement norm; - norm.name("parameter"); - norm.attribute("name", "Normalization"); - norm.attribute("value", "1"); - norm.attribute("min", "0"); - norm.attribute("max", "1000"); - - GXmlElement mjd; - mjd.name("parameter"); - mjd.attribute("name", "MJD"); - mjd.attribute("value", (*this)["mjd"].value()); - mjd.attribute("min", (*this)["mjd"].min()); - mjd.attribute("max", (*this)["mjd"].max()); - - GXmlElement p0; - p0.name("parameter"); - p0.attribute("name", "Phase"); - p0.attribute("value", (*this)["p0"].value()); - p0.attribute("min", (*this)["p0"].min()); - p0.attribute("max", (*this)["p0"].max()); - - GXmlElement f0; - f0.name("parameter"); - f0.attribute("name", "F0"); - f0.attribute("value", (*this)["f0"].value()); - f0.attribute("min", (*this)["f0"].min()); - f0.attribute("max", (*this)["f0"].max()); - - GXmlElement f1; - f1.name("parameter"); - f1.attribute("name", "F1"); - f1.attribute("value", (*this)["f1"].value()); - f1.attribute("min", (*this)["f1"].min()); - f1.attribute("max", (*this)["f1"].max()); - - GXmlElement f2; - f2.name("parameter"); - f2.attribute("name", "F2"); - f2.attribute("value", (*this)["f2"].value()); - f2.attribute("min", (*this)["f2"].min()); - f2.attribute("max", (*this)["f2"].max()); - - GXmlElement info; - info.append(norm); - info.append(mjd); - info.append(p0); - info.append(f0); - info.append(f1); - info.append(f2); - - info.attribute("file",""); - - m_phase = GModelTemporalPhaseCurve(); - m_phase.read(info); - } - else { - GException::invalid_value("ctphase::get_parameters()", - "Invalid parameter value for p0, mjd, f0, f1, or f2."); - } - - // Get other User parameters - m_chatter = static_cast((*this)["chatter"].integer()); - - // Optionally read ahead parameters so that they get correctly - // dumped into the log file - if (read_ahead()) { - m_outobs = (*this)["outobs"].filename(); - m_prefix = (*this)["prefix"].string(); - } - - // Write parameters into logger - log_parameters(TERSE); - - // Return - return; -} - -/************************************************************************* - * @brief Read phase information from XML file - * - * @exception GException::par_error - * Source does not have a temporal component with phase information. - ***************************************************************************/ -void ctphase::read_phase_info_from_xml(void) -{ - // Load the models from the xml file - GModels in_models((*this)["inmodel"].filename()); - - // Get the temporal component of this source - std::string src_name = (*this)["source"].string(); - GModelSky* model = dynamic_cast(in_models[src_name]); - GModelTemporalPhaseCurve* temporal_model = dynamic_cast(model->temporal()); - - // Make sure the model has a temporal component with phase info - if (temporal_model != NULL) { - - // Copy the phase model - m_phase = GModelTemporalPhaseCurve(*temporal_model); - - } else { - // Throw an error as there is no temporal component for this source - throw GException::par_error("ctphase::get_parameters(void)", - "source", - "\"" + src_name + "\" source model " + - "does not have a temporal component with "+ - "phase information!"); - } -} - - -/***********************************************************************//** - * @brief Define phase for events - * - * @param[in,out] obs CTA observation. - * @param[in] filename File name. - * @param[in] evtname Event extension name. - * @param[in] gtiname GTI extension name. - * - * @exception GException::invalid_value - * No events extension found in FITS file. - * - * Append phase information to events from a FITS file by adding a "PHASE" - * column to the FITS file. - ***************************************************************************/ -void ctphase::phase_events(GCTAObservation* obs, - const std::string& filename, - const std::string& evtname, - const std::string& gtiname) -{ - // Get CTA event list pointer - GCTAEventList* list = - static_cast(const_cast(obs->events())); - - GCTAEventAtom* event = (*list)[0]; - - // Print a warning into the log file if there are some precision concerns - double dt = std::abs(event->time() - m_phase.mjd()); - if (((m_phase.f1() * dt) > 1.0e8) || - (m_phase.f2() * dt * dt) > 1.0e8) { - log_string(NORMAL,"*** WARNING ***"); - log_string(NORMAL," Values supplied for reference MJD and/or f1 and/or f2"); - log_string(NORMAL," are large and may result in numerical precision issues."); - log_string(NORMAL," Consider using an ephemeris derived from a date closer"); - log_string(NORMAL," to the data being analyzed."); - log_value(NORMAL, " Observation ID", obs->obs_id()); - log_value(NORMAL, " Reference MJD", m_phase.mjd().mjd()); - log_value(NORMAL, " First event MJD", event->time().mjd()); - log_string(NORMAL,"*** WARNING ***"); - } - - for (int evnt=0; evntsize(); evnt++) { - // Get the next event - event = (*list)[evnt]; - - // Get the event time - GTime time = event->time(); - - // Set the event phase - event->phase(m_phase.phase(time)); - } - - // Make sure the file knows to save the phase information - if (!list->has_phase()) { - list->has_phase(true); - } - - // Return - return; -} - - -/***********************************************************************//** - * @brief Check input filename - * - * @param[in] filename File name. - * @param[in] evtname Event extension name. - * - * This method checks if the input FITS file is correct. - ***************************************************************************/ -std::string ctphase::check_infile(const std::string& filename, - const std::string& evtname) const -{ - // Initialise message string - std::string message = ""; - - // Open FITS file - GFits fits(filename); - - // Check for existence of events extensions - if (!fits.contains(evtname)) { - message = "No \""+evtname+"\" extension found in input file \""+ - filename + "\"."; - } - - // ... otherwise check column names - else { + // If a model definition XML file is provided the + GFilename inmodel = (*this)["inmodel"].filename(); + if (is_valid_filename(inmodel)) { - // Get pointer to FITS table - GFitsTable* table = fits.table(evtname); + // Load the models from the XML file + GModels models(inmodel); - // Initialise list of missing columns - std::vector missing; + // Get the source component + std::string source = (*this)["source"].string(); - // Check for existence of TIME column - if (!table->contains("TIME")) { - missing.push_back("TIME"); + // Throw an exception if source is not found in the model container + if (!models.contains(source)) { + std::string msg = "Source model \""+source+"\" not found in model " + "definition XML file \""+inmodel.url()+"\"."; + throw GException::invalid_value(G_GET_PARAMETERS, msg); } - // Check for existence of ENERGY column - if (!table->contains("ENERGY")) { - missing.push_back("ENERGY"); - } + // Get pointer to source model + GModelSky* sky = dynamic_cast(models[source]); - // Check for existence of RA column - if (!table->contains("RA")) { - missing.push_back("RA"); + // Throw an exception if model is not a sky model + if (sky == NULL) { + std::string msg = "Source model \""+source+"\" is not a sky model."; + throw GException::invalid_value(G_GET_PARAMETERS, msg); } - // Check for existence of DEC column - if (!table->contains("DEC")) { - missing.push_back("DEC"); - } + // Get pointer to temporal phase curve + GModelTemporalPhaseCurve* phasecurve = + dynamic_cast(sky->temporal()); - // Set error message for missing columns - if (!missing.empty()) { - message = "The following columns are missing in the " - "\""+evtname+"\" extension of input file \""+ - filename + "\": "; - for (int i = 0; i < missing.size(); ++i) { - message += "\"" + missing[i] + "\""; - if (i < missing.size()-1) { - message += ", "; - } - } + // Throw an exception if model has not a temporal phase curve + // component + if (phasecurve == NULL) { + std::string msg = "Source model \""+source+"\" does not have a " + "temporal component with phase information."; + throw GException::invalid_value(G_GET_PARAMETERS, msg); } - } // endelse: checked column names - - // Return - return message; -} - - -/***********************************************************************//** - * @brief Set output file name. - * - * @param[in] filename Input file name. - * - * Converts an input file name into an output filename by prepending the - * prefix stored in the member m_prefix to the input file name. Any path as - * well as extension will be stripped from the input file name. Also a - * trailing ".gz" will be stripped as one cannot write into gzipped files. - ***************************************************************************/ -std::string ctphase::set_outfile_name(const std::string& filename) const -{ - // Create filename - GFilename fname(filename); - - // Split input filename without any extensions into path elements - std::vector elements = gammalib::split(fname.url(), "/"); - - // The last path element is the filename - std::string outname = m_prefix + elements[elements.size()-1]; - - // Strip any ".gz" - outname = gammalib::strip_chars(outname, ".gz"); - - // Return output filename - return outname; -} - - -/***********************************************************************//** - * @brief Get Good Time Intervals extension name - * - * @param[in] filename Input file name. - * @param[in] evtname Events extension name. - * - * Extracts the Good Time Intervals extension name from the event file. We - * do this by loading the events and accessing the Good Time Intervals - * extension name using the GCTAEventList::gtiname() method. If the file name - * is empty, the method returns `GTI`. - ***************************************************************************/ -std::string ctphase::get_gtiname(const std::string& filename, - const std::string& evtname) const -{ - // Initialise GTI name - std::string gtiname = gammalib::extname_gti; - - // Continue only if the filename is not empty - if (!filename.empty()) { + // Set phase curve + m_phase = GModelTemporalPhaseCurve(*phasecurve); - // Load events - GCTAEventList events(filename+"["+evtname+"]"); + } // endif: model definition XML file provided - // Get GTI name - gtiname = events.gtiname(); - - } // endif: filename was not empty - - // Return GTI name - return (gtiname); -} - - -/***********************************************************************//** - * @brief Save event list in FITS format. - * - * Save the event list as a FITS file. The file name of the FITS file is - * specified by the "outobs" parameter. - ***************************************************************************/ -void ctphase::save_fits(void) -{ - // Save only if there are observations - if (m_obs.size() > 0) { - - // Get output filename - m_outobs = (*this)["outobs"].filename(); - - // Get CTA observation from observation container - GCTAObservation* obs = dynamic_cast(m_obs[0]); - - // Save only if it's a CTA observation - if (obs != NULL) { - - // Save only if file name is non-empty - if (m_infiles[0].length() > 0) { + // Otherwise read phase curve information from the parameter files + else { - // Create file name object - GFilename fname(m_outobs); + // Read phase curve user parameters + GFilename phasecurve = (*this)["phasecurve"].filename(); + double mjd_value = (*this)["mjd"].real(); + double phase = (*this)["phase"].real(); + double f0 = (*this)["f0"].real(); + double f1 = (*this)["f1"].real(); + double f2 = (*this)["f2"].real(); - // Extract filename and event extension name - std::string outfile = fname.url(); + // Set reference time + GTime mjd; + mjd.mjd(mjd_value); - // Append event extension name. We handle here the possibility - // to write the events into a different extension. - if (fname.has_extname()) { - outfile += "["+fname.extname()+"]"; - } - else { - outfile += "["+m_evtname[0]+"]"; - } + // Set phase curve + m_phase = GModelTemporalPhaseCurve(phasecurve, mjd, phase, f0, f1, f2); - // Log filename - log_value(NORMAL, "Event list file", outfile); + } // endelse: read phase curve information from User parameters - // Save event list - save_event_list(obs, m_infiles[0], m_evtname[0], m_gtiname[0], - outfile); - - } // endif: filename was non empty + // Get other User parameters + m_chatter = static_cast((*this)["chatter"].integer()); - } // endif: observation was CTA observation + // Optionally read ahead parameters so that they get correctly + // dumped into the log file + if (read_ahead()) { + (*this)["outobs"].filename(); + (*this)["prefix"].string(); + } - } // endif: there were observations + // Write parameters into logger + log_parameters(TERSE); // Return return; @@ -908,169 +498,68 @@ void ctphase::save_fits(void) /***********************************************************************//** - * @brief Save event list(s) in XML format. + * @brief Compute event phase for an observation * - * Save the event list(s) into FITS files and write the file path information - * into a XML file. The filename of the XML file is specified by the outfile - * parameter, the filename(s) of the event lists are built by prepending a - * prefix to the input event list filenames. Any path present in the input - * filename will be stripped, i.e. the event list(s) will be written in the - * local working directory (unless a path is specified in the prefix). + * @param[in,out] obs CTA observation. + * + * @exception GException::invalid_value + * No events extension found in FITS file. + * + * Append phase information to events from a FITS file by adding a "PHASE" + * column to the FITS file. ***************************************************************************/ -void ctphase::save_xml(void) +void ctphase::phase_events(GCTAObservation* obs) { - // Get output filename and prefix - m_outobs = (*this)["outobs"].filename(); - m_prefix = (*this)["prefix"].string(); - - // Issue warning if output filename has no .xml suffix - log_string(TERSE, warn_xml_suffix(m_outobs)); - - // Loop over all observation in the container - for (int i = 0; i < m_obs.size(); ++i) { - - // Get CTA observation - GCTAObservation* obs = dynamic_cast(m_obs[i]); - - // Skip observations that are no CTA observations - if (obs == NULL) { - continue; - } + // Write header into logger + log_header3(NORMAL, "Compute event phases"); + + // Make sure that the observation holds a CTA event list. If this + // is not the case then throw an exception. + GCTAEventList* events = dynamic_cast(const_cast + (obs->events())); + if (events == NULL) { + std::string msg = "CTA Observation does not contain an event list. An " + "event list is needed to compute the event phases."; + throw GException::invalid_value(G_PHASE_EVENTS, msg); + } - // Skip observations that have empty names - if (m_infiles[i].length() == 0) { - continue; + // Continue only if there are events + if (events->size() > 0) { + + // Print a warning into the log file if there are some precision concerns + GCTAEventAtom* event = (*events)[0]; + double dt = std::abs(event->time() - m_phase.mjd()); + if (((m_phase.f1() * dt) > 1.0e8) || + (m_phase.f2() * dt * dt) > 1.0e8) { + log_string(NORMAL,"*** WARNING ***"); + log_string(NORMAL," Values supplied for reference MJD and/or f1 and/or f2"); + log_string(NORMAL," are large and may result in numerical precision issues."); + log_string(NORMAL," Consider using an ephemeris derived from a date closer"); + log_string(NORMAL," to the data being analyzed."); + log_value(NORMAL, " Observation ID", obs->obs_id()); + log_value(NORMAL, " Reference MJD", m_phase.mjd().mjd()); + log_value(NORMAL, " First event MJD", event->time().mjd()); + log_string(NORMAL,"*** WARNING ***"); } - // Set event output file name - std::string outfile = set_outfile_name(m_infiles[i]); - - // Append event extension name - outfile += "["+m_evtname[i]+"]"; - - // Log filename - log_value(NORMAL, "Event list file", outfile); - - // Store output file name in observation - obs->eventfile(outfile); - - // Save event list - save_event_list(obs, m_infiles[i], m_evtname[i], m_gtiname[i], - outfile); - - } // endfor: looped over observations - - // Write observation definition XML file name into logger - log_value(NORMAL, "Obs. definition file", m_outobs); - - // Save observations in XML file - m_obs.save(m_outobs); - - // Return - return; -} - + // Compute the phase for all events + for (int i = 0; i < events->size(); ++i) { -/***********************************************************************//** - * @brief Save event list into FITS file - * - * @param[in] obs Pointer to CTA observation. - * @param[in] infile Input file name. - * @param[in] evtname Event extension name. - * @param[in] gtiname GTI extension name. - * @param[in] outfile Output file name. - * - * Saves an event list and the corresponding Good Time Intervals into a FITS - * file and copy all others extensions from the input file to the output - * file. - * - * If an extension name is specified in the @p outfile argument, the events - * and eventually also the Good Time Intervals will be extracted from the - * argument and used for writing the events. The format is - * - * [] - * - * where needs to be replaced by the name of the FITS file, - * and by the name of the events - * and Good Time Intervals extensions. For example - * - * myfits.fits[EVENTS1;GTI1] - * - * will write the updated events into the "EVENTS1" extension and the - * Good Time Intervals into the "GTI1" extension of the "myfits.fits" FITS - * file. If the Good Time Intervals extension name is skipped, e.g. - * - * myfits.fits[EVENTS1] - * - * the original extension name for the Good Time Intervals will be kept. - * Analogously, only the Good Time Intervals extension name can be changed - * by specifying - * - * myfits.fits[;GTI1] - * - * In none of the cases will the original events and Good Time Intervals be - * copied over to the output file. - ***************************************************************************/ -void ctphase::save_event_list(const GCTAObservation* obs, - const std::string& infile, - const std::string& evtname, - const std::string& gtiname, - const std::string& outfile) const -{ - // Save only if we have an event list - if (obs->eventtype() == "EventList") { - - // Set output FITS file event extension names - GFilename outname(outfile); - std::string outevt = evtname; - std::string outgti = gtiname; - if (outname.has_extname()) { - std::vector extnames = - gammalib::split(outname.extname(), ";"); - if (extnames.size() > 0) { - std::string extname = gammalib::strip_whitespace(extnames[0]); - if (!extname.empty()) { - outevt = extname; - } - } - if (extnames.size() > 1) { - std::string extname = gammalib::strip_whitespace(extnames[1]); - if (!extname.empty()) { - outgti = extname; - } - } - } + // Get pointer to the event + GCTAEventAtom* event = (*events)[i]; - // Create output FITS file - GFits outfits; - - // Write observation into FITS file - obs->write(outfits, outevt, outgti); - - // Copy all extensions other than evtname and gtiname extensions - // from the input to the output event list. The evtname and - // gtiname extensions are written by the save method, all others - // that may eventually be present have to be copied over - // explicitly. - GFits infits(infile); - for (int extno = 1; extno < infits.size(); ++extno) { - GFitsHDU* hdu = infits.at(extno); - if (hdu->extname() != evtname && - hdu->extname() != gtiname && - hdu->extname() != outevt && - hdu->extname() != outgti) { - outfits.append(*hdu); - } - } - - // Close input file - infits.close(); + // Get the event time + GTime time = event->time(); + + // Set the event phase + event->phase(m_phase.phase(time)); - // Save file to disk and close it (we need both operations) - outfits.saveto(outname.url(), clobber()); - outfits.close(); + } // endfor: looped over events - } // endif: observation was unbinned + } // endif: there were events + + // Signal that phase information is available + events->has_phase(true); // Return return; diff --git a/src/ctphase/ctphase.hpp b/src/ctphase/ctphase.hpp index 798125c..8387bb5 100644 --- a/src/ctphase/ctphase.hpp +++ b/src/ctphase/ctphase.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctphase - Append phase information to CTA events file * * ----------------------------------------------------------------------- * - * copyright (C) 2017 by Leonardo Di Venere * + * copyright (C) 2017 by Joshua Cardenzana * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -21,7 +21,7 @@ /** * @file ctphase.hpp * @brief Append phase information to CTA events file - * @author Leonardo Di Venere + * @author Joshua Cardenzana */ #ifndef CTPHASE_HPP @@ -35,8 +35,7 @@ #include "ctobservation.hpp" /* __Definitions _________________________________________________________ */ -#define CTPHASE_NAME "ctphase" -#define CTPHASE_VERSION "0.0.1" +#define CTPHASE_NAME "ctphase" /***********************************************************************//** @@ -65,38 +64,17 @@ public: protected: // Protected methods - void init_members(void); - void copy_members(const ctphase& app); - void free_members(void); - void get_parameters(void); - void read_phase_info_from_xml(void); - void phase_events(GCTAObservation* obs, - const std::string& filename, - const std::string& evtname, - const std::string& gtiname); - std::string check_infile(const std::string& filename, - const std::string& evtname) const; - std::string set_outfile_name(const std::string& filename) const; - std::string get_gtiname(const std::string& filename, - const std::string& evtname) const; - void save_fits(void); - void save_xml(void); - void save_event_list(const GCTAObservation* obs, - const std::string& infile, - const std::string& evtname, - const std::string& gtiname, - const std::string& outfile) const; + void init_members(void); + void copy_members(const ctphase& app); + void free_members(void); + void get_parameters(void); + void phase_events(GCTAObservation* obs); // User parameters - std::string m_outobs; //!< Output event list or XML file - std::string m_prefix; //!< Prefix for multiple event lists - GModelTemporalPhaseCurve m_phase; //!< Phase information object - GChatter m_chatter; //!< Chattiness + GChatter m_chatter; //!< Chattiness // Protected members - std::vector m_infiles; //!< Input event filenames - std::vector m_evtname; //!< Event extension names - std::vector m_gtiname; //!< GTI extension names + GModelTemporalPhaseCurve m_phase; //!< Phase curve model }; diff --git a/src/ctphase/ctphase.par b/src/ctphase/ctphase.par index 393dddb..49f1c83 100644 --- a/src/ctphase/ctphase.par +++ b/src/ctphase/ctphase.par @@ -5,11 +5,10 @@ # ------------------------------------------------------------------------- # # # # File: ctphase.par # -# Version: 0.0.1 # # # -# Author: Leonardo Di Venere # -# leonardo.divenere@ba.infn.it # -# INFN # +# Author: Joshua Cardenzana # +# jcardenzana@irap.omp.eu # +# IRAP # # # # Purpose: Parameter file for the ctphase application # # # @@ -17,20 +16,21 @@ # # General parameters #=================== -inobs, f, a,events.fits,,, "Input event list or observation definition XML file" -outobs, f, a,phased_events.fits,,, "Output event list or observation definition XML file" -prefix, s, h, "phased_",,, "Prefix for event lists in observation definition XML file" -inmodel, f, h, ,,, "Model file defining the phase information" -source, s, h, ,,, "Source in model file to use for phase information" +inobs, f, a, "events.fits",,, "Input event list or observation definition XML file" +outobs, f, a, "phased_events.fits",,, "Output event list or observation definition XML file" +prefix, s, h, "phased_",,, "Prefix for event lists in observation definition XML file" +inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Model file defining the phase information" +source, s, a, "Crab",,, "Source in model file to use for phase information" # # Phase information parameters #===================== -mjd, r, a, 51544.5,,, "Reference MJD for phase information" -p0, r, a, 0.0,0.0,1.0, "Phase at reference MJD" -f0, r, a, 1,,, "Frequency at reference MJD [Hz]" -f1, r, a, 0.0,,, "First freq. derivative at reference MJD [s^-1]" -f2, r, a, 0.0,,, "Second freq. derivative at reference MJD [s^-2]" +phasecurve, f, a, phase.fits,,, "Phase curve FITS file" +mjd, r, a, 51544.5,,, "Reference MJD for phase information" +phase, r, a, 0.0,0.0,1.0, "Phase at reference MJD" +f0, r, a, 1.0,,, "Frequency at reference MJD (Hz)" +f1, r, a, 0.1,,, "First freq. derivative at reference MJD (s^-1)" +f2, r, a, 0.01,,, "Second freq. derivative at reference MJD (s^-2)" # # Standard parameters diff --git a/src/ctphase/main.cpp b/src/ctphase/main.cpp index c0761b6..3644d39 100644 --- a/src/ctphase/main.cpp +++ b/src/ctphase/main.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctphase - Append phase information to CTA events file * * ----------------------------------------------------------------------- * - * copyright (C) 2017 by Leonardo Di Venere * + * copyright (C) 2017 by Joshua Cardenzana * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -21,7 +21,7 @@ /** * @file ctphase/main.cpp * @brief Append phase information to CTA events file - * @author Leonardo Di Venere + * @author Joshua Cardenzana */ /* __ Includes ___________________________________________________________ */ diff --git a/src/ctprob/Makefile.am b/src/ctprob/Makefile.am new file mode 100644 index 0000000..1ef279c --- /dev/null +++ b/src/ctprob/Makefile.am @@ -0,0 +1,47 @@ +############################################################################# +# Makefile.am for ctprob tool # +# ------------------------------------------------------------------------- # +# Copyright (C) 2017 Leonardo Di Venere # +# ------------------------------------------------------------------------- # +# # +# 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 # +# the Free Software Foundation, either version 3 of the License, or # +# (at your option) any later version. # +# # +# This program is distributed in the hope that it will be useful, # +# but WITHOUT ANY WARRANTY; without even the implied warranty of # +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # +# GNU General Public License for more details. # +# # +# You should have received a copy of the GNU General Public License # +# along with this program. If not, see . # +# # +############################################################################# + +# Set include directories +INCLUDES = -I$(top_srcdir)/src/support + +# Program name +bin_PROGRAMS = ctprob + +# Sources +sources = ctprob.cpp ctprob.hpp + +# Libtool library +noinst_LTLIBRARIES = libctprob.la +libctprob_la_SOURCES = $(sources) +libctprob_la_LDFLAGS = @OPENMP_LIB_LDFLAGS@ + +# Program definition +ctprob_SOURCES = main.cpp +ctprob_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@ +ctprob_CPPFLAGS = @CPPFLAGS@ +ctprob_LDADD = libctprob.la $(top_builddir)/src/support/libsupport.la + +# Parameter files +pfilesdir = $(prefix)/syspfiles +pfiles_DATA = $(srcdir)/ctprob.par + +# Extra files to be distributed in the tar file +EXTRA_DIST = $(pfiles_DATA) diff --git a/src/ctprob/ctprob.cpp b/src/ctprob/ctprob.cpp new file mode 100644 index 0000000..1f313be --- /dev/null +++ b/src/ctprob/ctprob.cpp @@ -0,0 +1,784 @@ +/*************************************************************************** + * ctprob - Computes event probability for a given model * + * ----------------------------------------------------------------------- * + * copyright (C) 2017 by Leonardo Di Venere * + * ----------------------------------------------------------------------- * + * * + * 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 * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file ctprob.cpp + * @brief Computes event probability for a given model + * @author Leonardo Di Venere + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#include "ctprob.hpp" +#include "GTools.hpp" + +/* __ Method name definitions ____________________________________________ */ +#define G_GET_OBS "ctprob::get_obs()" +#define G_EVALUATE_PROB "ctprob::evaluate_probability(GCTAObservation*)" + +/* __ Debug definitions __________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + ***************************************************************************/ +ctprob::ctprob(void) : ctobservation(CTPROB_NAME, VERSION) +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Observations constructor + * + * param[in] obs Observation container. + * + * Creates an instance of the class that is initialised using the information + * provided in an observation container. + ***************************************************************************/ +ctprob::ctprob(const GObservations& obs) : + ctobservation(CTPROB_NAME, VERSION, obs) +{ + // Initialise members + init_members(); + + // Return + return; +} + + + + +/***********************************************************************//** + * @brief Command line constructor + * + * @param[in] argc Number of arguments in command line. + * @param[in] argv Array of command line arguments. + ***************************************************************************/ +ctprob::ctprob(int argc, char *argv[]) : + ctobservation(CTPROB_NAME, VERSION, argc, argv) +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] app Application. + ***************************************************************************/ +ctprob::ctprob(const ctprob& app) : ctobservation(app) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(app); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +ctprob::~ctprob(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] app Application. + * @return Application. + ***************************************************************************/ +ctprob& ctprob::operator=(const ctprob& app) +{ + // Execute only if object is not identical + if (this != &app) { + + // Copy base class members + this->ctobservation::operator=(app); + + // Free members + free_members(); + + // Initialise members + init_members(); + + // Copy members + copy_members(app); + + } // endif: object was not identical + + // Return this object + return *this; +} + + +/*========================================================================== + = = + = Public methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Clear ctprob tool + * + * Clears ctprob tool. + ***************************************************************************/ +void ctprob::clear(void) +{ + // Free members + free_members(); + this->ctobservation::free_members(); + this->ctool::free_members(); + + // Clear base class (needed to conserve tool name and version) + this->GApplication::clear(); + + // Initialise members + this->ctool::init_members(); + this->ctobservation::init_members(); + init_members(); + + // Write header into logger + log_header(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Runs the ctprob tool + * + * This method reads in the application parameters and loops over all + * unbinned observations to compute the probability for each event that it + * comes from one of the model components. + ***************************************************************************/ +void ctprob::run(void) +{ + // Switch screen logging on in debug mode + if (logDebug()) { + log.cout(true); + } + + // Get parameters + get_parameters(); + + // Set energy dispersion flags of all CTA observations and save old + // values in save_edisp vector + std::vector save_edisp = set_edisp(m_obs, m_apply_edisp); + + // Write input observation container into logger + log_observations(NORMAL, m_obs, "Input observation"); + + // Write header into logger + log_header1(TERSE, "Compute event probabilities"); + + // Initialise counters + int n_observations = 0; + + // Loop over all unbinned CTA observations in the container + for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL; + obs = next_unbinned_observation()) { + + // Increment counter + n_observations++; + + // Evaluate event probabilities + evaluate_probability(obs); + + } // endfor: looped over observations + + // If more than a single observation has been handled then make sure that + // an XML file will be used for storage + if (n_observations > 1) { + m_use_xml = true; + } + + // Write observation(s) into logger + log_observations(NORMAL, m_obs, "Output observation"); + + // Restore energy dispersion flags of all CTA observations + restore_edisp(m_obs, save_edisp); + + // Optionally publish event list(s) + if ((*this)["publish"].boolean()) { + publish(); + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Save the selected event list(s) + * + * This method saves the selected event list(s) into FITS file(s). There are + * two modes, depending on the m_use_xml flag. + * + * If m_use_xml is true, all selected event list(s) will be saved into FITS + * files, where the output filenames are constructued from the input + * filenames by prepending the m_prefix string to name. Any path information + * will be stripped form the input name, hence event files will be written + * into the local working directory (unless some path information is present + * in the prefix). In addition, an XML file will be created that gathers + * the filename information for the selected event list(s). If an XML file + * was present on input, all metadata information will be copied from this + * input file. + * + * If m_use_xml is false, the selected event list will be saved into a FITS + * file. + ***************************************************************************/ +void ctprob::save(void) +{ + // Write header into logger + log_header1(TERSE, gammalib::number("Save event list", m_obs.size())); + + // Case A: Save event file(s) and XML metadata information + if (m_use_xml) { + save_xml(); + } + + // Case B: Save event file as FITS file + else { + save_fits(); + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Publish event lists + * + * @param[in] name Event list name. + ***************************************************************************/ +void ctprob::publish(const std::string& name) +{ + // Write header into logger + log_header1(TERSE, gammalib::number("Publish event list", m_obs.size())); + + // Loop over all observation in the container + for (int i = 0; i < m_obs.size(); ++i) { + + // Get CTA observation + GCTAObservation* obs = dynamic_cast(m_obs[i]); + + // Handle only CTA observations + if (obs != NULL) { + + // Continue only if there is an event list + if (obs->events()->size() != 0) { + + // Set default name if user name is empty + std::string user_name(name); + if (user_name.empty()) { + user_name = CTPROB_NAME; + } + + // If there are several event lists then add an index + if (m_use_xml) { + user_name += gammalib::str(i); + } + + // Write event list name into logger + log_value(NORMAL, "Event list name", user_name); + + // Write events into in-memory FITS file + GFits fits; + obs->write(fits); + + // Publish + fits.publish(gammalib::extname_cta_events, user_name); + + } // endif: there were events + + } // endif: observation was a CTA observation + + } // endfor: looped over observations + + // Return + return; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void ctprob::init_members(void) +{ + // Initialise parameters + m_outobs.clear(); + m_prefix.clear(); + m_chatter = static_cast(2); + + // Initialise protected members + m_infiles.clear(); + m_evtname.clear(); + m_gtiname.clear(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] app Application. + ***************************************************************************/ +void ctprob::copy_members(const ctprob& app) +{ + // Copy parameters + m_outobs = app.m_outobs; + m_prefix = app.m_prefix; + m_chatter = app.m_chatter; + + // Copy protected members + m_infiles = app.m_infiles; + m_evtname = app.m_evtname; + m_gtiname = app.m_gtiname; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void ctprob::free_members(void) +{ + // Return + return; +} + + +/***********************************************************************//** + * @brief Get application parameters + * + * Get all user parameters from parameter file or (if required) by querying + * the user. Times are assumed to be in the native CTA MJD format. + * + * This method also loads observations if no observations are yet allocated. + * Observations are either loaded from a single CTA even list, or from a + * XML file using the metadata information that is stored in that file. + ***************************************************************************/ +void ctprob::get_parameters(void) +{ + // Setup observations from "inobs" parameter. Do not request response + // information and do not accept counts cubes. + if (m_obs.size() == 0) { + get_obs(); + } + // ... otherwise add response information and energy boundaries in case + // that they are missing + else { + setup_observations(m_obs); + } + + // Read model definition file if required + if (m_obs.models().size() == 0) { + + // Get model filename + std::string inmodel = (*this)["inmodel"].filename(); + + // Load models from file + m_obs.models(inmodel); + + } // endif: there were no models + + // Get energy dispersion flag parameters + m_apply_edisp = (*this)["edisp"].boolean(); + + + // Get remaining parameters + m_publish = (*this)["publish"].boolean(); + m_chatter = static_cast((*this)["chatter"].integer()); + + // Optionally read ahead parameters so that they get correctly + // dumped into the log file + if (read_ahead()) { + m_outobs = (*this)["outobs"].filename(); + m_prefix = (*this)["prefix"].string(); + } + + // Write parameters into logger + log_parameters(TERSE); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Get observation container + * + * @exception GException::invalid_value + * Invalid FITS file. + * + * Get an observation container according to the user parameters. The method + * supports loading of a individual FITS file or an observation definition + * file in XML format. + ***************************************************************************/ +void ctprob::get_obs(void) +{ + // Get the filename from the input parameters + std::string filename = (*this)["inobs"].filename(); + + // If no observation definition file throw an exception + if (!is_valid_filename(filename)) { + throw GException::invalid_value(G_GET_OBS, "Input event list is not valid."); + } // endif: filename was not valid + + // ... otherwise we have a file name + else { + + // If file is a FITS file then create an empty CTA observation + // and load file into observation + if (GFilename(filename).is_fits()) { + + // Allocate empty CTA observation + GCTAObservation cta; + + // Load data + cta.load(filename); + + // Set response + set_obs_response(&cta); + + // Append observation to container + m_obs.append(cta); + + // Signal that no XML file should be used for storage + m_use_xml = false; + + } + + // ... otherwise load file into observation container + else { + + // Load observations from XML file + m_obs.load(filename); + + // For all observations that have no response, set the response + // from the task parameters + set_response(m_obs); + + // Set observation boundary parameters (emin, emax, rad) + set_obs_bounds(m_obs); + + // Signal that XML file should be used for storage + m_use_xml = true; + + } // endelse: file was an XML file + + } + + // Return + return; + +} + + + +/***********************************************************************//** + * @brief Evaluate probability for events + * + * @param[in,out] obs CTA observation. + * + * @exception GException::invalid_value + * Observation does not contain an event list. + * + * Computes for all events the probability that an event arises from a + * specific model component. The computation is done by by evaluating the + * differential event probability for a given model and normalizing these + * probabilities so that the sum of the probabilities for all model + * components is unity for each event. + * + * The method appends for each model component a single precision column to + * the event list that contains the event probabilities. The name of the + * column is build from the model name prefixed with `PROB_`. + ***************************************************************************/ +void ctprob::evaluate_probability(GCTAObservation* obs) +{ + // Write header into logger + log_header3(NORMAL, "Evaluate events probabilities"); + + // Make sure that the observation holds a CTA event list. If this + // is not the case then throw an exception. + GCTAEventList* events = dynamic_cast(const_cast + (obs->events())); + if (events == NULL) { + std::string msg = "CTA Observation does not contain an event list. An " + "event list is needed to compute the event " + "probabilities."; + throw GException::invalid_value(G_EVALUATE_PROB, msg); + } + + // Continue only if there are model components + int nmodels = m_obs.models().size(); + if (nmodels > 0) { + + // Allocate single precision FITS table columns for all model components + std::vector columns; + for (int k = 0; k < nmodels; ++k) { + const std::string colname = "PROB_" + m_obs.models()[k] ->name(); + GFitsTableFloatCol* col = + new GFitsTableFloatCol(colname, events->size()); + columns.push_back(col); + } + + // Loop over all events in the event list + for (int i = 0; i < events->size(); ++i) { + + // Initialise total probability + double total = 0.0; + + // Get point to event + const GEvent* event = (*events)[i]; + + // Initialise probability vector + std::vector probabilities(nmodels, 0.0); + + // Loop over models + for (int k = 0; k < nmodels; ++k) { + double value = m_obs.models()[k]->eval(*event, *obs); + probabilities[k] = value; + total += value; + } + + // Normalize probabilities to unity + if (total > 0.0) { + for (int k = 0; k < nmodels; ++k) { + probabilities[k] /= total; + } + } + + // Store probabilities in FITS columns + for (int k = 0; k < nmodels; ++k) { + (*(columns[k]))(i) = probabilities[k]; + } + + } // endfor: looped over all events + + // Append columns to event list + for (int k = 0; k < nmodels; ++k) { + events->append_column(*columns[k]); + } + + } // endif: there were model components + + // Return + return; +} + + +/***********************************************************************//** + * @brief Set output file name. + * + * @param[in] filename Input file name. + * + * Converts an input file name into an output filename by prepending the + * prefix stored in the member m_prefix to the input file name. Any path as + * well as extension will be stripped from the input file name. Also a + * trailing ".gz" will be stripped as one cannot write into gzipped files. + ***************************************************************************/ +std::string ctprob::set_outfile_name(const std::string& filename) const +{ + // Create filename + GFilename fname(filename); + + // Split input filename without any extensions into path elements + std::vector elements = gammalib::split(fname.url(), "/"); + + // The last path element is the filename + std::string outname = m_prefix + elements[elements.size()-1]; + + // Strip any ".gz" + outname = gammalib::strip_chars(outname, ".gz"); + + // Return output filename + return outname; +} + + +/***********************************************************************//** + * @brief Save event list in FITS format. + * + * Save the event list as a FITS file. The file name of the FITS file is + * specified by the `outobs` parameter. + ***************************************************************************/ +void ctprob::save_fits(void) +{ + // Save only if there are observations + if (m_obs.size() > 0) { + + // Get output filename + m_outobs = (*this)["outobs"].filename(); + + // Loop over all unbinned CTA observations in the container + for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL; + obs = next_unbinned_observation()) { + + // Get input file name and default extension names + std::string infile = obs->eventfile(); + std::string evtname = gammalib::extname_cta_events; + std::string gtiname = gammalib::extname_gti; + + // Create file name object + GFilename fname(m_outobs); + + // Extract filename and event extension name + std::string outfile = fname.url(); + + // Append event extension name. We handle here the possibility + // to write the events into a different extension. + if (fname.has_extname()) { + outfile += "["+fname.extname()+"]"; + } + else { + outfile += "["+gammalib::extname_cta_events+"]"; + } + + // Log filename + log_value(NORMAL, "Event list file", outfile); + + // Save event list + save_event_list(obs, infile, evtname, gtiname, outfile); + + // Exit the loop (if this method is called we should only + // have a single unbinned observation in the observation + // container + break; + + } // endfor: looped over all unbinned observations + + } // endif: there were observations + + // Return + return; +} + + +/***********************************************************************//** + * @brief Save event list(s) in XML format. + * + * Save the event list(s) into FITS files and write the file path information + * into a XML file. The filename of the XML file is specified by the `outobs` + * parameter, the filename(s) of the event lists are built by prepending a + * prefix to the input event list filenames. Any path present in the input + * filename will be stripped, i.e. the event list(s) will be written in the + * local working directory (unless a path is specified in the prefix). + ***************************************************************************/ +void ctprob::save_xml(void) +{ + // Get output filename and prefix + m_outobs = (*this)["outobs"].filename(); + m_prefix = (*this)["prefix"].string(); + + // Issue warning if output filename has no .xml suffix + log_string(TERSE, warn_xml_suffix(m_outobs)); + + // Loop over all unbinned CTA observations in the container + for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL; + obs = next_unbinned_observation()) { + + // Get input file name and default extension names + std::string infile = obs->eventfile(); + std::string evtname = gammalib::extname_cta_events; + std::string gtiname = gammalib::extname_gti; + + // Extract event and GTI extension names from input FITS file + GFilename fname(infile); + if (fname.has_extname()) { + evtname = fname.extname(); + } + gtiname = get_gtiname(fname.url(), evtname); + + // Set event output file name + std::string outfile = set_outfile_name(infile); + + // Append event extension name + outfile += "["+evtname+"]"; + + // Log filename + log_value(NORMAL, "Event list file", outfile); + + // Store output file name in observation + obs->eventfile(outfile); + + // Save event list + save_event_list(obs, infile, evtname, gtiname, outfile); + + } // endfor: looped over observations + + // Write observation definition XML file name into logger + log_value(NORMAL, "Obs. definition file", m_outobs); + + // Save observations in XML file + m_obs.save(m_outobs); + + // Return + return; +} diff --git a/src/ctprob/ctprob.hpp b/src/ctprob/ctprob.hpp new file mode 100644 index 0000000..607d2e8 --- /dev/null +++ b/src/ctprob/ctprob.hpp @@ -0,0 +1,90 @@ +/*************************************************************************** + * ctprob - Computes event probability for a given model * + * ----------------------------------------------------------------------- * + * copyright (C) 2017 by Leonardo Di Venere * + * ----------------------------------------------------------------------- * + * * + * 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 * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file ctprob.hpp + * @brief Computes event probability for a given model + * @author Leonardo Di Venere + */ + +#ifndef CTPROB_HPP +#define CTPROB_HPP + +/* __ Includes ___________________________________________________________ */ +#include +#include +#include "GammaLib.hpp" +#include "GCTALib.hpp" +#include "ctobservation.hpp" + +/* __Definitions _________________________________________________________ */ +#define CTPROB_NAME "ctprob" + + +/***********************************************************************//** + * @class ctprob + * + * @brief Data selection tool + ***************************************************************************/ +class ctprob : public ctobservation { + +public: + // Constructors and destructors + ctprob(void); + explicit ctprob(const GObservations& obs); + ctprob(int argc, char *argv[]); + ctprob(const ctprob& app); + virtual ~ctprob(void); + + // Operators + ctprob& operator=(const ctprob& app); + + // Methods + void clear(void); + void run(void); + void save(void); + void publish(const std::string& name = ""); + +protected: + // Protected methods + void init_members(void); + void copy_members(const ctprob& app); + void free_members(void); + void get_parameters(void); + void get_obs(void); + void evaluate_probability(GCTAObservation* obs); + std::string set_outfile_name(const std::string& filename) const; + void save_fits(void); + void save_xml(void); + + // User parameters + std::string m_outobs; //!< Output event list or XML file + std::string m_prefix; //!< Prefix for multiple event lists + bool m_apply_edisp; //!< Apply energy dispersion? + bool m_publish; //!< Publish event list? + GChatter m_chatter; //!< Chattiness + + // Protected members + std::vector m_infiles; //!< Input event filenames + std::vector m_evtname; //!< Event extension names + std::vector m_gtiname; //!< GTI extension names +}; + +#endif /* CTPROB_HPP */ diff --git a/src/ctprob/ctprob.par b/src/ctprob/ctprob.par new file mode 100644 index 0000000..39ee4ec --- /dev/null +++ b/src/ctprob/ctprob.par @@ -0,0 +1,35 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: ctprob.par # +# # +# Author: Leonardo Di Venere # +# leonardo.divenere@ba.infn.it # +# INFN # +# # +# Purpose: Parameter file for the ctprob application # +# # +############################################################################# +# +# General parameters +#=================== +inobs, f, a, "events.fits",,, "Input event list or observation definition XML file" +outobs, f, a, "prob_events.fits",,, "Output event list or observation definition XML file" +prefix, s, h, "prob_",,, "Prefix for event lists in observation definition XML file" +inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Input model definition XML file" +caldb, s, a, "prod2",,, "Calibration database" +irf, s, a, "South_0.5h",,, "Instrument response function" +edisp, b, h, no,,, "Apply energy dispersion?" + +# +# Standard parameters +#==================== +publish, b, h, no,,, "Publish event list on VO Hub?" +chatter, i, h, 2,0,4, "Chattiness of output" +clobber, b, h, yes,,, "Overwrite existing output files with new output files?" +debug, b, h, no,,, "Debugging mode activated" +mode, s, h, "ql",,, "Mode of automatic parameters" +logfile, f, h, "ctprob.log",,, "Log filename" diff --git a/src/ctprob/main.cpp b/src/ctprob/main.cpp new file mode 100644 index 0000000..05257c9 --- /dev/null +++ b/src/ctprob/main.cpp @@ -0,0 +1,77 @@ +/*************************************************************************** + * ctprob - Computes event probability for a given model * + * ----------------------------------------------------------------------- * + * copyright (C) 2017 by Leonardo Di Venere * + * ----------------------------------------------------------------------- * + * * + * 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 * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file main.cpp + * @brief Computes event probability for a given model + * @author Leonardo Di Venere + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include "support.hpp" +#include "ctprob.hpp" + + +/***********************************************************************//** + * @brief Main entry point of application + * + * @param[in] argc Number of command line arguments. + * @param[in] argv Command line arguments. + * + * This is the main entry point of the ctprob application. It allocates a + * ctprob object and executes the application. Any exceptions that occur + * will be catched and corresponding error messages written in the + * application logger and into the standard output. + ***************************************************************************/ +int main (int argc, char *argv[]) +{ + // Initialise return code + int rc = 1; + + // Initialise pointer on application + ctprob* application = NULL; + + // Execute application + try { + + // Create instance of application + application = new ctprob(argc, argv); + + // Execute ctool + rc = execute_ctool(application); + + } + + catch (std::exception &e) { + + // Report exception + report_ctool_failure("ctprob", e.what()); + + } + + // Delete application + delete application; + + // Return + return rc; +} diff --git a/src/ctpsfcube/ctpsfcube.cpp b/src/ctpsfcube/ctpsfcube.cpp index b751f38..5d69b4b 100644 --- a/src/ctpsfcube/ctpsfcube.cpp +++ b/src/ctpsfcube/ctpsfcube.cpp @@ -49,7 +49,7 @@ /***********************************************************************//** * @brief Void constructor ***************************************************************************/ -ctpsfcube::ctpsfcube(void) : ctobservation(CTPSFCUBE_NAME, CTPSFCUBE_VERSION) +ctpsfcube::ctpsfcube(void) : ctobservation(CTPSFCUBE_NAME, VERSION) { // Initialise members init_members(); @@ -68,7 +68,7 @@ ctpsfcube::ctpsfcube(void) : ctobservation(CTPSFCUBE_NAME, CTPSFCUBE_VERSION) * observations container. ***************************************************************************/ ctpsfcube::ctpsfcube(const GObservations& obs) : - ctobservation(CTPSFCUBE_NAME, CTPSFCUBE_VERSION, obs) + ctobservation(CTPSFCUBE_NAME, VERSION, obs) { // Initialise members init_members(); @@ -86,7 +86,7 @@ ctpsfcube::ctpsfcube(const GObservations& obs) : * @param[in] argv Array of command line arguments. ***************************************************************************/ ctpsfcube::ctpsfcube(int argc, char *argv[]) : - ctobservation(CTPSFCUBE_NAME, CTPSFCUBE_VERSION, argc, argv) + ctobservation(CTPSFCUBE_NAME, VERSION, argc, argv) { // Initialise members init_members(); diff --git a/src/ctpsfcube/ctpsfcube.hpp b/src/ctpsfcube/ctpsfcube.hpp index 7169ab4..76df42b 100644 --- a/src/ctpsfcube/ctpsfcube.hpp +++ b/src/ctpsfcube/ctpsfcube.hpp @@ -33,8 +33,7 @@ #include "ctobservation.hpp" /* __Definitions _________________________________________________________ */ -#define CTPSFCUBE_NAME "ctpsfcube" -#define CTPSFCUBE_VERSION "1.1.0" +#define CTPSFCUBE_NAME "ctpsfcube" /***********************************************************************//** diff --git a/src/ctselect/ctselect.cpp b/src/ctselect/ctselect.cpp index 870d8f0..42e0b37 100644 --- a/src/ctselect/ctselect.cpp +++ b/src/ctselect/ctselect.cpp @@ -53,7 +53,7 @@ /***********************************************************************//** * @brief Void constructor ***************************************************************************/ -ctselect::ctselect(void) : ctobservation(CTSELECT_NAME, CTSELECT_VERSION) +ctselect::ctselect(void) : ctobservation(CTSELECT_NAME, VERSION) { // Initialise members init_members(); @@ -72,7 +72,7 @@ ctselect::ctselect(void) : ctobservation(CTSELECT_NAME, CTSELECT_VERSION) * provided in an observation container. ***************************************************************************/ ctselect::ctselect(const GObservations& obs) : - ctobservation(CTSELECT_NAME, CTSELECT_VERSION, obs) + ctobservation(CTSELECT_NAME, VERSION, obs) { // Initialise members init_members(); @@ -91,7 +91,7 @@ ctselect::ctselect(const GObservations& obs) : * @param[in] argv Array of command line arguments. ***************************************************************************/ ctselect::ctselect(int argc, char *argv[]) : - ctobservation(CTSELECT_NAME, CTSELECT_VERSION, argc, argv) + ctobservation(CTSELECT_NAME, VERSION, argc, argv) { // Initialise members init_members(); @@ -1175,39 +1175,6 @@ std::string ctselect::set_outfile_name(const std::string& filename) const /***********************************************************************//** - * @brief Get Good Time Intervals extension name - * - * @param[in] filename Input file name. - * @param[in] evtname Events extension name. - * - * Extracts the Good Time Intervals extension name from the event file. We - * do this by loading the events and accessing the Good Time Intervals - * extension name using the GCTAEventList::gtiname() method. If the file name - * is empty, the method returns `GTI`. - ***************************************************************************/ -std::string ctselect::get_gtiname(const std::string& filename, - const std::string& evtname) const -{ - // Initialise GTI name - std::string gtiname = gammalib::extname_gti; - - // Continue only if the filename is not empty - if (!filename.empty()) { - - // Load events - GCTAEventList events(filename+"["+evtname+"]"); - - // Get GTI name - gtiname = events.gtiname(); - - } - - // Return GTI name - return (gtiname); -} - - -/***********************************************************************//** * @brief Save event list in FITS format. * * Save the event list as a FITS file. The file name of the FITS file is @@ -1325,109 +1292,3 @@ void ctselect::save_xml(void) // Return return; } - - -/***********************************************************************//** - * @brief Save event list into FITS file - * - * @param[in] obs Pointer to CTA observation. - * @param[in] infile Input file name. - * @param[in] evtname Event extension name. - * @param[in] gtiname GTI extension name. - * @param[in] outfile Output file name. - * - * Saves an event list and the corresponding Good Time Intervals into a FITS - * file and copy all others extensions from the input file to the output - * file. - * - * If an extension name is specified in the @p outfile argument, the events - * and eventually also the Good Time Intervals will be extracted from the - * argument and used for writing the events. The format is - * - * [] - * - * where needs to be replaced by the name of the FITS file, - * and by the name of the events - * and Good Time Intervals extensions. For example - * - * myfits.fits[EVENTS1;GTI1] - * - * will write the selected events into the "EVENTS1" extension and the - * Good Time Intervals into the "GTI1" extension of the "myfits.fits" FITS - * file. If the Good Time Intervals extension name is skipped, e.g. - * - * myfits.fits[EVENTS1] - * - * the original extension name for the Good Time Intervals will be kept. - * Analogously, only the Good Time Intervals extension name can be changed - * by specifying - * - * myfits.fits[;GTI1] - * - * In none of the cases will the original events and Good Time Intervals be - * copied over to the output file. - ***************************************************************************/ -void ctselect::save_event_list(const GCTAObservation* obs, - const std::string& infile, - const std::string& evtname, - const std::string& gtiname, - const std::string& outfile) const -{ - // Save only if we have an event list - if (obs->eventtype() == "EventList") { - - // Set output FITS file event extension names - GFilename outname(outfile); - std::string outevt = evtname; - std::string outgti = gtiname; - if (outname.has_extname()) { - std::vector extnames = - gammalib::split(outname.extname(), ";"); - if (extnames.size() > 0) { - std::string extname = gammalib::strip_whitespace(extnames[0]); - if (!extname.empty()) { - outevt = extname; - } - } - if (extnames.size() > 1) { - std::string extname = gammalib::strip_whitespace(extnames[1]); - if (!extname.empty()) { - outgti = extname; - } - } - } - - // Create output FITS file - GFits outfits; - - // Write observation into FITS file - obs->write(outfits, outevt, outgti); - - // Copy all extensions other than evtname and gtiname extensions - // from the input to the output event list. The evtname and - // gtiname extensions are written by the save method, all others - // that may eventually be present have to be copied over - // explicitly. - GFits infits(infile); - for (int extno = 1; extno < infits.size(); ++extno) { - GFitsHDU* hdu = infits.at(extno); - if (hdu->extname() != evtname && - hdu->extname() != gtiname && - hdu->extname() != outevt && - hdu->extname() != outgti) { - outfits.append(*hdu); - } - } - - // Close input file - infits.close(); - - // Save file to disk and close it (we need both operations) - outfits.saveto(outname.url(), clobber()); - outfits.close(); - - } // endif: observation was unbinned - - // Return - return; -} diff --git a/src/ctselect/ctselect.hpp b/src/ctselect/ctselect.hpp index f41e209..250b3c9 100644 --- a/src/ctselect/ctselect.hpp +++ b/src/ctselect/ctselect.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctselect - Data selection tool * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2016 by Juergen Knoedlseder * + * copyright (C) 2010-2017 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -35,8 +35,7 @@ #include "ctobservation.hpp" /* __Definitions _________________________________________________________ */ -#define CTSELECT_NAME "ctselect" -#define CTSELECT_VERSION "1.2.0" +#define CTSELECT_NAME "ctselect" /***********************************************************************//** @@ -78,15 +77,8 @@ protected: std::string check_infile(const std::string& filename, const std::string& evtname) const; std::string set_outfile_name(const std::string& filename) const; - std::string get_gtiname(const std::string& filename, - const std::string& evtname) const; void save_fits(void); void save_xml(void); - void save_event_list(const GCTAObservation* obs, - const std::string& infile, - const std::string& evtname, - const std::string& gtiname, - const std::string& outfile) const; // User parameters std::string m_outobs; //!< Output event list or XML file diff --git a/src/ctskymap/ctskymap.cpp b/src/ctskymap/ctskymap.cpp index 5372d86..33ea26f 100644 --- a/src/ctskymap/ctskymap.cpp +++ b/src/ctskymap/ctskymap.cpp @@ -52,7 +52,7 @@ * * Constructs an empty sky mapping tool. ***************************************************************************/ -ctskymap::ctskymap(void) : ctobservation(CTSKYMAP_NAME, CTSKYMAP_VERSION) +ctskymap::ctskymap(void) : ctobservation(CTSKYMAP_NAME, VERSION) { // Initialise members init_members(); @@ -70,7 +70,7 @@ ctskymap::ctskymap(void) : ctobservation(CTSKYMAP_NAME, CTSKYMAP_VERSION) * Constructs sky mapping tool from an observation container. ***************************************************************************/ ctskymap::ctskymap(const GObservations& obs) : - ctobservation(CTSKYMAP_NAME, CTSKYMAP_VERSION, obs) + ctobservation(CTSKYMAP_NAME, VERSION, obs) { // Initialise members init_members(); @@ -91,7 +91,7 @@ ctskymap::ctskymap(const GObservations& obs) : * parameter setting. ***************************************************************************/ ctskymap::ctskymap(int argc, char *argv[]) : - ctobservation(CTSKYMAP_NAME, CTSKYMAP_VERSION, argc, argv) + ctobservation(CTSKYMAP_NAME, VERSION, argc, argv) { // Initialise members init_members(); diff --git a/src/ctskymap/ctskymap.hpp b/src/ctskymap/ctskymap.hpp index 67f0985..cf16a85 100644 --- a/src/ctskymap/ctskymap.hpp +++ b/src/ctskymap/ctskymap.hpp @@ -32,8 +32,7 @@ #include "ctobservation.hpp" /* __Definitions _________________________________________________________ */ -#define CTSKYMAP_NAME "ctskymap" -#define CTSKYMAP_VERSION "1.3.0" +#define CTSKYMAP_NAME "ctskymap" /***********************************************************************//** diff --git a/src/cttsmap/cttsmap.cpp b/src/cttsmap/cttsmap.cpp index 7372bab..e932dca 100644 --- a/src/cttsmap/cttsmap.cpp +++ b/src/cttsmap/cttsmap.cpp @@ -51,7 +51,7 @@ * * Constructs an empty cttsmap tool. ***************************************************************************/ -cttsmap::cttsmap(void) : ctlikelihood(CTTSMAP_NAME, CTTSMAP_VERSION) +cttsmap::cttsmap(void) : ctlikelihood(CTTSMAP_NAME, VERSION) { // Initialise members init_members(); @@ -69,7 +69,7 @@ cttsmap::cttsmap(void) : ctlikelihood(CTTSMAP_NAME, CTTSMAP_VERSION) * Constructs cttsmap tool from an observations container. ***************************************************************************/ cttsmap::cttsmap(const GObservations& obs) : - ctlikelihood(CTTSMAP_NAME, CTTSMAP_VERSION, obs) + ctlikelihood(CTTSMAP_NAME, VERSION, obs) { // Initialise members init_members(); @@ -89,7 +89,7 @@ cttsmap::cttsmap(const GObservations& obs) : * that are provided as command line arguments. ***************************************************************************/ cttsmap::cttsmap(int argc, char *argv[]) : - ctlikelihood(CTTSMAP_NAME, CTTSMAP_VERSION, argc, argv) + ctlikelihood(CTTSMAP_NAME, VERSION, argc, argv) { // Initialise members init_members(); diff --git a/src/cttsmap/cttsmap.hpp b/src/cttsmap/cttsmap.hpp index c4995b7..a916fd0 100644 --- a/src/cttsmap/cttsmap.hpp +++ b/src/cttsmap/cttsmap.hpp @@ -31,8 +31,7 @@ #include "ctlikelihood.hpp" /* __Definitions _________________________________________________________ */ -#define CTTSMAP_NAME "cttsmap" -#define CTTSMAP_VERSION "1.2.0" +#define CTTSMAP_NAME "cttsmap" /***********************************************************************//** diff --git a/src/ctulimit/ctulimit.cpp b/src/ctulimit/ctulimit.cpp index 93790b5..ce6f95b 100644 --- a/src/ctulimit/ctulimit.cpp +++ b/src/ctulimit/ctulimit.cpp @@ -53,7 +53,7 @@ * * Constructs an empty ctulimit tool. ***************************************************************************/ -ctulimit::ctulimit(void) : ctlikelihood(CTULIMIT_NAME, CTULIMIT_VERSION) +ctulimit::ctulimit(void) : ctlikelihood(CTULIMIT_NAME, VERSION) { // Initialise members init_members(); @@ -71,7 +71,7 @@ ctulimit::ctulimit(void) : ctlikelihood(CTULIMIT_NAME, CTULIMIT_VERSION) * Constructs ctulimit tool from an observations container. ***************************************************************************/ ctulimit::ctulimit(const GObservations& obs) : - ctlikelihood(CTULIMIT_NAME, CTULIMIT_VERSION, obs) + ctlikelihood(CTULIMIT_NAME, VERSION, obs) { // Initialise members init_members(); @@ -92,7 +92,7 @@ ctulimit::ctulimit(const GObservations& obs) : * parameters that are provided as command line arguments. ***************************************************************************/ ctulimit::ctulimit(int argc, char *argv[]) : - ctlikelihood(CTULIMIT_NAME, CTULIMIT_VERSION, argc, argv) + ctlikelihood(CTULIMIT_NAME, VERSION, argc, argv) { // Initialise members init_members(); diff --git a/src/ctulimit/ctulimit.hpp b/src/ctulimit/ctulimit.hpp index cb16e5f..f942418 100644 --- a/src/ctulimit/ctulimit.hpp +++ b/src/ctulimit/ctulimit.hpp @@ -31,8 +31,7 @@ #include "ctlikelihood.hpp" /* __Definitions _________________________________________________________ */ -#define CTULIMIT_NAME "ctulimit" -#define CTULIMIT_VERSION "1.2.0" +#define CTULIMIT_NAME "ctulimit" /***********************************************************************//** diff --git a/src/support/ctobservation.cpp b/src/support/ctobservation.cpp index ea89be4..4006a9b 100644 --- a/src/support/ctobservation.cpp +++ b/src/support/ctobservation.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctobservation - Base class for observation tools * * ----------------------------------------------------------------------- * - * copyright (C) 2016 by Juergen Knoedlseder * + * copyright (C) 2016-2017 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -326,3 +326,123 @@ const GCTAObservation* ctobservation::next_unbinned_observation(void) const // Return next CTA observation return obs; } + + +/***********************************************************************//** + * @brief Save event list in FITS format. + * + * Save the event list as a FITS file. The file name of the FITS file is + * specified by the `outobs` parameter. + ***************************************************************************/ +void ctobservation::save_events_fits(void) +{ + // Save only if there are observations + if (m_obs.size() > 0) { + + // Get output filename + std::string outobs = (*this)["outobs"].filename(); + + // Loop over all unbinned CTA observations in the container + for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL; + obs = next_unbinned_observation()) { + + // Get input file name and default extension names + std::string infile = obs->eventfile(); + std::string evtname = gammalib::extname_cta_events; + std::string gtiname = gammalib::extname_gti; + + // Create file name object + GFilename fname(outobs); + + // Extract filename and event extension name + std::string outfile = fname.url(); + + // Append event extension name. We handle here the possibility + // to write the events into a different extension. + if (fname.has_extname()) { + outfile += "["+fname.extname()+"]"; + } + else { + outfile += "["+gammalib::extname_cta_events+"]"; + } + + // Log filename + log_value(NORMAL, "Event list file", outfile); + + // Save event list + save_event_list(obs, infile, evtname, gtiname, outfile); + + // Exit the loop (if this method is called we should only + // have a single unbinned observation in the observation + // container + break; + + } // endfor: looped over all unbinned observations + + } // endif: there were observations + + // Return + return; +} + + +/***********************************************************************//** + * @brief Save event list(s) in XML format. + * + * Save the event list(s) into FITS files and write the file path information + * into a XML file. The filename of the XML file is specified by the `outobs` + * parameter, the filename(s) of the event lists are built by prepending a + * prefix to the input event list filenames. Any path present in the input + * filename will be stripped, i.e. the event list(s) will be written in the + * local working directory (unless a path is specified in the prefix). + ***************************************************************************/ +void ctobservation::save_events_xml(void) +{ + // Get output filename and prefix + std::string outobs = (*this)["outobs"].filename(); + + // Issue warning if output filename has no .xml suffix + log_string(TERSE, warn_xml_suffix(outobs)); + + // Loop over all unbinned CTA observations in the container + for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL; + obs = next_unbinned_observation()) { + + // Get input file name and default extension names + std::string infile = obs->eventfile(); + std::string evtname = gammalib::extname_cta_events; + std::string gtiname = gammalib::extname_gti; + + // Extract event and GTI extension names from input FITS file + GFilename fname(infile); + if (fname.has_extname()) { + evtname = fname.extname(); + } + gtiname = get_gtiname(fname.url(), evtname); + + // Set event output file name + std::string outfile = set_outfile_name(infile); + + // Append event extension name + outfile += "["+evtname+"]"; + + // Log filename + log_value(NORMAL, "Event list file", outfile); + + // Store output file name in observation + obs->eventfile(outfile); + + // Save event list + save_event_list(obs, infile, evtname, gtiname, outfile); + + } // endfor: looped over observations + + // Write observation definition XML file name into logger + log_value(NORMAL, "Obs. definition file", outobs); + + // Save observations in XML file + m_obs.save(outobs); + + // Return + return; +} diff --git a/src/support/ctobservation.hpp b/src/support/ctobservation.hpp index 7d5f9b2..3e6dc7a 100644 --- a/src/support/ctobservation.hpp +++ b/src/support/ctobservation.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctobservation - Base class for observation tools * * ----------------------------------------------------------------------- * - * copyright (C) 2016 by Juergen Knoedlseder * + * copyright (C) 2016-2017 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -78,6 +78,8 @@ protected: void init_members(void); void copy_members(const ctobservation& app); void free_members(void); + void save_events_fits(void); + void save_events_xml(void); // Protected members GObservations m_obs; //!< Observation container diff --git a/src/support/ctool.cpp b/src/support/ctool.cpp index d71b914..a8b5488 100644 --- a/src/support/ctool.cpp +++ b/src/support/ctool.cpp @@ -1028,6 +1028,41 @@ GCTARoi ctool::get_roi(void) /***********************************************************************//** + * @brief Set output file name. + * + * @param[in] filename Input file name. + * + * Required user parameters: + * prefix + * + * Converts an input file name into an output filename by prepending the + * prefix specified by the `prefix` parameter to the input file name. Any + * path as well as extension will be stripped from the input file name. Also + * a trailing `.gz` will be stripped as one cannot write into gzipped files. + ***************************************************************************/ +std::string ctool::set_outfile_name(const std::string& filename) +{ + // Get prefix User parameter + std::string prefix = (*this)["prefix"].string(); + + // Create filename + GFilename fname(filename); + + // Split input filename without any extensions into path elements + std::vector elements = gammalib::split(fname.url(), "/"); + + // The last path element is the filename + std::string outname = prefix + elements[elements.size()-1]; + + // Strip any ".gz" + outname = gammalib::strip_chars(outname, ".gz"); + + // Return output filename + return outname; +} + + +/***********************************************************************//** * @brief Set response for all CTA observations in container * * @param[in,out] obs Observation container @@ -1730,6 +1765,145 @@ std::vector ctool::cube_layer_usage(const GEbounds& cube_ebounds, /***********************************************************************//** + * @brief Save event list into FITS file + * + * @param[in] obs Pointer to CTA observation. + * @param[in] infile Input file name. + * @param[in] evtname Event extension name. + * @param[in] gtiname GTI extension name. + * @param[in] outfile Output file name. + * + * Saves an event list and the corresponding Good Time Intervals into a FITS + * file and copy all others extensions from the input file to the output + * file. + * + * If an extension name is specified in the @p outfile argument, the events + * and eventually also the Good Time Intervals will be extracted from the + * argument and used for writing the events. The format is + * + * [] + * + * where needs to be replaced by the name of the FITS file, + * and by the name of the events + * and Good Time Intervals extensions. For example + * + * myfits.fits[EVENTS1;GTI1] + * + * will write the selected events into the `EVENTS1` extension and the + * Good Time Intervals into the `GTI1` extension of the `myfits.fits` FITS + * file. If the Good Time Intervals extension name is skipped, e.g. + * + * myfits.fits[EVENTS1] + * + * the original extension name for the Good Time Intervals will be kept. + * Analogously, only the Good Time Intervals extension name can be changed + * by specifying + * + * myfits.fits[;GTI1] + * + * In none of the cases will the original events and Good Time Intervals be + * copied over to the output file. + ***************************************************************************/ +void ctool::save_event_list(const GCTAObservation* obs, + const std::string& infile, + const std::string& evtname, + const std::string& gtiname, + const std::string& outfile) const +{ + // Save only if we have an event list + if (obs->eventtype() == "EventList") { + + // Set output FITS file event extension names + GFilename outname(outfile); + std::string outevt = evtname; + std::string outgti = gtiname; + if (outname.has_extname()) { + std::vector extnames = + gammalib::split(outname.extname(), ";"); + if (extnames.size() > 0) { + std::string extname = gammalib::strip_whitespace(extnames[0]); + if (!extname.empty()) { + outevt = extname; + } + } + if (extnames.size() > 1) { + std::string extname = gammalib::strip_whitespace(extnames[1]); + if (!extname.empty()) { + outgti = extname; + } + } + } + + // Create output FITS file + GFits outfits; + + // Write observation into FITS file + obs->write(outfits, outevt, outgti); + + // Copy all extensions other than evtname and gtiname extensions + // from the input to the output event list. The evtname and + // gtiname extensions are written by the save method, all others + // that may eventually be present have to be copied over + // explicitly. + GFits infits(infile); + for (int extno = 1; extno < infits.size(); ++extno) { + GFitsHDU* hdu = infits.at(extno); + if (hdu->extname() != evtname && + hdu->extname() != gtiname && + hdu->extname() != outevt && + hdu->extname() != outgti) { + outfits.append(*hdu); + } + } + + // Close input file + infits.close(); + + // Save file to disk and close it (we need both operations) + outfits.saveto(outname.url(), clobber()); + outfits.close(); + + } // endif: observation was unbinned + + // Return + return; +} + + +/***********************************************************************//** + * @brief Get Good Time Intervals extension name + * + * @param[in] filename Input file name. + * @param[in] evtname Events extension name. + * + * Extracts the Good Time Intervals extension name from the event file. We + * do this by loading the events and accessing the Good Time Intervals + * extension name using the GCTAEventList::gtiname() method. If the file name + * is empty, the method returns `GTI`. + ***************************************************************************/ +std::string ctool::get_gtiname(const std::string& filename, + const std::string& evtname) const +{ + // Initialise GTI name + std::string gtiname = gammalib::extname_gti; + + // Continue only if the filename is not empty + if (!filename.empty()) { + + // Load events + GCTAEventList events(filename+"["+evtname+"]"); + + // Get GTI name + gtiname = events.gtiname(); + + } + + // Return GTI name + return (gtiname); +} + + +/***********************************************************************//** * @brief Set warning string if there are too few energies * * @param[in] energies Energies. diff --git a/src/support/ctool.hpp b/src/support/ctool.hpp index 5a6dde5..43876e8 100644 --- a/src/support/ctool.hpp +++ b/src/support/ctool.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctool - ctool base class * * ----------------------------------------------------------------------- * - * copyright (C) 2014-2016 by Juergen Knoedlseder * + * copyright (C) 2014-2017 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -104,6 +104,7 @@ protected: // Protected methods that extract user parameters GCTARoi get_roi(void); + std::string set_outfile_name(const std::string& filename); // Protected support methods void set_response(GObservations& obs); @@ -121,6 +122,13 @@ protected: std::vector cube_layer_usage(const GEbounds& cube_ebounds, const GEbounds& list_ebounds) const; bool is_valid_filename(const GFilename& filename) const; + std::string get_gtiname(const std::string& filename, + const std::string& evtname) const; + void save_event_list(const GCTAObservation* obs, + const std::string& infile, + const std::string& evtname, + const std::string& gtiname, + const std::string& outfile) const; // Protected warning strings std::string warn_too_few_energies(const GEnergies& energies) const; diff --git a/test/data/irf.root b/test/data/irf.root index 58a5941..e34aa8e 100644 --- a/test/data/irf.root +++ b/test/data/irf.root diff --git a/test/data/model_temporal_phasecurve.fits b/test/data/model_temporal_phasecurve.fits new file mode 100644 index 0000000..05491b5 --- /dev/null +++ b/test/data/model_temporal_phasecurve.fits diff --git a/test/data/model_temporal_phasecurve.xml b/test/data/model_temporal_phasecurve.xml new file mode 100644 index 0000000..2ad10b5 --- /dev/null +++ b/test/data/model_temporal_phasecurve.xml @@ -0,0 +1,22 @@ + + + + + + + + + + + + + + + + + + + + + + diff --git a/test/test_cssens.py b/test/test_cssens.py index 2de7474..bdd09db 100755 --- a/test/test_cssens.py +++ b/test/test_cssens.py @@ -2,7 +2,7 @@ # ========================================================================== # This scripts performs unit tests for the cssens script # -# Copyright (C) 2016 Juergen Knoedlseder +# Copyright (C) 2016-2017 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 @@ -168,7 +168,7 @@ class Test(test): # Check dimensions self.test_value(results.nrows(), 2, 'Check number of rows in sensitivity file') - self.test_value(results.ncols(), 7, + self.test_value(results.ncols(), 10, 'Check number of columns in sensitivity file') # Return diff --git a/test/test_ctbutterfly.py b/test/test_ctbutterfly.py index ba7fe16..0e7562a 100755 --- a/test/test_ctbutterfly.py +++ b/test/test_ctbutterfly.py @@ -168,6 +168,9 @@ class Test(test): cpy_butterfly['outfile'] = 'ctbutterfly_py3.dat' cpy_butterfly['logfile'] = 'ctbutterfly_py3.log' cpy_butterfly['chatter'] = 4 + + # Execute ctbutterfly tool + cpy_butterfly.logFileOpen() # Make sure we get a log file cpy_butterfly.execute() # Check result file diff --git a/test/test_ctphase.py b/test/test_ctphase.py index 2412ff4..3cf8a89 100644 --- a/test/test_ctphase.py +++ b/test/test_ctphase.py @@ -1,8 +1,8 @@ #! /usr/bin/env python # ========================================================================== -# This scripts performs unit tests for the ctselect tool. +# This scripts performs unit tests for the ctphase tool. # -# Copyright (C) 2014-2016 Juergen Knoedlseder +# Copyright (C) 2017 by Joshua Cardenzana # # 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 @@ -40,6 +40,8 @@ class Test(test): test.__init__(self) # Set test data + self._model_file = self._datadir + '/model_temporal_phasecurve.xml' + self._phase_file = self._datadir + '/model_temporal_phasecurve.fits' self._invalid_events = self._datadir + '/invalid_event_list.fits' # Return @@ -70,10 +72,9 @@ class Test(test): # Setup ctphase command cmd = ctphase+' inobs="'+self._events+'"'+ \ - ' outobs="ctphase_cmd1.fits"'+ \ - ' p0=0.1 mjd=51544.5'+ \ - ' f0=0.1 f1=0. f2=0.'+ \ - ' logfile="ctphase_cmd1.log" chatter=1' + ' outobs="ctphase_cmd1.fits"'+ \ + ' inmodel="'+self._model_file+'" source="Crab"'+ \ + ' logfile="ctphase_cmd1.log" chatter=1' # Check if execution of wrong command fails self.test_assert(self._execute('command_that_does_not_exist') != 0, @@ -88,10 +89,9 @@ class Test(test): # Setup ctphase command cmd = ctphase+' inobs="event_file_that_does_not_exist.fits"'+ \ - ' outobs="ctphase_cmd2.fits"'+ \ - ' p0=0.1 mjd=51544.5'+ \ - ' f0=0.1 f1=0. f2=0.'+ \ - ' logfile="ctphase_cmd2.log" chatter=1' + ' outobs="ctphase_cmd2.fits"'+ \ + ' inmodel="'+self._model_file+'" source="Crab"'+ \ + ' logfile="ctphase_cmd2.log" chatter=1' # Check if execution failed self.test_assert(self._execute(cmd) != 0, @@ -124,6 +124,8 @@ class Test(test): self._check_obs(phase.obs(), nobs=0) # Check that saving does nothing + phase['inmodel'] = self._model_file + phase['source'] = 'Crab' phase['outobs'] = 'ctphase_py0.fits' phase['logfile'] = 'ctphase_py0.log' phase.logFileOpen() @@ -136,11 +138,8 @@ class Test(test): # Now set ctphase parameters phase['inobs'] = self._events - phase['p0'] = 0. - phase['mjd'] = 51544.5 - phase['f0'] = 1. - phase['f1'] = 0.0 - phase['f2'] = 0. + phase['inmodel'] = self._model_file + phase['source'] = 'Crab' phase['outobs'] = 'ctphase_py1.fits' phase['logfile'] = 'ctphase_py1.log' phase['chatter'] = 2 @@ -151,34 +150,50 @@ class Test(test): phase.save() # Check result file - self._check_result_file('ctphase_py1.fits', nevents=6141) + self._check_result_file('ctphase_py1.fits') + + # Now run ctphase tool without a model file + phase['inmodel'] = 'NONE' + phase['source'] = 'NONE' + phase['phasecurve'] = self._phase_file + phase['mjd'] = 51544.5 + phase['phase'] = 0.0 + phase['f0'] = 1.0 + phase['f1'] = 0.1 + phase['f2'] = 0.01 + phase['outobs'] = 'ctphase_py2.fits' + phase['logfile'] = 'ctphase_py2.log' + phase['chatter'] = 3 + + # Execute ctphase tool + phase.logFileOpen() # Make sure we get a log file + phase.execute() + + # Check result file + self._check_result_file('ctphase_py2.fits') # Copy ctphase tool cpy_phase = phase.copy() - # Execute copy of ctphase tool again, now with a negative value - # of initial phase to test for invalid parameter value - self.test_try('Test invalid parameter value') - try: - cpy_phase['mjd'] = 51544.5 - cpy_phase['f0'] = 1. - cpy_phase['f1'] = 0.2 - cpy_phase['f2'] = 0. - cpy_phase['outobs'] = 'ctphase_py2.fits' - cpy_phase['logfile'] = 'ctphase_py2.log' - cpy_phase['chatter'] = 3 - cpy_phase['publish'] = True - cpy_phase.logFileOpen() # Needed to get a new log file - cpy_phase['p0'] = -0.5 - cpy_phase.execute() - self.test_try_failure('Exception not thrown') - except ValueError: - self.test_try_success() + # Clear phase tool + phase.clear() + + # Execute copy + cpy_phase['outobs'] = 'ctphase_py3.fits' + cpy_phase['logfile'] = 'ctphase_py3.log' + cpy_phase['chatter'] = 4 + cpy_phase['publish'] = True + cpy_phase.logFileOpen() # Needed to get a new log file + cpy_phase.execute() + # Check result file + self._check_result_file('ctphase_py3.fits') + + # Return return # Check result file - def _check_result_file(self, filename, nevents=591, dec=22.01, rad=3.0): + def _check_result_file(self, filename, nevents=6141): """ Check result file @@ -188,14 +203,10 @@ class Test(test): Event list file name nevents : int, optional Expected number of events - dec : float, optional - Expected Declination (deg) - rad : float, optional - Expected radius (deg) """ # Open result file - events = gammalib.GCTAEventList(filename) - phases = [ events[i].phase() for i in range(events.size())] + events = gammalib.GCTAEventList(filename) + phases = [events[i].phase() for i in range(events.size())] test_result = int(sum(phases)>0) self.test_value(test_result, 1, 'Check number of events') @@ -213,8 +224,6 @@ class Test(test): Models nobs : int, optional Expected number of observations - nevents : int, optional - Expected number of events """ # Check number of observations self.test_value(obs.size(), nobs, 'Check number of observations') diff --git a/test/test_ctprob.py b/test/test_ctprob.py new file mode 100755 index 0000000..0ad8306 --- /dev/null +++ b/test/test_ctprob.py @@ -0,0 +1,333 @@ +#! /usr/bin/env python +# ========================================================================== +# This scripts performs unit tests for the ctprob tool. +# +# Copyright (C) 2017 by Leonardo Di Venere +# +# 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 +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# ========================================================================== +import os +import gammalib +import ctools +from testing import test + + +# ============================ # +# Test class for ctprob tool # +# ============================ # +class Test(test): + """ + Test class for ctprob tool + """ + # Constructor + def __init__(self): + """ + Constructor + """ + # Call base class constructor + test.__init__(self) + + # Set test data + self._invalid_events = self._datadir + '/invalid_event_list.fits' + + # Return + return + + # Set test functions + def set(self): + """ + Set all test functions + """ + # Set test name + self.name('ctprob') + + # Append tests + self.append(self._test_cmd, 'Test ctprob on command line') + self.append(self._test_python, 'Test ctprob from Python') + + # Return + return + + # Test ctprob on command line + def _test_cmd(self): + """ + Test ctprob on the command line + """ + # Set tool name + ctprob = self._tool('ctprob') + + # Setup ctprob command + cmd = ctprob+' inobs="'+self._events+'"'+ \ + ' outobs="ctprob_cmd1.fits"'+ \ + ' inmodel='+self._model+ \ + ' caldb='+self._caldb+' irf='+self._irf+ \ + ' edisp=no'+ \ + ' logfile="ctprob_cmd1.log" chatter=1' + + # Check if execution was successful + self.test_assert(self._execute(cmd) == 0, + 'Check successful execution from command line') + + # Check result file + self._check_result_file('ctprob_cmd1.fits') + + # Setup ctprob command + cmd = ctprob+' inobs="event_file_that_does_not_exist.fits"'+ \ + ' outobs="ctprob_cmd2.fits"'+ \ + ' inmodel='+self._model+ \ + ' caldb='+self._caldb+' irf='+self._irf+ \ + ' edisp=no'+ \ + ' logfile="ctprob_cmd2.log" chatter=1' + + # Check if execution failed + self.test_assert(self._execute(cmd) != 0, + 'Check invalid input file when executed from command line') + + # Setup ctprob --help option + cmd = ctprob+' --help' + + # Check if execution was successful in case that the CTOOLS + # environment variable was set or failed otherwise + if 'CTOOLS' in os.environ: + self.test_value(self._execute(cmd), 0, + 'Check successful execution with --help option') + else: + self.test_assert(self._execute(cmd) != 0, + 'Check execution failure with --help option') + + # Return + return + + # Test ctprob from Python + def _test_python(self): + """ + Test ctprob from Python + """ + # Allocate empty ctprob tool + prob = ctools.ctprob() + + # Check that empty ctprob tool holds an empty observation + self._check_obs(prob.obs(), nobs=0) + + # Check that saving does nothing + prob['outobs'] = 'ctprob_py0.fits' + prob['logfile'] = 'ctprob_py0.log' + prob.logFileOpen() + prob.save() + self.test_assert(not os.path.isfile('ctprob_py0.fits'), + 'Check that no event list has been created') + + # Check that clearing does not lead to an exception or segfault + prob.clear() + + # Now set ctprob parameters + prob['inobs'] = self._events + prob['outobs'] = 'ctprob_py1.fits' + prob['inmodel'] = self._model + prob['caldb'] = self._caldb + prob['irf'] = self._irf + prob['edisp'] = False + prob['logfile'] = 'ctprob_py1.log' + prob['chatter'] = 2 + + # Run ctprob tool + prob.logFileOpen() # Make sure we get a log file + prob.run() + prob.save() + + # Check result file + self._check_result_file('ctprob_py1.fits') + + # Setup ctprob tool for an invalid event file + prob = ctools.ctprob() + prob['inobs'] = self._invalid_events + prob['outobs'] = 'ctprob_py2.fits' + prob['inmodel'] = self._model + prob['caldb'] = self._caldb + prob['irf'] = self._irf + prob['edisp'] = False + prob['logfile'] = 'ctprob_py2.log' + prob['chatter'] = 3 + + # Execute tool + prob.logFileOpen() # Needed to get a new log file + self.test_try('Test invalid event file') + try: + prob.execute() + self.test_try_failure('Exception not thrown') + except RuntimeError: + self.test_try_success() + + # Setup ctprob with a list of observations and test prefix + prob = ctools.ctprob() + prob['inobs'] = self._datadir + '/obs_unbinned.xml' + prob['prefix'] = 'ctprob_py3_' + prob['outobs'] = 'ctprob_py3_obs.xml' + prob['inmodel'] = self._model + prob['caldb'] = self._caldb + prob['irf'] = self._irf + prob['edisp'] = True + prob['logfile'] = 'ctprob_py3.log' + prob['chatter'] = 4 + + # Run ctprob tool + prob.logFileOpen() # Make sure we get a log file + prob.execute() + + # Check result file + self._check_result_file('ctprob_py3_crab_events.fits') + + # Return + return + + # Check result file + def _check_result_file(self, filename, nevents=6141): + """ + Check result file + + Parameters + ---------- + filename : str + Event list file name + nevents : int, optional + Expected number of events + """ + # Open result file + events = gammalib.GCTAEventList(filename) + + # Check event list + self._check_events(events, nevents=nevents) + + # Check that probability columns exist + colname1 = 'PROB_Background' + colname2 = 'PROB_Crab' + self._check_column(filename, colname1) + self._check_column(filename, colname2) + + # Check that the probability values are correctly normalized + self._check_normalization(filename, [colname1,colname2]) + + # Return + return + + # Check observation and event list + def _check_obs(self, obs, nobs=1, nevents=6141): + """ + Check observation and event list + + Parameters + ---------- + obs : `~gammalib.GObservations` + Models + nobs : int, optional + Expected number of observations + nevents : int, optional + Expected number of events + """ + # Check number of observations + self.test_value(obs.size(), nobs, 'Check number of observations') + + # If there is an observation then check the event list of the first + # one + if obs.size() > 0: + self._check_events(obs[0].events(), nevents=nevents) + + # Return + return + + # Check events + def _check_events(self, events, nevents=6141): + """ + Check event list + + Parameters + ---------- + events : `~gammalib.GCTAEventList` + Event list + nevents : int, optional + Expected number of events + """ + # Check event list + self.test_value(events.size(), nevents, 'Check number of events') + + # Return + return + + + def _check_column(self, filename, colname): + """ + Check that column colname exists. + + Parameters + ---------- + filename : str + Event list file name + colname : str + Name of the column to check + """ + # Check that column exists + gfits = gammalib.GFits(filename) + gtable = gfits[1] + check = int(gtable.contains(colname)) + self.test_value(check, 1, 'Check column '+colname) + + # Return + return + + + def _check_normalization(self, filename, colnames): + """ + Check that the probability values in columns identified by colnames are + correctly normalized. + + Parameters + ---------- + filename : str + Event list file name + colnames : list + Names of the columns to check + """ + # Check that the probability values are correctly normalized + + # Open fits table + gfits = gammalib.GFits(filename) + gtable = gfits[1] + ncols = gtable.ncols() + nevt = gtable.nrows() + srcs = [] + + # Loop over colnames + for colname in colnames: + + # Get the correct column corresponding to colname from the table + for icol in range(ncols): + gcol = gtable[icol] + if gcol.name()==colname: + break + + # Read and save column content + src1 = [] + for ievt in range(nevt): + src1.append(gcol[ievt]) + srcs.extend(src1) + + # Compute sum + tot = sum(srcs) + self.test_value(int(tot+0.5), nevt, 'Check that probability columns ' + 'are normalized') + + # Return + return + diff --git a/test/test_python_ctools.py b/test/test_python_ctools.py index 0ba05e0..18c98d5 100755 --- a/test/test_python_ctools.py +++ b/test/test_python_ctools.py @@ -40,6 +40,7 @@ import test_ctcubemask import test_ctbutterfly import test_ctulimit import test_cterror +import test_ctprob import test_pipelines @@ -118,6 +119,7 @@ def test(installed=False): test_ctbutterfly.Test(), test_ctulimit.Test(), test_cterror.Test(), + test_ctprob.Test(), test_pipelines.Test()] # Allocate test suite container