Modify csspec to work on binned input, too. Added test case and slighlty updated…

… ctcubemask for this purpose

Modify csspec to work on binned input, too. Added test case and slighlty updated…
… ctcubemask for this purpose
Michael Mayer
1 parent 794678a8
Showing 4 changed files with 144 additions and 183 deletions
cscripts/csspec.py
src/ctcubemask/ctcubemask.cpp
src/ctcubemask/ctcubemask.hpp
test/test_cscripts.sh.in
@@ -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
@@ -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&amp; 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;
@@ -92,6 +92,8 @@ protected:
     // Protected members
     GObservations            m_obs;        //!< Observations container
     std::vector<std::string> m_infiles;    //!< Input event filenames
+    bool                     m_select_energy; //!< Perform energy selection
+    bool                     m_select_roi;    //!< Perform ROI selection
 };
@@ -36,7 +36,7 @@ $ECHO &quot;*****************&quot;
 #
 # 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 &quot; ok&quot;
 $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"