added phase selection functionality

ldivenere
1 parent 1a11125e
Showing 3 changed files with 45 additions and 0 deletions
src/ctselect/ctselect.cpp
src/ctselect/ctselect.hpp
src/ctselect/ctselect.par
@@ -484,6 +484,7 @@ void ctselect::init_members(void)
     m_select_energy = true;
     m_select_roi    = true;
     m_select_time   = true;
+    m_select_phase  = true;
  
     // Return
     return;
@@ -519,6 +520,7 @@ void ctselect::copy_members(const ctselect&amp; app)
     m_select_energy = app.m_select_energy;
     m_select_roi    = app.m_select_roi;
     m_select_time   = app.m_select_time;
+    m_select_phase  = app.m_select_phase;
  
     // Return
     return;
@@ -551,6 +553,7 @@ void ctselect::get_parameters(void)
     m_select_energy = true;
     m_select_roi    = true;
     m_select_time   = true;
+    m_select_phase  = true;
  
     // Setup observations from "inobs" parameter. Do not request response
     // information and do not accept counts cubes.
@@ -592,6 +595,16 @@ void ctselect::get_parameters(void)
         m_select_energy = false;
     }
  
+    // Check for sanity of phase selection parameters
+    if ((*this)["phasemin"].is_valid() && (*this)["phasemax"].is_valid()) {
+        m_phasemin     = (*this)["phasemin"].real();
+        m_phasemax     = (*this)["phasemax"].real();
+        m_select_phase = true;
+    }
+    else {
+        m_select_phase = false;
+    }
+
     // Get other User parameters
     m_expr     = (*this)["expr"].string();
     m_usethres = (*this)["usethres"].string();
@@ -722,6 +735,33 @@ void ctselect::select_events(GCTAObservation*   obs,
  
     } // endif: made time selection
  
+    // Make phase selection
+    if (m_select_phase) {
+        // Check if event_list has phase (m_has_phase)
+        if (list->has_phase()) {  // gammalib version from branch 1982-add-phase-info necessary
+            // Check if phasemax is larger than phasemin
+	    if (m_phasemax > m_phasemin) {
+	        // Format phase with sufficient accuracy and add to selection string
+	        char cmin[80];
+		char cmax[80];
+		sprintf(cmin, "%.8f", m_phasemin);
+		sprintf(cmax, "%.8f", m_phasemax);
+		selection += add + "PHASE >= "+std::string(cmin)+" && PHASE <= "+std::string(cmax);
+		add       = " && ";
+		log_value(NORMAL, "Phase range",
+			  gammalib::str(m_phasemin)+" - "+gammalib::str(m_phasemax));
+	    }
+	    else {
+	        log << "Invalid values for phasemin and phasemax. ";
+		log << "Phase selection skipped." << std::endl;
+	    }
+	}
+	else {
+	    log << "Event list has no PHASE column. ";
+	    log << "Phase selection skipped." << std::endl;
+	}
+    } // endif: made phase selection
+
     // Make energy selection
     if (m_select_energy) {
  
@@ -97,6 +97,8 @@ protected:
     double      m_tmax;       //!< Stop time
     double      m_emin;       //!< Lower energy
     double      m_emax;       //!< Upper energy
+    double      m_phasemin;   //!< Lower phase
+    double      m_phasemax;   //!< Upper phase
     std::string m_expr;       //!< Selection expression
     std::string m_usethres;   //!< Energy threshold type
     GChatter    m_chatter;    //!< Chattiness
@@ -110,6 +112,7 @@ protected:
     bool                     m_select_energy; //!< Perform energy selection
     bool                     m_select_roi;    //!< Perform RoI selection
     bool                     m_select_time;   //!< Perform time selection
+    bool                     m_select_phase;  //!< Perform phase selection
 };
  
  
@@ -32,6 +32,8 @@ tmin,    r, a, 0.0,,, &quot;Start time (CTA MET in seconds)&quot;
 tmax,    r, a, 0.0,,, "End time (CTA MET in seconds)"
 emin,    r, a, 0.1,,, "Lower energy limit (TeV)"
 emax,    r, a, 100.0,,, "Upper energy limit (TeV)"
+phasemin,r, h, 0.0,,, "Lower phase limit"
+phasemax,r, h, 1.0,,, "Upper phase limit"
 expr,    s, h, "",,, "Additional selection expression"
 usethres,s, h, "NONE",NONE|DEFAULT|USER,, "Energy threshold type"