Olga Sergijenko / ctools

Commit 4d7363b7a9e2c62dd1b4da0647ddff6d3b2a948a

Authored by knodlseder
1 parent 8f12da81

initial version

Inline Side-by-side
Showing 1 changed file with 358 additions and 0 deletions
scripts/cspull.py 0 → 100755
  View file @4d7363b
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script generates the pull distribution for all model parameters.
  4 +#
  5 +# Copyright (C) 2011 Jurgen Knodlseder
  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 +from ctatools import *
  22 +from gammalib import *
  23 +import obsutils
  24 +import sys
  25 +import csv
  26 +import math
  27 +
  28 +
  29 +# ============ #
  30 +# cspull class #
  31 +# ============ #
  32 +class cspull(GApplication):
  33 + """
  34 + This class implements the pull distribution generation script. It derives
  35 + from the GammaLib::GApplication class which provides support for parameter
  36 + files, command line arguments, and logging. In that way the Python
  37 + script behaves just as a regular ctool.
  38 + """
  39 + def __init__(self, *argv):
  40 + """
  41 + Constructor.
  42 + """
  43 + # Set name
  44 + self.name = "cspull"
  45 + self.version = "0.1.0"
  46 +
  47 + # Initialise some members
  48 + self.obs = None
  49 + self.model = None
  50 + self.m_srcmdl = None
  51 +
  52 + # Make sure that parfile exists
  53 + file = self.parfile()
  54 +
  55 + # Initialise application
  56 + if len(argv) == 0:
  57 + GApplication.__init__(self, self.name, self.version)
  58 + elif len(argv) ==1:
  59 + GApplication.__init__(self, self.name, self.version, *argv)
  60 + else:
  61 + raise TypeError("Invalid number of arguments given.")
  62 +
  63 + # Set logger properties
  64 + self.log_header()
  65 + self.log.date(True)
  66 +
  67 + # Return
  68 + return
  69 +
  70 + def __del__(self):
  71 + """
  72 + Destructor.
  73 + """
  74 + # Write separator into logger
  75 + if self.logTerse():
  76 + self.log("\n")
  77 +
  78 + # Return
  79 + return
  80 +
  81 + def parfile(self):
  82 + """
  83 + Check if parfile exists. If parfile does not exist then create a
  84 + default parfile. This kluge avoids shipping the cscript with a parfile.
  85 + """
  86 + # Set parfile name
  87 + parfile = self.name+".par"
  88 +
  89 + try:
  90 + pars = GPars(parfile)
  91 + except:
  92 + # Signal if parfile was not found
  93 + print "Parfile "+parfile+" not found. Create default parfile."
  94 +
  95 + # Create default parfile
  96 + pars = GPars()
  97 + pars.append(GPar("srcmdl","f","a","$GAMMALIB/share/models/crab.xml","","","Source model"))
  98 + pars.append(GPar("outfile","f","a","pull.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("ra","r","a","83.6331","0","360","RA of pointing (deg)"))
  103 + pars.append(GPar("dec","r","a","22.0145","-90","90","Dec of pointing (deg)"))
  104 + pars.append(GPar("emin","r","a","0.1","0.0","","Lower energy limit (TeV)"))
  105 + pars.append(GPar("emax","r","a","100.0","0.0","","Upper energy limit (TeV)"))
  106 + pars.append(GPar("enumbins","i","a","0","","","Number of energy bins (0=unbinned)"))
  107 + pars.append(GPar("duration","r","a","180000.0","","","Effective exposure time (s)"))
  108 + pars.append(GPar("rad","r","h","5.0","","","Radius of ROI (deg)"))
  109 + pars.append(GPar("npix","i","h","200","","","Number of pixels for binned"))
  110 + pars.append(GPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)"))
  111 + pars.append_standard()
  112 + pars.save(parfile)
  113 +
  114 + # Return
  115 + return
  116 +
  117 + def get_parameters(self):
  118 + """
  119 + Get parameters from parfile and setup the observation.
  120 + """
  121 + # Get parameters
  122 + if self.model == None:
  123 + self.m_srcmdl = self["srcmdl"].filename()
  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_ra = self["ra"].real()
  129 + self.m_dec = self["dec"].real()
  130 + self.m_emin = self["emin"].real()
  131 + self.m_emax = self["emax"].real()
  132 + self.m_enumbins = self["enumbins"].integer()
  133 + self.m_duration = self["duration"].real()
  134 + self.m_rad = self["rad"].real()
  135 + self.m_npix = self["npix"].integer()
  136 + self.m_binsz = self["binsz"].real()
  137 +
  138 + # Set some fixed parameters
  139 + self.m_log = False # Logging in client tools
  140 + self.m_debug = False # Debugging in client tools
  141 +
  142 + # Setup observation
  143 + self.obs = self.set_obs(ra=self.m_ra, dec=self.m_dec, \
  144 + emin=self.m_emin, emax=self.m_emax)
  145 +
  146 + # Load source model
  147 + if self.m_srcmdl != None:
  148 + self.model = GModels(self.m_srcmdl)
  149 +
  150 + # Append source model to observation
  151 + self.obs.models(self.model)
  152 +
  153 + # Return
  154 + return
  155 +
  156 + def models(self, models):
  157 + """
  158 + Set model.
  159 + """
  160 + # Copy models
  161 + self.model = models.copy()
  162 +
  163 + # Return
  164 + return
  165 +
  166 + def execute(self):
  167 + """
  168 + Execute the script.
  169 + """
  170 + # Run the script
  171 + self.run()
  172 +
  173 + # Return
  174 + return
  175 +
  176 + def run(self):
  177 + """
  178 + Run the script.
  179 + """
  180 + # Switch screen logging on in debug mode
  181 + if self.logDebug():
  182 + self.log.cout(True)
  183 +
  184 + # Get parameters
  185 + self.get_parameters()
  186 +
  187 + # Write input parameters into logger
  188 + if self.logTerse():
  189 + self.log_parameters()
  190 + self.log("\n")
  191 +
  192 + # Write observation into logger
  193 + if self.logTerse():
  194 + self.log("\n")
  195 + self.log.header1("Observation")
  196 + self.log(str(self.obs))
  197 + self.log("\n")
  198 +
  199 + # Write header
  200 + if self.logTerse():
  201 + self.log("\n")
  202 + self.log.header1("Generate pull distribution")
  203 +
  204 + # Loop over trials
  205 + for seed in range(self.m_ntrials):
  206 +
  207 + # Make a trial
  208 + result = self.trial(seed)
  209 +
  210 + # Write out result immediately
  211 + if seed == 0:
  212 + file = open(self.m_outfile, 'w')
  213 + writer = csv.DictWriter(file, result['colnames'])
  214 + writer.writerow(dict((_,_) for _ in result['colnames']))
  215 + else:
  216 + file = open(self.m_outfile, 'a')
  217 + writer = csv.DictWriter(file, result['colnames'])
  218 + writer.writerow(result['pulls'])
  219 + file.close()
  220 +
  221 + # Return
  222 + return
  223 +
  224 + def set_obs(self, ra=83.6331, dec=22.0145, emin=0.1, emax=100.0):
  225 + """
  226 + Returns an observation container with a single CTA observation.
  227 +
  228 + Keywords:
  229 + ra - Right Ascension of pointing [deg] (default: 83.6331)
  230 + dec - Declination of pointing [deg] (default: 22.0145)
  231 + emin - Minimum energy [TeV] (default: 0.1)
  232 + emax - Maximum energy [TeV] (default: 100.0)
  233 + """
  234 + # Allocate observation container
  235 + obs = GObservations()
  236 +
  237 + # Set single pointing
  238 + pntdir = GSkyDir()
  239 + pntdir.radec_deg(ra, dec)
  240 +
  241 + # Create CTA observation
  242 + run = obsutils.set(pntdir, caldb=self.m_caldb, irf=self.m_irf, \
  243 + duration=self.m_duration, \
  244 + emin=emin, emax=emax, rad=self.m_rad)
  245 +
  246 + # Append observation to container
  247 + obs.append(run)
  248 +
  249 + # Return observation container
  250 + return obs
  251 +
  252 + def trial(self, seed):
  253 + """
  254 + Create the pull for a single trial.
  255 +
  256 + Parameters:
  257 + seed - Random number generator seed
  258 + """
  259 + # Write header
  260 + if self.logExplicit():
  261 + self.log.header2("Trial "+str(seed+1))
  262 +
  263 + # Simulate events
  264 + obs = obsutils.sim(self.obs, \
  265 + nbins=self.m_enumbins, \
  266 + seed=seed, \
  267 + binsz=self.m_binsz, \
  268 + npix=self.m_npix, \
  269 + log=self.m_log, debug=self.m_debug)
  270 +
  271 + # Determine number of events in simulation
  272 + nevents = 0.0
  273 + for run in obs:
  274 + nevents += run.events().number()
  275 +
  276 + # Write simulation results
  277 + if self.logExplicit():
  278 + self.log.header3("Simulation")
  279 + self.log.parformat("Number of simulated events")
  280 + self.log(nevents)
  281 + self.log("\n")
  282 +
  283 + # Fit model
  284 + like = obsutils.fit(obs, log=self.m_log, debug=self.m_debug)
  285 +
  286 + # Store results
  287 + logL = like.opt().value()
  288 + npred = like.obs().npred()
  289 + models = like.obs().models()
  290 +
  291 + # Write result header
  292 + if self.logExplicit():
  293 + self.log.header3("Pulls")
  294 +
  295 + # Gather results
  296 + colnames = []
  297 + values = {}
  298 + errors = {}
  299 + pulls = {}
  300 + for i in range(models.size()):
  301 + model = models[i]
  302 + model_name = model.name()
  303 + for k in range(model.size()):
  304 + par = model[k]
  305 + if par.isfree():
  306 +
  307 + # Set parameter name
  308 + name = model_name+"_"+par.name()
  309 + colnames.append(name)
  310 +
  311 + # Compute pull
  312 + fitted_value = par.value()
  313 + real_value = self.model[i][k].value()
  314 + error = par.error()
  315 + if error != 0.0:
  316 + pull = (fitted_value - real_value) / error
  317 + else:
  318 + pull = 99.0
  319 +
  320 + # Store pull
  321 + values[name] = fitted_value
  322 + errors[name] = error
  323 + pulls[name] = pull
  324 +
  325 + # Write result
  326 + if self.logExplicit():
  327 + self.log.parformat(name)
  328 + self.log(pull)
  329 + self.log(" (")
  330 + self.log(fitted_value)
  331 + self.log(" +/- ")
  332 + self.log(error)
  333 + self.log(")\n")
  334 +
  335 + # Bundle together results
  336 + result = {'colnames': colnames, 'values': values, \
  337 + 'errors': errors, 'pulls': pulls}
  338 +
  339 + # Return
  340 + return result
  341 +
  342 +
  343 +# ======================== #
  344 +# Main routine entry point #
  345 +# ======================== #
  346 +if __name__ == '__main__':
  347 + """
  348 + Generates pull distribution for a source model.
  349 + """
  350 + # Create instance of application
  351 + app = cspull(sys.argv)
  352 +
  353 + # Open logfile
  354 + app.logFileOpen()
  355 +
  356 + # Execute application
  357 + app.execute()
  358 +
0 \ No newline at end of file 359 \ No newline at end of file