Modified ctmodel following the discussions on #1439. In addition, csresmap was a…

…lso 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

Modified ctmodel following the discussions on #1439. In addition, csresmap was a…
…lso 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
Michael Mayer
1 parent 51c98949
Showing 4 changed files with 179 additions and 66 deletions
scripts/csresmap.py
src/ctmodel/ctmodel.cpp
src/ctmodel/ctmodel.hpp
test/test_cscripts.sh.in
@@ -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
@@ -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&amp; 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<GCTAObservation*>(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<GCTAEventCube*>(const_cast<GEvents*>(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;
 }
+
@@ -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
 };
@@ -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