From 9afc4de28fcd628c59acf45ccf57b727dfaa1133 Mon Sep 17 00:00:00 2001 From: Michael Mayer Date: Fri, 6 Nov 2015 14:14:03 +0100 Subject: [PATCH] Modify csspec to work on binned input, too. Added test case and slighlty updated ctcubemask for this purpose --- cscripts/csspec.py | 207 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------------------------------- src/ctcubemask/ctcubemask.cpp | 81 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------- src/ctcubemask/ctcubemask.hpp | 2 ++ test/test_cscripts.sh.in | 37 ++++++++++++++----------------------- 4 files changed, 144 insertions(+), 183 deletions(-) diff --git a/cscripts/csspec.py b/cscripts/csspec.py index 6b3b06e..2e037a7 100755 --- a/cscripts/csspec.py +++ b/cscripts/csspec.py @@ -43,7 +43,7 @@ class csspec(ctools.cscript): self.version = "1.0.0" # Initialise some members - self.obs = None + self.obs = gammalib.GObservations() self.fits = None # Initialise some members @@ -106,7 +106,7 @@ class csspec(ctools.cscript): pars.append(gammalib.GApplicationPar("outfile","f","a","spectrum.fits","","","Output file name")) pars.append(gammalib.GApplicationPar("expcube","f","a","NONE","","","Exposure cube file (only needed for stacked analysis)")) pars.append(gammalib.GApplicationPar("psfcube","f","a","NONE","","","PSF cube file (only needed for stacked analysis)")) - pars.append(gammalib.GApplicationPar("bkgcube","s","a","NONE","","","Background cube file (only needed for stacked analysis)")) + pars.append(gammalib.GApplicationPar("bkgcube","f","a","NONE","","","Background cube file (only needed for stacked analysis)")) pars.append(gammalib.GApplicationPar("caldb","s","a","prod2","","","Calibration database")) pars.append(gammalib.GApplicationPar("irf","s","a","South_50h","","","Instrument response function")) pars.append(gammalib.GApplicationPar("edisp","b","h","no","","","Apply energy dispersion?")) @@ -114,16 +114,6 @@ class csspec(ctools.cscript): 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("anumbins","i","h","200","","","Angular separation bins for PSF cube in binned mode")) 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")) @@ -139,10 +129,18 @@ class csspec(ctools.cscript): """ Get parameters from parfile and setup the observation. """ + # Set observation if not done before if self.obs == None or self.obs.size() == 0: self.require_inobs("csspec::get_parameters()") self.obs = self.get_observations() + + # Check if we have one binned cta observation, i.e. if we are in binned mode + self.m_binned_mode = False + if self.obs.size() == 1: + if self.obs[0].classname() == "GCTAObservation": + if self.obs[0].eventtype() == "CountsCube": + self.m_binned_mode = True # Set models if we have none if self.obs.models().size() == 0: @@ -151,25 +149,44 @@ class csspec(ctools.cscript): # Get source name self.m_srcname = self["srcname"].string() - # Get ebounds - self.m_ebounds = self.create_ebounds() - - # Get binning flag - self.m_binned = self["binned"].boolean() - + # Get ebounds + if self.m_binned_mode: + + # Extract cube ebounds + cube_ebounds = self.obs[0].events().ebounds() + + # Read user parameters + self.m_enumbins = self["enumbins"].integer() + self.m_emin = self["emin"].real() - 1e-6 #rounding tolerance + self.m_emax = self["emax"].real() + 1e-6 #rounding tolerance + + # Compute number of energy layers used for one spectral bin + n_layers = cube_ebounds.size() / self.m_enumbins + + # Don't allow zero layers + if n_layers == 0: + n_layers = 1 + + # Create new ebounds object and fill it + self.m_ebounds = gammalib.GEbounds() + for i in range(cube_ebounds.size())[::n_layers]: + if i + n_layers > cube_ebounds.size(): + break + + # Get new emin and emax of bin + emin = cube_ebounds.emin(i) + emax = cube_ebounds.emax(i + n_layers - 1) + + # Append energy bin if within user range + if self.m_emin <= emin.TeV() and self.m_emax >= emax.TeV(): + self.m_ebounds.append(emin, emax) + + # Unbinned mode + else: + self.m_ebounds = self.create_ebounds() + # Get edisp flag self.m_edisp = self["edisp"].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() - self.m_anumbins = self["anumbins"].integer() # Read other parameters self.m_outfile = self["outfile"].filename() @@ -286,7 +303,8 @@ class csspec(ctools.cscript): # Write header if self.logTerse(): self.log("\n") - self.log.header1("Generate spectrum") + self.log.header1("Generate spectrum") + self.log(str(self.m_ebounds)) # Initialise FITS Table with extension "SPECTRUM" table = gammalib.GFitsBinTable(self.m_ebounds.size()) @@ -334,8 +352,8 @@ class csspec(ctools.cscript): energy_low[i] = (elogmean - emin).TeV() energy_high[i] = (emax - elogmean).TeV() - # unbinned analysis - if not self.m_binned: + # use ctselect for unbinned analysis + if not self.m_binned_mode: # Log information if self.logExplicit(): @@ -355,121 +373,32 @@ class csspec(ctools.cscript): # Retrieve observation obs = select.obs() - # Binned analysis + # use ctcubemask for binned analysis else: # Header if self.logTerse(): - self.log.header3("Binning events") - - # Bin events - bin = ctools.ctbin(self.obs) - bin["usepnt"] = False - bin["ebinalg"] = "LOG" - bin["xref"] = self.m_xref - bin["yref"] = self.m_yref - bin["binsz"] = self.m_binsz - bin["nxpix"] = self.m_nxpix - bin["nypix"] = self.m_nypix - bin["enumbins"] = self.m_ebins - bin["emin"] = emin.TeV() - bin["emax"] = emax.TeV() - bin["coordsys"] = self.m_coordsys - bin["proj"] = self.m_proj - bin.run() - - # Header - if self.logTerse(): - self.log.header3("Creating exposure cube") - - # Create exposure cube - expcube = ctools.ctexpcube(self.obs) - expcube["incube"] = "NONE" - expcube["usepnt"] = False - expcube["ebinalg"] = "LOG" - expcube["edisp"] = self.m_edisp - expcube["xref"] = self.m_xref - expcube["yref"] = self.m_yref - expcube["binsz"] = self.m_binsz - expcube["nxpix"] = self.m_nxpix - expcube["nypix"] = self.m_nypix - expcube["enumbins"] = self.m_ebins - expcube["emin"] = emin.TeV() - expcube["emax"] = emax.TeV() - expcube["coordsys"] = self.m_coordsys - expcube["proj"] = self.m_proj - expcube.run() - - # Header - if self.logTerse(): - self.log.header3("Creating PSF cube") - - # Create psf cube - psfcube = ctools.ctpsfcube(self.obs) - psfcube["incube"] = "NONE" - psfcube["usepnt"] = False - psfcube["ebinalg"] = "LOG" - psfcube["edisp"] = self.m_edisp - psfcube["xref"] = self.m_xref - psfcube["yref"] = self.m_yref - psfcube["binsz"] = self.m_binsz*20.0 # Slow Psf variation - psfcube["nxpix"] = self.m_nxpix/20+1 - psfcube["nypix"] = self.m_nypix/20+1 - psfcube["enumbins"] = self.m_ebins - psfcube["anumbins"] = self.m_anumbins - psfcube["emin"] = emin.TeV() - psfcube["emax"] = emax.TeV() - psfcube["coordsys"] = self.m_coordsys - psfcube["proj"] = self.m_proj - psfcube.run() - - # Header - if self.logTerse(): - self.log.header3("Creating background cube") - - # Create background cube - bkgcube = ctools.ctbkgcube(self.obs) - bkgcube["incube"] = "NONE" - bkgcube["usepnt"] = False - bkgcube["ebinalg"] = "LOG" - bkgcube["xref"] = self.m_xref - bkgcube["yref"] = self.m_yref - bkgcube["binsz"] = self.m_binsz - bkgcube["nxpix"] = self.m_nxpix - bkgcube["nypix"] = self.m_nypix - bkgcube["enumbins"] = self.m_ebins - bkgcube["emin"] = emin.TeV() - bkgcube["emax"] = emax.TeV() - bkgcube["coordsys"] = self.m_coordsys - bkgcube["proj"] = self.m_proj - bkgcube.run() + self.log.header3("Selecting events") + cubemask = ctools.ctcubemask(self.obs) + cubemask["regfile"] = "NONE" + cubemask["ra"] = "UNDEFINED" + cubemask["dec"] = "UNDEFINED" + cubemask["rad"] = "UNDEFINED" + cubemask["emin"] = emin.TeV() + cubemask["emax"] = emax.TeV() + cubemask.run() + # Set new binned observation - obs = bin.obs() - - # Set precomputed binned response - obs[0].response(expcube.expcube(), psfcube.psfcube(), bkgcube.bkgcube()) - - # 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) + obs = cubemask.obs() # Header if self.logTerse(): self.log.header3("Performing fit") # Likelihood - like = ctools.ctlike(obs) - like["edisp"].boolean(self.m_edisp) + like = ctools.ctlike(obs) + like["edisp"] = self.m_edisp like.run() # Skip bin if no event was present @@ -518,9 +447,9 @@ class csspec(ctools.cscript): if self.m_calc_ts: TS = source.ts() - # Compute Npred value + # Compute Npred value (only works for unbinned analysis) Npred = 0.0 - if not self.m_binned: + if not self.m_binned_mode: for observation in like.obs(): Npred += observation.npred(source) @@ -529,8 +458,8 @@ class csspec(ctools.cscript): # Compute flux error parvalue = source.spectral()[0].value() - rel_error = source.spectral()[0].error()/parvalue - e_flux = fitted_flux*rel_error + rel_error = source.spectral()[0].error() / parvalue + e_flux = fitted_flux * rel_error # Set values for storage TSvalues[i] = TS diff --git a/src/ctcubemask/ctcubemask.cpp b/src/ctcubemask/ctcubemask.cpp index 50ac4a0..3cd4c13 100644 --- a/src/ctcubemask/ctcubemask.cpp +++ b/src/ctcubemask/ctcubemask.cpp @@ -349,6 +349,10 @@ void ctcubemask::save(void) ***************************************************************************/ void ctcubemask::get_parameters(void) { + // Initialise selection flags + m_select_energy = true; + m_select_roi = true; + // If there are no observations in container then load them via user // parameters if (m_obs.size() == 0) { @@ -373,13 +377,36 @@ void ctcubemask::get_parameters(void) m_regfile = (*this)["regfile"].filename(); m_usepnt = (*this)["usepnt"].boolean(); if (!m_usepnt) { - m_ra = (*this)["ra"].real(); - m_dec = (*this)["dec"].real(); + + // Check RA/DEC parameters for validity to read + if ((*this)["ra"].is_valid() && (*this)["dec"].is_valid()) { + m_ra = (*this)["ra"].real(); + m_dec = (*this)["dec"].real(); + m_select_roi = true; + } + else { + m_select_roi = false; + } + } + + // Check if radius is vaild for a RoI selection + if (m_select_roi && (*this)["rad"].is_valid()) { + m_rad = (*this)["rad"].real(); } - m_rad = (*this)["rad"].real(); - m_emin = (*this)["emin"].real(); - m_emax = (*this)["emax"].real(); - + else { + m_select_roi = false; + } + + // Check for sanity of energy selection parameters + if ((*this)["emin"].is_valid() && (*this)["emax"].is_valid()) { + m_emin = (*this)["emin"].real(); + m_emax = (*this)["emax"].real(); + m_select_energy = true; + } + else { + m_select_energy = false; + } + // Optionally read ahead parameters so that they get correctly // dumped into the log file if (read_ahead()) { @@ -410,6 +437,12 @@ void ctcubemask::apply_mask(GCTAObservation* obs) GSkyMap map = cube->map(); const GEbounds& ebounds = cube->ebounds(); + // If no energy selection is required set energy boundaries to cube boundaries + if (!m_select_energy) { + m_emin = ebounds.emin().TeV(); + m_emax = ebounds.emax().TeV(); + } + // Initialise energy selection int n_ebin = ebounds.size(); int e_idx1 = 0; @@ -455,24 +488,26 @@ void ctcubemask::apply_mask(GCTAObservation* obs) } // Set all pixels inside selected energy bands but outside ROI - // to -1.0 - GSkyRegionCircle roi(m_ra, m_dec, m_rad); - for (int i = e_idx1; i <= e_idx2; ++i) { - for (int pixel = 0; pixel < npix; ++pixel) { - GSkyDir dir = map.inx2dir(pixel); - if (!roi.contains(dir)) { - map(pixel,i) = -1.0; + // to -1.0 if requested + if (m_select_roi) { + GSkyRegionCircle roi(m_ra, m_dec, m_rad); + for (int i = e_idx1; i <= e_idx2; ++i) { + for (int pixel = 0; pixel < npix; ++pixel) { + GSkyDir dir = map.inx2dir(pixel); + if (!roi.contains(dir)) { + map(pixel,i) = -1.0; + } } } - } - // Log selected energy band - if (logTerse()) { - log << gammalib::parformat("Selected ROI"); - log << "RA=" << m_ra << " deg, "; - log << "DEC=" << m_dec << " deg, "; - log << "Radius=" << m_rad << " deg"; - log << std::endl; + // Log selected energy band + if (logTerse()) { + log << gammalib::parformat("Selected ROI"); + log << "RA=" << m_ra << " deg, "; + log << "DEC=" << m_dec << " deg, "; + log << "Radius=" << m_rad << " deg"; + log << std::endl; + } } // Set all pixels inside selected energy bands and inside exclusion @@ -538,6 +573,8 @@ void ctcubemask::init_members(void) // Initialise protected members m_obs.clear(); m_infiles.clear(); + m_select_energy = true; + m_select_roi = true; // Return return; @@ -565,6 +602,8 @@ void ctcubemask::copy_members(const ctcubemask& app) // Copy protected members m_obs = app.m_obs; m_infiles = app.m_infiles; + m_select_energy = app.m_select_energy; + m_select_roi = app.m_select_roi; // Return return; diff --git a/src/ctcubemask/ctcubemask.hpp b/src/ctcubemask/ctcubemask.hpp index 6db740b..06c8812 100644 --- a/src/ctcubemask/ctcubemask.hpp +++ b/src/ctcubemask/ctcubemask.hpp @@ -92,6 +92,8 @@ protected: // Protected members GObservations m_obs; //!< Observations container std::vector m_infiles; //!< Input event filenames + bool m_select_energy; //!< Perform energy selection + bool m_select_roi; //!< Perform ROI selection }; diff --git a/test/test_cscripts.sh.in b/test/test_cscripts.sh.in index 48bf410..252659e 100755 --- a/test/test_cscripts.sh.in +++ b/test/test_cscripts.sh.in @@ -36,7 +36,7 @@ $ECHO "*****************" # # Remove any existing result files # ================================ -rm -rf *.fits *.xml *.dat +rm -rf spectrum*.fits lightcurve*.fits resmap*.fits obs.xml *.dat # @@ -209,7 +209,6 @@ chmod +x @abs_top_srcdir@/cscripts/csspec.py srcname="Crab" \ emin=0.1 \ emax=100.0 \ - binned=no \ outfile="spectrum.fits" > csspec.log $ECHO -n "." if [ -s "spectrum.fits" ] @@ -228,30 +227,22 @@ $ECHO " ok" $ECHO -n "Test csspec (binned): " chmod +x @abs_top_srcdir@/cscripts/csspec.py @abs_top_srcdir@/cscripts/csspec.py \ - inobs="data/crab_events.fits" \ - inmodel="data/crab.xml" \ - caldb="irf" \ - irf="cta_dummy_irf" \ - enumbins=3 \ - srcname="Crab" \ + inobs="cntmap2.fits" \ + expcube="expcube2.fits" \ + psfcube="psfcube2.fits" \ + bkgcube="bkgcube2.fits" \ + inmodel="bkgcube2.xml" \ + enumbins=2 \ 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 + srcname="Crab" \ + outfile="spectrum_binned.fits" > csspec.log $ECHO -n "." -if [ -s "spectrum.fits" ] +if [ -s "spectrum_binned.fits" ] then $ECHO -n "." else - $ECHO " spectrum.fits file is not valid" + $ECHO " spectrum_binned.fits file is not valid" exit 1 fi $ECHO " ok" @@ -309,13 +300,13 @@ chmod +x @abs_top_srcdir@/cscripts/cslightcrv.py nxpix=10 \ nypix=10 \ binsz=0.1 \ - outfile="lightcurve.fits" > cslightcrv.log + outfile="lightcurve_binned.fits" > cslightcrv.log $ECHO -n "." -if [ -s "lightcurve.fits" ] +if [ -s "lightcurve_binned.fits" ] then $ECHO -n "." else - $ECHO " lightcurve.fits file is not valid" + $ECHO " lightcurve_binned.fits file is not valid" exit 1 fi $ECHO " ok" -- libgit2 0.21.4