diff --git a/src/ctmodel/ctmodel.cpp b/src/ctmodel/ctmodel.cpp index b3010d7..c4fb162 100644 --- a/src/ctmodel/ctmodel.cpp +++ b/src/ctmodel/ctmodel.cpp @@ -235,48 +235,55 @@ void ctmodel::run(void) log_header1(TERSE, "Generate model cube"); // Loop over all observations in the container - #pragma omp parallel for - for (int i = 0; i < m_obs.size(); ++i) { + #pragma omp parallel + { + // Create a models object for this core + GModels models(m_obs.models()); - // Get CTA observation - GCTAObservation* obs = dynamic_cast(m_obs[i]); - - // Skip observation if it's not CTA - if (obs == NULL) { - log_header3(TERSE, get_obs_header(m_obs[i])); - std::string msg = " Skipping "+m_obs[i]->instrument()+ - " observation"; - log_string(NORMAL, msg); - } - - // Fill cube and leave loop if we are binned mode (meaning we - // only have one binned observation) - else if (m_binned) { - fill_cube(obs); - i = m_obs.size(); - } - - // Skip observation if we have a binned observation - else if (obs->eventtype() == "CountsCube") { - log_header3(TERSE, get_obs_header(m_obs[i])); - std::string msg = " Skipping binned "+m_obs[i]->instrument()+ - " observation"; - log_string(NORMAL, msg); - } - - // Otherwise, everything seems to be fine, so fill the cube from obs - else { - // Fill the cube - fill_cube(obs); + #pragma omp for schedule(static,1) + for (int i = 0; i < m_obs.size(); ++i) { + + // Get CTA observation + GCTAObservation* obs = dynamic_cast(m_obs[i]); - // Dispose events to free memory if event file exists on disk - if (obs->eventfile().length() > 0 && obs->eventfile().exists()) { - obs->dispose_events(); + // Skip observation if it's not CTA + if (obs == NULL) { + log_header3(TERSE, get_obs_header(m_obs[i])); + std::string msg = " Skipping "+m_obs[i]->instrument()+ + " observation"; + log_string(NORMAL, msg); } - } - - } // endfor: looped over observations + + // Fill cube and leave loop if we are binned mode (meaning we + // only have one binned observation) + else if (m_binned) { + fill_cube(obs, models); + i = m_obs.size(); + } + + // Skip observation if we have a binned observation + else if (obs->eventtype() == "CountsCube") { + log_header3(TERSE, get_obs_header(m_obs[i])); + std::string msg = " Skipping binned "+m_obs[i]->instrument()+ + " observation"; + log_string(NORMAL, msg); + } + + // Otherwise, everything seems to be fine, so fill the cube from obs + else { + // Fill the cube + fill_cube(obs, models); + + // Dispose events to free memory if event file exists on disk + if (obs->eventfile().length() > 0 && obs->eventfile().exists()) { + obs->dispose_events(); + } + } + + } // endfor: looped over observations + } // end openmp parallel section + // Write model cube into header log_header1(NORMAL, "Model cube"); log_string(NORMAL, m_cube.print()); @@ -789,7 +796,7 @@ void ctmodel::extract_cube_properties(void) * GTI of the model cube so that cube GTI is a list of the GTIs of all * observations that were used to generate the model cube. ***************************************************************************/ -void ctmodel::fill_cube(const GCTAObservation* obs) +void ctmodel::fill_cube(const GCTAObservation* obs, GModels& models) { // Get number of spatial and spectral bins in counts cube int npix = m_cube.npix(); @@ -819,7 +826,7 @@ void ctmodel::fill_cube(const GCTAObservation* obs) int num_outside_roi = 0; // Extract copy of models from observation container - GModels models = trim_models(m_obs.models(), roi); + GModels ex_models = trim_models(models, roi); // Get pointer to event cube pixels double* pixels = const_cast(m_cube.counts().pixels()); @@ -877,6 +884,9 @@ void ctmodel::fill_cube(const GCTAObservation* obs) } // endfor: looped over all cube bins + // Re-append the models to the main model + models.extend(ex_models); + // Update GTIs #pragma omp critical(ctmodel_fill_cube) { @@ -918,27 +928,24 @@ void ctmodel::fill_cube(const GCTAObservation* obs) * point sources (which have radius=0) are appropriately included if they * fall near the edge of the observation. ***************************************************************************/ -GModels ctmodel::trim_models(const GModels& all_models, const GCTARoi& roi) +GModels ctmodel::trim_models(GModels& all_models, const GCTARoi& roi) { + // Create the model object to be returned + GModels excluded_models; + // Do nothing if roi is not valid. This will be the case for binned data // sets as they have the default roi radius of 0. In this case we would // remove ALL models, which is obviously incorrect. So we keep all models // when filling based on a counts cube. if (!roi.is_valid()) { - return all_models; + return excluded_models; } // Remove all models that dont overlap with the region of interest. Note // that an extra factor is used since point sources have regions of // radius 0, which is a problem when they fall just barely outside the // ROI. This ensures point sources are appropriately included. - GSkyRegionCircle obsreg( - roi.centre().dir().ra_deg(), - roi.centre().dir().dec_deg(), - roi.radius()+0.5); // Note the buffer - - // Create the model object to be returned - GModels contained_models; + GSkyRegionCircle obsreg(roi.centre().dir(), roi.radius()+0.5); // Loop over the models in the passed container for (int mod=0; modclone()); - - } else if (model->spatial()->region()->overlaps(obsreg)) { - // Model represents a source, that overlaps the observation - // and should be kept. - contained_models.append(*model); - - // Delete model since it isnt needed anymore - delete model; + continue; + } + // Otherwise, exclude the model if it isn't a diffuse cube and it + // doesn't overlap with the observation. + else if (!(model->spatial()->code() == GMODEL_SPATIAL_DIFFUSE) && + (!model->spatial()->region()->overlaps(obsreg))) { + // Model represents a source, that doesn't overlap the observation + // and should be removed for now. + excluded_models.append(*model); + all_models.remove(mod); + mod--; // decrement to prevent skipping models } + delete model; } - return contained_models; + return excluded_models; } diff --git a/src/ctmodel/ctmodel.hpp b/src/ctmodel/ctmodel.hpp index df58c7c..587061a 100644 --- a/src/ctmodel/ctmodel.hpp +++ b/src/ctmodel/ctmodel.hpp @@ -80,9 +80,9 @@ protected: void get_parameters(void); void get_obs(void); void extract_cube_properties(void); - void fill_cube(const GCTAObservation* obs); + void fill_cube(const GCTAObservation* obs, GModels& models); bool has_cube(void) const; - GModels trim_models(const GModels& all_models, const GCTARoi& roi); + GModels trim_models(GModels& all_models, const GCTARoi& roi); // User parameters GFilename m_outcube; //!< Output model cube