From 706e318f0551e239951821e28ddd02a442680bf9 Mon Sep 17 00:00:00 2001 From: Michael Mayer Date: Wed, 1 Jul 2015 12:21:23 +0200 Subject: [PATCH] Modified ctmodel following the discussions on #1439. In addition, csresmap was also modified consistently, csresmap can now also take two maps (cnt and model) as input and compute the residuals directly. The algorithm parameter has been made unhidden --- scripts/csresmap.py | 183 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------------------------- src/ctmodel/ctmodel.cpp | 58 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/ctmodel/ctmodel.hpp | 1 + test/test_cscripts.sh.in | 3 ++- 4 files changed, 179 insertions(+), 66 deletions(-) diff --git a/scripts/csresmap.py b/scripts/csresmap.py index 1ab36fe..ae4f30b 100755 --- a/scripts/csresmap.py +++ b/scripts/csresmap.py @@ -101,6 +101,8 @@ class csresmap(ctools.cscript): # Create default parfile pars = gammalib.GApplicationPars() pars.append(gammalib.GApplicationPar("inobs","f","a","events.fits","","","Event list, counts cube, or observation definition file")) + pars.append(gammalib.GApplicationPar("cntcube","f","h","NONE","","","Input model cube (generated with ctmodel)")) + pars.append(gammalib.GApplicationPar("modcube","f","h","NONE","","","Input counts cube (generated with ctbin)")) pars.append(gammalib.GApplicationPar("inmodel","f","a","$CTOOLS/share/models/crab.xml","","","Source model")) pars.append(gammalib.GApplicationPar("outmap","f","a","resmap.fits","","","Output residual map")) pars.append(gammalib.GApplicationPar("caldb","s","a","prod2","","","Calibration database")) @@ -117,7 +119,7 @@ class csresmap(ctools.cscript): 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("algorithm","s","h","SUBDIV","SUB|SUBDIV|SUBDIVSQRT","","Residual map computation algorithm")) + pars.append(gammalib.GApplicationPar("algorithm","s","a","SUBDIV","SUB|SUBDIV|SUBDIVSQRT","","Residual map computation algorithm")) pars.append_standard() pars.append(gammalib.GApplicationPar("logfile","f","h","csresmap.log","","","Log filename")) pars.save(parfile) @@ -131,30 +133,58 @@ class csresmap(ctools.cscript): """ # Get parameters - # Set observation if not done before - if self.obs.size() == 0: - self.require_inobs("csresmap.get_parameters()") - self.obs = self.get_observations() - - # Set models if we have none - if self.obs.models().size() == 0: - self.obs.models(self["inmodel"].filename()) - - # Read other parameters - self.m_outfile = self["outmap"].filename() - self.m_xref = self["xref"].real() - self.m_yref = self["yref"].real() - self.m_emin = self["emin"].real() - self.m_emax = self["emax"].real() - self.m_enumbins = self["enumbins"].integer() - self.m_ebinalg = self["ebinalg"].string() - self.m_coordsys = self["coordsys"].string() - self.m_proj = self["proj"].string() - self.m_nxpix = self["nxpix"].integer() - self.m_nypix = self["nypix"].integer() - self.m_binsz = self["binsz"].real() + # Initialise some flags + self.m_use_maps = False + self.m_skip_binning = False + + # First check for cntcube and modcube parameters + self.m_cntcube = self["cntcube"].filename() + self.m_modcube = self["modcube"].filename() + if not (self.m_cntcube == "NONE" or self.m_modcube == "NONE"): + + # Signal that we can use the maps directly without ctbin and ctmodel + self.m_use_maps = True + + # If not two maps are given, proceed to set up observation + if not self.m_use_maps: + + # Set observation if not done before + if self.obs.size() == 0: + self.require_inobs("csresmap.get_parameters()") + self.obs = self.get_observations() + + # Check if we have exactly one binned CTA observation + if self.obs.size() == 1: + + if self.obs[0].classname() == "GCTAObservation": + if self.obs[0].eventtype() == "CountsCube": + # Skip ctbin step later on + self.m_skip_binning = True + + + # Set models if we have none + if self.obs.models().size() == 0: + self.obs.models(self["inmodel"].filename()) + + # Skip query for spatial parameters if a binning is provided in the observation + if not self.m_skip_binning: + # Read other parameters + self.m_xref = self["xref"].real() + self.m_yref = self["yref"].real() + self.m_emin = self["emin"].real() + self.m_emax = self["emax"].real() + self.m_enumbins = self["enumbins"].integer() + self.m_ebinalg = self["ebinalg"].string() + self.m_coordsys = self["coordsys"].string() + self.m_proj = self["proj"].string() + self.m_nxpix = self["nxpix"].integer() + self.m_nypix = self["nypix"].integer() + self.m_binsz = self["binsz"].real() + + # Read necessary parameters + self.m_outfile = self["outmap"].filename() self.m_algorithm = self["algorithm"].string() - + # Set some fixed parameters self.m_log = False # Logging in client tools self.m_chatter = self["chatter"].integer() @@ -210,53 +240,76 @@ class csresmap(ctools.cscript): self.log(str(self.obs)) self.log("\n") - # Write header - if self.logTerse(): - self.log("\n") - self.log.header1("Generate binned map (ctbin)") - - # Create countsmap - bin = ctools.ctbin(self.obs) - bin["nxpix"].integer(self.m_nxpix) - bin["nypix"].integer(self.m_nypix) - bin["proj"].string(self.m_proj) - bin["coordsys"].string(self.m_coordsys) - bin["xref"].real(self.m_xref) - bin["yref"].real(self.m_yref) - bin["enumbins"].integer(self.m_enumbins) - bin["ebinalg"].string(self.m_ebinalg) - bin["emin"].real(self.m_emin) - bin["emax"].real(self.m_emax) - bin["binsz"].real(self.m_binsz) - bin["chatter"].integer(self.m_chatter) - bin["clobber"].boolean(self.m_clobber) - bin["debug"].boolean(self.m_debug) - bin.run() + countmap = gammalib.GSkymap() + modelmap = gammalib.GSkymap() + + # Use input file directly if given + if self.m_use_maps: + countmap = gammalib.GSkymap(self["cntcube"].filename()) + modelmap = gammalib.GSkymap(self["modcube"].filename()) + + else: + + cta_counts_cube = gammalib.GCTAEventCube() + + if self.m_skip_binning: + cta_counts_cube = gammalib.GCTAEventCube(self.obs[0].events().clone()) + else: + + # Write header + if self.logTerse(): + self.log("\n") + self.log.header1("Generate binned map (ctbin)") + + # Create countsmap + bin = ctools.ctbin(self.obs) + bin["nxpix"].integer(self.m_nxpix) + bin["nypix"].integer(self.m_nypix) + bin["proj"].string(self.m_proj) + bin["coordsys"].string(self.m_coordsys) + bin["xref"].real(self.m_xref) + bin["yref"].real(self.m_yref) + bin["enumbins"].integer(self.m_enumbins) + bin["ebinalg"].string(self.m_ebinalg) + bin["emin"].real(self.m_emin) + bin["emax"].real(self.m_emax) + bin["binsz"].real(self.m_binsz) + bin["chatter"].integer(self.m_chatter) + bin["clobber"].boolean(self.m_clobber) + bin["debug"].boolean(self.m_debug) + bin.run() + + # Retrieve counts cube + cta_counts_cube = bin.cube() + + # Assign GCTAEventCube to skymap + countmap = cta_counts_cube.map() + + # Write header + if self.logTerse(): + self.log("\n") + self.log.header1("Generate model map (ctmodel)") + + # Create model map + model = ctools.ctmodel(self.obs) + model.cube(cta_counts_cube) + model["chatter"].integer(self.m_chatter) + model["clobber"].boolean(self.m_clobber) + model["debug"].boolean(self.m_debug) + model.run() + + # Get model map into GSkymap object + modelmap = model.cube().map().copy() + + # Store counts map as residual map. Note that we need a # special construct here to avoid memory leaks. This seems # to be a SWIG feature as SWIG creates a new object when # calling bin.cube() - #residualmap = bin.cube().map() - counts_cube = bin.cube() - self.resmap = counts_cube.map().copy() + #residualmap = bin.cube().map() + self.resmap = countmap.copy() self.resmap.stack_maps() - - # Write header - if self.logTerse(): - self.log("\n") - self.log.header1("Generate model map (ctmodel)") - - # Create model map - model = ctools.ctmodel(self.obs) - model.cube(bin.cube()) - model["chatter"].integer(self.m_chatter) - model["clobber"].boolean(self.m_clobber) - model["debug"].boolean(self.m_debug) - model.run() - - # Get model map into GSkymap object - modelmap = model.cube().map() modelmap.stack_maps() # Continue calculations depending on given algorithm diff --git a/src/ctmodel/ctmodel.cpp b/src/ctmodel/ctmodel.cpp index 862c0fa..7cafb99 100644 --- a/src/ctmodel/ctmodel.cpp +++ b/src/ctmodel/ctmodel.cpp @@ -270,9 +270,34 @@ void ctmodel::run(void) } } + // Fill cube and leave loop if we are binned mode (meaning we only have + // one binned observation) + if (m_binned && m_obs.size() == 1) { + fill_cube(obs); + break; + } + + // Skip observation if we don't have an unbinned observation + if (obs->eventtype() != "EventList") { + + // Log that we skip the this observation + if (logTerse()) { + log << "Warning: Skipping binned observation \""+obs->name()+"\""; + log << std::endl; + } + continue; + } + + // Fill the cube fill_cube(obs); + // Leave loop if we are binned mode (meaning we only have + // one binned observation + if (m_binned && m_obs.size() == 1) { + break; + } + // Dispose events to free memory if event file exists on disk if (obs->eventfile().length() > 0 && gammalib::file_exists(obs->eventfile())) { @@ -374,6 +399,7 @@ void ctmodel::init_members(void) m_gti.clear(); m_has_cube = false; m_append_cube = false; + m_binned = false; // Return return; @@ -397,6 +423,7 @@ void ctmodel::copy_members(const ctmodel& app) m_gti = app.m_gti; m_has_cube = app.m_has_cube; m_append_cube = app.m_append_cube; + m_binned = app.m_binned; // Return return; @@ -430,6 +457,36 @@ void ctmodel::get_parameters(void) m_obs = get_observations(); } + // Check if we got excactly one binned CTA observation + if (m_obs.size() == 1) { + + // Get CTA observation + GCTAObservation* obs = dynamic_cast(m_obs[0]); + + // Continue only if observation is a CTA observation + if (obs != NULL) { + + // Check for binned observation + if (obs->eventtype() == "CountsCube") { + + // Set cube from binned observation + GCTAEventCube* evtcube = dynamic_cast(const_cast(obs->events())); + + cube(*evtcube); + + // Signal that cube has been set + m_has_cube = true; + + // Signal that we are in binned mode + m_binned = true; + + } // endif: observation was binned + + } // endif: observation was CTA + + } // endif: had exactly one observation + + // Read model definition file if required if (m_obs.models().size() == 0) { @@ -583,3 +640,4 @@ void ctmodel::fill_cube(const GCTAObservation* obs) // Return return; } + diff --git a/src/ctmodel/ctmodel.hpp b/src/ctmodel/ctmodel.hpp index 9c37587..5d600cc 100644 --- a/src/ctmodel/ctmodel.hpp +++ b/src/ctmodel/ctmodel.hpp @@ -92,6 +92,7 @@ protected: GGti m_gti; //!< Model cube GTIs bool m_has_cube; //!< Signal if cube has been set or loaded bool m_append_cube; //!< Signal that cube should be appended + bool m_binned; //!< Signals that we are in binned mode }; diff --git a/test/test_cscripts.sh.in b/test/test_cscripts.sh.in index 2010c40..f1d0d8e 100755 --- a/test/test_cscripts.sh.in +++ b/test/test_cscripts.sh.in @@ -182,7 +182,8 @@ chmod +x @abs_top_srcdir@/scripts/csresmap.py xref=83.63 \ yref=22.01 \ proj="CAR" \ - binsz=0.05 > csresmap.log + binsz=0.05 \ + algorithm="SUBDIV" > csresmap.log $ECHO -n "." if [ -s "resmap.fits" ] then -- libgit2 0.21.4