Dmitry Malyshev / ctools

Commit e6c90a83bf3dfaa7479e1588f446501fc7bfc985

Authored by Jürgen Knödlseder
2 parents 473d708f b3dc3387

Merge branch '1365-write-cslightcrv' into integration

Inline Side-by-side
Showing 3 changed files with 668 additions and 8 deletions
scripts/Makefile.am
  View file @e6c90a8
@@ -6,13 +6,13 @@ binscripts = cta_kb_bkgrates.py cta_root2irf.py cta_make_nodes.py \ @@ -6,13 +6,13 @@ binscripts = cta_kb_bkgrates.py cta_root2irf.py cta_make_nodes.py \
6 6
7 # Scripts to be distributed in pythondir 7 # Scripts to be distributed in pythondir
8 pythonscripts = cspull.py cssens.py cstsdist.py cscaldb.py csresmap.py cshessobs.py \ 8 pythonscripts = cspull.py cssens.py cstsdist.py cscaldb.py csresmap.py cshessobs.py \
9 - csobsdef.py csspec.py \ 9 + csobsdef.py csspec.py cslightcrv.py \
10 lsmodel.py \ 10 lsmodel.py \
11 obsutils.py 11 obsutils.py
12 12
13 # Scripts for which symbolic links will be created in bindir 13 # Scripts for which symbolic links will be created in bindir
14 linkscripts = cspull cssens cstsdist cscaldb csresmap \ 14 linkscripts = cspull cssens cstsdist cscaldb csresmap \
15 - csobsdef csspec cshessobs\ 15 + csobsdef csspec cshessobs cslightcrv \
16 lsmodel 16 lsmodel
17 17
18 # Distribute scripts 18 # Distribute scripts
scripts/cslightcrv.py 0 → 100755
  View file @e6c90a8
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# Light curve generation script.
  4 +#
  5 +# Copyright (C) 2014-2015 Michael Mayer
  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 gammalib
  22 +import ctools
  23 +import sys
  24 +
  25 +# ============ #
  26 +# csspec class #
  27 +# ============ #
  28 +class cslightcrv(ctools.cscript):
  29 + """
  30 + This class implements the creation of a light curve. It derives from
  31 + the ctools.cscript class which provides support for parameter files,
  32 + command line arguments, and logging. In that way the Python script
  33 + behaves just as a regular ctool.
  34 + """
  35 + def __init__(self, *argv):
  36 + """
  37 + Constructor.
  38 + """
  39 +
  40 + # Set name
  41 + self.name = "cslightcrv"
  42 + self.version = "1.0.0"
  43 +
  44 + # Initialise some members
  45 + self.obs = None
  46 +
  47 + # Initialise some members
  48 + if isinstance(argv[0],gammalib.GObservations):
  49 + self.obs = argv[0]
  50 + argv = argv[1:]
  51 + else:
  52 + self.obs = gammalib.GObservations()
  53 + self.obs.clear()
  54 + self.m_outfile = ""
  55 +
  56 + # Make sure that parfile exists
  57 + file = self.parfile()
  58 +
  59 + # Initialise application
  60 + if len(argv) == 0:
  61 + ctools.cscript.__init__(self, self.name, self.version)
  62 + elif len(argv) ==1:
  63 + ctools.cscript.__init__(self, self.name, self.version, *argv)
  64 + else:
  65 + raise TypeError("Invalid number of arguments given.")
  66 +
  67 + # Set logger properties
  68 + self.log_header()
  69 + self.log.date(True)
  70 +
  71 + # Return
  72 + return
  73 +
  74 + def __del__(self):
  75 + """
  76 + Destructor.
  77 + """
  78 + # Write separator into logger
  79 + if self.logTerse():
  80 + self.log("\n")
  81 +
  82 + # Return
  83 + return
  84 +
  85 + def parfile(self):
  86 + """
  87 + Check if parfile exists. If parfile does not exist then create a
  88 + default parfile. This kluge avoids shipping the cscript with a parfile.
  89 + """
  90 + # Set parfile name
  91 + parfile = self.name+".par"
  92 +
  93 + try:
  94 + pars = gammalib.GApplicationPars(parfile)
  95 + except:
  96 + # Signal if parfile was not found
  97 + sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")
  98 +
  99 + # Create default parfile
  100 + pars = gammalib.GApplicationPars()
  101 + pars.append(gammalib.GApplicationPar("inobs","f","a","events.fits","","","Event list, counts cube, or observation definition file"))
  102 + pars.append(gammalib.GApplicationPar("inmodel","f","a","$CTOOLS/share/models/crab.xml","","","Source model"))
  103 + pars.append(gammalib.GApplicationPar("srcname","s","a","Crab","","","Source name"))
  104 + pars.append(gammalib.GApplicationPar("outfile","f","a","lightcurve.fits","","","Output file name"))
  105 + pars.append(gammalib.GApplicationPar("expcube","f","a","NONE","","","Exposure cube file (only needed for stacked analysis)"))
  106 + pars.append(gammalib.GApplicationPar("psfcube","f","a","NONE","","","PSF cube file (only needed for stacked analysis)"))
  107 + pars.append(gammalib.GApplicationPar("bkgcube","s","a","NONE","","","Background cube file (only needed for stacked analysis)"))
  108 + pars.append(gammalib.GApplicationPar("caldb","s","a","dummy","","","Calibration database"))
  109 + pars.append(gammalib.GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
  110 + pars.append(gammalib.GApplicationPar("tmin","r","a","51544.5","","days","Lightcurve Time start [MJD]"))
  111 + pars.append(gammalib.GApplicationPar("tmax","r","a","51544.6","","days","Lightcurve Time stop [MJD]"))
  112 + pars.append(gammalib.GApplicationPar("tbins","i","a","5","","","Number of time bins points"))
  113 + pars.append(gammalib.GApplicationPar("tbinalg","s", "a","GTI","FILE|LIN|GTI","", "Algorithm for defining energy bins"))
  114 + pars.append(gammalib.GApplicationPar("tbinfile","s", "a","tbins.fits","","", "File defining the time binning"))
  115 + pars.append(gammalib.GApplicationPar("binned","b","a","no","yes|no","","Use binned analysis in each energy bin"))
  116 + pars.append(gammalib.GApplicationPar("nebins","i","h","20","","","Number of energy bins per light curve point"))
  117 + pars.append(gammalib.GApplicationPar("coordsys","s","a","CEL","CEL|GAL","","Coordinate System"))
  118 + pars.append(gammalib.GApplicationPar("proj","s","a","CAR","AIT|AZP|CAR|MER|MOL|STG|TAN","","Projection method"))
  119 + pars.append(gammalib.GApplicationPar("xref","r","a","83.63","0","360","First coordinate of image center in degrees (RA or galactic l)"))
  120 + pars.append(gammalib.GApplicationPar("yref","r","a","22.01","-90","90","Second coordinate of image center in degrees (DEC or galactic b)"))
  121 + pars.append(gammalib.GApplicationPar("nxpix","i","a","200","","","Size of the X axis in pixels"))
  122 + pars.append(gammalib.GApplicationPar("nypix","i","a","200","","","Size of the Y axis in pixels"))
  123 + pars.append(gammalib.GApplicationPar("binsz","r","a","0.02","","","Pixel size (deg/pixel)"))
  124 + pars.append(gammalib.GApplicationPar("emin","r","a","0.1","","","Lower energy limit for binned analysis (TeV)"))
  125 + pars.append(gammalib.GApplicationPar("emax","r","a","100.0","","","Upper energy limit for binned analysis (TeV)"))
  126 + pars.append(gammalib.GApplicationPar("calc_ts","b","h","yes","yes|no","","Compute TS value in each bin"))
  127 + pars.append(gammalib.GApplicationPar("calc_ulim","b","h","yes","yes|no","","Compute upper limit in each bin"))
  128 + pars.append(gammalib.GApplicationPar("fix_srcs","b","h","yes","yes|no","","Fix other skymodel parameters"))
  129 + pars.append(gammalib.GApplicationPar("fix_bkg","b","h","no","yes|no","","Fix background parameters"))
  130 + pars.append_standard()
  131 + pars.append(gammalib.GApplicationPar("logfile","f","h","cslightcrv.log","","","Log filename"))
  132 + pars.save(parfile)
  133 +
  134 + # Return
  135 + return
  136 +
  137 + def create_tbounds(self):
  138 +
  139 + # Create time bin container
  140 + self.m_tbins = gammalib.GGti()
  141 +
  142 + # check for binning algorithm
  143 + if self["tbinalg"].string() == "FILE":
  144 + try:
  145 + self.m_tbins.load(self["tbinfile"].filename())
  146 + except:
  147 + csv = gammalib.GCsv(self["tbinfile"].filename())
  148 + for i in range(csv.nrows()):
  149 + tmin = gammalib.GTime()
  150 + tmax = gammalib.GTime()
  151 + tmin.mjd(csv.real(i,0))
  152 + tmax.mjd(csv.real(i,1))
  153 + self.m_tbins.append(tmin,tmax)
  154 + elif self["tbinalg"].string() == "LIN":
  155 + # Use linear time binning
  156 + time_min = self["tmin"].real()
  157 + time_max = self["tmax"].real()
  158 + nbins = self["tbins"].integer()
  159 + step = (time_max - time_min) / float(nbins)
  160 + for i in range(nbins):
  161 + tmin = gammalib.GTime()
  162 + tmin.mjd(time_min + i*step)# ref)
  163 + tmax = gammalib.GTime()
  164 + tmax.mjd(time_min + (i+1)*step)#ref)
  165 + self.m_tbins.append(tmin,tmax)
  166 + elif self["tbinalg"].string() == "GTI":
  167 + # Use GTIs of observations
  168 + for obs in self.obs:
  169 + for i in range(obs.events().gti().size()):
  170 + self.m_tbins.append(obs.events().gti().tstart(i),obs.events().gti().tstop(i))
  171 + else:
  172 + raise AttributeError("tbinalg=\""+self["tbinalg"].string()+"\" unkown. Must be one of \"FILE\", \"LIN\" or \"GTI\"")
  173 +
  174 + # Return
  175 + return
  176 +
  177 + def get_parameters(self):
  178 + """
  179 + Get parameters from parfile and setup the observation.
  180 + """
  181 + # Set observation if not done before
  182 + if self.obs == None or self.obs.size() == 0:
  183 + self.require_inobs("cslightcrv::get_parameters()")
  184 + self.obs = self.get_observations()
  185 +
  186 + # Set models if we have none
  187 + if self.obs.models().size() == 0:
  188 + self.obs.models(self["inmodel"].filename())
  189 +
  190 + # Get source name
  191 + self.m_srcname = self["srcname"].string()
  192 +
  193 + # Get ebounds
  194 + self.create_tbounds()
  195 +
  196 + # Get binning flag
  197 + self.m_binned = self["binned"].boolean()
  198 + if self.m_binned:
  199 + self.m_xref = self["xref"].real()
  200 + self.m_yref = self["yref"].real()
  201 + self.m_nxpix = self["nxpix"].integer()
  202 + self.m_nypix = self["nypix"].integer()
  203 + self.m_binsz = self["binsz"].real()
  204 + self.m_coordsys = self["coordsys"].string()
  205 + self.m_proj = self["proj"].string()
  206 + self.m_ebins = self["nebins"].integer()
  207 + self.m_emin = self["emin"].real()
  208 + self.m_emax = self["emax"].real()
  209 +
  210 + # Read other parameters
  211 + self.m_outfile = self["outfile"].filename()
  212 +
  213 + # Get other parameeters
  214 + self.m_calc_ulimit = self["calc_ulim"].boolean()
  215 + self.m_calc_ts = self["calc_ts"].boolean()
  216 + self.m_fix_bkg = self["fix_bkg"].boolean()
  217 + self.m_fix_srcs = self["fix_srcs"].boolean()
  218 +
  219 + # Set some fixed parameters
  220 + self.m_log = False # Logging in client tools
  221 + self.m_chatter = self["chatter"].integer()
  222 + self.m_clobber = self["clobber"].boolean()
  223 + self.m_debug = self["debug"].boolean()
  224 +
  225 + # Return
  226 + return
  227 +
  228 + def models(self, models):
  229 + """
  230 + Set model.
  231 + """
  232 + # Copy models
  233 + self.obs.models(models.clone())
  234 +
  235 + # Return
  236 + return
  237 +
  238 + def execute(self):
  239 + """
  240 + Execute the script.
  241 + """
  242 + # Run the script
  243 + self.run()
  244 +
  245 + # Save residual map
  246 + self.fits.saveto(self.m_outfile, self.m_clobber)
  247 +
  248 + # Return
  249 + return
  250 +
  251 + def run(self):
  252 + """
  253 + Run the script.
  254 + """
  255 + # Switch screen logging on in debug mode
  256 + if self.logDebug():
  257 + self.log.cout(True)
  258 +
  259 + # Get parameters
  260 + self.get_parameters()
  261 +
  262 + # Write input parameters into logger
  263 + if self.logTerse():
  264 + self.log_parameters()
  265 + self.log("\n")
  266 +
  267 + # Write observation into logger
  268 + if self.logTerse():
  269 + self.log("\n")
  270 + self.log.header1("Observation")
  271 + self.log(str(self.obs))
  272 + self.log("\n")
  273 +
  274 + # Write header
  275 + if self.logTerse():
  276 + self.log("\n")
  277 + self.log.header1("Adjust model parameters")
  278 +
  279 + # Adjust model parameters dependent on input user parameters
  280 + for model in self.obs.models():
  281 +
  282 + # Set TS flag for all models to false.
  283 + # Source of interest will be set to true later
  284 + model.tscalc(False)
  285 +
  286 + # Log model name
  287 + if self.logExplicit():
  288 + self.log.header3(model.name())
  289 +
  290 + # Deal with the source of interest
  291 + if model.name() == self.m_srcname:
  292 + if self.m_calc_ts:
  293 + model.tscalc(True)
  294 +
  295 + elif self.m_fix_bkg and not model.classname() == "GModelSky":
  296 + for par in model:
  297 + if par.is_free() and self.logExplicit():
  298 + self.log(" Fixing \""+par.name()+"\"\n")
  299 + par.fix()
  300 +
  301 + elif self.m_fix_srcs and model.classname() == "GModelSky":
  302 + for par in model:
  303 + if par.is_free() and self.logExplicit():
  304 + self.log(" Fixing \""+par.name()+"\"\n")
  305 + par.fix()
  306 +
  307 + # Write header
  308 + if self.logTerse():
  309 + self.log("\n")
  310 + self.log.header1("Generate lightcurve")
  311 +
  312 + # Initialise FITS Table with extension "LIGHTCURVE"
  313 + table = gammalib.GFitsBinTable(self.m_tbins.size())
  314 + table.extname("LIGHTCURVE")
  315 +
  316 + # Add Header for compatibility with gammalib.GMWLSpectrum
  317 + table.card("INSTRUME", "CTA", "Name of Instrument")
  318 + table.card("TELESCOP", "CTA", "Name of Telescope")
  319 +
  320 + # Create FITS table columns
  321 + MJD = gammalib.GFitsTableDoubleCol("MJD", self.m_tbins.size())
  322 + MJD.unit("days")
  323 + e_MJD = gammalib.GFitsTableDoubleCol("e_MJD", self.m_tbins.size())
  324 + e_MJD.unit("days")
  325 +
  326 + # Create a FITS column for every free parameter
  327 + columns = []
  328 + for par in self.obs.models()[self.m_srcname]:
  329 + if par.is_free():
  330 + col = gammalib.GFitsTableDoubleCol(par.name(), self.m_tbins.size())
  331 + col.unit(par.unit())
  332 + columns.append(col)
  333 + e_col = gammalib.GFitsTableDoubleCol("e_"+par.name(), self.m_tbins.size())
  334 + e_col.unit(par.unit())
  335 + columns.append(e_col)
  336 +
  337 + # Create TS and upper limit columns
  338 + TSvalues = gammalib.GFitsTableDoubleCol("TS", self.m_tbins.size())
  339 + ulim_values = gammalib.GFitsTableDoubleCol("UpperLimit", self.m_tbins.size())
  340 + ulim_values.unit("ph/cm2/s")
  341 +
  342 + # Loop over energy bins
  343 + for i in range(self.m_tbins.size()):
  344 +
  345 + # Log information
  346 + if self.logTerse():
  347 + self.log("\n")
  348 + self.log.header2("Time bin "+str(i))
  349 +
  350 + # Get time boundaries
  351 + tmin = self.m_tbins.tstart(i)
  352 + tmax = self.m_tbins.tstop(i)
  353 +
  354 + # Compute time bin center and time width
  355 + tmean = (tmin + tmax)
  356 + tmean *= 0.5
  357 + twidth = (tmax - tmin)
  358 + twidth *= 0.5
  359 +
  360 + # Store time as MJD
  361 + MJD[i] = tmean.mjd()
  362 + e_MJD[i] = twidth.days()
  363 +
  364 + # Log information
  365 + if self.logExplicit():
  366 + self.log.header3("Selecting events")
  367 +
  368 + # Select events
  369 + select = ctools.ctselect(self.obs)
  370 + select["emin"].value("UNDEFINED")
  371 + select["emax"].value("UNDEFINED")
  372 + select["tmin"].real(tmin.convert(select.time_reference()))
  373 + select["tmax"].real(tmax.convert(select.time_reference()))
  374 + select["rad"].value("UNDEFINED")
  375 + select["ra"].value("UNDEFINED")
  376 + select["dec"].value("UNDEFINED")
  377 + select.run()
  378 +
  379 + # Retrieve observation
  380 + obs = select.obs()
  381 +
  382 + # Binned analysis
  383 + if self.m_binned:
  384 +
  385 + # Header
  386 + if self.logTerse():
  387 + self.log.header3("Binning events")
  388 +
  389 + # Bin events
  390 + bin = ctools.ctbin(select.obs())
  391 + bin["usepnt"].boolean(False)
  392 + bin["ebinalg"].string("LOG")
  393 + bin["xref"].real(self.m_xref)
  394 + bin["yref"].real(self.m_yref)
  395 + bin["binsz"].real(self.m_binsz)
  396 + bin["nxpix"].integer(self.m_nxpix)
  397 + bin["nypix"].integer(self.m_nypix)
  398 + bin["enumbins"].integer(self.m_ebins)
  399 + bin["emin"].real(self.m_emin)
  400 + bin["emax"].real(self.m_emax)
  401 + bin["coordsys"].string(self.m_coordsys)
  402 + bin["proj"].string(self.m_proj)
  403 + bin.run()
  404 +
  405 + # Header
  406 + if self.logTerse():
  407 + self.log.header3("Creating exposure cube")
  408 +
  409 + # Create exposure cube
  410 + expcube = ctools.ctexpcube(select.obs())
  411 + expcube["incube"].filename("NONE")
  412 + expcube["usepnt"].boolean(False)
  413 + expcube["ebinalg"].string("LOG")
  414 + expcube["xref"].real(self.m_xref)
  415 + expcube["yref"].real(self.m_yref)
  416 + expcube["binsz"].real(self.m_binsz)
  417 + expcube["nxpix"].integer(self.m_nxpix)
  418 + expcube["nypix"].integer(self.m_nypix)
  419 + expcube["enumbins"].integer(self.m_ebins)
  420 + expcube["emin"].real(self.m_emin)
  421 + expcube["emax"].real(self.m_emax)
  422 + expcube["coordsys"].string(self.m_coordsys)
  423 + expcube["proj"].string(self.m_proj)
  424 + expcube.run()
  425 +
  426 + # Header
  427 + if self.logTerse():
  428 + self.log.header3("Creating PSF cube")
  429 +
  430 + # Create psf cube
  431 + psfcube = ctools.ctpsfcube(select.obs())
  432 + psfcube["incube"].filename("NONE")
  433 + psfcube["usepnt"].boolean(False)
  434 + psfcube["ebinalg"].string("LOG")
  435 + psfcube["xref"].real(self.m_xref)
  436 + psfcube["yref"].real(self.m_yref)
  437 + psfcube["binsz"].real(self.m_binsz)
  438 + psfcube["nxpix"].integer(self.m_nxpix)
  439 + psfcube["nypix"].integer(self.m_nypix)
  440 + psfcube["enumbins"].integer(self.m_ebins)
  441 + psfcube["emin"].real(self.m_emin)
  442 + psfcube["emax"].real(self.m_emax)
  443 + psfcube["coordsys"].string(self.m_coordsys)
  444 + psfcube["proj"].string(self.m_proj)
  445 + psfcube.run()
  446 +
  447 + # Header
  448 + if self.logTerse():
  449 + self.log.header3("Creating background cube")
  450 +
  451 + # Create background cube
  452 + bkgcube = ctools.ctbkgcube(select.obs())
  453 + bkgcube["incube"].filename("NONE")
  454 + bkgcube["usepnt"].boolean(False)
  455 + bkgcube["ebinalg"].string("LOG")
  456 + bkgcube["xref"].real(self.m_xref)
  457 + bkgcube["yref"].real(self.m_yref)
  458 + bkgcube["binsz"].real(self.m_binsz)
  459 + bkgcube["nxpix"].integer(self.m_nxpix)
  460 + bkgcube["nypix"].integer(self.m_nypix)
  461 + bkgcube["enumbins"].integer(self.m_ebins)
  462 + bkgcube["emin"].real(self.m_emin)
  463 + bkgcube["emax"].real(self.m_emax)
  464 + bkgcube["coordsys"].string(self.m_coordsys)
  465 + bkgcube["proj"].string(self.m_proj)
  466 + bkgcube.run()
  467 +
  468 + # Set new binned observation
  469 + obs = bin.obs()
  470 +
  471 + # Set precomputed binned response
  472 + obs[0].response(expcube.expcube(), psfcube.psfcube(), bkgcube.cube())
  473 +
  474 + # Get new models
  475 + models = bkgcube.models()
  476 +
  477 + # Fix background models if required
  478 + if self.m_fix_bkg:
  479 + for model in models:
  480 + if not model.classname() == "GModelSky":
  481 + for par in model:
  482 + par.fix()
  483 +
  484 + # Set new models to binned observation
  485 + obs.models(models)
  486 +
  487 + # Header
  488 + if self.logTerse():
  489 + self.log.header3("Performing fit")
  490 +
  491 + # Likelihood
  492 + like = ctools.ctlike(obs)
  493 + like.run()
  494 +
  495 + # Skip bin if no event was present
  496 + if like.obs().logL() == 0.0:
  497 +
  498 + # Log information
  499 + if self.logTerse():
  500 + self.log("No event in this time bin. Bin is skipped\n")
  501 +
  502 + # Set all values to 0
  503 + for col in columns:
  504 + col[i] = 0.0
  505 + TSvalues[i] = 0.0
  506 + ulim_values[i] = 0.0
  507 + continue
  508 +
  509 + # Get results
  510 + fitted_models = like.obs().models()
  511 + source = fitted_models[self.m_srcname]
  512 +
  513 + # Calculate Upper Limit
  514 + ulimit_value = -1.0
  515 + if self.m_calc_ulimit:
  516 +
  517 + # Logging information
  518 + if self.logTerse():
  519 + self.log.header3("Computing upper limit")
  520 +
  521 + # Create upper limit object
  522 + ulimit = ctools.ctulimit(like.obs())
  523 + ulimit["srcname"].string(self.m_srcname)
  524 + ulimit["eref"].real(1.0)
  525 +
  526 + # Try to run upper limit and catch exceptions
  527 + try:
  528 + ulimit.run()
  529 + ulimit_value = ulimit.flux_ulimit()
  530 + except:
  531 + if self.logTerse():
  532 + self.log("Upper limit calculation failed\n")
  533 + ulimit_value = -1.0
  534 +
  535 + # Get TS value
  536 + TS = -1.0
  537 + if self.m_calc_ts:
  538 + TS = source.ts()
  539 +
  540 + # Set values for storage
  541 + TSvalues[i] = TS
  542 +
  543 + # Set FITS column values
  544 + for col in columns:
  545 + if "e_" == col.name()[:2]:
  546 + col[i] = source.spectral()[col.name()[2:]].error()
  547 + else:
  548 + col[i] = source.spectral()[col.name()].value()
  549 +
  550 + # Store upper limit value if available
  551 + if ulimit_value > 0.0:
  552 + ulim_values[i] = ulimit_value
  553 +
  554 + # Log information
  555 + if self.logExplicit():
  556 + self.log.header3("Results of bin "+str(i)+": MJD "+str(tmin.mjd())+"-"+str(tmax.mjd()))
  557 + for col in columns:
  558 + if "e_" == col.name()[:2]:
  559 + continue
  560 + value = source.spectral()[col.name()].value()
  561 + error = source.spectral()[col.name()].error()
  562 + unit = source.spectral()[col.name()].unit()
  563 + self.log(" > "+col.name()+": "+str(value)+" +- "+str(error)+" "+unit+"\n")
  564 + if self.m_calc_ts and TSvalues[i] > 0.0:
  565 + self.log(" > TS = "+str(TS)+" \n")
  566 + if self.m_calc_ulimit and ulim_values[i] > 0.0:
  567 + self.log(" > UL = "+str(ulim_values[i])+" [ph/cm2/s]")
  568 + self.log("\n")
  569 +
  570 + # Append filles columns to fits table
  571 + table.append(MJD)
  572 + table.append(e_MJD)
  573 + for col in columns:
  574 + table.append(col)
  575 + table.append(TSvalues)
  576 + table.append(ulim_values)
  577 +
  578 + # Create the FITS file now
  579 + self.fits = gammalib.GFits()
  580 + self.fits.append(table)
  581 +
  582 + # Return
  583 + return
  584 +
  585 +
  586 +# ======================== #
  587 +# Main routine entry point #
  588 +# ======================== #
  589 +if __name__ == '__main__':
  590 + """
  591 + Generates light curve.
  592 + """
  593 + # Create instance of application
  594 + app = cslightcrv(sys.argv)
  595 +
  596 + # Open logfile
  597 + app.logFileOpen()
  598 +
  599 + # Execute application
  600 + app.execute()
test/test_cscripts.sh.in
  View file @e6c90a8
@@ -190,9 +190,9 @@ fi @@ -190,9 +190,9 @@ fi
190 $ECHO " ok" 190 $ECHO " ok"
191 191
192 # 192 #
193 -# Test csspec  
194 -# ===========  
195 -$ECHO -n "Test csspec: " 193 +# Test csspec (unbinned)
  194 +# =====================
  195 +$ECHO -n "Test csspec (unbinned): "
196 chmod +x @abs_top_srcdir@/scripts/csspec.py 196 chmod +x @abs_top_srcdir@/scripts/csspec.py
197 @abs_top_srcdir@/scripts/csspec.py \ 197 @abs_top_srcdir@/scripts/csspec.py \
198 inobs="data/crab_events.fits" \ 198 inobs="data/crab_events.fits" \
@@ -216,9 +216,9 @@ fi @@ -216,9 +216,9 @@ fi
216 $ECHO " ok" 216 $ECHO " ok"
217 217
218 # 218 #
219 -# Test csspec Run2  
220 -# ================  
221 -$ECHO -n "Test csspec: " 219 +# Test csspec (binned)
  220 +# ====================
  221 +$ECHO -n "Test csspec (binned): "
222 chmod +x @abs_top_srcdir@/scripts/csspec.py 222 chmod +x @abs_top_srcdir@/scripts/csspec.py
223 @abs_top_srcdir@/scripts/csspec.py \ 223 @abs_top_srcdir@/scripts/csspec.py \
224 inobs="data/crab_events.fits" \ 224 inobs="data/crab_events.fits" \
@@ -249,3 +249,63 @@ else @@ -249,3 +249,63 @@ else
249 fi 249 fi
250 $ECHO " ok" 250 $ECHO " ok"
251 251
  252 +#
  253 +# Test cslightcrv (unbinned)
  254 +# ==========================
  255 +$ECHO -n "Test cslightcrv (unbinned): "
  256 +chmod +x @abs_top_srcdir@/scripts/cslightcrv.py
  257 +@abs_top_srcdir@/scripts/cslightcrv.py \
  258 + inobs="data/crab_events.fits" \
  259 + inmodel="data/crab.xml" \
  260 + caldb="irf" \
  261 + irf="cta_dummy_irf" \
  262 + tbins=3 \
  263 + tmin=51544.5 \
  264 + tmax=51544.53 \
  265 + srcname="Crab" \
  266 + tbinalg="LIN" \
  267 + binned=no \
  268 + outfile="lightcurve.fits" > cslightcrv.log
  269 +$ECHO -n "."
  270 +if [ -s "lightcurve.fits" ]
  271 +then
  272 + $ECHO -n "."
  273 +else
  274 + $ECHO " lightcurve.fits file is not valid"
  275 + exit 1
  276 +fi
  277 +$ECHO " ok"
  278 +
  279 +#
  280 +# Test cslightcrv (binned)
  281 +# ========================
  282 +$ECHO -n "Test cslightcrv (binned): "
  283 +chmod +x @abs_top_srcdir@/scripts/cslightcrv.py
  284 +@abs_top_srcdir@/scripts/cslightcrv.py \
  285 + inobs="data/crab_events.fits" \
  286 + inmodel="data/crab.xml" \
  287 + caldb="irf" \
  288 + irf="cta_dummy_irf" \
  289 + tbinalg="GTI" \
  290 + srcname="Crab" \
  291 + emin=0.1 \
  292 + emax=100.0 \
  293 + binned=yes \
  294 + nxpix=10 \
  295 + nypix=10 \
  296 + binsz=0.1 \
  297 + nebins=2 \
  298 + coordsys="CEL" \
  299 + xref=83.63 \
  300 + yref=22.01 \
  301 + proj="CAR" \
  302 + outfile="lightcurve.fits" > cslightcrv.log
  303 +$ECHO -n "."
  304 +if [ -s "lightcurve.fits" ]
  305 +then
  306 + $ECHO -n "."
  307 +else
  308 + $ECHO " lightcurve.fits file is not valid"
  309 + exit 1
  310 +fi
  311 +$ECHO " ok"