Michael Mayer / ctools

Commit 0645558ac45e2dac1ab141bf4001cb943f01b295

Authored by knodlseder
1 parent 78547611

initial version

Inline Side-by-side
Showing 1 changed file with 518 additions and 0 deletions
scripts/cstsdist.py 0 → 100755
  View file @0645558
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script generates the TS distribution for a particular model based
  4 +# on Monte-Carlo simulations.
  5 +#
  6 +# Copyright (C) 2011 Jurgen Knodlseder
  7 +#
  8 +# This program is free software: you can redistribute it and/or modify
  9 +# it under the terms of the GNU General Public License as published by
  10 +# the Free Software Foundation, either version 3 of the License, or
  11 +# (at your option) any later version.
  12 +#
  13 +# This program is distributed in the hope that it will be useful,
  14 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  15 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  16 +# GNU General Public License for more details.
  17 +#
  18 +# You should have received a copy of the GNU General Public License
  19 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  20 +#
  21 +# ==========================================================================
  22 +from ctatools import *
  23 +from gammalib import *
  24 +import obsutils
  25 +import sys
  26 +import csv
  27 +import math
  28 +
  29 +
  30 +# ============== #
  31 +# cstsdist class #
  32 +# ============== #
  33 +class cstsdist(GApplication):
  34 + """
  35 + This class implements the TS distribution generation script. It derives
  36 + from the GammaLib::GApplication class which provides support for parameter
  37 + files, command line arguments, and logging. In that way the Python
  38 + script behaves just as a regular ctool.
  39 + """
  40 + def __init__(self, *argv):
  41 + """
  42 + Constructor.
  43 + """
  44 + # Set name
  45 + self.name = "cstsdist"
  46 + self.version = "0.1.0"
  47 +
  48 + # Initialise some members
  49 + self.obs = None
  50 + self.bkg_model = None
  51 + self.full_model = None
  52 +
  53 + # Make sure that parfile exists
  54 + file = self.parfile()
  55 +
  56 + # Initialise application
  57 + if len(argv) == 0:
  58 + GApplication.__init__(self, self.name, self.version)
  59 + elif len(argv) ==1:
  60 + GApplication.__init__(self, self.name, self.version, *argv)
  61 + else:
  62 + raise TypeError("Invalid number of arguments given.")
  63 +
  64 + # Set logger properties
  65 + self.log_header()
  66 + self.log.date(True)
  67 +
  68 + # Return
  69 + return
  70 +
  71 + def __del__(self):
  72 + """
  73 + Destructor.
  74 + """
  75 + # Write separator into logger
  76 + if self.logTerse():
  77 + self.log("\n")
  78 +
  79 + # Return
  80 + return
  81 +
  82 + def parfile(self):
  83 + """
  84 + Check if parfile exists. If parfile does not exist then create a
  85 + default parfile. This kluge avoids shipping the cscript with a parfile.
  86 + """
  87 + # Set parfile name
  88 + parfile = self.name+".par"
  89 +
  90 + try:
  91 + pars = GPars(parfile)
  92 + except:
  93 + # Signal if parfile was not found
  94 + print "Parfile "+parfile+" not found. Create default parfile."
  95 +
  96 + # Create default parfile
  97 + pars = GPars()
  98 + pars.append(GPar("outfile","f","h","ts.dat","","","Output file name"))
  99 + pars.append(GPar("ntrials","i","a","10","","","Number of trials"))
  100 + pars.append(GPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database"))
  101 + pars.append(GPar("irf","s","a","kb_E_50h_v3","","","Instrument response function"))
  102 + pars.append(GPar("type","s","a","point","","","Source model type (point/gauss/shell/disk)"))
  103 + pars.append(GPar("index","r","h","-2.48","","","Spectral index"))
  104 + pars.append(GPar("offset","r","a","0.0","0.0","","Source offset angle (deg)"))
  105 + pars.append(GPar("bkg","s","a","$GAMMALIB/share/models/bkg_kb_E_50h_v3.txt","","","Background model file function (none=power law for E)"))
  106 + pars.append(GPar("emin","r","a","0.1","0.0","","Lower energy limit (TeV)"))
  107 + pars.append(GPar("emax","r","a","100.0","0.0","","Upper energy limit (TeV)"))
  108 + pars.append(GPar("enumbins","i","a","0","","","Number of energy bins (0=unbinned)"))
  109 + pars.append(GPar("duration","r","a","180000.0","","","Effective exposure time (s)"))
  110 + pars.append(GPar("rad","r","h","5.0","","","Radius of ROI (deg)"))
  111 + pars.append(GPar("npix","i","h","200","","","Number of pixels for binned"))
  112 + pars.append(GPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)"))
  113 + pars.append_standard()
  114 + pars.save(parfile)
  115 +
  116 + # Return
  117 + return
  118 +
  119 + def get_parameters(self):
  120 + """
  121 + Get parameters from parfile and setup the observation.
  122 + """
  123 + # Get parameters
  124 + self.m_outfile = self["outfile"].filename()
  125 + self.m_ntrials = self["ntrials"].integer()
  126 + self.m_caldb = self["caldb"].string()
  127 + self.m_irf = self["irf"].string()
  128 + self.m_type = self["type"].string()
  129 + self.m_index = self["index"].real()
  130 + self.m_offset = self["offset"].real()
  131 + self.m_bkg = self["bkg"].string()
  132 + self.m_emin = self["emin"].real()
  133 + self.m_emax = self["emax"].real()
  134 + self.m_enumbins = self["enumbins"].integer()
  135 + self.m_duration = self["duration"].real()
  136 + self.m_rad = self["rad"].real()
  137 + self.m_npix = self["npix"].integer()
  138 + self.m_binsz = self["binsz"].real()
  139 +
  140 + # Set some fixed parameters
  141 + self.m_log = False # Logging in client tools
  142 + self.m_debug = False # Debugging in client tools
  143 +
  144 + # Setup observation
  145 + self.obs = self.set_obs(emin=self.m_emin, emax=self.m_emax)
  146 +
  147 + # Initialise models. Note that we centre the point source at the Galactic
  148 + # center as our observation is also centred at the Galactic centre, so
  149 + # we're onaxis.
  150 + self.bkg_model = GModels()
  151 + self.full_model = GModels()
  152 + self.bkg_model.append(self.set_bkg_model())
  153 + self.full_model.append(self.set_bkg_model())
  154 + self.full_model.append(self.set_src_model(0.0, self.m_offset, \
  155 + flux=0.010, \
  156 + type=self.m_type, \
  157 + index=self.m_index))
  158 +
  159 + # Attach background model to observation container
  160 + self.obs.models(self.bkg_model)
  161 +
  162 + # Return
  163 + return
  164 +
  165 + def models(self, models):
  166 + """
  167 + Set model.
  168 + """
  169 + # Copy models
  170 + self.model = models.copy()
  171 +
  172 + # Return
  173 + return
  174 +
  175 + def execute(self):
  176 + """
  177 + Execute the script.
  178 + """
  179 + # Run the script
  180 + self.run()
  181 +
  182 + # Return
  183 + return
  184 +
  185 + def run(self):
  186 + """
  187 + Run the script.
  188 + """
  189 + # Switch screen logging on in debug mode
  190 + if self.logDebug():
  191 + self.log.cout(True)
  192 +
  193 + # Get parameters
  194 + self.get_parameters()
  195 +
  196 + # Write input parameters into logger
  197 + if self.logTerse():
  198 + self.log_parameters()
  199 + self.log("\n")
  200 +
  201 + # Write observation into logger
  202 + if self.logTerse():
  203 + self.log("\n")
  204 + self.log.header1("Observation")
  205 + self.log(str(self.obs))
  206 + self.log("\n")
  207 +
  208 + # Write models into logger
  209 + if self.logTerse():
  210 + self.log("\n")
  211 + self.log.header1("Test model")
  212 + self.log(str(self.full_model))
  213 + self.log("\n")
  214 +
  215 + # Write header
  216 + if self.logTerse():
  217 + self.log("\n")
  218 + self.log.header1("Generate TS distribution")
  219 +
  220 + # Loop over trials
  221 + for seed in range(self.m_ntrials):
  222 +
  223 + # Make a trial
  224 + result = self.trial(seed)
  225 +
  226 + # Write out result immediately
  227 + if seed == 0:
  228 + file = open(self.m_outfile, 'w')
  229 + writer = csv.DictWriter(file, result['colnames'])
  230 + writer.writerow(dict((_,_) for _ in result['colnames']))
  231 + else:
  232 + file = open(self.m_outfile, 'a')
  233 + writer = csv.DictWriter(file, result['colnames'])
  234 + writer.writerow(result['values'])
  235 + file.close()
  236 +
  237 + # Return
  238 + return
  239 +
  240 + def set_obs(self, lpnt=0.0, bpnt=0.0, emin=0.1, emax=100.0):
  241 + """
  242 + Returns an observation container with a single CTA observation.
  243 +
  244 + Keywords:
  245 + lpnt - Galactic longitude of pointing [deg] (default: 0.0)
  246 + bpnt - Galactic latitude of pointing [deg] (default: 0.0)
  247 + emin - Minimum energy [TeV] (default: 0.1)
  248 + emax - Maximum energy [TeV] (default: 100.0)
  249 + """
  250 + # Allocate observation container
  251 + obs = GObservations()
  252 +
  253 + # Set single pointing
  254 + pntdir = GSkyDir()
  255 + pntdir.lb_deg(lpnt, bpnt)
  256 +
  257 + # Create CTA observation
  258 + run = obsutils.set(pntdir, caldb=self.m_caldb, irf=self.m_irf, \
  259 + duration=self.m_duration, \
  260 + emin=emin, emax=emax, rad=self.m_rad)
  261 +
  262 + # Append observation to container
  263 + obs.append(run)
  264 +
  265 + # Return observation container
  266 + return obs
  267 +
  268 + def set_bkg_model(self):
  269 + """
  270 + Setup CTA background model.
  271 + """
  272 + # Define radial component
  273 + radial = GCTAModelRadialGauss(3.0)
  274 +
  275 + # Define spectral component
  276 + spectrum = GModelSpectralFunc(self.m_bkg)
  277 +
  278 + # Create background model
  279 + model = GCTAModelRadialAcceptance(radial, spectrum)
  280 + model.name("Background")
  281 + model.instruments("CTA")
  282 +
  283 + # Return background model
  284 + return model
  285 +
  286 + def set_src_model(self, l, b, flux=1.0, index=-2.48, \
  287 + type="point", sigma=1.0, radius=1.0, width=0.1, \
  288 + fitpos=False, fitidx=False):
  289 + """
  290 + Returns a single source with Crab-like spectrum. The source flux
  291 + can be scaled in Crab units. The Crab spectrum is based on MAGIC
  292 + observations (Albert et al. 2008, ApJ, 674, 1037).
  293 +
  294 + Parameters:
  295 + l - Galactic longitude of source location [deg]
  296 + b - Galactic latitude of source location [deg]
  297 + Keywords:
  298 + flux - Source flux [Crabs]
  299 + index - Spectral index
  300 + type - Source type ("point", "gauss", "disk", "shell")
  301 + sigma - Gaussian sigma (for type="gauss")
  302 + radius - Disk or shell inner radius [deg] (for type="disk" and type="shell")
  303 + width - Shell width [deg] (for type="shell")
  304 + fitpos - Fit position? (default: True)
  305 + fitidx - Fit index? (default: True)
  306 + """
  307 + # Set source location
  308 + location = GSkyDir()
  309 + location.lb_deg(l, b)
  310 +
  311 + # Set source spectrum
  312 + spectrum = GModelSpectralPlaw(flux, index)
  313 + spectrum["Prefactor"].scale(5.7e-16)
  314 + spectrum["PivotEnergy"].value(0.3)
  315 + spectrum["PivotEnergy"].scale(1.0e6)
  316 + if fitidx:
  317 + spectrum["Index"].free()
  318 + else:
  319 + spectrum["Index"].fix()
  320 +
  321 + # Set source
  322 + if type == "point":
  323 + spatial = GModelSpatialPtsrc(location)
  324 + if fitpos:
  325 + spatial[0].free()
  326 + spatial[1].free()
  327 + source = GModelPointSource(spatial, spectrum)
  328 + elif type == "gauss":
  329 + radial = GModelRadialGauss(location, sigma)
  330 + if fitpos:
  331 + radial[0].free()
  332 + radial[1].free()
  333 + source = GModelExtendedSource(radial, spectrum)
  334 + elif type == "disk":
  335 + radial = GModelRadialDisk(location, radius)
  336 + if fitpos:
  337 + radial[0].free()
  338 + radial[1].free()
  339 + source = GModelExtendedSource(radial, spectrum)
  340 + elif type == "shell":
  341 + radial = GModelRadialShell(location, radius, width)
  342 + if fitpos:
  343 + radial[0].free()
  344 + radial[1].free()
  345 + source = GModelExtendedSource(radial, spectrum)
  346 + else:
  347 + self.log("ERROR: Unknown source type '"+type+"'.\n")
  348 + return None
  349 +
  350 + # Set source name
  351 + source.name("Test")
  352 +
  353 + # Return source
  354 + return source
  355 +
  356 + def trial(self, seed):
  357 + """
  358 + Create the TS for a single trial.
  359 +
  360 + Parameters:
  361 + seed - Random number generator seed
  362 + """
  363 + # Write header
  364 + if self.logExplicit():
  365 + self.log.header2("Trial "+str(seed+1))
  366 +
  367 + # Simulate events
  368 + sim = obsutils.sim(self.obs, \
  369 + nbins=self.m_enumbins, \
  370 + seed=seed, \
  371 + binsz=self.m_binsz, \
  372 + npix=self.m_npix, \
  373 + log=self.m_log, debug=self.m_debug)
  374 +
  375 + # Determine number of events in simulation
  376 + nevents = 0.0
  377 + for run in sim:
  378 + nevents += run.events().number()
  379 +
  380 + # Write simulation results
  381 + if self.logExplicit():
  382 + self.log.header3("Simulation")
  383 + self.log.parformat("Number of simulated events")
  384 + self.log(nevents)
  385 + self.log("\n")
  386 +
  387 + # Fit background only
  388 + sim.models(self.bkg_model)
  389 + like = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
  390 + result_bgm = like.obs().models().copy()
  391 + LogL_bgm = like.opt().value()
  392 + npred_bgm = like.obs().npred()
  393 +
  394 + # Write background fit results
  395 + if self.logExplicit():
  396 + self.log.header3("Background model fit")
  397 + self.log.parformat("log likelihood")
  398 + self.log(LogL_bgm)
  399 + self.log("\n")
  400 + self.log.parformat("Number of predicted events")
  401 + self.log(npred_bgm)
  402 + self.log("\n")
  403 + for model in result_bgm:
  404 + self.log.parformat("Model")
  405 + self.log(model.name())
  406 + self.log("\n")
  407 + for par in model:
  408 + self.log(str(par)+"\n")
  409 +
  410 + # Fit background and test source
  411 + sim.models(self.full_model)
  412 + like = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
  413 + result_all = like.obs().models().copy()
  414 + LogL_all = like.opt().value()
  415 + npred_all = like.obs().npred()
  416 + ts = 2.0*(LogL_bgm-LogL_all)
  417 +
  418 + # Write background and test source fit results
  419 + if self.logExplicit():
  420 + self.log.header3("Background and test source model fit")
  421 + self.log.parformat("Test statistics")
  422 + self.log(ts)
  423 + self.log("\n")
  424 + self.log.parformat("log likelihood")
  425 + self.log(LogL_all)
  426 + self.log("\n")
  427 + self.log.parformat("Number of predicted events")
  428 + self.log(npred_all)
  429 + self.log("\n")
  430 + for model in result_all:
  431 + self.log.parformat("Model")
  432 + self.log(model.name())
  433 + self.log("\n")
  434 + for par in model:
  435 + self.log(str(par)+"\n")
  436 +
  437 + # Write result
  438 + elif self.logTerse():
  439 + self.log.parformat("Trial "+str(seed))
  440 + self.log("TS=")
  441 + self.log(ts)
  442 + self.log(" Prefactor=")
  443 + self.log(result_all["Test"]["Prefactor"].value())
  444 + self.log("+/-")
  445 + self.log(result_all["Test"]["Prefactor"].error())
  446 + self.log("\n")
  447 +
  448 + # Initialise results
  449 + colnames = []
  450 + values = {}
  451 +
  452 + # Set TS value
  453 + colnames.append("TS")
  454 + values["TS"] = ts
  455 +
  456 + # Set logL for background fit
  457 + colnames.append("LogL_bgm")
  458 + values["LogL_bgm"] = LogL_bgm
  459 +
  460 + # Set logL for full fit
  461 + colnames.append("LogL_all")
  462 + values["LogL_all"] = LogL_all
  463 +
  464 + # Set Nevents
  465 + colnames.append("Nevents")
  466 + values["Nevents"] = nevents
  467 +
  468 + # Set Npred for background fit
  469 + colnames.append("Npred_bkg")
  470 + values["Npred_bkg"] = npred_bgm
  471 +
  472 + # Set Npred for full fit
  473 + colnames.append("Npred_all")
  474 + values["Npred_all"] = npred_all
  475 +
  476 + # Gather free full fit parameters
  477 + for i in range(result_all.size()):
  478 + model = result_all[i]
  479 + model_name = model.name()
  480 + for k in range(model.size()):
  481 + par = model[k]
  482 + if par.isfree():
  483 +
  484 + # Set parameter name
  485 + name = model_name+"_"+par.name()
  486 +
  487 + # Append value
  488 + colnames.append(name)
  489 + values[name] = par.value()
  490 +
  491 + # Append error
  492 + name = "Unc_"+name
  493 + colnames.append(name)
  494 + values[name] = par.error()
  495 +
  496 + # Bundle together results
  497 + result = {'colnames': colnames, 'values': values}
  498 +
  499 + # Return
  500 + return result
  501 +
  502 +
  503 +# ======================== #
  504 +# Main routine entry point #
  505 +# ======================== #
  506 +if __name__ == '__main__':
  507 + """
  508 + Generates TS distribution.
  509 + """
  510 + # Create instance of application
  511 + app = cstsdist(sys.argv)
  512 +
  513 + # Open logfile
  514 + app.logFileOpen()
  515 +
  516 + # Execute application
  517 + app.execute()
  518 +
0 519 \ No newline at end of file
... ...