diff --git a/src/ctmodel/ctmodel.cpp b/src/ctmodel/ctmodel.cpp index f902157..7b5b5f7 100644 --- a/src/ctmodel/ctmodel.cpp +++ b/src/ctmodel/ctmodel.cpp @@ -365,7 +365,8 @@ void ctmodel::init_members(void) m_obs.clear(); m_cube.clear(); m_gti.clear(); - m_has_cube = false; + m_has_cube = false; + m_append_cube = false; // Return return; @@ -384,10 +385,11 @@ void ctmodel::copy_members(const ctmodel& app) m_apply_edisp = app.m_apply_edisp; // Copy protected members - m_obs = app.m_obs; - m_cube = app.m_cube; - m_gti = app.m_gti; - m_has_cube = app.m_has_cube; + m_obs = app.m_obs; + m_cube = app.m_cube; + m_gti = app.m_gti; + m_has_cube = app.m_has_cube; + m_append_cube = app.m_append_cube; // Return return; @@ -412,10 +414,13 @@ void ctmodel::free_members(void) ***************************************************************************/ void ctmodel::get_parameters(void) { + // Reset cube append flag + m_append_cube = false; + // If there are no observations in container then load them via user // parameters if (m_obs.size() == 0) { - m_obs = get_observations(); + get_obs(); } // Read model definition file if required @@ -468,10 +473,145 @@ void ctmodel::get_parameters(void) m_outcube = (*this)["outcube"].filename(); } + // If cube should be appended to first observation then do that now. + // This is a kluge that makes sure that the cube is passed as part + // of the observation in case that a cube response is used. The kluge + // is needed because the GCTACubeSourceDiffuse::set method needs to + // get the full event cube from the observation. + if (m_append_cube) { + + //TODO: Check that energy boundaries are compatible + + // Attach GTI of observations to model cube + m_cube.gti(m_obs[0]->events()->gti()); + + // Attach model cube to observations + m_obs[0]->events(m_cube); + + } // endif: cube was scheduled for appending + // Return return; } + +/***********************************************************************//** + * @brief Get observation container + * + * Get an observation container according to the user parameters. The method + * supports loading of a individual FITS file or an observation definition + * file in XML format. If the input filename is empty, parameters are read + * to build a CTA observation from scratch. + ***************************************************************************/ +void ctmodel::get_obs(void) +{ + // Get the filename from the input parameters + std::string filename = (*this)["inobs"].filename(); + + // If no observation definition file has been specified then read all + // parameters that are necessary to create an observation from scratch + if ((filename == "NONE") || (gammalib::strip_whitespace(filename) == "")) { + + // Get response cube filenames + std::string expcube = (*this)["expcube"].filename(); + std::string psfcube = (*this)["psfcube"].filename(); + + // If the filenames are valid then build an observation from cube + // response information + if ((expcube != "NONE") && (psfcube != "NONE") && + (gammalib::strip_whitespace(expcube) != "") && + (gammalib::strip_whitespace(psfcube) != "")) { + + // Get exposure and PSF cubes + GCTACubeExposure exposure(expcube); + GCTACubePsf psf(psfcube); + + // Create energy boundaries + GEbounds ebounds = create_ebounds(); + + // Create dummy sky map cube + GSkymap map("CAR","GAL",0.0,0.0,1.0,1.0,1,1,ebounds.size()); + m_append_cube = true; + + // Create event cube + GCTAEventCube cube(map, ebounds, exposure.gti()); + + // Create CTA observation + GCTAObservation cta; + cta.events(cube); + cta.response(exposure, psf); + + // Append observation to container + m_obs.append(cta); + + } // endif: cube response information was available + + // ... otherwise build an observation from IRF response information + else { + + // Create CTA observation + GCTAObservation cta = create_cta_obs(); + + // Set response + set_obs_response(&cta); + + // Append observation to container + m_obs.append(cta); + + } + + } // endif: filename was "NONE" or "" + + // ... otherwise we have a file name + else { + + // If file is a FITS file then create an empty CTA observation + // and load file into observation + if (gammalib::is_fits(filename)) { + + // Allocate empty CTA observation + GCTAObservation cta; + + // Load data + cta.load(filename); + + // Set response + set_obs_response(&cta); + + // Append observation to container + m_obs.append(cta); + + // Signal that no XML file should be used for storage + m_use_xml = false; + + } + + // ... otherwise load file into observation container + else { + + // Load observations from XML file + m_obs.load(filename); + + // For all observations that have no response, set the response + // from the task parameters + set_response(m_obs); + + // Set observation boundary parameters (emin, emax, rad) + set_obs_bounds(m_obs); + + // Signal that XML file should be used for storage + m_use_xml = true; + + } // endelse: file was an XML file + + } + + // Return + return; + +} + + /***********************************************************************//** * @brief Fill model into model cube * @@ -515,6 +655,7 @@ void ctmodel::fill_cube(const GCTAObservation* obs) // Compute model value for cube bin double model = m_obs.models().eval(*bin, *obs) * bin->size(); + //double model = m_obs.models().eval(*bin, *ptr) * bin->size(); // Add model to actual value value += model; diff --git a/src/ctmodel/ctmodel.hpp b/src/ctmodel/ctmodel.hpp index 9bcc51c..fe6e9a8 100644 --- a/src/ctmodel/ctmodel.hpp +++ b/src/ctmodel/ctmodel.hpp @@ -78,6 +78,7 @@ protected: void copy_members(const ctmodel& app); void free_members(void); void get_parameters(void); + void get_obs(void); void fill_cube(const GCTAObservation* obs); bool has_cube(void) const; @@ -86,10 +87,11 @@ protected: bool m_apply_edisp; //!< Apply energy dispersion? // Protected members - GObservations m_obs; //!< Observation container - GCTAEventCube m_cube; //!< Model cube - GGti m_gti; //!< Model cube GTIs - bool m_has_cube; //!< Signal if cube has been set or loaded + GObservations m_obs; //!< Observation container + GCTAEventCube m_cube; //!< Model cube + 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 }; diff --git a/src/support/ctool.cpp b/src/support/ctool.cpp index b067e97..391bfc1 100644 --- a/src/support/ctool.cpp +++ b/src/support/ctool.cpp @@ -304,7 +304,7 @@ GObservations ctool::get_observations(const bool& get_response) if (get_response) { // Set response - set_obs_response(&cta_obs); + set_obs_response(&cta_obs); } // endif: response was required diff --git a/test/test_ctmodel.py b/test/test_ctmodel.py index 41ce959..92b2e4f 100755 --- a/test/test_ctmodel.py +++ b/test/test_ctmodel.py @@ -70,6 +70,8 @@ class Test(gammalib.GPythonTestSuite): model["outcube"].filename("modmap.fits") model["inmodel"].filename(self.model_name) model["inobs"].filename("NONE") + model["expcube"].filename("NONE") + model["psfcube"].filename("NONE") model["caldb"].string(self.caldb) model["irf"].string(self.irf) model["rad"].real(5.0) diff --git a/test/test_ctools.sh b/test/test_ctools.sh index aa51ca1..6303fa4 100755 --- a/test/test_ctools.sh +++ b/test/test_ctools.sh @@ -218,10 +218,12 @@ else exit 1 fi # -# Run 3 +# Run 2 $ctmodel inobs="NONE" \ incube="NONE" \ outcube="modmap3.fits" \ + expcube="NONE" \ + psfcube="NONE" \ caldb="irf" \ irf="cta_dummy_irf" \ inmodel="data/crab.xml" \ @@ -249,7 +251,7 @@ else exit 1 fi # -# Run 2 +# Run 3 $ctmodel inobs="cntmap2.fits" \ incube="cntmap2.fits" \ outcube="modmap2.fits" \