From 50d6ad8850eb912ac8f5d7405bdb7faa952b18df Mon Sep 17 00:00:00 2001 From: Michael Mayer Date: Mon, 30 Mar 2015 13:30:17 +0200 Subject: [PATCH] csspec now supports binned analysis in each energy bin --- scripts/csspec.py | 242 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------- test/test_cscripts.sh.in | 34 ++++++++++++++++++++++++++++++++++ 2 files changed, 247 insertions(+), 29 deletions(-) diff --git a/scripts/csspec.py b/scripts/csspec.py index 00e32c8..fb6b34b 100755 --- a/scripts/csspec.py +++ b/scripts/csspec.py @@ -21,6 +21,7 @@ import gammalib import ctools import sys +import tempfile # ============== # # cstsdist class # @@ -106,10 +107,19 @@ class csspec(ctools.cscript): pars.append(gammalib.GApplicationPar("psfcube","f","a","NONE","","","PSF cube file (only needed for stacked analysis)")) pars.append(gammalib.GApplicationPar("caldb","s","a","dummy","","","Calibration database")) pars.append(gammalib.GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function")) - pars.append(gammalib.GApplicationPar("emin","r","h","0.1","","","Lower energy limit (TeV)")) - pars.append(gammalib.GApplicationPar("emax","r","h","100.0","","","Upper energy limit (TeV)")) - pars.append(gammalib.GApplicationPar("enumbins","i","a","20","","","Number of energy bins")) - pars.append(gammalib.GApplicationPar("ebinalg","s","h","LOG","LIN|LOG|FILE","","Binning algorithm")) + pars.append(gammalib.GApplicationPar("emin","r","h","0.1","","","Lower energy limit for spectral points(TeV)")) + pars.append(gammalib.GApplicationPar("emax","r","h","100.0","","","Upper energy limit for spectral points(TeV)")) + pars.append(gammalib.GApplicationPar("enumbins","i","a","20","","","Number of spectral points")) + pars.append(gammalib.GApplicationPar("ebinalg","s", "h","LOG","FILE|LIN|LOG","", "Algorithm for defining energy bins")) + pars.append(gammalib.GApplicationPar("binned","b","a","no","yes|no","","Use binned analysis in each energy bin")) + pars.append(gammalib.GApplicationPar("nebins","i","h","5","","","Number of energy bins per spectral point")) + pars.append(gammalib.GApplicationPar("coordsys","s","a","CEL","CEL|GAL","","Coordinate System")) + pars.append(gammalib.GApplicationPar("proj","s","a","CAR","AIT|AZP|CAR|MER|MOL|STG|TAN","","Projection method")) + pars.append(gammalib.GApplicationPar("xref","r","a","83.63","0","360","First coordinate of image center in degrees (RA or galactic l)")) + pars.append(gammalib.GApplicationPar("yref","r","a","22.01","-90","90","Second coordinate of image center in degrees (DEC or galactic b)")) + pars.append(gammalib.GApplicationPar("nxpix","i","a","200","","","Size of the X axis in pixels")) + pars.append(gammalib.GApplicationPar("nypix","i","a","200","","","Size of the Y axis in pixels")) + pars.append(gammalib.GApplicationPar("binsz","r","a","0.02","","","Pixel size (deg/pixel)")) pars.append(gammalib.GApplicationPar("calc_ts","b","h","yes","yes|no","","Compute TS value in each bin")) pars.append(gammalib.GApplicationPar("calc_ulim","b","h","yes","yes|no","","Compute upper limit in each bin")) pars.append(gammalib.GApplicationPar("fix_srcs","b","h","yes","yes|no","","Fix other skymodel parameters")) @@ -138,9 +148,21 @@ class csspec(ctools.cscript): # Get source name self.m_srcname = self["srcname"].string() - # Get ebounds + # Get ebounds self.m_ebounds = self.create_ebounds() + # Get binning flag + self.m_binned = self["binned"].boolean() + if self.m_binned: + self.m_xref = self["xref"].real() + self.m_yref = self["yref"].real() + self.m_nxpix = self["nxpix"].integer() + self.m_nypix = self["nypix"].integer() + self.m_binsz = self["binsz"].real() + self.m_coordsys = self["coordsys"].string() + self.m_proj = self["proj"].string() + self.m_ebins = self["nebins"].integer() + # Read other parameters self.m_outfile = self["outfile"].filename() @@ -210,21 +232,39 @@ class csspec(ctools.cscript): self.log("\n") self.log.header1("Adjust model parameters") + # Adjust model parameters dependent on input user parameters for model in self.obs.models(): + + # Set TS flag for all models to false. + # Source of interest will be set to true later model.tscalc(False) + + # Log model name + if self.logExplicit(): + self.log.header3(model.name()) + + # Deal with the source of interest if model.name() == self.m_srcname: for par in model: + if par.is_free() and self.logExplicit(): + self.log(" Fixing \""+par.name()+"\"\n") par.fix() + if par.is_fixed() and self.logExplicit(): + self.log(" Freeing \""+par.name()+"\"\n") model.spectral()[0].free() if self.m_calc_ts: model.tscalc(True) - elif self.m_fix_bkg and not model.classname() == "GModelData": + elif self.m_fix_bkg and not model.classname() == "GModelSky": for par in model: + if par.is_free() and self.logExplicit(): + self.log(" Fixing \""+par.name()+"\"\n") par.fix() elif self.m_fix_srcs and model.classname() == "GModelSky": for par in model: + if par.is_free() and self.logExplicit(): + self.log(" Fixing \""+par.name()+"\"\n") par.fix() # Write header @@ -243,8 +283,10 @@ class csspec(ctools.cscript): # Create FITS table columns energy = gammalib.GFitsTableDoubleCol("Energy", self.m_ebounds.size()) energy.unit("TeV") - energy_err = gammalib.GFitsTableDoubleCol("e_Energy", self.m_ebounds.size()) - energy_err.unit("TeV") + energy_low = gammalib.GFitsTableDoubleCol("ed_Energy", self.m_ebounds.size()) + energy_low.unit("TeV") + energy_high = gammalib.GFitsTableDoubleCol("eu_Energy", self.m_ebounds.size()) + energy_high.unit("TeV") flux = gammalib.GFitsTableDoubleCol("Flux", self.m_ebounds.size()) flux.unit("erg/cm2/s") flux_err = gammalib.GFitsTableDoubleCol("e_Flux", self.m_ebounds.size()) @@ -252,16 +294,31 @@ class csspec(ctools.cscript): TSvalues = gammalib.GFitsTableDoubleCol("TS", self.m_ebounds.size()) ulim_values = gammalib.GFitsTableDoubleCol("UpperLimit", self.m_ebounds.size()) ulim_values.unit("erg/cm2/s") - + # Loop over energy bins for i in range(self.m_ebounds.size()): + # Log information + if self.logTerse(): + self.log("\n") + self.log.header2("Energy bin "+str(i)) + # Get energy boundaries emin = self.m_ebounds.emin(i) emax = self.m_ebounds.emax(i) - emean = self.m_ebounds.emean(i) elogmean = self.m_ebounds.elogmean(i) - elogmean2 = elogmean.MeV() * elogmean.MeV() + elogmean2 = elogmean.MeV() * elogmean.MeV() + + # Store energy as TeV + energy[i] = elogmean.TeV() + + # Store energy errors + energy_low[i] = (elogmean - emin).TeV() + energy_high[i] = (emax - elogmean).TeV() + + # Log information + if self.logExplicit(): + self.log.header3("Selecting events") # Select events select = ctools.ctselect(self.obs) @@ -272,23 +329,156 @@ class csspec(ctools.cscript): select["rad"].value("UNDEFINED") select["ra"].value("UNDEFINED") select["dec"].value("UNDEFINED") - select.run() + select.run() + obs = select.obs() + + if self.m_binned: + + if self.logTerse(): + self.log.header3("Binning events") + + # bin events + bin = ctools.ctbin(select.obs()) + bin["usepnt"].boolean(False) + bin["ebinalg"].string("LOG") + bin["xref"].real(self.m_xref) + bin["yref"].real(self.m_yref) + bin["binsz"].real(self.m_binsz) + bin["nxpix"].integer(self.m_nxpix) + bin["nypix"].integer(self.m_nypix) + bin["enumbins"].integer(self.m_ebins) + bin["emin"].real(emin.TeV()) + bin["emax"].real(emax.TeV()) + bin["coordsys"].string(self.m_coordsys) + bin["proj"].string(self.m_proj) + bin.run() + + if self.logTerse(): + self.log.header3("Creating exposure cube") + + # create exposure cube + expcube = ctools.ctexpcube(select.obs()) + expcube["incube"].filename("NONE") + expcube["usepnt"].boolean(False) + expcube["ebinalg"].string("LOG") + expcube["xref"].real(self.m_xref) + expcube["yref"].real(self.m_yref) + expcube["binsz"].real(self.m_binsz) + expcube["nxpix"].integer(self.m_nxpix) + expcube["nypix"].integer(self.m_nypix) + expcube["enumbins"].integer(self.m_ebins) + expcube["emin"].real(emin.TeV()) + expcube["emax"].real(emax.TeV()) + expcube["coordsys"].string(self.m_coordsys) + expcube["proj"].string(self.m_proj) + expcube.run() + + if self.logTerse(): + self.log.header3("Creating PSF cube") + + # create psf cube + psfcube = ctools.ctpsfcube(select.obs()) + psfcube["incube"].filename("NONE") + psfcube["usepnt"].boolean(False) + psfcube["ebinalg"].string("LOG") + psfcube["xref"].real(self.m_xref) + psfcube["yref"].real(self.m_yref) + psfcube["binsz"].real(self.m_binsz) + psfcube["nxpix"].integer(self.m_nxpix) + psfcube["nypix"].integer(self.m_nypix) + psfcube["enumbins"].integer(self.m_ebins) + psfcube["emin"].real(emin.TeV()) + psfcube["emax"].real(emax.TeV()) + psfcube["coordsys"].string(self.m_coordsys) + psfcube["proj"].string(self.m_proj) + psfcube.run() + + if self.logTerse(): + self.log.header3("Creating background cube") + + # create background cube + bkgcube = ctools.ctbkgcube(select.obs()) + bkgcube["incube"].filename("NONE") + bkgcube["usepnt"].boolean(False) + bkgcube["ebinalg"].string("LOG") + bkgcube["xref"].real(self.m_xref) + bkgcube["yref"].real(self.m_yref) + bkgcube["binsz"].real(self.m_binsz) + bkgcube["nxpix"].integer(self.m_nxpix) + bkgcube["nypix"].integer(self.m_nypix) + bkgcube["enumbins"].integer(self.m_ebins) + bkgcube["emin"].real(emin.TeV()) + bkgcube["emax"].real(emax.TeV()) + bkgcube["coordsys"].string(self.m_coordsys) + bkgcube["proj"].string(self.m_proj) + bkgcube.run() + + # Set new binned observation + obs = bin.obs() + + # Set precomputed binned response + obs[0].response(expcube.expcube(), psfcube.psfcube(), bkgcube.cube()) + + # Get new models + models = bkgcube.models() + + # Fix background models if required + if self.m_fix_bkg: + for model in models: + if not model.classname() == "GModelSky": + for par in model: + par.fix() + + # Set new models to binned observation + obs.models(models) + + if self.logTerse(): + self.log.header3("Performing fit") + # likelihood - like = ctools.ctlike(select.obs()) + like = ctools.ctlike(obs) like.run() + # Skip bin if no event was present + if like.obs().logL() == 0.0: + + # Log information + if self.logTerse(): + self.log("No event in this bin. Bin is skipped\n") + + flux[i] = 0.0 + flux_err[i] = 0.0 + TSvalues[i] = 0.0 + ulim_values[i] = 0.0 + continue + # Get results fitted_models = like.obs().models() source = fitted_models[self.m_srcname] - + # Calculate Upper Limit ulimit_value = -1.0 - #if self.m_calc_ulimit: - # ulimit = ctools.ctulimit(like.obs()) - # ulimit["eref"].real(elogmean.TeV()) - # ulimit.run() - # ulimit_value.ulimit.diff_ulimit() + if self.m_calc_ulimit: + + # Logging information + if self.logTerse(): + self.log.header3("Computing upper limit") + + # Create upper limit object + ulimit = ctools.ctulimit(like.obs()) + ulimit["srcname"].string(self.m_srcname) + ulimit["eref"].real(elogmean.TeV()) + + # Try to run upper limit and catch exceptions + try: + ulimit.run() + ulimit_value = ulimit.diff_ulimit() + except: + + if self.logTerse(): + self.log("Upper limit calculation failed\n") + ulimit_value = -1.0 # Get TS value TS = -1.0 @@ -306,19 +496,12 @@ class csspec(ctools.cscript): # Set values for storage TSvalues[i] = TS - # Store energy as TeV - energy[i] = elogmean.TeV() - - # Store energy error - energy_err[i] = (emax - emean).TeV() - # Convert fluxes to nuFnu flux[i] = fitted_flux * elogmean2 * gammalib.MeV2erg flux_err[i] = e_flux * elogmean2 * gammalib.MeV2erg if ulimit_value > 0.0: ulim_values[i] = ulimit_value * elogmean2 * gammalib.MeV2erg - - + # Log information if self.logExplicit(): self.log("Bin "+str(i)+" ["+str(emin.TeV())+"-"+str(emax.TeV())+"] TeV: ") @@ -329,10 +512,11 @@ class csspec(ctools.cscript): if self.m_calc_ulimit and ulim_values[i] > 0.0: self.log(", UL = "+str(ulim_values[i])+" [erg/cm2/s]") self.log("\n") - + # Append filles columns to fits table table.append(energy) - table.append(energy_err) + table.append(energy_low) + table.append(energy_high) table.append(flux) table.append(flux_err) table.append(TSvalues) diff --git a/test/test_cscripts.sh.in b/test/test_cscripts.sh.in index eb326b1..e66f5cd 100755 --- a/test/test_cscripts.sh.in +++ b/test/test_cscripts.sh.in @@ -203,6 +203,7 @@ chmod +x @abs_top_srcdir@/scripts/csspec.py srcname="Crab" \ emin=0.1 \ emax=100.0 \ + binned=no \ outfile="spectrum.fits" > csspec.log $ECHO -n "." if [ -s "spectrum.fits" ] @@ -214,4 +215,37 @@ else fi $ECHO " ok" +# +# Test csspec Run2 +# ============= +$ECHO -n "Test csspec: " +chmod +x @abs_top_srcdir@/scripts/csspec.py +@abs_top_srcdir@/scripts/csspec.py \ + inobs="data/crab_events.fits" \ + inmodel="data/crab.xml" \ + caldb="irf" \ + irf="cta_dummy_irf" \ + enumbins=3 \ + srcname="Crab" \ + emin=0.1 \ + emax=100.0 \ + binned=yes \ + nxpix=10 \ + nypix=10 \ + binsz=0.1 \ + nebins=2 \ + coordsys="CEL" \ + xref=83.63 \ + yref=22.01 \ + proj="CAR" \ + outfile="spectrum.fits" > csspec.log +$ECHO -n "." +if [ -s "spectrum.fits" ] +then + $ECHO -n "." +else + $ECHO " spectrum.fits file is not valid" + exit 1 +fi +$ECHO " ok" -- libgit2 0.21.4