diff --git a/scripts/cssens.py b/scripts/cssens.py index 725bf50..242bee8 100755 --- a/scripts/cssens.py +++ b/scripts/cssens.py @@ -17,7 +17,7 @@ # # The source is modelled as a point source # -# Copyright (C) 2011 Jurgen Knodlseder +# Copyright (C) 2011-2014 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 @@ -45,682 +45,682 @@ import math # cssens class # # ============ # class cssens(GApplication): - """ - This class implements the sensitivity computation script. It derives - from the GammaLib::GApplication class which provides support for parameter - files, command line arguments, and logging. In that way the Python - script behaves just as a regular ctool. - """ - def __init__(self, *argv): - """ - Constructor. - """ - # Set name - self.name = "cssens" - self.version = "0.3.0" - - # Make sure that parfile exists - file = self.parfile() - - # Initialise application - if len(argv) == 0: - GApplication.__init__(self, self.name, self.version) - elif len(argv) ==1: - GApplication.__init__(self, self.name, self.version, *argv) - else: - raise TypeError("Invalid number of arguments given.") - - # Set logger properties - self.logFileOpen() - self.log_header() - self.log.date(True) - - # Return - return - - def __del__(self): - """ - Destructor. - """ - # Write separator into logger - if self.logTerse(): - self.log("\n") - - # Return - return - - def parfile(self): - """ - Check if parfile exists. If parfile does not exist then create a - default parfile. This kluge avoids shipping the cscript with a parfile. - """ - # Set parfile name - parfile = self.name+".par" - - try: - pars = GApplicationPars(parfile) - except: - # Signal if parfile was not found - sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n") - - # Create default parfile - pars = GApplicationPars() - pars.append(GApplicationPar("outfile","f","h","sensitivity.dat","","","Output file name")) - pars.append(GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database")) - pars.append(GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function")) - pars.append(GApplicationPar("type","s","a","point","","","Source model type (point/gauss/shell/disk)")) - pars.append(GApplicationPar("offset","r","a","0.0","0.0","","Source offset angle (deg)")) - pars.append(GApplicationPar("bkg","s","a","$GAMMALIB/share/models/bkg_dummy.txt","","","Background model file function (none=power law for E)")) - pars.append(GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)")) - pars.append(GApplicationPar("rad","r","a","5.0","","","Radius of ROI (deg)")) - pars.append(GApplicationPar("enumbins","i","h","0","","","Number of energy bins (0=unbinned)")) - pars.append(GApplicationPar("npix","i","h","200","","","Number of pixels for binned")) - pars.append(GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)")) - pars.append(GApplicationPar("sigma","r","h","5.0","","","Significance threshold")) - pars.append(GApplicationPar("ts_use","b","h","yes","","","Use TS?")) - pars.append(GApplicationPar("index","r","h","-2.48","","","Assumed spectral index")) - pars.append(GApplicationPar("radius","r","h","0.1","","","Extended source model radius")) - pars.append(GApplicationPar("width","r","h","0.05","","","Extended source model width")) - pars.append(GApplicationPar("max_iter","i","h","50","","","Maximum number of iterations")) - pars.append(GApplicationPar("num_avg","i","h","3","","","Number of iterations for sliding average")) - pars.append_standard() - pars.save(parfile) - - # Return - return - - def get_parameters(self): - """ - Get parameters from parfile. - """ - # Get parameters - self.m_outfile = self["outfile"].filename() - self.m_caldb = self["caldb"].string() - self.m_irf = self["irf"].string() - self.m_type = self["type"].string() - self.m_offset = self["offset"].real() - self.m_bkg = self["bkg"].string() - self.m_duration = self["duration"].real() - self.m_roi = self["rad"].real() - self.m_nbins = self["enumbins"].integer() - self.m_npix = self["npix"].integer() - self.m_binsz = self["binsz"].real() - self.m_ts_thres = self["sigma"].real()*self["sigma"].real() - self.m_use_ts = self["ts_use"].boolean() - self.m_index = self["index"].real() - self.m_radius = self["radius"].real() - self.m_width = self["width"].real() - self.m_max_iter = self["max_iter"].integer() - self.m_num_avg = self["num_avg"].integer() - - # Set some fixed parameters - self.m_log = False # Logging in client tools - self.m_debug = False # Debugging in client tools - - # Return - return - - def execute(self): - """ - Execute the script. - """ - # Run the script - self.run() - - # Return - return - - def run(self): - """ - Run the script. - """ - # Switch screen logging on in debug mode - if self.logDebug(): - self.log.cout(True) - - # Get parameters - self.get_parameters() - - # Write input parameters into logger - if self.logTerse(): - self.log_parameters() - self.log("\n") - - # Initialise script - colnames = ['loge', 'emin', 'emax', 'crab', 'pflux', 'eflux', 'diffSens'] - results = [] - - # Initialise models. Note that we centre the point source at the Galactic - # center as our observation is also centred at the Galactic centre, so - # we're onaxis. - bkg_model = GModels() - full_model = GModels() - bkg_model.append(self.set_bkg_model()) - full_model.append(self.set_bkg_model()) - full_model.append(self.set_src_model(0.0, \ - self.m_offset, \ - type=self.m_type, \ - index=self.m_index)) - - # Write models into logger - if self.logTerse(): - self.log("\n") - self.log.header1("Models") - self.log.header2("Background model") - self.log(str(bkg_model)) - self.log("\n\n") - self.log.header2("Full model") - self.log(str(full_model)) - self.log("\n") - - # Write heder - if self.logTerse(): - self.log("\n") - self.log.header1("Sensitivity determination") - - # Loop over energy bands. The number of energy band is still - # fixed ... - for ieng in range(21): - - # Set energies - loge = -1.7 + ieng * 0.2 - emean = pow(10.0, loge) - emin = pow(10.0, loge-0.1) - emax = pow(10.0, loge+0.1) - - # Setup observation(s) - obs = self.set_obs(emin=emin, emax=emax) - - # Determine sensitivity - result = self.get_sensitivity(obs, bkg_model, full_model) - - # Write results - if ieng == 0: - f = open(self.m_outfile, 'w') - writer = csv.DictWriter(f, colnames) - writer.writerow(dict((_,_) for _ in colnames)) - writer.writerow(result) - f.close() - else: - f = open(self.m_outfile, 'a') - writer = csv.DictWriter(f, colnames) - writer.writerow(result) - f.close() - - # Append results - results.append(result) - - # Return - return - - def set_obs(self, lpnt=0.0, bpnt=0.0, emin=0.1, emax=100.0): - """ - Returns an observation container with a single CTA observation. - - Keywords: - lpnt - Galactic longitude of pointing [deg] (default: 0.0) - bpnt - Galactic latitude of pointing [deg] (default: 0.0) - emin - Minimum energy [TeV] (default: 0.1) - emax - Maximum energy [TeV] (default: 100.0) - """ - # Allocate observation container - obs = GObservations() - - # Set single pointing - pntdir = GSkyDir() - pntdir.lb_deg(lpnt, bpnt) - - # Create CTA observation - run = obsutils.set(pntdir, caldb=self.m_caldb, irf=self.m_irf, \ - duration=self.m_duration, \ + """ + This class implements the sensitivity computation script. It derives + from the GammaLib::GApplication class which provides support for parameter + files, command line arguments, and logging. In that way the Python + script behaves just as a regular ctool. + """ + def __init__(self, *argv): + """ + Constructor. + """ + # Set name + self.name = "cssens" + self.version = "0.3.0" + + # Make sure that parfile exists + file = self.parfile() + + # Initialise application + if len(argv) == 0: + GApplication.__init__(self, self.name, self.version) + elif len(argv) ==1: + GApplication.__init__(self, self.name, self.version, *argv) + else: + raise TypeError("Invalid number of arguments given.") + + # Set logger properties + self.logFileOpen() + self.log_header() + self.log.date(True) + + # Return + return + + def __del__(self): + """ + Destructor. + """ + # Write separator into logger + if self.logTerse(): + self.log("\n") + + # Return + return + + def parfile(self): + """ + Check if parfile exists. If parfile does not exist then create a + default parfile. This kluge avoids shipping the cscript with a parfile. + """ + # Set parfile name + parfile = self.name+".par" + + try: + pars = GApplicationPars(parfile) + except: + # Signal if parfile was not found + sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n") + + # Create default parfile + pars = GApplicationPars() + pars.append(GApplicationPar("outfile","f","h","sensitivity.dat","","","Output file name")) + pars.append(GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database")) + pars.append(GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function")) + pars.append(GApplicationPar("type","s","a","point","","","Source model type (point/gauss/shell/disk)")) + pars.append(GApplicationPar("offset","r","a","0.0","0.0","","Source offset angle (deg)")) + pars.append(GApplicationPar("bkg","s","a","$GAMMALIB/share/models/bkg_dummy.txt","","","Background model file function (none=power law for E)")) + pars.append(GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)")) + pars.append(GApplicationPar("rad","r","a","5.0","","","Radius of ROI (deg)")) + pars.append(GApplicationPar("enumbins","i","h","0","","","Number of energy bins (0=unbinned)")) + pars.append(GApplicationPar("npix","i","h","200","","","Number of pixels for binned")) + pars.append(GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)")) + pars.append(GApplicationPar("sigma","r","h","5.0","","","Significance threshold")) + pars.append(GApplicationPar("ts_use","b","h","yes","","","Use TS?")) + pars.append(GApplicationPar("index","r","h","-2.48","","","Assumed spectral index")) + pars.append(GApplicationPar("radius","r","h","0.1","","","Extended source model radius")) + pars.append(GApplicationPar("width","r","h","0.05","","","Extended source model width")) + pars.append(GApplicationPar("max_iter","i","h","50","","","Maximum number of iterations")) + pars.append(GApplicationPar("num_avg","i","h","3","","","Number of iterations for sliding average")) + pars.append_standard() + pars.save(parfile) + + # Return + return + + def get_parameters(self): + """ + Get parameters from parfile. + """ + # Get parameters + self.m_outfile = self["outfile"].filename() + self.m_caldb = self["caldb"].string() + self.m_irf = self["irf"].string() + self.m_type = self["type"].string() + self.m_offset = self["offset"].real() + self.m_bkg = self["bkg"].string() + self.m_duration = self["duration"].real() + self.m_roi = self["rad"].real() + self.m_nbins = self["enumbins"].integer() + self.m_npix = self["npix"].integer() + self.m_binsz = self["binsz"].real() + self.m_ts_thres = self["sigma"].real()*self["sigma"].real() + self.m_use_ts = self["ts_use"].boolean() + self.m_index = self["index"].real() + self.m_radius = self["radius"].real() + self.m_width = self["width"].real() + self.m_max_iter = self["max_iter"].integer() + self.m_num_avg = self["num_avg"].integer() + + # Set some fixed parameters + self.m_log = False # Logging in client tools + self.m_debug = False # Debugging in client tools + + # Return + return + + def execute(self): + """ + Execute the script. + """ + # Run the script + self.run() + + # Return + return + + def run(self): + """ + Run the script. + """ + # Switch screen logging on in debug mode + if self.logDebug(): + self.log.cout(True) + + # Get parameters + self.get_parameters() + + # Write input parameters into logger + if self.logTerse(): + self.log_parameters() + self.log("\n") + + # Initialise script + colnames = ['loge', 'emin', 'emax', 'crab', 'pflux', 'eflux', 'diffSens'] + results = [] + + # Initialise models. Note that we centre the point source at the Galactic + # center as our observation is also centred at the Galactic centre, so + # we're onaxis. + bkg_model = GModels() + full_model = GModels() + bkg_model.append(self.set_bkg_model()) + full_model.append(self.set_bkg_model()) + full_model.append(self.set_src_model(0.0, \ + self.m_offset, \ + type=self.m_type, \ + index=self.m_index)) + + # Write models into logger + if self.logTerse(): + self.log("\n") + self.log.header1("Models") + self.log.header2("Background model") + self.log(str(bkg_model)) + self.log("\n\n") + self.log.header2("Full model") + self.log(str(full_model)) + self.log("\n") + + # Write heder + if self.logTerse(): + self.log("\n") + self.log.header1("Sensitivity determination") + + # Loop over energy bands. The number of energy band is still + # fixed ... + for ieng in range(21): + + # Set energies + loge = -1.7 + ieng * 0.2 + emean = pow(10.0, loge) + emin = pow(10.0, loge-0.1) + emax = pow(10.0, loge+0.1) + + # Setup observation(s) + obs = self.set_obs(emin=emin, emax=emax) + + # Determine sensitivity + result = self.get_sensitivity(obs, bkg_model, full_model) + + # Write results + if ieng == 0: + f = open(self.m_outfile, 'w') + writer = csv.DictWriter(f, colnames) + writer.writerow(dict((_,_) for _ in colnames)) + writer.writerow(result) + f.close() + else: + f = open(self.m_outfile, 'a') + writer = csv.DictWriter(f, colnames) + writer.writerow(result) + f.close() + + # Append results + results.append(result) + + # Return + return + + def set_obs(self, lpnt=0.0, bpnt=0.0, emin=0.1, emax=100.0): + """ + Returns an observation container with a single CTA observation. + + Keywords: + lpnt - Galactic longitude of pointing [deg] (default: 0.0) + bpnt - Galactic latitude of pointing [deg] (default: 0.0) + emin - Minimum energy [TeV] (default: 0.1) + emax - Maximum energy [TeV] (default: 100.0) + """ + # Allocate observation container + obs = GObservations() + + # Set single pointing + pntdir = GSkyDir() + pntdir.lb_deg(lpnt, bpnt) + + # Create CTA observation + run = obsutils.set(pntdir, caldb=self.m_caldb, irf=self.m_irf, \ + duration=self.m_duration, \ emin=emin, emax=emax, rad=self.m_roi) - - # Append observation to container - obs.append(run) - - # Return observation container - return obs - - def set_bkg_model(self, fitidx=False): - """ - Setup standard CTA background model. We use a simple power law here, - scaled to Konrad's E configuration performance table. - - Keywords: - fitidx - Fit spectral index (default: False) - """ - # Define radial component - radial = GCTAModelRadialGauss(3.0) - - # Define spectral component. If a background model is given then - # use a file function. Otherwise use a static power law. - if len(self.m_bkg) > 0: - spectrum = GModelSpectralFunc(self.m_bkg, 1.0) - else: - spectrum = GModelSpectralPlaw(61.8e-6, -1.85, GEnergy(1.0, "TeV")) - if fitidx: - spectrum["Index"].free() - else: - spectrum["Index"].fix() - - # Create background model - model = GCTAModelRadialAcceptance(radial, spectrum) - model.name("Background") - model.instruments("CTA") - - # Return background model - return model - - def set_src_model(self, l, b, flux=1.0, index=-2.48, \ + + # Append observation to container + obs.append(run) + + # Return observation container + return obs + + def set_bkg_model(self, fitidx=False): + """ + Setup standard CTA background model. We use a simple power law here, + scaled to Konrad's E configuration performance table. + + Keywords: + fitidx - Fit spectral index (default: False) + """ + # Define radial component + radial = GCTAModelRadialGauss(3.0) + + # Define spectral component. If a background model is given then + # use a file function. Otherwise use a static power law. + if len(self.m_bkg) > 0: + spectrum = GModelSpectralFunc(self.m_bkg, 1.0) + else: + spectrum = GModelSpectralPlaw(61.8e-6, -1.85, GEnergy(1.0, "TeV")) + if fitidx: + spectrum["Index"].free() + else: + spectrum["Index"].fix() + + # Create background model + model = GCTAModelRadialAcceptance(radial, spectrum) + model.name("Background") + model.instruments("CTA") + + # Return background model + return model + + def set_src_model(self, l, b, flux=1.0, index=-2.48, \ type="point", fitpos=False, fitidx=False): - """ - Returns a single source with Crab-like spectrum. The source flux - can be scaled in Crab units. The Crab spectrum is based on MAGIC - observations (Albert et al. 2008, ApJ, 674, 1037). - - Parameters: - l - Galactic longitude of source location [deg] - b - Galactic latitude of source location [deg] - Keywords: - flux - Source flux [Crabs] - index - Spectral index - type - Source type ("point", "gauss", "disk", "shell") - fitpos - Fit position and size? (default: False) - fitidx - Fit index? (default: False) - """ - # Set source location - location = GSkyDir() - location.lb_deg(l, b) - - # Set source spectrum - spectrum = GModelSpectralPlaw(flux*5.7e-16, index, GEnergy(0.3, "TeV") - if fitidx: - spectrum["Index"].free() - else: - spectrum["Index"].fix() - - # Set source - if type == "point": - spatial = GModelSpatialPointSource(location) - if fitpos: - spatial["RA"].free() - spatial["DEC"].free() - else: - spatial["RA"].fix() - spatial["DEC"].fix() - elif type == "gauss": - spatial = GModelSpatialRadialGauss(location, self.m_radius) - if fitpos: - spatial["RA"].free() - spatial["DEC"].free() - spatial["Sigma"].free() - else: - spatial["RA"].fix() - spatial["DEC"].fix() - spatial["Sigma"].fix() - elif type == "disk": - spatial = GModelSpatialRadialDisk(location, self.m_radius) - if fitpos: - spatial["RA"].free() - spatial["DEC"].free() - spatial["Radius"].free() - else: - spatial["RA"].fix() - spatial["DEC"].fix() - spatial["Radius"].fix() - elif type == "shell": - spatial = GModelSpatialRadialShell(location, self.m_radius, self.m_width) - if fitpos: - spatial["RA"].free() - spatial["DEC"].free() - spatial["Radius"].free() - spatial["Width"].free() - else: - spatial["RA"].fix() - spatial["DEC"].fix() - spatial["Radius"].fix() - spatial["Width"].fix() - else: - self.log("ERROR: Unknown source type '"+type+"'.\n") - return None + """ + Returns a single source with Crab-like spectrum. The source flux + can be scaled in Crab units. The Crab spectrum is based on MAGIC + observations (Albert et al. 2008, ApJ, 674, 1037). + + Parameters: + l - Galactic longitude of source location [deg] + b - Galactic latitude of source location [deg] + Keywords: + flux - Source flux [Crabs] + index - Spectral index + type - Source type ("point", "gauss", "disk", "shell") + fitpos - Fit position and size? (default: False) + fitidx - Fit index? (default: False) + """ + # Set source location + location = GSkyDir() + location.lb_deg(l, b) + + # Set source spectrum + spectrum = GModelSpectralPlaw(flux*5.7e-16, index, GEnergy(0.3, "TeV") + if fitidx: + spectrum["Index"].free() + else: + spectrum["Index"].fix() + + # Set source + if type == "point": + spatial = GModelSpatialPointSource(location) + if fitpos: + spatial["RA"].free() + spatial["DEC"].free() + else: + spatial["RA"].fix() + spatial["DEC"].fix() + elif type == "gauss": + spatial = GModelSpatialRadialGauss(location, self.m_radius) + if fitpos: + spatial["RA"].free() + spatial["DEC"].free() + spatial["Sigma"].free() + else: + spatial["RA"].fix() + spatial["DEC"].fix() + spatial["Sigma"].fix() + elif type == "disk": + spatial = GModelSpatialRadialDisk(location, self.m_radius) + if fitpos: + spatial["RA"].free() + spatial["DEC"].free() + spatial["Radius"].free() + else: + spatial["RA"].fix() + spatial["DEC"].fix() + spatial["Radius"].fix() + elif type == "shell": + spatial = GModelSpatialRadialShell(location, self.m_radius, self.m_width) + if fitpos: + spatial["RA"].free() + spatial["DEC"].free() + spatial["Radius"].free() + spatial["Width"].free() + else: + spatial["RA"].fix() + spatial["DEC"].fix() + spatial["Radius"].fix() + spatial["Width"].fix() + else: + self.log("ERROR: Unknown source type '"+type+"'.\n") + return None source = GModelSky(spatial, spectrum) - - # Set source name - source.name("Test") - - # Return source - return source - - def get_sensitivity(self, obs, bkg_model, full_model): - """ - Determine sensitivity for a given observations. - - Parameters: - obs - Observation container - bkg_model - Background model - full_model - Source model - """ - # Set TeV->erg conversion factor - tev2erg = 1.6021764 - - # Determine energy boundaries from first observation in - # the container - for run in obs: - emin = run.events().ebounds().emin() - emax = run.events().ebounds().emax() - loge = math.log10(math.sqrt(emin.TeV()*emax.TeV())) - erg_mean = math.pow(10.0, loge) * tev2erg - erg_width = (emax.TeV()-emin.TeV()) * tev2erg - break - - # Write header - if self.logTerse(): - self.log("\n") - self.log.header2("Energies: "+str(emin)+" - "+str(emax)) - - # Initialise loop - flux_value = [] - pflux_value = [] - eflux_value = [] - diffSens_value = [] - iter = 0 - test_flux = 0.1 # This is the initial test flux in Crab units - - # Loop until we break - while True: - - # Update iteration counter - iter += 1 - - # Write header - if self.logExplicit(): - self.log.header2("Iteration "+str(iter)) - - # Set test flux - full_model["Test"]['Prefactor'].value(test_flux) - obs.models(full_model) - - # Simulate events - sim = obsutils.sim(obs, nbins=self.m_nbins, seed=iter, \ - binsz=self.m_binsz, npix=self.m_npix, \ - log=self.m_log, debug=self.m_debug) - - # Determine number of events in simulation - nevents = 0.0 - for run in sim: - nevents += run.events().number() - - # Write simulation results - if self.logExplicit(): - self.log.header3("Simulation") - self.log.parformat("Number of simulated events") - self.log(nevents) - self.log("\n") - - # Fit background only - sim.models(bkg_model) - like = obsutils.fit(sim, log=self.m_log, debug=self.m_debug) - result_bgm = like.obs().models().copy() - LogL_bgm = like.opt().value() - npred_bgm = like.obs().npred() - - # Assess quality based on a comparison between Npred and Nevents - quality_bgm = npred_bgm-nevents - - # Write background fit results - if self.logExplicit(): - self.log.header3("Background model fit") - self.log.parformat("log likelihood") - self.log(LogL_bgm) - self.log("\n") - self.log.parformat("Number of predicted events") - self.log(npred_bgm) - self.log("\n") - self.log.parformat("Fit quality") - self.log(quality_bgm) - self.log("\n") - - # Start over if the fit quality was bad - if abs(quality_bgm) > 3.0: - if self.logExplicit(): - self.log("Fit quality outside required range. Start over.\n") - continue - - # Write model fit results - if self.logExplicit(): - for model in result_bgm: - self.log.parformat("Model") - self.log(model.name()) - self.log("\n") - for par in model: - self.log(str(par)+"\n") - - # Fit background and test source - sim.models(full_model) - like = obsutils.fit(sim, log=self.m_log, debug=self.m_debug) - result_all = like.obs().models().copy() - LogL_all = like.opt().value() - npred_all = like.obs().npred() - ts = 2.0*(LogL_bgm-LogL_all) - - # Assess quality based on a comparison between Npred and Nevents - quality_all = npred_all-nevents - - # Write background and test source fit results - if self.logExplicit(): - self.log.header3("Background and test source model fit") - self.log.parformat("Test statistics") - self.log(ts) - self.log("\n") - self.log.parformat("log likelihood") - self.log(LogL_all) - self.log("\n") - self.log.parformat("Number of predicted events") - self.log(npred_all) - self.log("\n") - self.log.parformat("Fit quality") - self.log(quality_all) - self.log("\n") - # - for model in result_all: - self.log.parformat("Model") - self.log(model.name()) - self.log("\n") - for par in model: - self.log(str(par)+"\n") - - # Start over if the fit quality was bad - if abs(quality_all) > 3.0: - if self.logExplicit(): - self.log("Fit quality outside required range. Start over.\n") - continue - - # Start over if TS was non-positive - if ts <= 0.0: - if self.logExplicit(): - self.log("Non positive TS. Start over.\n") - continue - - # Get fitted Crab, photon and energy fluxes - cflux = result_all["Test"]['Prefactor'].value() - cflux_err = result_all["Test"]['Prefactor'].error() - pflux = result_all["Test"].spectral().flux(emin, emax) - eflux = result_all["Test"].spectral().eflux(emin, emax) - - # Compute differential sensitivity - diffSens = pflux / erg_width * erg_mean*erg_mean - - # Compute flux correction factor based on average TS - correct = 1.0 - if self.m_use_ts: - if ts > 0: - correct = math.sqrt(self.m_ts_thres/ts) - else: - if cflux > 0 and cflux_err > 0: - correct = math.sqrt(self.m_ts_thres)/(cflux/cflux_err) - - # Compute extrapolated fluxes - flux = correct * cflux - pflux = correct * pflux - eflux = correct * eflux - diffSens = correct * diffSens - flux_value.append(flux) - pflux_value.append(pflux) - eflux_value.append(eflux) - diffSens_value.append(diffSens) - - # Write background and test source fit results - if self.logExplicit(): - self.log.parformat("Photon flux") - self.log(pflux) - self.log(" ph/cm2/s\n") - self.log.parformat("Energy flux") - self.log(eflux) - self.log(" erg/cm2/s\n") - self.log.parformat("Crab flux") - self.log(cflux*1000.0) - self.log(" mCrab\n") - self.log.parformat("Differential sensitivity") - self.log(diffSens) - self.log(" erg/cm2/s\n") - for model in result_all: - self.log.parformat("Model") - self.log(model.name()) - self.log("\n") - for par in model: - self.log(str(par)+"\n") - elif self.logTerse(): - self.log.parformat("Iteration "+str(iter)) - self.log("TS=") - self.log(ts) - self.log(" ") - self.log("corr=") - self.log(correct) - self.log(" ") - self.log(pflux) - self.log(" ph/cm2/s = ") - self.log(eflux) - self.log(" erg/cm2/s = ") - self.log(cflux*1000.0) - self.log(" mCrab = ") - self.log(diffSens) - self.log(" erg/cm2/s\n") - - # Compute sliding average of extrapolated fitted prefactor, - # photon and energy flux. This damps out fluctuations and - # improves convergence - flux = 0.0 - pflux = 0.0 - eflux = 0.0 - diffSens = 0.0 - num = 0.0 - for k in range(self.m_num_avg): - inx = len(flux_value) - k - 1 - if inx >= 0: - flux += flux_value[inx] - pflux += pflux_value[inx] - eflux += eflux_value[inx] - diffSens += diffSens_value[inx] - num += 1.0 - flux /= num - pflux /= num - eflux /= num - diffSens /= num - - # Compare average flux to last average - if iter > self.m_num_avg: - if test_flux > 0: - ratio = flux/test_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") - break - else: - if self.logTerse(): - self.log(" Flux is zero.\n") - break - - # Use average for next iteration - test_flux = flux - - # Exit loop if number of trials exhausted - if (iter >= self.m_max_iter): - if self.logTerse(): - self.log(" Test ended after "+str(self.m_max_iter)+" iterations.\n") - break - - # Write fit results - if self.logTerse(): - self.log.header3("Fit results") - self.log.parformat("Test statistics") - self.log(ts) - self.log("\n") - self.log.parformat("Photon flux") - self.log(pflux) - self.log(" ph/cm2/s\n") - self.log.parformat("Energy flux") - self.log(eflux) - self.log(" erg/cm2/s\n") - self.log.parformat("Crab flux") - self.log(cflux*1000.0) - self.log(" mCrab\n") - self.log.parformat("Differential sensitivity") - self.log(diffSens) - self.log(" erg/cm2/s\n") - self.log.parformat("Number of simulated events") - self.log(nevents) - self.log("\n") - self.log.header3("Background and test source model fitting") - self.log.parformat("log likelihood") - self.log(LogL_all) - self.log("\n") - self.log.parformat("Number of predicted events") - self.log(npred_all) - self.log("\n") - for model in result_all: - self.log.parformat("Model") - self.log(model.name()) - self.log("\n") - for par in model: - self.log(str(par)+"\n") - self.log.header3("Background model fit") - self.log.parformat("log likelihood") - self.log(LogL_bgm) - self.log("\n") - self.log.parformat("Number of predicted events") - self.log(npred_bgm) - self.log("\n") - for model in result_bgm: - self.log.parformat("Model") - self.log(model.name()) - self.log("\n") - for par in model: - self.log(str(par)+"\n") - - # Store result - result = {'loge': loge, 'emin': emin.TeV(), 'emax': emax.TeV(), \ - 'crab': flux, 'pflux': pflux, 'eflux': eflux, \ - 'diffSens': diffSens} - - # Return result - return result + + # Set source name + source.name("Test") + + # Return source + return source + + def get_sensitivity(self, obs, bkg_model, full_model): + """ + Determine sensitivity for a given observations. + + Parameters: + obs - Observation container + bkg_model - Background model + full_model - Source model + """ + # Set TeV->erg conversion factor + tev2erg = 1.6021764 + + # Determine energy boundaries from first observation in + # the container + for run in obs: + emin = run.events().ebounds().emin() + emax = run.events().ebounds().emax() + loge = math.log10(math.sqrt(emin.TeV()*emax.TeV())) + erg_mean = math.pow(10.0, loge) * tev2erg + erg_width = (emax.TeV()-emin.TeV()) * tev2erg + break + + # Write header + if self.logTerse(): + self.log("\n") + self.log.header2("Energies: "+str(emin)+" - "+str(emax)) + + # Initialise loop + flux_value = [] + pflux_value = [] + eflux_value = [] + diffSens_value = [] + iter = 0 + test_flux = 0.1 # This is the initial test flux in Crab units + + # Loop until we break + while True: + + # Update iteration counter + iter += 1 + + # Write header + if self.logExplicit(): + self.log.header2("Iteration "+str(iter)) + + # Set test flux + full_model["Test"]['Prefactor'].value(test_flux) + obs.models(full_model) + + # Simulate events + sim = obsutils.sim(obs, nbins=self.m_nbins, seed=iter, \ + binsz=self.m_binsz, npix=self.m_npix, \ + log=self.m_log, debug=self.m_debug) + + # Determine number of events in simulation + nevents = 0.0 + for run in sim: + nevents += run.events().number() + + # Write simulation results + if self.logExplicit(): + self.log.header3("Simulation") + self.log.parformat("Number of simulated events") + self.log(nevents) + self.log("\n") + + # Fit background only + sim.models(bkg_model) + like = obsutils.fit(sim, log=self.m_log, debug=self.m_debug) + result_bgm = like.obs().models().copy() + LogL_bgm = like.opt().value() + npred_bgm = like.obs().npred() + + # Assess quality based on a comparison between Npred and Nevents + quality_bgm = npred_bgm-nevents + + # Write background fit results + if self.logExplicit(): + self.log.header3("Background model fit") + self.log.parformat("log likelihood") + self.log(LogL_bgm) + self.log("\n") + self.log.parformat("Number of predicted events") + self.log(npred_bgm) + self.log("\n") + self.log.parformat("Fit quality") + self.log(quality_bgm) + self.log("\n") + + # Start over if the fit quality was bad + if abs(quality_bgm) > 3.0: + if self.logExplicit(): + self.log("Fit quality outside required range. Start over.\n") + continue + + # Write model fit results + if self.logExplicit(): + for model in result_bgm: + self.log.parformat("Model") + self.log(model.name()) + self.log("\n") + for par in model: + self.log(str(par)+"\n") + + # Fit background and test source + sim.models(full_model) + like = obsutils.fit(sim, log=self.m_log, debug=self.m_debug) + result_all = like.obs().models().copy() + LogL_all = like.opt().value() + npred_all = like.obs().npred() + ts = 2.0*(LogL_bgm-LogL_all) + + # Assess quality based on a comparison between Npred and Nevents + quality_all = npred_all-nevents + + # Write background and test source fit results + if self.logExplicit(): + self.log.header3("Background and test source model fit") + self.log.parformat("Test statistics") + self.log(ts) + self.log("\n") + self.log.parformat("log likelihood") + self.log(LogL_all) + self.log("\n") + self.log.parformat("Number of predicted events") + self.log(npred_all) + self.log("\n") + self.log.parformat("Fit quality") + self.log(quality_all) + self.log("\n") + # + for model in result_all: + self.log.parformat("Model") + self.log(model.name()) + self.log("\n") + for par in model: + self.log(str(par)+"\n") + + # Start over if the fit quality was bad + if abs(quality_all) > 3.0: + if self.logExplicit(): + self.log("Fit quality outside required range. Start over.\n") + continue + + # Start over if TS was non-positive + if ts <= 0.0: + if self.logExplicit(): + self.log("Non positive TS. Start over.\n") + continue + + # Get fitted Crab, photon and energy fluxes + cflux = result_all["Test"]['Prefactor'].value() + cflux_err = result_all["Test"]['Prefactor'].error() + pflux = result_all["Test"].spectral().flux(emin, emax) + eflux = result_all["Test"].spectral().eflux(emin, emax) + + # Compute differential sensitivity + diffSens = pflux / erg_width * erg_mean*erg_mean + + # Compute flux correction factor based on average TS + correct = 1.0 + if self.m_use_ts: + if ts > 0: + correct = math.sqrt(self.m_ts_thres/ts) + else: + if cflux > 0 and cflux_err > 0: + correct = math.sqrt(self.m_ts_thres)/(cflux/cflux_err) + + # Compute extrapolated fluxes + flux = correct * cflux + pflux = correct * pflux + eflux = correct * eflux + diffSens = correct * diffSens + flux_value.append(flux) + pflux_value.append(pflux) + eflux_value.append(eflux) + diffSens_value.append(diffSens) + + # Write background and test source fit results + if self.logExplicit(): + self.log.parformat("Photon flux") + self.log(pflux) + self.log(" ph/cm2/s\n") + self.log.parformat("Energy flux") + self.log(eflux) + self.log(" erg/cm2/s\n") + self.log.parformat("Crab flux") + self.log(cflux*1000.0) + self.log(" mCrab\n") + self.log.parformat("Differential sensitivity") + self.log(diffSens) + self.log(" erg/cm2/s\n") + for model in result_all: + self.log.parformat("Model") + self.log(model.name()) + self.log("\n") + for par in model: + self.log(str(par)+"\n") + elif self.logTerse(): + self.log.parformat("Iteration "+str(iter)) + self.log("TS=") + self.log(ts) + self.log(" ") + self.log("corr=") + self.log(correct) + self.log(" ") + self.log(pflux) + self.log(" ph/cm2/s = ") + self.log(eflux) + self.log(" erg/cm2/s = ") + self.log(cflux*1000.0) + self.log(" mCrab = ") + self.log(diffSens) + self.log(" erg/cm2/s\n") + + # Compute sliding average of extrapolated fitted prefactor, + # photon and energy flux. This damps out fluctuations and + # improves convergence + flux = 0.0 + pflux = 0.0 + eflux = 0.0 + diffSens = 0.0 + num = 0.0 + for k in range(self.m_num_avg): + inx = len(flux_value) - k - 1 + if inx >= 0: + flux += flux_value[inx] + pflux += pflux_value[inx] + eflux += eflux_value[inx] + diffSens += diffSens_value[inx] + num += 1.0 + flux /= num + pflux /= num + eflux /= num + diffSens /= num + + # Compare average flux to last average + if iter > self.m_num_avg: + if test_flux > 0: + ratio = flux/test_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") + break + else: + if self.logTerse(): + self.log(" Flux is zero.\n") + break + + # Use average for next iteration + test_flux = flux + + # Exit loop if number of trials exhausted + if (iter >= self.m_max_iter): + if self.logTerse(): + self.log(" Test ended after "+str(self.m_max_iter)+" iterations.\n") + break + + # Write fit results + if self.logTerse(): + self.log.header3("Fit results") + self.log.parformat("Test statistics") + self.log(ts) + self.log("\n") + self.log.parformat("Photon flux") + self.log(pflux) + self.log(" ph/cm2/s\n") + self.log.parformat("Energy flux") + self.log(eflux) + self.log(" erg/cm2/s\n") + self.log.parformat("Crab flux") + self.log(cflux*1000.0) + self.log(" mCrab\n") + self.log.parformat("Differential sensitivity") + self.log(diffSens) + self.log(" erg/cm2/s\n") + self.log.parformat("Number of simulated events") + self.log(nevents) + self.log("\n") + self.log.header3("Background and test source model fitting") + self.log.parformat("log likelihood") + self.log(LogL_all) + self.log("\n") + self.log.parformat("Number of predicted events") + self.log(npred_all) + self.log("\n") + for model in result_all: + self.log.parformat("Model") + self.log(model.name()) + self.log("\n") + for par in model: + self.log(str(par)+"\n") + self.log.header3("Background model fit") + self.log.parformat("log likelihood") + self.log(LogL_bgm) + self.log("\n") + self.log.parformat("Number of predicted events") + self.log(npred_bgm) + self.log("\n") + for model in result_bgm: + self.log.parformat("Model") + self.log(model.name()) + self.log("\n") + for par in model: + self.log(str(par)+"\n") + + # Store result + result = {'loge': loge, 'emin': emin.TeV(), 'emax': emax.TeV(), \ + 'crab': flux, 'pflux': pflux, 'eflux': eflux, \ + 'diffSens': diffSens} + + # Return result + return result # ======================== # # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - CTA ctlike-based sensitivity estimator. - """ - # Create instance of application - app = cssens(sys.argv) - - # Run application - app.execute() - \ No newline at end of file + """ + CTA ctlike-based sensitivity estimator. + """ + # Create instance of application + app = cssens(sys.argv) + + # Run application + app.execute() + \ No newline at end of file diff --git a/scripts/cstsdist.py b/scripts/cstsdist.py index c5a0c36..95ec225 100755 --- a/scripts/cstsdist.py +++ b/scripts/cstsdist.py @@ -3,7 +3,7 @@ # This script generates the TS distribution for a particular model based # on Monte-Carlo simulations. # -# Copyright (C) 2011-2013 Jurgen Knodlseder +# Copyright (C) 2011-2014 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 @@ -31,486 +31,486 @@ import math # cstsdist class # # ============== # class cstsdist(GApplication): - """ - This class implements the TS distribution generation script. It derives - from the GammaLib::GApplication class which provides support for parameter - files, command line arguments, and logging. In that way the Python - script behaves just as a regular ctool. - """ - def __init__(self, *argv): - """ - Constructor. - """ - # Set name - self.name = "cstsdist" - self.version = "0.1.0" - - # Initialise some members - self.obs = None - self.bkg_model = None - self.full_model = None - - # Make sure that parfile exists - file = self.parfile() + """ + This class implements the TS distribution generation script. It derives + from the GammaLib::GApplication class which provides support for parameter + files, command line arguments, and logging. In that way the Python + script behaves just as a regular ctool. + """ + def __init__(self, *argv): + """ + Constructor. + """ + # Set name + self.name = "cstsdist" + self.version = "0.1.0" + + # Initialise some members + self.obs = None + self.bkg_model = None + self.full_model = None + + # Make sure that parfile exists + file = self.parfile() - # Initialise application - if len(argv) == 0: - GApplication.__init__(self, self.name, self.version) - elif len(argv) ==1: - GApplication.__init__(self, self.name, self.version, *argv) - else: - raise TypeError("Invalid number of arguments given.") + # Initialise application + if len(argv) == 0: + GApplication.__init__(self, self.name, self.version) + elif len(argv) ==1: + GApplication.__init__(self, self.name, self.version, *argv) + else: + raise TypeError("Invalid number of arguments given.") - # Set logger properties - self.log_header() - self.log.date(True) + # Set logger properties + self.log_header() + self.log.date(True) - # Return - return - - def __del__(self): - """ - Destructor. - """ - # Write separator into logger - if self.logTerse(): - self.log("\n") - - # Return - return + # Return + return + + def __del__(self): + """ + Destructor. + """ + # Write separator into logger + if self.logTerse(): + self.log("\n") + + # Return + return - def parfile(self): - """ - Check if parfile exists. If parfile does not exist then create a - default parfile. This kluge avoids shipping the cscript with a parfile. - """ - # Set parfile name - parfile = self.name+".par" - - try: - pars = GApplicationPars(parfile) - except: - # Signal if parfile was not found - sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n") - - # Create default parfile - pars = GApplicationPars() - pars.append(GApplicationPar("outfile","f","h","ts.dat","","","Output file name")) - pars.append(GApplicationPar("ntrials","i","a","10","","","Number of trials")) - pars.append(GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database")) - pars.append(GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function")) - pars.append(GApplicationPar("type","s","a","point","","","Source model type (point/gauss/shell/disk)")) - pars.append(GApplicationPar("index","r","h","-2.48","","","Spectral index")) - pars.append(GApplicationPar("offset","r","a","0.0","0.0","","Source offset angle (deg)")) - pars.append(GApplicationPar("bkg","s","a","$GAMMALIB/share/models/bkg_dummy.txt","","","Background model file function (none=power law for E)")) - pars.append(GApplicationPar("emin","r","a","0.1","0.0","","Lower energy limit (TeV)")) - pars.append(GApplicationPar("emax","r","a","100.0","0.0","","Upper energy limit (TeV)")) - pars.append(GApplicationPar("enumbins","i","a","0","","","Number of energy bins (0=unbinned)")) - pars.append(GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)")) - pars.append(GApplicationPar("rad","r","h","5.0","","","Radius of ROI (deg)")) - pars.append(GApplicationPar("npix","i","h","200","","","Number of pixels for binned")) - pars.append(GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)")) - pars.append_standard() - pars.save(parfile) - - # Return - return - - def get_parameters(self): - """ - Get parameters from parfile and setup the observation. - """ - # Get parameters - self.m_outfile = self["outfile"].filename() - self.m_ntrials = self["ntrials"].integer() - self.m_caldb = self["caldb"].string() - self.m_irf = self["irf"].string() - self.m_type = self["type"].string() - self.m_index = self["index"].real() - self.m_offset = self["offset"].real() - self.m_bkg = self["bkg"].string() - self.m_emin = self["emin"].real() - self.m_emax = self["emax"].real() - self.m_enumbins = self["enumbins"].integer() - self.m_duration = self["duration"].real() - self.m_rad = self["rad"].real() - self.m_npix = self["npix"].integer() - self.m_binsz = self["binsz"].real() + def parfile(self): + """ + Check if parfile exists. If parfile does not exist then create a + default parfile. This kluge avoids shipping the cscript with a parfile. + """ + # Set parfile name + parfile = self.name+".par" + + try: + pars = GApplicationPars(parfile) + except: + # Signal if parfile was not found + sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n") + + # Create default parfile + pars = GApplicationPars() + pars.append(GApplicationPar("outfile","f","h","ts.dat","","","Output file name")) + pars.append(GApplicationPar("ntrials","i","a","10","","","Number of trials")) + pars.append(GApplicationPar("caldb","s","a","$GAMMALIB/share/caldb/cta","","","Calibration database")) + pars.append(GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function")) + pars.append(GApplicationPar("type","s","a","point","","","Source model type (point/gauss/shell/disk)")) + pars.append(GApplicationPar("index","r","h","-2.48","","","Spectral index")) + pars.append(GApplicationPar("offset","r","a","0.0","0.0","","Source offset angle (deg)")) + pars.append(GApplicationPar("bkg","s","a","$GAMMALIB/share/models/bkg_dummy.txt","","","Background model file function (none=power law for E)")) + pars.append(GApplicationPar("emin","r","a","0.1","0.0","","Lower energy limit (TeV)")) + pars.append(GApplicationPar("emax","r","a","100.0","0.0","","Upper energy limit (TeV)")) + pars.append(GApplicationPar("enumbins","i","a","0","","","Number of energy bins (0=unbinned)")) + pars.append(GApplicationPar("duration","r","a","180000.0","","","Effective exposure time (s)")) + pars.append(GApplicationPar("rad","r","h","5.0","","","Radius of ROI (deg)")) + pars.append(GApplicationPar("npix","i","h","200","","","Number of pixels for binned")) + pars.append(GApplicationPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)")) + pars.append_standard() + pars.save(parfile) + + # Return + return + + def get_parameters(self): + """ + Get parameters from parfile and setup the observation. + """ + # Get parameters + self.m_outfile = self["outfile"].filename() + self.m_ntrials = self["ntrials"].integer() + self.m_caldb = self["caldb"].string() + self.m_irf = self["irf"].string() + self.m_type = self["type"].string() + self.m_index = self["index"].real() + self.m_offset = self["offset"].real() + self.m_bkg = self["bkg"].string() + self.m_emin = self["emin"].real() + self.m_emax = self["emax"].real() + self.m_enumbins = self["enumbins"].integer() + self.m_duration = self["duration"].real() + self.m_rad = self["rad"].real() + self.m_npix = self["npix"].integer() + self.m_binsz = self["binsz"].real() - # Set some fixed parameters - self.m_log = False # Logging in client tools - self.m_debug = False # Debugging in client tools - - # Setup observation - self.obs = self.set_obs(emin=self.m_emin, emax=self.m_emax) - - # Initialise models. Note that we centre the point source at the Galactic - # center as our observation is also centred at the Galactic centre, so - # we're onaxis. - self.bkg_model = GModels() - self.full_model = GModels() - self.bkg_model.append(self.set_bkg_model()) - self.full_model.append(self.set_bkg_model()) - self.full_model.append(self.set_src_model(0.0, self.m_offset, \ - flux=0.010, \ - type=self.m_type, \ - index=self.m_index)) + # Set some fixed parameters + self.m_log = False # Logging in client tools + self.m_debug = False # Debugging in client tools + + # Setup observation + self.obs = self.set_obs(emin=self.m_emin, emax=self.m_emax) + + # Initialise models. Note that we centre the point source at the Galactic + # center as our observation is also centred at the Galactic centre, so + # we're onaxis. + self.bkg_model = GModels() + self.full_model = GModels() + self.bkg_model.append(self.set_bkg_model()) + self.full_model.append(self.set_bkg_model()) + self.full_model.append(self.set_src_model(0.0, self.m_offset, \ + flux=0.010, \ + type=self.m_type, \ + index=self.m_index)) - # Attach background model to observation container - self.obs.models(self.bkg_model) - - # Return - return - - def models(self, models): - """ - Set model. - """ - # Copy models - self.model = models.copy() - - # Return - return - - def execute(self): - """ - Execute the script. - """ - # Run the script - self.run() - - # Return - return + # Attach background model to observation container + self.obs.models(self.bkg_model) + + # Return + return + + def models(self, models): + """ + Set model. + """ + # Copy models + self.model = models.copy() + + # Return + return + + def execute(self): + """ + Execute the script. + """ + # Run the script + self.run() + + # Return + return - def run(self): - """ - Run the script. - """ - # Switch screen logging on in debug mode - if self.logDebug(): - self.log.cout(True) + def run(self): + """ + Run the script. + """ + # Switch screen logging on in debug mode + if self.logDebug(): + self.log.cout(True) - # Get parameters - self.get_parameters() - - # Write input parameters into logger - if self.logTerse(): - self.log_parameters() - self.log("\n") - - # Write observation into logger - if self.logTerse(): - self.log("\n") - self.log.header1("Observation") - self.log(str(self.obs)) - self.log("\n") + # Get parameters + self.get_parameters() + + # Write input parameters into logger + if self.logTerse(): + self.log_parameters() + self.log("\n") + + # Write observation into logger + if self.logTerse(): + self.log("\n") + self.log.header1("Observation") + self.log(str(self.obs)) + self.log("\n") - # Write models into logger - if self.logTerse(): - self.log("\n") - self.log.header1("Test model") - self.log(str(self.full_model)) - self.log("\n") + # Write models into logger + if self.logTerse(): + self.log("\n") + self.log.header1("Test model") + self.log(str(self.full_model)) + self.log("\n") - # Write header - if self.logTerse(): - self.log("\n") - self.log.header1("Generate TS distribution") + # Write header + if self.logTerse(): + self.log("\n") + self.log.header1("Generate TS distribution") - # Loop over trials - for seed in range(self.m_ntrials): - - # Make a trial - result = self.trial(seed) - - # Write out result immediately - if seed == 0: - file = open(self.m_outfile, 'w') - writer = csv.DictWriter(file, result['colnames']) - writer.writerow(dict((_,_) for _ in result['colnames'])) - else: - file = open(self.m_outfile, 'a') - writer = csv.DictWriter(file, result['colnames']) - writer.writerow(result['values']) - file.close() - - # Return - return - - def set_obs(self, lpnt=0.0, bpnt=0.0, emin=0.1, emax=100.0): - """ - Returns an observation container with a single CTA observation. - - Keywords: - lpnt - Galactic longitude of pointing [deg] (default: 0.0) - bpnt - Galactic latitude of pointing [deg] (default: 0.0) - emin - Minimum energy [TeV] (default: 0.1) - emax - Maximum energy [TeV] (default: 100.0) - """ - # Allocate observation container - obs = GObservations() - - # Set single pointing - pntdir = GSkyDir() - pntdir.lb_deg(lpnt, bpnt) - - # Create CTA observation - run = obsutils.set(pntdir, caldb=self.m_caldb, irf=self.m_irf, \ - duration=self.m_duration, \ + # Loop over trials + for seed in range(self.m_ntrials): + + # Make a trial + result = self.trial(seed) + + # Write out result immediately + if seed == 0: + file = open(self.m_outfile, 'w') + writer = csv.DictWriter(file, result['colnames']) + writer.writerow(dict((_,_) for _ in result['colnames'])) + else: + file = open(self.m_outfile, 'a') + writer = csv.DictWriter(file, result['colnames']) + writer.writerow(result['values']) + file.close() + + # Return + return + + def set_obs(self, lpnt=0.0, bpnt=0.0, emin=0.1, emax=100.0): + """ + Returns an observation container with a single CTA observation. + + Keywords: + lpnt - Galactic longitude of pointing [deg] (default: 0.0) + bpnt - Galactic latitude of pointing [deg] (default: 0.0) + emin - Minimum energy [TeV] (default: 0.1) + emax - Maximum energy [TeV] (default: 100.0) + """ + # Allocate observation container + obs = GObservations() + + # Set single pointing + pntdir = GSkyDir() + pntdir.lb_deg(lpnt, bpnt) + + # Create CTA observation + run = obsutils.set(pntdir, caldb=self.m_caldb, irf=self.m_irf, \ + duration=self.m_duration, \ emin=emin, emax=emax, rad=self.m_rad) - - # Append observation to container - obs.append(run) - - # Return observation container - return obs - - def set_bkg_model(self, fitsigma=False): - """ - Setup CTA background model. - """ - # Define radial component - radial = GCTAModelRadialGauss(3.0) - if fitsigma: - radial["Sigma"].free() - else: - radial["Sigma"].fix() - - # Define spectral component - spectrum = GModelSpectralFunc(self.m_bkg, 1.0) - - # Create background model - model = GCTAModelRadialAcceptance(radial, spectrum) - model.name("Background") - model.instruments("CTA") - - # Return background model - return model - - def set_src_model(self, l, b, flux=1.0, index=-2.48, \ + + # Append observation to container + obs.append(run) + + # Return observation container + return obs + + def set_bkg_model(self, fitsigma=False): + """ + Setup CTA background model. + """ + # Define radial component + radial = GCTAModelRadialGauss(3.0) + if fitsigma: + radial["Sigma"].free() + else: + radial["Sigma"].fix() + + # Define spectral component + spectrum = GModelSpectralFunc(self.m_bkg, 1.0) + + # Create background model + model = GCTAModelRadialAcceptance(radial, spectrum) + model.name("Background") + model.instruments("CTA") + + # Return background model + return model + + def set_src_model(self, l, b, flux=1.0, index=-2.48, \ type="point", sigma=1.0, radius=1.0, width=0.1, \ fitpos=False, fitidx=False): - """ - Returns a single source with Crab-like spectrum. The source flux - can be scaled in Crab units. The Crab spectrum is based on MAGIC - observations (Albert et al. 2008, ApJ, 674, 1037). + """ + Returns a single source with Crab-like spectrum. The source flux + can be scaled in Crab units. The Crab spectrum is based on MAGIC + observations (Albert et al. 2008, ApJ, 674, 1037). - Parameters: - l - Galactic longitude of source location [deg] - b - Galactic latitude of source location [deg] - Keywords: - flux - Source flux [Crabs] - index - Spectral index - type - Source type ("point", "gauss", "disk", "shell") - sigma - Gaussian sigma (for type="gauss") - radius - Disk or shell inner radius [deg] (for type="disk" and type="shell") - width - Shell width [deg] (for type="shell") - fitpos - Fit position? (default: True) - fitidx - Fit index? (default: True) - """ - # Set source location - location = GSkyDir() - location.lb_deg(l, b) - - # Set source spectrum - spectrum = GModelSpectralPlaw(flux*5.7e-16, index, GEnergy(0.3, "TeV")) - if fitidx: - spectrum["Index"].free() - else: - spectrum["Index"].fix() + Parameters: + l - Galactic longitude of source location [deg] + b - Galactic latitude of source location [deg] + Keywords: + flux - Source flux [Crabs] + index - Spectral index + type - Source type ("point", "gauss", "disk", "shell") + sigma - Gaussian sigma (for type="gauss") + radius - Disk or shell inner radius [deg] (for type="disk" and type="shell") + width - Shell width [deg] (for type="shell") + fitpos - Fit position? (default: True) + fitidx - Fit index? (default: True) + """ + # Set source location + location = GSkyDir() + location.lb_deg(l, b) + + # Set source spectrum + spectrum = GModelSpectralPlaw(flux*5.7e-16, index, GEnergy(0.3, "TeV")) + if fitidx: + spectrum["Index"].free() + else: + spectrum["Index"].fix() - # Set source - if type == "point": - spatial = GModelSpatialPointSource(location) - if fitpos: - spatial[0].free() - spatial[1].free() - elif type == "gauss": - spatial = GModelSpatialRadialGauss(location, sigma) - if fitpos: - spatial[0].free() - spatial[1].free() - elif type == "disk": - spatial = GModelSpatialRadialDisk(location, radius) - if fitpos: - spatial[0].free() - spatial[1].free() - elif type == "shell": - spatial = GModelSpatialRadialShell(location, radius, width) - if fitpos: - spatial[0].free() - spatial[1].free() - else: - self.log("ERROR: Unknown source type '"+type+"'.\n") - return None - source = GModelSky(spatial, spectrum) - - # Set source name - source.name("Test") - - # Return source - return source - - def trial(self, seed): - """ - Create the TS for a single trial. - - Parameters: - seed - Random number generator seed - """ - # Write header - if self.logExplicit(): - self.log.header2("Trial "+str(seed+1)) + # Set source + if type == "point": + spatial = GModelSpatialPointSource(location) + if fitpos: + spatial[0].free() + spatial[1].free() + elif type == "gauss": + spatial = GModelSpatialRadialGauss(location, sigma) + if fitpos: + spatial[0].free() + spatial[1].free() + elif type == "disk": + spatial = GModelSpatialRadialDisk(location, radius) + if fitpos: + spatial[0].free() + spatial[1].free() + elif type == "shell": + spatial = GModelSpatialRadialShell(location, radius, width) + if fitpos: + spatial[0].free() + spatial[1].free() + else: + self.log("ERROR: Unknown source type '"+type+"'.\n") + return None + source = GModelSky(spatial, spectrum) + + # Set source name + source.name("Test") + + # Return source + return source + + def trial(self, seed): + """ + Create the TS for a single trial. + + Parameters: + seed - Random number generator seed + """ + # Write header + if self.logExplicit(): + self.log.header2("Trial "+str(seed+1)) - # Simulate events - sim = obsutils.sim(self.obs, \ - nbins=self.m_enumbins, \ - seed=seed, \ - binsz=self.m_binsz, \ - npix=self.m_npix, \ - log=self.m_log, debug=self.m_debug) + # Simulate events + sim = obsutils.sim(self.obs, \ + nbins=self.m_enumbins, \ + seed=seed, \ + binsz=self.m_binsz, \ + npix=self.m_npix, \ + log=self.m_log, debug=self.m_debug) - # Determine number of events in simulation - nevents = 0.0 - for run in sim: - nevents += run.events().number() + # Determine number of events in simulation + nevents = 0.0 + for run in sim: + nevents += run.events().number() - # Write simulation results - if self.logExplicit(): - self.log.header3("Simulation") - self.log.parformat("Number of simulated events") - self.log(nevents) - self.log("\n") + # Write simulation results + if self.logExplicit(): + self.log.header3("Simulation") + self.log.parformat("Number of simulated events") + self.log(nevents) + self.log("\n") - # Fit background only - sim.models(self.bkg_model) - like_bgm = obsutils.fit(sim, log=self.m_log, debug=self.m_debug) - result_bgm = like_bgm.obs().models() - LogL_bgm = like_bgm.opt().value() - npred_bgm = like_bgm.obs().npred() + # Fit background only + sim.models(self.bkg_model) + like_bgm = obsutils.fit(sim, log=self.m_log, debug=self.m_debug) + result_bgm = like_bgm.obs().models() + LogL_bgm = like_bgm.opt().value() + npred_bgm = like_bgm.obs().npred() - # Write background fit results - if self.logExplicit(): - self.log.header3("Background model fit") - self.log.parformat("log likelihood") - self.log(LogL_bgm) - self.log("\n") - self.log.parformat("Number of predicted events") - self.log(npred_bgm) - self.log("\n") - for model in result_bgm: - self.log.parformat("Model") - self.log(model.name()) - self.log("\n") - for par in model: - self.log(str(par)+"\n") + # Write background fit results + if self.logExplicit(): + self.log.header3("Background model fit") + self.log.parformat("log likelihood") + self.log(LogL_bgm) + self.log("\n") + self.log.parformat("Number of predicted events") + self.log(npred_bgm) + self.log("\n") + for model in result_bgm: + self.log.parformat("Model") + self.log(model.name()) + self.log("\n") + for par in model: + self.log(str(par)+"\n") - # Fit background and test source - sim.models(self.full_model) - like_all = obsutils.fit(sim, log=self.m_log, debug=self.m_debug) - result_all = like_all.obs().models() - LogL_all = like_all.opt().value() - npred_all = like_all.obs().npred() - ts = 2.0*(LogL_bgm-LogL_all) - - # Write background and test source fit results - if self.logExplicit(): - self.log.header3("Background and test source model fit") - self.log.parformat("Test statistics") - self.log(ts) - self.log("\n") - self.log.parformat("log likelihood") - self.log(LogL_all) - self.log("\n") - self.log.parformat("Number of predicted events") - self.log(npred_all) - self.log("\n") - for model in result_all: - self.log.parformat("Model") - self.log(model.name()) - self.log("\n") - for par in model: - self.log(str(par)+"\n") - - # Write result - elif self.logTerse(): - self.log.parformat("Trial "+str(seed)) - self.log("TS=") - self.log(ts) - self.log(" Prefactor=") - self.log(result_all["Test"]["Prefactor"].value()) - self.log("+/-") - self.log(result_all["Test"]["Prefactor"].error()) - self.log("\n") - - # Initialise results - colnames = [] - values = {} - - # Set TS value - colnames.append("TS") - values["TS"] = ts + # Fit background and test source + sim.models(self.full_model) + like_all = obsutils.fit(sim, log=self.m_log, debug=self.m_debug) + result_all = like_all.obs().models() + LogL_all = like_all.opt().value() + npred_all = like_all.obs().npred() + ts = 2.0*(LogL_bgm-LogL_all) + + # Write background and test source fit results + if self.logExplicit(): + self.log.header3("Background and test source model fit") + self.log.parformat("Test statistics") + self.log(ts) + self.log("\n") + self.log.parformat("log likelihood") + self.log(LogL_all) + self.log("\n") + self.log.parformat("Number of predicted events") + self.log(npred_all) + self.log("\n") + for model in result_all: + self.log.parformat("Model") + self.log(model.name()) + self.log("\n") + for par in model: + self.log(str(par)+"\n") + + # Write result + elif self.logTerse(): + self.log.parformat("Trial "+str(seed)) + self.log("TS=") + self.log(ts) + self.log(" Prefactor=") + self.log(result_all["Test"]["Prefactor"].value()) + self.log("+/-") + self.log(result_all["Test"]["Prefactor"].error()) + self.log("\n") + + # Initialise results + colnames = [] + values = {} + + # Set TS value + colnames.append("TS") + values["TS"] = ts - # Set logL for background fit - colnames.append("LogL_bgm") - values["LogL_bgm"] = LogL_bgm + # Set logL for background fit + colnames.append("LogL_bgm") + values["LogL_bgm"] = LogL_bgm - # Set logL for full fit - colnames.append("LogL_all") - values["LogL_all"] = LogL_all + # Set logL for full fit + colnames.append("LogL_all") + values["LogL_all"] = LogL_all - # Set Nevents - colnames.append("Nevents") - values["Nevents"] = nevents + # Set Nevents + colnames.append("Nevents") + values["Nevents"] = nevents - # Set Npred for background fit - colnames.append("Npred_bkg") - values["Npred_bkg"] = npred_bgm + # Set Npred for background fit + colnames.append("Npred_bkg") + values["Npred_bkg"] = npred_bgm - # Set Npred for full fit - colnames.append("Npred_all") - values["Npred_all"] = npred_all - - # Gather free full fit parameters - for i in range(result_all.size()): - model = result_all[i] - model_name = model.name() - for k in range(model.size()): - par = model[k] - if par.is_free(): - - # Set parameter name - name = model_name+"_"+par.name() - - # Append value - colnames.append(name) - values[name] = par.value() - - # Append error - name = "Unc_"+name - colnames.append(name) - values[name] = par.error() - - # Bundle together results - result = {'colnames': colnames, 'values': values} - - # Return - return result + # Set Npred for full fit + colnames.append("Npred_all") + values["Npred_all"] = npred_all + + # Gather free full fit parameters + for i in range(result_all.size()): + model = result_all[i] + model_name = model.name() + for k in range(model.size()): + par = model[k] + if par.is_free(): + + # Set parameter name + name = model_name+"_"+par.name() + + # Append value + colnames.append(name) + values[name] = par.value() + + # Append error + name = "Unc_"+name + colnames.append(name) + values[name] = par.error() + + # Bundle together results + result = {'colnames': colnames, 'values': values} + + # Return + return result # ======================== # # Main routine entry point # # ======================== # if __name__ == '__main__': - """ - Generates TS distribution. - """ - # Create instance of application - app = cstsdist(sys.argv) - - # Open logfile - app.logFileOpen() - - # Execute application - app.execute() - \ No newline at end of file + """ + Generates TS distribution. + """ + # Create instance of application + app = cstsdist(sys.argv) + + # Open logfile + app.logFileOpen() + + # Execute application + app.execute() + \ No newline at end of file diff --git a/scripts/cta_kb_bkgrates.py b/scripts/cta_kb_bkgrates.py index bea8e07..75386e3 100755 --- a/scripts/cta_kb_bkgrates.py +++ b/scripts/cta_kb_bkgrates.py @@ -4,7 +4,7 @@ # into a column separated file function table. The table is automatically # named "bkg_XXX.txt" where "XXX" is for example "kb_E_50h_v3". # -# Copyright (C) 2011 Jurgen Knodlseder +# Copyright (C) 2011-2014 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 @@ -29,78 +29,78 @@ import sys # Generate file function from model # # ================================= # def make_file_function(irfname, filename): - """ - Generate file function from model. - - Parameters: - irfname - IRF filename - filename - File function filename - """ - # Open performance file and file function - irf = open(irfname, "r") - file = open(filename, "w") - - # Read data - for row in irf: - - # Skip header and trailer - if row.find("log(E)") != -1: - continue - if row.find("--------") != -1: - break - - # Extract data - fields = row.replace(" ",",").replace(",,",",").rstrip("\n").split(",") - logE = float(fields[0]) - aeff = float(fields[1]) - r68 = float(fields[2]) - r80 = float(fields[3]) - eres = float(fields[4]) - bgd = float(fields[5]) - sens = float(fields[6]) - - # Compute energy (in MeV) - energy = math.pow(10.0, logE)*1.0e6 - emin = math.pow(10.0, logE-0.1)*1.0e6 - emax = math.pow(10.0, logE+0.1)*1.0e6 - ewidth = emax-emin - - # Compute solid angle of r80 region in sr - omega = 2.0 * math.pi * (1.0 - math.cos(math.radians(r80))) - - # Compute background rate per steradian and MeV - bkg_rate = bgd/omega/ewidth + """ + Generate file function from model. + + Parameters: + irfname - IRF filename + filename - File function filename + """ + # Open performance file and file function + irf = open(irfname, "r") + file = open(filename, "w") + + # Read data + for row in irf: + + # Skip header and trailer + if row.find("log(E)") != -1: + continue + if row.find("--------") != -1: + break + + # Extract data + fields = row.replace(" ",",").replace(",,",",").rstrip("\n").split(",") + logE = float(fields[0]) + aeff = float(fields[1]) + r68 = float(fields[2]) + r80 = float(fields[3]) + eres = float(fields[4]) + bgd = float(fields[5]) + sens = float(fields[6]) + + # Compute energy (in MeV) + energy = math.pow(10.0, logE)*1.0e6 + emin = math.pow(10.0, logE-0.1)*1.0e6 + emax = math.pow(10.0, logE+0.1)*1.0e6 + ewidth = emax-emin + + # Compute solid angle of r80 region in sr + omega = 2.0 * math.pi * (1.0 - math.cos(math.radians(r80))) + + # Compute background rate per steradian and MeV + bkg_rate = bgd/omega/ewidth - # Write - file.write(str(energy)+" "+str(bkg_rate)+"\n") - - # Close files - irf.close() - file.close() - - # Return - return + # Write + file.write(str(energy)+" "+str(bkg_rate)+"\n") + + # Close files + irf.close() + file.close() + + # Return + return # ================ # # Main entry point # # ================ # if __name__ == '__main__': - """ - Main entry point - """ - # Get IRF from command line argument - if (len(sys.argv) > 1): - irf = sys.argv[1] - else: - irf = "kb_E_50h_v3" + """ + Main entry point + """ + # Get IRF from command line argument + if (len(sys.argv) > 1): + irf = sys.argv[1] + else: + irf = "kb_E_50h_v3" - # Set parameters - caldb = "/usr/local/gamma/share/caldb/cta" - - # Set filenames - irfname = caldb+"/"+irf+".dat" - filename = "bkg_"+irf+".txt" - - # Read - make_file_function(irfname, filename) + # Set parameters + caldb = "/usr/local/gamma/share/caldb/cta" + + # Set filenames + irfname = caldb+"/"+irf+".dat" + filename = "bkg_"+irf+".txt" + + # Read + make_file_function(irfname, filename) diff --git a/scripts/cta_make_nodes.py b/scripts/cta_make_nodes.py index 6100304..9555a15 100755 --- a/scripts/cta_make_nodes.py +++ b/scripts/cta_make_nodes.py @@ -4,7 +4,7 @@ # logarithmically spaced energy bins. The intensity scale is set to the # HESS Crab intensity (assuming a power law). # -# Copyright (C) 2012 Juergen Knoedlseder +# Copyright (C) 2012-2014 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 @@ -29,91 +29,91 @@ import math # Write one node # # ============== # def write_node(file, energy, scale): - """ - Writes one node to XML file. - """ - # Convert to strings - str_energy = str(energy) - str_scale = str(scale) + """ + Writes one node to XML file. + """ + # Convert to strings + str_energy = str(energy) + str_scale = str(scale) - # Write start tag - file.write(' \n') + # Write start tag + file.write(' \n') - # Write energy - file.write(' \n') + # Write energy + file.write(' \n') - # Write intensity - file.write(' \n') + # Write intensity + file.write(' \n') - # Write end tag - file.write(' \n') + # Write end tag + file.write(' \n') - # Return - return - + # Return + return + # ============ # # Create nodes # # ============ # def create_nodes(emin, emax, enumbins): - """ - Create nodes (energies in TeV). - """ - # Open file - file = open("nodes.xml", "w") + """ + Create nodes (energies in TeV). + """ + # Open file + file = open("nodes.xml", "w") - # Set node boundaries - elogmin = math.log10(float(emin)) - elogmax = math.log10(float(emax)) - elogbin = (elogmax - elogmin)/(float(enumbins)-1.0) - - # Fill arrays - for i in range(int(enumbins)): - - # Compute energy - energy = math.pow(10.0, i*elogbin+elogmin) - - # Compute scale (HESS Crab spectrum) - scale = 3.45e-17 * math.pow(energy, -2.63) - - # Write node - write_node(file, energy, scale) - - # Debug - #sys.stdout.write(energy) + # Set node boundaries + elogmin = math.log10(float(emin)) + elogmax = math.log10(float(emax)) + elogbin = (elogmax - elogmin)/(float(enumbins)-1.0) + + # Fill arrays + for i in range(int(enumbins)): + + # Compute energy + energy = math.pow(10.0, i*elogbin+elogmin) + + # Compute scale (HESS Crab spectrum) + scale = 3.45e-17 * math.pow(energy, -2.63) + + # Write node + write_node(file, energy, scale) + + # Debug + #sys.stdout.write(energy) #sys.stdout.write(math.pow(10.0, int(math.log10(scale))-1.0)±"\n") - - # Close file - file.close() - - # Return - return - + + # Close file + file.close() + + # Return + return + # ======================= # # Main script entry point # # ======================= # if __name__ == '__main__': - """ - Compute nodes. - """ - # Check command line - usage = "Usage: cta_make_nodes emin emax enumbins" - if len(sys.argv) < 3: - sys.stdout.write(usage+"\n") - sys.exit() - - # Extract parameters - emin = sys.argv[1] - emax = sys.argv[2] - enumbins = sys.argv[3] + """ + Compute nodes. + """ + # Check command line + usage = "Usage: cta_make_nodes emin emax enumbins" + if len(sys.argv) < 3: + sys.stdout.write(usage+"\n") + sys.exit() + + # Extract parameters + emin = sys.argv[1] + emax = sys.argv[2] + enumbins = sys.argv[3] - # Create nodes - create_nodes(emin, emax, enumbins) - + # Create nodes + create_nodes(emin, emax, enumbins) + diff --git a/scripts/cta_root2irf.py b/scripts/cta_root2irf.py index 5afebab..c12f7b4 100755 --- a/scripts/cta_root2irf.py +++ b/scripts/cta_root2irf.py @@ -11,7 +11,7 @@ # area (i.e. the area before applying any theta cut). This is required by # GammaLib as the maximum likelihood analysis does not imply any theta cut. # -# Copyright (C) 2011 Jurgen Knodlseder +# Copyright (C) 2011-2014 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