diff --git a/pyext/ctbkgcube.i b/pyext/ctbkgcube.i index bc950b9..45e18eb 100644 --- a/pyext/ctbkgcube.i +++ b/pyext/ctbkgcube.i @@ -47,8 +47,7 @@ public: void clear(void); void run(void); void save(void); - const GCTACubeBackground& bgcube(void) const; - const GCTAEventCube& cube(void) const; + const GCTACubeBackground& cube(void) const; const GModels& models(void) const; }; diff --git a/src/ctbkgcube/ctbkgcube.cpp b/src/ctbkgcube/ctbkgcube.cpp index 43718b4..e1cec89 100644 --- a/src/ctbkgcube/ctbkgcube.cpp +++ b/src/ctbkgcube/ctbkgcube.cpp @@ -209,9 +209,6 @@ void ctbkgcube::run(void) log.cout(true); } - // Initialise run parameters - m_livetime = 0.0; - // Get task parameters get_parameters(); @@ -248,9 +245,9 @@ void ctbkgcube::run(void) // Copy models from observation container and reset output model // container + GModels models_orig = m_obs.models(); m_bkgmdl = m_obs.models(); m_outmdl.clear(); - m_cube_model = -1; // Remove all models that are not CTA background models from the // container and put all removed components in the output @@ -309,59 +306,11 @@ void ctbkgcube::run(void) log.header1("Generate background cube"); } - // Initialise observation counter - int n_observations = 0; - - // Loop over all observations in the container - for (int i = 0; i < m_obs.size(); ++i) { - - // Get CTA observation - GCTAObservation* obs = dynamic_cast(m_obs[i]); - - // Continue only if observation is a CTA observation - if (obs != NULL) { - - // Write header for observation - if (logTerse()) { - if (obs->name().length() > 1) { - log.header3("Observation "+obs->name()); - } - else { - log.header3("Observation "+gammalib::str(i)); - } - } - - // Increment number of observations - n_observations++; - - // Fill the cube - fill_cube(obs); - - } // endif: CTA observation found - - } // endfor: looped over observations - - // Re-normalize cube to get units of counts/MeV/s/sr - if (m_livetime > 0.0) { - - // Loop over all bins in background cube and divide the content - // by the total livetime. - for (int i = 0; i < m_bkgcube.size(); ++i) { - GCTAEventBin* bin = m_bkgcube[i]; - double rate = bin->counts() / m_livetime; - bin->counts(rate); - } - - } // endif: livetime was positive - - // Log results - if (logTerse()) { - log << gammalib::parformat("Number of observations"); - log << n_observations << std::endl; - } + // Assign background models to container + m_obs.models(m_bkgmdl); - // Set cube background response - m_background.set(m_bkgcube.map(), m_bkgcube.ebounds()); + // Fill background cube from observations + m_background.fill(m_obs); // Create a background model for the output background cube and append // that model to the input model in place of the original @@ -391,6 +340,9 @@ void ctbkgcube::run(void) log << m_outmdl << std::endl; } + // Recover original models + m_obs.models(models_orig); + // Return return; } @@ -444,13 +396,10 @@ void ctbkgcube::init_members(void) m_outcube.clear(); m_outmodel.clear(); m_obs.clear(); - m_bkgcube.clear(); m_background.clear(); m_bkgmdl.clear(); m_outmdl.clear(); m_ebounds.clear(); - m_livetime = 0.0; - m_cube_model = -1; // Return return; @@ -468,13 +417,10 @@ void ctbkgcube::copy_members(const ctbkgcube& app) m_outmodel = app.m_outmodel; m_outcube = app.m_outcube; m_obs = app.m_obs; - m_bkgcube = app.m_bkgcube; m_background = app.m_background; m_bkgmdl = app.m_bkgmdl; m_outmdl = app.m_outmdl; m_ebounds = app.m_ebounds; - m_livetime = app.m_livetime; - m_cube_model = app.m_cube_model; // Return return; @@ -522,23 +468,18 @@ void ctbkgcube::get_parameters(void) (gammalib::strip_whitespace(incube) == "")) { // Create an event cube based on task parameters - m_bkgcube = create_cube(m_obs); + GCTAEventCube cube = create_cube(m_obs); + + // Define background cube + m_background = GCTACubeBackground(cube); } // endif: filename was not valid // ... otherwise setup the background cube from the counts map else { - - // Load event cube from filename - m_bkgcube.load(incube); - // Loop over all bins in background cube - for (int i = 0; i < m_bkgcube.size(); ++i) { - // Get event bin - GCTAEventBin* bin = m_bkgcube[i]; - bin->counts(0.0); - - } + // Define background cube from file + m_background = GCTACubeBackground(incube); } // endelse: cube was loaded from file @@ -557,8 +498,8 @@ void ctbkgcube::get_parameters(void) m_obs.models(models); } - // Get energy definition - m_ebounds = m_bkgcube.ebounds(); +// // Get energy definition +// m_ebounds = m_background.ebounds(); // Read output filenames (if needed) if (read_ahead()) { @@ -569,92 +510,3 @@ void ctbkgcube::get_parameters(void) // Return return; } - -/***********************************************************************//** - * @brief Generate background cube - * - * @param[in] obs Pointer to CTA observation. - * - * Fills the background cube with the model value for a CTA observation. - * This method will not need to load the event file as it gets the event - * ROI and GTI directly from the observation. - ***************************************************************************/ -void ctbkgcube::fill_cube(GCTAObservation* obs) -{ - - // Continue only if observation pointer is valid - if (obs != NULL) { - - // Initialise statistics - double sum = 0.0; - int n_bins_inside = 0; - int n_bins_outside = 0; - - // Extract region of interest from CTA observation - GCTARoi roi = obs->roi(); - - // Set GTI of actual observations as the GTI of the event cube - m_bkgcube.gti(obs->gti()); - - // Get observation livetime - double livetime = obs->livetime(); - - // Loop over all bins in background cube - for (int i = 0; i < m_bkgcube.size(); ++i) { - - // Get event bin - GCTAEventBin* bin = m_bkgcube[i]; - - // Continue only if binned in contained in ROI - if (roi.contains(*bin)) { - - // Compute model value for event bin. The model value is - // given in counts/MeV/s/sr. - double model = m_bkgmdl.eval(*bin, *obs); - - // Compute number events - sum += model * bin->size(); - - // Multiply by livetime to get the correct weighting for - // each observation. We divide by the total livetime later - // to get the background model in units of counts/MeV/s/sr. - model *= livetime; - - // Add existing number of counts - model += bin->counts(); - - // Store cumulated value (units: counts/MeV/sr) - bin->counts(model); - - // Increment contained bin counter - n_bins_inside++; - - } // endif: bin was contained in RoI - - // ... otherwise count bin as an outsider - else { - n_bins_outside++; - } - - } // endfor: looped over all bins - - // Accumulate livetime - m_livetime += livetime; - - // Log results - if (logTerse()) { - log << gammalib::parformat("Bins within RoI"); - log << n_bins_inside << std::endl; - log << gammalib::parformat("Bins outside RoI"); - log << n_bins_outside << std::endl; - log << gammalib::parformat("Background events in cube"); - log << sum << std::endl; - log << gammalib::parformat("Cube livetime"); - log << livetime << " sec" << std::endl; - } - - } // endif: observation pointer was not valid - - // Return - return; -} diff --git a/src/ctbkgcube/ctbkgcube.hpp b/src/ctbkgcube/ctbkgcube.hpp index 1a1306e..eb1b30c 100644 --- a/src/ctbkgcube/ctbkgcube.hpp +++ b/src/ctbkgcube/ctbkgcube.hpp @@ -59,8 +59,7 @@ public: void clear(void); void run(void); void save(void); - const GCTACubeBackground& bgcube(void) const; - const GCTAEventCube& cube(void) const; + const GCTACubeBackground& cube(void) const; const GModels& models(void) const; protected: @@ -75,27 +74,13 @@ protected: std::string m_outcube; //!< Filename of output cube std::string m_outmodel; //!< Filename of output XML model GObservations m_obs; //!< Observation container - GCTAEventCube m_bkgcube; //!< Background event cube GCTACubeBackground m_background; //!< Background cube response GModels m_bkgmdl; //!< CTA background models GModels m_outmdl; //!< Output models GEbounds m_ebounds; //!< Energy boundaries - double m_livetime; //!< Total livetime for normalization - int m_cube_model; //!< Cube model index -}; +}; -/***********************************************************************//** - * @brief Return event cube containing background rate - * - * @return Background event cube. - ***************************************************************************/ -inline -const GCTAEventCube& ctbkgcube::cube(void) const -{ - // Return background cube - return (m_bkgcube); -} /***********************************************************************//** * @brief Return background response cube containing background rate @@ -103,7 +88,7 @@ const GCTAEventCube& ctbkgcube::cube(void) const * @return Background response cube. ***************************************************************************/ inline -const GCTACubeBackground& ctbkgcube::bgcube(void) const +const GCTACubeBackground& ctbkgcube::cube(void) const { // Return background cube return (m_background);