Commit 1110a3531438088c055e47893f245839ec1458bd
Merge branch 'pull-request-1981-add-csphasecrv' into integration
Showing
4 changed files
with
1273 additions
and
0 deletions
cscripts/csphasecrv.par
0 → 100644
| 1 | +############################################################################# | |
| 2 | +# # | |
| 3 | +# CTA Science Analysis Tools # | |
| 4 | +# # | |
| 5 | +# ------------------------------------------------------------------------- # | |
| 6 | +# # | |
| 7 | +# File: csphasecrv.par # | |
| 8 | +# # | |
| 9 | +# Author: Rolf Buehler # | |
| 10 | +# rolf.buehler@desy.de # | |
| 11 | +# DESY # | |
| 12 | +# # | |
| 13 | +# Purpose: Parameter file for the csphasecrv script # | |
| 14 | +# # | |
| 15 | +############################################################################# | |
| 16 | +# | |
| 17 | +# Input and output parameters | |
| 18 | +#============================ | |
| 19 | +inobs, f, a, phased_events.fits,,, "Input event list or observation definition XML file" | |
| 20 | +inmodel, f, a, $CTOOLS/share/models/crab.xml,,, "Input model definition XML file" | |
| 21 | +srcname, s, a, Crab,,, "Source name" | |
| 22 | +caldb, s, a, prod2,,, "Calibration database" | |
| 23 | +irf, s, a, South_0.5h,,, "Instrument response function" | |
| 24 | +edisp, b, h, no,,, "Apply energy dispersion?" | |
| 25 | +outfile, f, a, phasecrv.fits,,, "Output fits file name." | |
| 26 | + | |
| 27 | +# | |
| 28 | +# Script parameters | |
| 29 | +#================== | |
| 30 | +tbinalg, s, a, LIN,LIN|FILE,, "Algorithm for defining time bins" | |
| 31 | +phbins, i, a, 5,1,10000, "Number of phase bins" | |
| 32 | +phbinfile, f, a, NONE,,, "File defining the phase binning" | |
| 33 | +emin, r, a, 0.1,,, "Lower energy limit of events (TeV)" | |
| 34 | +emax, r, a, 100.0,,, "Upper energy limit of events (TeV)" | |
| 35 | +enumbins, i, a, 0,0,100, "Number of energy bins for binned (0=unbinned)" | |
| 36 | +coordsys, s, a, CEL,CEL|GAL,, "Coordinate System" | |
| 37 | +proj, s, a, CAR,AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method" | |
| 38 | +xref, r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)" | |
| 39 | +yref, r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)" | |
| 40 | +nxpix, i, a, 200,,, "Size of the X axis in pixels" | |
| 41 | +nypix, i, a, 200,,, "Size of the Y axis in pixels" | |
| 42 | +binsz, r, a, 0.02,,, "Pixel size (deg/pixel)" | |
| 43 | + | |
| 44 | +# | |
| 45 | +# Standard parameters | |
| 46 | +#==================== | |
| 47 | +chatter, i, h, 2,0,4, "Chattiness of output" | |
| 48 | +clobber, b, h, yes,,, "Overwrite existing output files with new output files?" | |
| 49 | +debug, b, h, no,,, "Debugging mode activated" | |
| 50 | +mode, s, h, ql,,, "Mode of automatic parameters" | |
| 51 | +logfile, f, h, csphasecrv.log,,, "Log filename" | |
| 52 | + | |
| 53 | + | |
| 54 | + | ... | ... |
cscripts/csphasecrv.py
0 → 100755
| 1 | +#! /usr/bin/env python | |
| 2 | +# ========================================================================== | |
| 3 | +# Generates pulsar spectra as a function of phase | |
| 4 | +# | |
| 5 | +# Copyright (C) 2014-2017 Rolf Buehler | |
| 6 | +# | |
| 7 | +# This program is free software: you can redistribute it and/or modify | |
| 8 | +# it under the terms of the GNU General Public License as published by | |
| 9 | +# the Free Software Foundation, either version 3 of the License, or | |
| 10 | +# (at your option) any later version. | |
| 11 | +# | |
| 12 | +# This program is distributed in the hope that it will be useful, | |
| 13 | +# but WITHOUT ANY WARRANTY; without even the implied warranty of | |
| 14 | +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
| 15 | +# GNU General Public License for more details. | |
| 16 | +# | |
| 17 | +# You should have received a copy of the GNU General Public License | |
| 18 | +# along with this program. If not, see <http://www.gnu.org/licenses/>. | |
| 19 | +# | |
| 20 | +# ========================================================================== | |
| 21 | +import sys | |
| 22 | +import gammalib | |
| 23 | +import ctools | |
| 24 | +from cscripts import obsutils | |
| 25 | + | |
| 26 | + | |
| 27 | +# ================ # | |
| 28 | +# csphasecrv class # | |
| 29 | +# ================ # | |
| 30 | +class csphasecrv(ctools.cscript): | |
| 31 | + """ | |
| 32 | + Generates a spectra in phase bins | |
| 33 | + | |
| 34 | + This script computes spectra by performing a maximum likelihood fit | |
| 35 | + using :doc:`ctlike` in a series of phase bins for pulsars. | |
| 36 | + The phase bins can be either specified in an ASCII file, as an interval | |
| 37 | + divided into equally sized phase bins. | |
| 38 | + | |
| 39 | + Examples: | |
| 40 | + >>> phrv = csphasecrv() | |
| 41 | + >>> phrv.run() | |
| 42 | + >>> ... (querying for parameters) ... | |
| 43 | + >>> phrv = phrv.phasecrv() | |
| 44 | + Generates phase fits and retrieves dictionary with best fit models. | |
| 45 | + | |
| 46 | + >>> lcrv = csphasecrv(obs) | |
| 47 | + >>> lcrv.execute() | |
| 48 | + >>> ... (querying for parameters) ... | |
| 49 | + Generates phase fits from the observations | |
| 50 | + and saves results in XML files. | |
| 51 | + """ | |
| 52 | + | |
| 53 | + # Constructor | |
| 54 | + def __init__(self, *argv): | |
| 55 | + """ | |
| 56 | + Constructor | |
| 57 | + """ | |
| 58 | + # Set name | |
| 59 | + self._name = 'csphasecrv' | |
| 60 | + self._version = '1.0.0' | |
| 61 | + | |
| 62 | + # Initialise some members | |
| 63 | + self._srcname = '' | |
| 64 | + # Phases are stored in a nested list [[ph1min,ph1max], [ph2min,ph2max],..] | |
| 65 | + self._phbins = [[0.0,1.0]] | |
| 66 | + self._stacked = False | |
| 67 | + | |
| 68 | + # Initialise observation container from constructor arguments. | |
| 69 | + self._obs, argv = self._set_input_obs(argv) | |
| 70 | + | |
| 71 | + # Initialise script by calling the appropriate class constructor. | |
| 72 | + self._init_cscript(argv) | |
| 73 | + | |
| 74 | + # Return | |
| 75 | + return | |
| 76 | + | |
| 77 | + # Private methods | |
| 78 | + def _get_parameters(self): | |
| 79 | + """ | |
| 80 | + Get parameters from parfile | |
| 81 | + """ | |
| 82 | + # Setup observations (require response and allow event list, don't | |
| 83 | + # allow counts cube) | |
| 84 | + self._setup_observations(self._obs, True, True, False) | |
| 85 | + | |
| 86 | + # Set models if there are none in the container | |
| 87 | + if self._obs.models().size() == 0: | |
| 88 | + self._obs.models(self['inmodel'].filename()) | |
| 89 | + | |
| 90 | + # Get source name | |
| 91 | + self._srcname = self['srcname'].string() | |
| 92 | + | |
| 93 | + # Get time boundaries | |
| 94 | + self._phbins = self._create_tbounds() | |
| 95 | + | |
| 96 | + # Set stacked analysis flag to True if the requested number of | |
| 97 | + # energy bins is positive. Otherwise an unbinned analysis will | |
| 98 | + # be done and the stacked analysis flag will be set to False. | |
| 99 | + if self['enumbins'].integer() > 0: | |
| 100 | + self._stacked = True | |
| 101 | + else: | |
| 102 | + self._stacked = False | |
| 103 | + | |
| 104 | + # Make sure that remaining user parameters are queried now. We | |
| 105 | + # do not store the actual parameter values as we do not want | |
| 106 | + # too many instance attributes with enhances the maintenance | |
| 107 | + # costs. | |
| 108 | + self['emin'].real() | |
| 109 | + self['emax'].real() | |
| 110 | + if self._stacked: | |
| 111 | + self['coordsys'].string() | |
| 112 | + self['proj'].string() | |
| 113 | + self['xref'].real() | |
| 114 | + self['yref'].real() | |
| 115 | + self['nxpix'].integer() | |
| 116 | + self['nypix'].integer() | |
| 117 | + self['binsz'].real() | |
| 118 | + | |
| 119 | + # Do the same for the hidden parameters, just in case | |
| 120 | + self['edisp'].boolean() | |
| 121 | + | |
| 122 | + # Read ahead output parameters | |
| 123 | + if self._read_ahead(): | |
| 124 | + self['outfile'].filename() | |
| 125 | + | |
| 126 | + # Write input parameters into logger | |
| 127 | + self._log_parameters(gammalib.TERSE) | |
| 128 | + | |
| 129 | + # Return | |
| 130 | + return | |
| 131 | + | |
| 132 | + def _create_tbounds(self): | |
| 133 | + """ | |
| 134 | + Creates phase bins | |
| 135 | + | |
| 136 | + Returns | |
| 137 | + ------- | |
| 138 | + phbins : `[]` | |
| 139 | + List of phase bins | |
| 140 | + """ | |
| 141 | + # Initialise Good Time Intervals | |
| 142 | + phbins = [] | |
| 143 | + | |
| 144 | + # Get algorithm to use for defining the time intervals. Possible | |
| 145 | + # values are "FILE" or "LIN". This is enforced at | |
| 146 | + # parameter file level, hence no checking is needed. | |
| 147 | + algorithm = self['tbinalg'].string() | |
| 148 | + | |
| 149 | + # If the algorithm is "FILE" then handle a FITS or an ASCII file for | |
| 150 | + # the time bin definition | |
| 151 | + if algorithm == 'FILE': | |
| 152 | + | |
| 153 | + # Get the filename | |
| 154 | + filename = self['phbinfile'].filename() | |
| 155 | + | |
| 156 | + # Load ASCII file as CSV file and construct | |
| 157 | + # the GTIs from the rows of the CSV file. It is expected that the | |
| 158 | + # CSV file has two columns containing the "START" and "STOP" | |
| 159 | + # values of the phase bins. No header row is expected. | |
| 160 | + csv = gammalib.GCsv(filename) | |
| 161 | + for i in range(csv.nrows()): | |
| 162 | + phmin = csv.real(i,0) | |
| 163 | + phmax = csv.real(i,1) | |
| 164 | + phbins.append([phmin,phmax]) | |
| 165 | + | |
| 166 | + # If the algorithm is "LIN" then use a linear time binning, defined by | |
| 167 | + # the "phbins" user parameters | |
| 168 | + elif algorithm == 'LIN': | |
| 169 | + | |
| 170 | + # Get start and stop time and number of time bins | |
| 171 | + nbins = self['phbins'].integer() | |
| 172 | + | |
| 173 | + # Compute time step and setup the GTIs | |
| 174 | + ph_step = 1.0 / float(nbins) | |
| 175 | + for i in range(nbins): | |
| 176 | + phmin = i *ph_step | |
| 177 | + phmax = (i+1)*ph_step | |
| 178 | + phbins.append([phmin,phmax]) | |
| 179 | + | |
| 180 | + # Return Good Time Intervals | |
| 181 | + return phbins | |
| 182 | + | |
| 183 | + def _bin_observation(self, obs): | |
| 184 | + """ | |
| 185 | + Bin an observation if a binned analysis was requested | |
| 186 | + | |
| 187 | + Parameters | |
| 188 | + ---------- | |
| 189 | + obs : `~gammalib.GObservations` | |
| 190 | + Observation container | |
| 191 | + | |
| 192 | + Returns | |
| 193 | + ------- | |
| 194 | + obs : `~gammalib.GObservations` | |
| 195 | + Observation container where the first observation is a binned observation | |
| 196 | + """ | |
| 197 | + # Header | |
| 198 | + if self._logExplicit(): | |
| 199 | + self._log.header3('Binning events') | |
| 200 | + | |
| 201 | + # Bin events | |
| 202 | + cntcube = ctools.ctbin(obs) | |
| 203 | + cntcube['usepnt'] = False | |
| 204 | + cntcube['ebinalg'] = 'LOG' | |
| 205 | + cntcube['xref'] = self['xref'].real() | |
| 206 | + cntcube['yref'] = self['yref'].real() | |
| 207 | + cntcube['binsz'] = self['binsz'].real() | |
| 208 | + cntcube['nxpix'] = self['nxpix'].integer() | |
| 209 | + cntcube['nypix'] = self['nypix'].integer() | |
| 210 | + cntcube['enumbins'] = self['enumbins'].integer() | |
| 211 | + cntcube['emin'] = self['emin'].real() | |
| 212 | + cntcube['emax'] = self['emax'].real() | |
| 213 | + cntcube['coordsys'] = self['coordsys'].string() | |
| 214 | + cntcube['proj'] = self['proj'].string() | |
| 215 | + cntcube.run() | |
| 216 | + | |
| 217 | + # Header | |
| 218 | + if self._logExplicit(): | |
| 219 | + self._log.header3('Creating stacked response') | |
| 220 | + | |
| 221 | + # Get stacked response | |
| 222 | + response = obsutils.get_stacked_response(obs, | |
| 223 | + self['xref'].real(), | |
| 224 | + self['yref'].real(), | |
| 225 | + binsz=self['binsz'].real(), | |
| 226 | + nxpix=self['nxpix'].integer(), | |
| 227 | + nypix=self['nypix'].integer(), | |
| 228 | + emin=self['emin'].real(), | |
| 229 | + emax=self['emax'].real(), | |
| 230 | + enumbins=self['enumbins'].integer(), | |
| 231 | + edisp=self['edisp'].boolean(), | |
| 232 | + coordsys=self['coordsys'].string(), | |
| 233 | + proj=self['proj'].string()) | |
| 234 | + | |
| 235 | + # Retrieve a new oberservation container | |
| 236 | + new_obs = cntcube.obs().copy() | |
| 237 | + | |
| 238 | + # Set stacked response | |
| 239 | + if self['edisp'].boolean(): | |
| 240 | + new_obs[0].response(response['expcube'], response['psfcube'], | |
| 241 | + response['edispcube'], response['bkgcube']) | |
| 242 | + else: | |
| 243 | + new_obs[0].response(response['expcube'], response['psfcube'], | |
| 244 | + response['bkgcube']) | |
| 245 | + | |
| 246 | + # Get new models | |
| 247 | + models = response['models'] | |
| 248 | + | |
| 249 | + # Set models for new oberservation container | |
| 250 | + new_obs.models(models) | |
| 251 | + | |
| 252 | + # Return new oberservation container | |
| 253 | + return new_obs | |
| 254 | + | |
| 255 | + | |
| 256 | + # Public methods | |
| 257 | + def run(self): | |
| 258 | + """ | |
| 259 | + Run the script | |
| 260 | + """ | |
| 261 | + # Switch screen logging on in debug mode | |
| 262 | + if self._logDebug(): | |
| 263 | + self._log.cout(True) | |
| 264 | + | |
| 265 | + # Get parameters | |
| 266 | + self._get_parameters() | |
| 267 | + | |
| 268 | + # Write observation into logger | |
| 269 | + if self._logTerse(): | |
| 270 | + self._log('\n') | |
| 271 | + self._log.header1(gammalib.number('Input observation',len(self._obs))) | |
| 272 | + self._log(str(self._obs)) | |
| 273 | + self._log('\n') | |
| 274 | + | |
| 275 | + #Store and clear observations to split them in phase | |
| 276 | + orig_obs = self._obs.copy() | |
| 277 | + #self._obs.clear() | |
| 278 | + | |
| 279 | + #Dictionary to save phase fitted models | |
| 280 | + self._fitmodels = {} | |
| 281 | + | |
| 282 | + #Loop over all phases | |
| 283 | + for phbin in self._phbins: | |
| 284 | + | |
| 285 | + # Select events | |
| 286 | + select = ctools.ctselect(orig_obs) | |
| 287 | + select['emin'] = self['emin'].real() | |
| 288 | + select['emax'] = self['emax'].real() | |
| 289 | + select['tmin'] = 'UNDEFINED' | |
| 290 | + select['tmax'] = 'UNDEFINED' | |
| 291 | + select['rad'] = 'UNDEFINED' | |
| 292 | + select['ra'] = 'UNDEFINED' | |
| 293 | + select['dec'] = 'UNDEFINED' | |
| 294 | + select['expr'] = "PHASE>"+str(phbin[0])+" && PHASE<"+str(phbin[1]) | |
| 295 | + select.run() | |
| 296 | + | |
| 297 | + #Add phase to observation id | |
| 298 | + for i in range(0,select.obs().size()): | |
| 299 | + oldid = select.obs()[i].id() | |
| 300 | + select.obs()[i].id(oldid+"_"+str(phbin[0])+"-"+str(phbin[1])) | |
| 301 | + obs = select.obs() | |
| 302 | + | |
| 303 | + # If a stacked analysis is requested then bin the events | |
| 304 | + # and compute the stacked response functions and setup | |
| 305 | + # an observation container with a single stacked observation. | |
| 306 | + if self._stacked: | |
| 307 | + obs = self._bin_observation(select.obs()) | |
| 308 | + | |
| 309 | + # Header | |
| 310 | + if self._logExplicit(): | |
| 311 | + self._log.header3('Fitting the data') | |
| 312 | + | |
| 313 | + # Do maximum likelihood model fitting | |
| 314 | + like = ctools.ctlike(obs) | |
| 315 | + like['edisp'] = self['edisp'].boolean() | |
| 316 | + like.run() | |
| 317 | + | |
| 318 | + # Renormalize models to phase selection | |
| 319 | + # TODO move the scaling from the temporal to the spectral component | |
| 320 | + for model in like.obs().models(): | |
| 321 | + scaled_norm = model["Normalization"].value()/(phbin[1] - phbin[0]) | |
| 322 | + model["Normalization"].value(scaled_norm) | |
| 323 | + | |
| 324 | + #Store fit model | |
| 325 | + self._fitmodels[str(phbin[0])+"-"+str(phbin[1])] = like.obs().models().copy() | |
| 326 | + | |
| 327 | + ## Store observations | |
| 328 | + #for ob in select.obs(): | |
| 329 | + #self._obs.append(ob.copy()) | |
| 330 | + | |
| 331 | + # Return | |
| 332 | + return | |
| 333 | + | |
| 334 | + def execute(self): | |
| 335 | + """ | |
| 336 | + Execute the script | |
| 337 | + """ | |
| 338 | + # Open logfile | |
| 339 | + self.logFileOpen() | |
| 340 | + | |
| 341 | + # Read ahead output parameters | |
| 342 | + self._read_ahead(True) | |
| 343 | + | |
| 344 | + # Run the script | |
| 345 | + self.run() | |
| 346 | + | |
| 347 | + # Save fitted models | |
| 348 | + self.save() | |
| 349 | + | |
| 350 | + # Return | |
| 351 | + return | |
| 352 | + | |
| 353 | + def saveXML(self): | |
| 354 | + """ | |
| 355 | + Save phase fits into one XML model each | |
| 356 | + """ | |
| 357 | + # Write header | |
| 358 | + self._log_header1(gammalib.TERSE, 'Save phase XML') | |
| 359 | + | |
| 360 | + # Get out filename | |
| 361 | + outfile = str() | |
| 362 | + | |
| 363 | + # Save XMLs | |
| 364 | + for phbin in self._phbins: | |
| 365 | + phname = str(phbin[0])+"-"+str(phbin[1]) | |
| 366 | + outname = self['outfile'].filename() | |
| 367 | + outname_ph = (outname.file()).replace(".fits","_ph_"+phname+".xml") | |
| 368 | + outfile = outname.path()+outname_ph | |
| 369 | + print outfile | |
| 370 | + self._fitmodels[phname].save(outfile) | |
| 371 | + | |
| 372 | + # Return | |
| 373 | + return | |
| 374 | + | |
| 375 | + def save(self): | |
| 376 | + """ | |
| 377 | + Saves results to fits and XML | |
| 378 | + """ | |
| 379 | + print "SAVING" | |
| 380 | + self.savefits() | |
| 381 | + self.saveXML() | |
| 382 | + | |
| 383 | + # Return | |
| 384 | + return | |
| 385 | + | |
| 386 | + def phaserv(self): | |
| 387 | + """ | |
| 388 | + Return dictionary with best fit models | |
| 389 | + """ | |
| 390 | + # Return | |
| 391 | + return self._fitmodels | |
| 392 | + | |
| 393 | + def _get_free_par_names(self): | |
| 394 | + """ | |
| 395 | + Return list of free parameter names | |
| 396 | + | |
| 397 | + Returns | |
| 398 | + ------- | |
| 399 | + names : list of str | |
| 400 | + List of free parameter names. | |
| 401 | + """ | |
| 402 | + # Initialise empty list of free parameter names | |
| 403 | + names = [] | |
| 404 | + | |
| 405 | + # Collect list of free parameter names | |
| 406 | + for par in self._obs.models()[self._srcname]: | |
| 407 | + if par.is_free(): | |
| 408 | + names.append(par.name()) | |
| 409 | + | |
| 410 | + # Return names | |
| 411 | + return names | |
| 412 | + | |
| 413 | + def savefits(self): | |
| 414 | + """ | |
| 415 | + Saved phase fit results to a fits file | |
| 416 | + """ | |
| 417 | + # Initialise list of result dictionaries | |
| 418 | + results = [] | |
| 419 | + | |
| 420 | + # Get source parameters | |
| 421 | + pars = self._get_free_par_names() | |
| 422 | + | |
| 423 | + # Loop over time bins | |
| 424 | + for i in range(len(self._phbins)): | |
| 425 | + | |
| 426 | + # Get time boundaries | |
| 427 | + phmin = self._phbins[i][0] | |
| 428 | + phmax = self._phbins[i][1] | |
| 429 | + | |
| 430 | + # Write time bin into header | |
| 431 | + if self._logTerse(): | |
| 432 | + self._log.header2('PHASE '+str(phmin)+'-'+str(phmax)) | |
| 433 | + | |
| 434 | + # Initialise result dictionary | |
| 435 | + result = {'phmin': phmin, | |
| 436 | + 'phmax': phmax, | |
| 437 | + 'pars': pars, | |
| 438 | + 'values': {}} | |
| 439 | + # Extract parameter values | |
| 440 | + | |
| 441 | + #Store fit results | |
| 442 | + phname = str(phmin)+"-"+str(phmax) | |
| 443 | + source = self._fitmodels[phname][self._srcname] | |
| 444 | + for par in pars: | |
| 445 | + result['values'][par] = source[par].value() | |
| 446 | + result['values']['e_'+par] = source[par].error() | |
| 447 | + | |
| 448 | + # Append result to list of dictionaries | |
| 449 | + results.append(result) | |
| 450 | + | |
| 451 | + # Create FITS table from results | |
| 452 | + table = self._create_fits_table(results) | |
| 453 | + | |
| 454 | + # Create FITS file and append FITS table to FITS file | |
| 455 | + self._fits = gammalib.GFits() | |
| 456 | + self._fits.append(table) | |
| 457 | + | |
| 458 | + # Save to fits | |
| 459 | + outfile = self['outfile'].filename() | |
| 460 | + print outfile | |
| 461 | + self._fits.saveto(outfile, self._clobber()) | |
| 462 | + | |
| 463 | + #Return | |
| 464 | + return | |
| 465 | + | |
| 466 | + def _create_fits_table(self, results): | |
| 467 | + """ | |
| 468 | + Creates FITS binary table containing light curve results | |
| 469 | + | |
| 470 | + Parameters | |
| 471 | + ---------- | |
| 472 | + results : list of dict | |
| 473 | + List of result dictionaries | |
| 474 | + | |
| 475 | + Returns | |
| 476 | + ------- | |
| 477 | + table : `~gammalib.GFitsBinTable` | |
| 478 | + FITS binary table containing light curve | |
| 479 | + """ | |
| 480 | + # Determine number of rows in FITS table | |
| 481 | + nrows = len(results) | |
| 482 | + | |
| 483 | + # Create FITS Table with extension "LIGHTCURVE" | |
| 484 | + table = gammalib.GFitsBinTable(nrows) | |
| 485 | + table.extname('PHASESPEC') | |
| 486 | + | |
| 487 | + # Append time columns "MJD" and "e_MJD" | |
| 488 | + phmin = gammalib.GFitsTableDoubleCol('PHMIN', nrows) | |
| 489 | + phmax = gammalib.GFitsTableDoubleCol('PHMAX', nrows) | |
| 490 | + #mjd.unit('days') | |
| 491 | + #e_mjd.unit('days') | |
| 492 | + for i, result in enumerate(results): | |
| 493 | + phmin[i] = result['phmin'] | |
| 494 | + phmax[i] = result['phmax'] | |
| 495 | + table.append(phmin) | |
| 496 | + table.append(phmax) | |
| 497 | + | |
| 498 | + # Create parameter columns | |
| 499 | + for par in self._obs.models()[self._srcname]: | |
| 500 | + if par.is_free(): | |
| 501 | + name = par.name() | |
| 502 | + e_name = 'e_'+par.name() | |
| 503 | + col = gammalib.GFitsTableDoubleCol(name, nrows) | |
| 504 | + e_col = gammalib.GFitsTableDoubleCol(e_name, nrows) | |
| 505 | + col.unit(par.unit()) | |
| 506 | + e_col.unit(par.unit()) | |
| 507 | + for i, result in enumerate(results): | |
| 508 | + col[i] = result['values'][name] | |
| 509 | + e_col[i] = result['values'][e_name] | |
| 510 | + table.append(col) | |
| 511 | + table.append(e_col) | |
| 512 | + | |
| 513 | + # Return table | |
| 514 | + return table | |
| 515 | + | |
| 516 | + | |
| 517 | + | |
| 518 | +# ======================== # | |
| 519 | +# Main routine entry point # | |
| 520 | +# ======================== # | |
| 521 | +if __name__ == '__main__': | |
| 522 | + | |
| 523 | + # Create instance of application | |
| 524 | + app = csphasecrv(sys.argv) | |
| 525 | + | |
| 526 | + # Execute application | |
| 527 | + app.execute() | ... | ... |
doc/source/users/reference_manual/csphasecrv.rst
0 → 100644
| 1 | +.. _csphasecrv: | |
| 2 | + | |
| 3 | +csphasecrv | |
| 4 | +========== | |
| 5 | + | |
| 6 | +Computes phase dependent spectra for a given source. | |
| 7 | + | |
| 8 | + | |
| 9 | +Synopsis | |
| 10 | +-------- | |
| 11 | + | |
| 12 | +This script computes spectra by performing a maximum likelihood fit | |
| 13 | +using :doc:`ctlike` in a series of phase bins for pulsars. | |
| 14 | +The phase bins can be either specified in an ASCII file, as an interval | |
| 15 | +divided into equally sized time bins. | |
| 16 | +The format of the ASCII file is one row per time bin, each specifying the | |
| 17 | +start of stop value of the phase bin, separated by a whitespace. The phase | |
| 18 | +goes from 0.0 to 1.0 . | |
| 19 | + | |
| 20 | + | |
| 21 | +General parameters | |
| 22 | +------------------ | |
| 23 | + | |
| 24 | +``inobs [file]`` | |
| 25 | + Input event list or observation definition XML file. | |
| 26 | + | |
| 27 | +``inmodel [file]`` | |
| 28 | + Input model XML file. | |
| 29 | + | |
| 30 | +``srcname [string]`` | |
| 31 | + Name of the periodic source in the source model XML file | |
| 32 | + | |
| 33 | +``caldb [string]`` | |
| 34 | + Calibration database. | |
| 35 | + | |
| 36 | +``irf [string]`` | |
| 37 | + Instrumental response function. | |
| 38 | + | |
| 39 | +``(edisp = no) [boolean]`` | |
| 40 | + Applies energy dispersion to response computation. | |
| 41 | + | |
| 42 | +``outfile [file]`` | |
| 43 | + Name of the XML output file. The phase intervall will be automatically appended to the name. | |
| 44 | + | |
| 45 | +``(tbinalg = LIN) <FILE|LIN> [string]`` | |
| 46 | + Algorithm for defining time bins. | |
| 47 | + | |
| 48 | +``phbins [integer]`` | |
| 49 | + Number of phase bins. | |
| 50 | + | |
| 51 | +``phbinfile [file]`` | |
| 52 | + File defining the phase binning. | |
| 53 | + | |
| 54 | +``emin [real]`` | |
| 55 | + Lower energy limit of events (in TeV). | |
| 56 | + | |
| 57 | +``emax [real]`` | |
| 58 | + Upper energy limit of events (in TeV). | |
| 59 | + | |
| 60 | +``enumbins [integer]`` | |
| 61 | + Number of energy bins per phase bin (0=unbinned). | |
| 62 | + | |
| 63 | +``coordsys <CEL|GAL> [string]`` | |
| 64 | + Coordinate system (CEL - celestial, GAL - galactic). | |
| 65 | + | |
| 66 | +``proj <AIT|AZP|CAR|MER|MOL|STG|TAN> [string]`` | |
| 67 | + Projection method. | |
| 68 | + | |
| 69 | +``xref [real]`` | |
| 70 | + Right Ascension / Galactic longitude of image centre (J2000, in degrees). | |
| 71 | + | |
| 72 | +``yref [real]`` | |
| 73 | + Declination / Galactic latitude of image centre (J2000, in degrees). | |
| 74 | + | |
| 75 | +``nxpix [integer]`` | |
| 76 | + Size of the Right Ascension / Galactic longitude axis (in pixels). | |
| 77 | + | |
| 78 | +``nypix [integer]`` | |
| 79 | + Size of the Declination / Galactic latitude axis (in pixels). | |
| 80 | + | |
| 81 | +``binsz [real]`` | |
| 82 | + Pixel size (in degrees/pixel). | |
| 83 | + | |
| 84 | +Standard parameters | |
| 85 | +------------------- | |
| 86 | + | |
| 87 | +``(chatter = 2) [integer]`` | |
| 88 | + Verbosity of the executable: | |
| 89 | + chatter = 0: no information will be logged | |
| 90 | + | |
| 91 | + chatter = 1: only errors will be logged | |
| 92 | + | |
| 93 | + chatter = 2: errors and actions will be logged | |
| 94 | + | |
| 95 | + chatter = 3: report about the task execution | |
| 96 | + | |
| 97 | + chatter = 4: detailed report about the task execution | |
| 98 | + | |
| 99 | +``(clobber = yes) [boolean]`` | |
| 100 | + Specifies whether an existing light curve output file should be overwritten. | |
| 101 | + | |
| 102 | +``(debug = no) [boolean]`` | |
| 103 | + Enables debug mode. In debug mode the executable will dump any log file output to the console. | |
| 104 | + | |
| 105 | +``(mode = ql) [string]`` | |
| 106 | + Mode of automatic parameters (default is "ql", i.e. "query and learn"). | |
| 107 | + | |
| 108 | +``(logfile = csphasecrv.log) [filename]`` | |
| 109 | + Log filename. | |
| 110 | + | |
| 111 | + | |
| 112 | +Related tools or scripts | |
| 113 | +------------------------ | |
| 114 | + | |
| 115 | +:doc:`ctlike` | ... | ... |
test/data/phased_events.fits
0 → 100644
No preview for this file type