From 89cecd4a8dc4bbe61de2d555f061ade20cec803b Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Fri, 31 Jan 2014 01:01:36 +0100 Subject: [PATCH] Adapt to new GCaldb interface. --- scripts/cspull.py | 4 ++-- scripts/obsutils.py | 551 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- src/ctlike/ctlike.cpp | 16 ++++++++++++++-- src/ctlike/ctlike.par | 2 +- src/ctmodel/ctmodel.cpp | 16 +++++++++++++++- src/ctmodel/ctmodel.par | 2 +- src/ctobssim/ctobssim.cpp | 14 +++++++++++++- 7 files changed, 325 insertions(+), 280 deletions(-) diff --git a/scripts/cspull.py b/scripts/cspull.py index cf61755..cea3a80 100755 --- a/scripts/cspull.py +++ b/scripts/cspull.py @@ -2,7 +2,7 @@ # ========================================================================== # This script generates the pull distribution for all model parameters. # -# Copyright (C) 2011-2013 Juergen Knoedlseder +# 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 @@ -97,7 +97,7 @@ class cspull(GApplication): pars.append(GApplicationPar("srcmdl","f","a","$GAMMALIB/share/models/crab.xml","","","Source model")) pars.append(GApplicationPar("outfile","f","a","pull.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("caldb","s","h","","","","Calibration database")) pars.append(GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function")) pars.append(GApplicationPar("ra","r","a","83.6331","0","360","RA of pointing (deg)")) pars.append(GApplicationPar("dec","r","a","22.0145","-90","90","Dec of pointing (deg)")) diff --git a/scripts/obsutils.py b/scripts/obsutils.py index 37b19a2..d9cdca0 100755 --- a/scripts/obsutils.py +++ b/scripts/obsutils.py @@ -2,7 +2,7 @@ # This script provides a number of functions that are useful for handling # CTA observations. # -# Copyright (C) 2011-2013 Juergen Knoedlseder +# 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 @@ -18,136 +18,136 @@ # along with this program. If not, see . # # ========================================================================== -from ctools import * -from gammalib import * +import ctools +import gammalib # ===================== # # Simulate observations # # ===================== # def sim(obs, log=False, debug=False, seed=0, nbins=0, binsz=0.05, npix=200): - """ - Simulate events for all observations in the container. - - Parameters: - obs - Observation container - Keywords: - log - Create log file(s) - debug - Create screen dump - seed - Seed value for simulations (default: 0) - nbins - Number of energy bins (default: 0=unbinned) - binsz - Pixel size for binned simulation (deg/pixel) - npix - Number of pixels in X and Y for binned simulation - """ - # Allocate ctobssim application and set parameters - sim = ctobssim(obs) - sim['seed'].integer(seed) - - # Optionally open the log file - if log: - sim.logFileOpen() + """ + Simulate events for all observations in the container. + + Parameters: + obs - Observation container + Keywords: + log - Create log file(s) + debug - Create screen dump + seed - Seed value for simulations (default: 0) + nbins - Number of energy bins (default: 0=unbinned) + binsz - Pixel size for binned simulation (deg/pixel) + npix - Number of pixels in X and Y for binned simulation + """ + # Allocate ctobssim application and set parameters + sim = ctools.ctobssim(obs) + sim['seed'].integer(seed) + + # Optionally open the log file + if log: + sim.logFileOpen() - # Optionally switch-on debugging model - if debug: - sim["debug"].boolean(True) - - # Run ctobssim application. This will loop over all observations in the - # container and simulation the events for each observation. Note that - # events are not added together, they still apply to each observation - # separately. - sim.run() - - # Binned option? - if nbins > 0: - - # Determine common energy boundaries for observations - emin = None - emax = None - for run in sim.obs(): - run_emin = run.events().ebounds().emin().TeV() - run_emax = run.events().ebounds().emax().TeV() - if emin == None: - emin = run_emin - elif run_emin > emin: - emin = run_emin - if emax == None: - emax = run_emax - elif run_emax > emax: - emax = run_emax - - # Allocate ctbin application and set parameters - bin = ctbin(sim.obs()) - bin["ebinalg"].string("LOG") - bin["emin"].real(emin) - bin["emax"].real(emax) - bin["enumbins"].integer(nbins) - bin["usepnt"].boolean(True) # Use pointing for map centre - bin["nxpix"].integer(npix) - bin["nypix"].integer(npix) - bin["binsz"].real(binsz) - bin["coordsys"].string("GAL") - bin["proj"].string("TAN") - - # Optionally open the log file - if log: - bin.logFileOpen() + # Optionally switch-on debugging model + if debug: + sim["debug"].boolean(True) + + # Run ctobssim application. This will loop over all observations in the + # container and simulation the events for each observation. Note that + # events are not added together, they still apply to each observation + # separately. + sim.run() + + # Binned option? + if nbins > 0: + + # Determine common energy boundaries for observations + emin = None + emax = None + for run in sim.obs(): + run_emin = run.events().ebounds().emin().TeV() + run_emax = run.events().ebounds().emax().TeV() + if emin == None: + emin = run_emin + elif run_emin > emin: + emin = run_emin + if emax == None: + emax = run_emax + elif run_emax > emax: + emax = run_emax + + # Allocate ctbin application and set parameters + bin = ctools.ctbin(sim.obs()) + bin["ebinalg"].string("LOG") + bin["emin"].real(emin) + bin["emax"].real(emax) + bin["enumbins"].integer(nbins) + bin["usepnt"].boolean(True) # Use pointing for map centre + bin["nxpix"].integer(npix) + bin["nypix"].integer(npix) + bin["binsz"].real(binsz) + bin["coordsys"].string("GAL") + bin["proj"].string("TAN") + + # Optionally open the log file + if log: + bin.logFileOpen() - # Optionally switch-on debugging model - if debug: - bin["debug"].boolean(True) + # Optionally switch-on debugging model + if debug: + bin["debug"].boolean(True) - # Run ctbin application. This will loop over all observations in - # the container and bin the events in counts maps - bin.run() - - # Make a deep copy of the observation that will be returned - # (the ctbin object will go out of scope one the function is - # left) - obs = bin.obs().copy() - - else: - - # Make a deep copy of the observation that will be returned - # (the ctobssim object will go out of scope one the function is - # left) - obs = sim.obs().copy() + # Run ctbin application. This will loop over all observations in + # the container and bin the events in counts maps + bin.run() + + # Make a deep copy of the observation that will be returned + # (the ctbin object will go out of scope one the function is + # left) + obs = bin.obs().copy() + + else: + + # Make a deep copy of the observation that will be returned + # (the ctobssim object will go out of scope one the function is + # left) + obs = sim.obs().copy() - # Delete the simulation - del sim - - # Return observation container - return obs + # Delete the simulation + del sim + + # Return observation container + return obs # ================ # # Fit observations # # ================ # def fit(obs, log=False, debug=False): - """ - Perform maximum likelihood fitting of observations in the container. - - Parameters: - obs - Observation container - Keywords: - log - Create log file(s) - debug - Create screen dump - """ - # Allocate ctlike application - like = ctlike(obs) - - # Optionally open the log file - if log: - like.logFileOpen() - - # Optionally switch-on debugging model - if debug: - like["debug"].boolean(True) + """ + Perform maximum likelihood fitting of observations in the container. + + Parameters: + obs - Observation container + Keywords: + log - Create log file(s) + debug - Create screen dump + """ + # Allocate ctlike application + like = ctools.ctlike(obs) + + # Optionally open the log file + if log: + like.logFileOpen() + + # Optionally switch-on debugging model + if debug: + like["debug"].boolean(True) - # Run ctlike application. - like.run() - - # Return observations - return like + # Run ctlike application. + like.run() + + # Return observations + return like # ================= # @@ -155,45 +155,45 @@ def fit(obs, log=False, debug=False): # ================= # def cntmap(obs, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \ binsz=0.05, nxpix=200, nypix=200, \ - outname="cntmap.fits"): - """ - Creates a counts map by combining the events of all observations. - The counts map will be a summed map over all energies. - - Parameters: - obs - Observation container - Keywords: - proj - Projection type (e.g. TAN, CAR, STG, ...) (default: TAN) - coord - Coordinate type (GAL, CEL) (default: GAL) - xval - Reference longitude value [deg] (default: 0.0) - yval - Reference latitude value [deg] (default: 0.0) - binsz - Pixel size [deg/pixel] (default: 0.05) - nxpix - Number of pixels in X direction (default: 200) - nypix - Number of pixels in Y direction (default: 200) - outname - Counts map FITS filename (default: cntmap.fits) - """ - # Allocate counts map - map = GSkymap(proj, coord, xval, yval, -binsz, binsz, nxpix, nypix, 1) - - # Fill all observations - for run in obs: - - # Loop over all events - for event in run.events(): - - # Determine sky pixel - skydir = GCTAInstDir(event.dir()).dir() - pixel = map.dir2inx(skydir) - - # Set pixel - map[pixel] += 1.0 - - # Save sky map. The clobber flag is set to True, so any existing FITS - # file will be overwritten. - map.save(outname, True) - - # Return counts map - return map + outname="cntmap.fits"): + """ + Creates a counts map by combining the events of all observations. + The counts map will be a summed map over all energies. + + Parameters: + obs - Observation container + Keywords: + proj - Projection type (e.g. TAN, CAR, STG, ...) (default: TAN) + coord - Coordinate type (GAL, CEL) (default: GAL) + xval - Reference longitude value [deg] (default: 0.0) + yval - Reference latitude value [deg] (default: 0.0) + binsz - Pixel size [deg/pixel] (default: 0.05) + nxpix - Number of pixels in X direction (default: 200) + nypix - Number of pixels in Y direction (default: 200) + outname - Counts map FITS filename (default: cntmap.fits) + """ + # Allocate counts map + map = gammalib.GSkymap(proj, coord, xval, yval, -binsz, binsz, nxpix, nypix, 1) + + # Fill all observations + for run in obs: + + # Loop over all events + for event in run.events(): + + # Determine sky pixel + skydir = gammalib.GCTAInstDir(event.dir()).dir() + pixel = map.dir2inx(skydir) + + # Set pixel + map[pixel] += 1.0 + + # Save sky map. The clobber flag is set to True, so any existing FITS + # file will be overwritten. + map.save(outname, True) + + # Return counts map + return map # ================ # @@ -202,63 +202,63 @@ def cntmap(obs, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \ def modmap(obs, eref=0.1, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \ binsz=0.05, nxpix=200, nypix=200, \ outname="modmap.fits"): - """ - Make model map for a given reference energy by combining all observations. - The model map will be evaluated for a given reference energy 'eref' and will - be given in units of [counts/(sr MeV s)]. - - Parameters: - obs - Observation container - Keywords: - eref - Reference energy for which model is created [TeV] (default: 0.1) - proj - Projection type (e.g. TAN, CAR, STG, ...) (default: TAN) - coord - Coordinate type (GAL, CEL) (default: GAL) - xval - Reference longitude value [deg] (default: 0.0) - yval - Reference latitude value [deg] (default: 0.0) - binsz - Pixel size [deg/pixel] (default: 0.05) - nxpix - Number of pixels in X direction (default: 200) - nypix - Number of pixels in Y direction (default: 200) - outname - Model map FITS filename (default: modmap.fits) - """ - # Allocate model map - map = GSkymap(proj, coord, xval, yval, -binsz, binsz, nxpix, nypix, 1) - - # Set reference energy, time and direction. The time is not initialised and is - # in fact not used (as the IRF is assumed to be time independent for now). - # The sky direction is set later using the pixel values. - energy = GEnergy() - time = GTime() - instdir = GCTAInstDir() - energy.TeV(eref) - - # Loop over all map pixels - for pixel in range(map.npix()): - - # Get sky direction - skydir = map.inx2dir(pixel) - instdir.dir(skydir) - - # Create event atom for map pixel - atom = GCTAEventAtom() - atom.dir(instdir) - atom.energy(energy) - atom.time(time) - - # Initialise model value - value = 0.0 - - # Loop over all observations - for run in obs: - value += obs.models().eval(atom, run) - - # Set map value - map[pixel] = value - - # Save sky map - map.save(outname, True) - - # Return model map - return map + """ + Make model map for a given reference energy by combining all observations. + The model map will be evaluated for a given reference energy 'eref' and will + be given in units of [counts/(sr MeV s)]. + + Parameters: + obs - Observation container + Keywords: + eref - Reference energy for which model is created [TeV] (default: 0.1) + proj - Projection type (e.g. TAN, CAR, STG, ...) (default: TAN) + coord - Coordinate type (GAL, CEL) (default: GAL) + xval - Reference longitude value [deg] (default: 0.0) + yval - Reference latitude value [deg] (default: 0.0) + binsz - Pixel size [deg/pixel] (default: 0.05) + nxpix - Number of pixels in X direction (default: 200) + nypix - Number of pixels in Y direction (default: 200) + outname - Model map FITS filename (default: modmap.fits) + """ + # Allocate model map + map = gammalib.GSkymap(proj, coord, xval, yval, -binsz, binsz, nxpix, nypix, 1) + + # Set reference energy, time and direction. The time is not initialised and is + # in fact not used (as the IRF is assumed to be time independent for now). + # The sky direction is set later using the pixel values. + energy = gammalib.GEnergy() + time = gammalib.GTime() + instdir = gammalib.GCTAInstDir() + energy.TeV(eref) + + # Loop over all map pixels + for pixel in range(map.npix()): + + # Get sky direction + skydir = map.inx2dir(pixel) + instdir.dir(skydir) + + # Create event atom for map pixel + atom = gammalib.GCTAEventAtom() + atom.dir(instdir) + atom.energy(energy) + atom.time(time) + + # Initialise model value + value = 0.0 + + # Loop over all observations + for run in obs: + value += obs.models().eval(atom, run) + + # Set map value + map[pixel] = value + + # Save sky map + map.save(outname, True) + + # Return model map + return map # ======================= # @@ -266,66 +266,73 @@ def modmap(obs, eref=0.1, proj="TAN", coord="GAL", xval=0.0, yval=0.0, \ # ======================= # def set(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \ emin=0.1, emax=100.0, rad=5.0, \ - irf="cta_dummy_irf", caldb="$GAMMALIB/share/caldb/cta"): - """ - Returns a single CTA observation. By looping over this function we can - add CTA observations to the observation container. - - Parameters: - pntdir - Pointing direction - Keywords: - tstart - Start time [seconds] (default: 0.0) - duration - Duration of observation [seconds] (default: 1800.0) - deadc - Deadtime correction factor (default: 0.95) - emin - Minimum event energy [TeV] (default: 0.1) - emax - Maximum event energy [TeV] (default: 100.0) - rad - ROI radius used for analysis [deg] (default: 5.0) - irf - Instrument response function (default: cta_dummy_irf) - caldb - Calibration database path (default: $GAMMALIB/share/caldb/cta) - """ - # Allocate CTA observation - obs = GCTAObservation() - - # Set pointing direction - pnt = GCTAPointing() - pnt.dir(pntdir) - obs.pointing(pnt) - - # Set ROI - roi = GCTARoi() - instdir = GCTAInstDir() - instdir.dir(pntdir) - roi.centre(instdir) - roi.radius(rad) - - # Set GTI - gti = GGti() - start = GTime(tstart) - stop = GTime(tstart+duration) - gti.append(start, stop) - - # Set energy boundaries - ebounds = GEbounds() - e_min = GEnergy() - e_max = GEnergy() - e_min.TeV(emin) - e_max.TeV(emax) - ebounds.append(e_min, e_max) + irf="cta_dummy_irf", caldb=""): + """ + Returns a single CTA observation. By looping over this function we can + add CTA observations to the observation container. + + Parameters: + pntdir - Pointing direction + Keywords: + tstart - Start time [seconds] (default: 0.0) + duration - Duration of observation [seconds] (default: 1800.0) + deadc - Deadtime correction factor (default: 0.95) + emin - Minimum event energy [TeV] (default: 0.1) + emax - Maximum event energy [TeV] (default: 100.0) + rad - ROI radius used for analysis [deg] (default: 5.0) + irf - Instrument response function (default: cta_dummy_irf) + caldb - Calibration database path (default: "") + """ + # Allocate CTA observation + obs = gammalib.GCTAObservation() - # Allocate event list - events = GCTAEventList() - events.roi(roi) - events.gti(gti) - events.ebounds(ebounds) - obs.events(events) - - # Set instrument response - obs.response(irf, caldb) - - # Set ontime, livetime, and deadtime correction factor - obs.ontime(duration) - obs.livetime(duration*deadc) - obs.deadc(deadc) - - # Return observation - return obs + # Set calibration database + db = gammalib.GCaldb() + if (gammalib.dir_exists(caldb)): + db.dir(caldb) + else: + db.open("cta", caldb) + + # Set pointing direction + pnt = gammalib.GCTAPointing() + pnt.dir(pntdir) + obs.pointing(pnt) + + # Set ROI + roi = gammalib.GCTARoi() + instdir = gammalib.GCTAInstDir() + instdir.dir(pntdir) + roi.centre(instdir) + roi.radius(rad) + + # Set GTI + gti = gammalib.GGti() + start = gammalib.GTime(tstart) + stop = gammalib.GTime(tstart+duration) + gti.append(start, stop) + + # Set energy boundaries + ebounds = gammalib.GEbounds() + e_min = gammalib.GEnergy() + e_max = gammalib.GEnergy() + e_min.TeV(emin) + e_max.TeV(emax) + ebounds.append(e_min, e_max) + + # Allocate event list + events = gammalib.GCTAEventList() + events.roi(roi) + events.gti(gti) + events.ebounds(ebounds) + obs.events(events) + + # Set instrument response + obs.response(irf, db) + + # Set ontime, livetime, and deadtime correction factor + obs.ontime(duration) + obs.livetime(duration*deadc) + obs.deadc(deadc) + + # Return observation + return obs diff --git a/src/ctlike/ctlike.cpp b/src/ctlike/ctlike.cpp index 5314df0..d05699a 100644 --- a/src/ctlike/ctlike.cpp +++ b/src/ctlike/ctlike.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctlike - CTA maximum likelihood tool * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2013 by Juergen Knoedlseder * + * copyright (C) 2010-2014 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -361,8 +361,20 @@ void ctlike::get_parameters(void) // Set statistics obs.statistics(m_stat); + // Set calibration database. If specified parameter is a + // directory then use this as the pathname to the calibration + // database. Otherwise interpret this as the instrument name, + // the mission being "cta" + GCaldb caldb; + if (gammalib::dir_exists(m_caldb)) { + caldb.dir(m_caldb); + } + else { + caldb.open("cta", m_caldb); + } + // Set reponse - obs.response(m_irf, m_caldb); + obs.response(m_irf, caldb); // Append observation to container m_obs.append(obs); diff --git a/src/ctlike/ctlike.par b/src/ctlike/ctlike.par index 8c83dfb..bd8b77d 100644 --- a/src/ctlike/ctlike.par +++ b/src/ctlike/ctlike.par @@ -20,7 +20,7 @@ infile, f, a, "test/data/crab_events.fits.gz",,,"Event list, counts map or observation definition file" stat, s, h, "POISSON",,, "Optimization statistics" refit, b, h, no,,, "Do refitting?" -caldb, s, a, "test/irf",,, "Calibration database" +caldb, s, h, "",,, "Calibration database" irf, s, a, "cta_dummy_irf",,, "Instrument response function" srcmdl, f, a, "test/data/crab.xml",,, "Source model" outmdl, f, a, "crab_results.xml",,, "Source model output file" diff --git a/src/ctmodel/ctmodel.cpp b/src/ctmodel/ctmodel.cpp index 89f19f2..c02e0d8 100644 --- a/src/ctmodel/ctmodel.cpp +++ b/src/ctmodel/ctmodel.cpp @@ -581,7 +581,21 @@ void ctmodel::setup_obs(void) // Set response if it isn't set already if (obs->response().aeff() == NULL) { - obs->response(m_irf, m_caldb); + + // Set calibration database. If specified parameter is a + // directory then use this as the pathname to the calibration + // database. Otherwise interpret this as the instrument name, + // the mission being "cta" + GCaldb caldb; + if (gammalib::dir_exists(m_caldb)) { + caldb.dir(m_caldb); + } + else { + caldb.open("cta", m_caldb); + } + + // Set reponse + obs->response(m_irf, caldb); } // endif: observation already has a response diff --git a/src/ctmodel/ctmodel.par b/src/ctmodel/ctmodel.par index a87e677..a48ba56 100644 --- a/src/ctmodel/ctmodel.par +++ b/src/ctmodel/ctmodel.par @@ -20,7 +20,7 @@ infile, f, a, "NONE",,, "Input counts map or observation definition file" outfile, f, a, "modmap.fits",,, "Output counts map or observation definition file" prefix, s, h, "modmap_",,, "Prefix for counts maps in observation definition file" -caldb, s, a, "test/irf",,, "Calibration database" +caldb, s, h, "",,, "Calibration database" irf, s, a, "cta_dummy_irf",,, "Instrument response function" srcmdl, f, a, "test/data/crab.xml",,, "Source model" # diff --git a/src/ctobssim/ctobssim.cpp b/src/ctobssim/ctobssim.cpp index 6163a60..a65f1ad 100644 --- a/src/ctobssim/ctobssim.cpp +++ b/src/ctobssim/ctobssim.cpp @@ -432,9 +432,21 @@ void ctobssim::get_parameters(void) // Allocate CTA observation GCTAObservation obs; + // Set calibration database. If specified parameter is a directory + // then use this as the pathname to the calibration database. Other- + // wise interpret this as the instrument name, the mission being + // "cta" + GCaldb caldb; + if (gammalib::dir_exists(m_caldb)) { + caldb.dir(m_caldb); + } + else { + caldb.open("cta", m_caldb); + } + // Set CTA observation attributes obs.pointing(pnt); - obs.response(m_irf, m_caldb); + obs.response(m_irf, caldb); // Set event list (queries remaining parameters) set_list(&obs); -- libgit2 0.21.4