diff --git a/src/ctobssim/ctobssim.cpp b/src/ctobssim/ctobssim.cpp index f46bca4..27c1c22 100644 --- a/src/ctobssim/ctobssim.cpp +++ b/src/ctobssim/ctobssim.cpp @@ -256,6 +256,18 @@ void ctobssim::run(void) } } + // Set energy dispersion flag for all CTA observations and save old + // values in save_edisp vector + std::vector save_edisp; + save_edisp.assign(m_obs.size(), false); + for (int i = 0; i < m_obs.size(); ++i) { + GCTAObservation* obs = dynamic_cast(m_obs[i]); + if (obs != NULL) { + save_edisp[i] = obs->response().apply_edisp(); + obs->response().apply_edisp(m_apply_edisp); + } + } + // Write observation(s) into logger if (logTerse()) { log << std::endl; @@ -349,6 +361,14 @@ void ctobssim::run(void) m_use_xml = true; } + // Restore energy dispersion flag for all CTA observations + for (int i = 0; i < m_obs.size(); ++i) { + GCTAObservation* obs = dynamic_cast(m_obs[i]); + if (obs != NULL) { + obs->response().apply_edisp(save_edisp[i]); + } + } + // Return return; } @@ -468,7 +488,8 @@ void ctobssim::get_parameters(void) } // Get other parameters - m_seed = (*this)["seed"].integer(); + m_seed = (*this)["seed"].integer(); + m_apply_edisp = (*this)["edisp"].boolean(); // Optionally read ahead parameters so that they get correctly // dumped into the log file @@ -544,15 +565,14 @@ void ctobssim::set_list(GCTAObservation* obs) if (obs != NULL) { // Get CTA observation parameters - m_ra = (*this)["ra"].real(); - m_dec = (*this)["dec"].real(); - m_rad = (*this)["rad"].real(); - m_tmin = (*this)["tmin"].real(); - m_tmax = (*this)["tmax"].real(); - m_emin = (*this)["emin"].real(); - m_emax = (*this)["emax"].real(); - m_apply_edisp = (*this)["edisp"].boolean(); - m_deadc = (*this)["deadc"].real(); + m_ra = (*this)["ra"].real(); + m_dec = (*this)["dec"].real(); + m_rad = (*this)["rad"].real(); + m_tmin = (*this)["tmin"].real(); + m_tmax = (*this)["tmax"].real(); + m_emin = (*this)["emin"].real(); + m_emax = (*this)["emax"].real(); + m_deadc = (*this)["deadc"].real(); // Allocate CTA event list GCTAEventList events; @@ -628,10 +648,8 @@ void ctobssim::simulate_source(GCTAObservation* obs, const GModels& models, wrklog = &log; } - // Get pointer on CTA response. Throw an exception if the response - // is not defined. - GCTAResponse& rsp = const_cast(obs->response()); - rsp.apply_edisp(m_apply_edisp); + // Get CTA response + const GCTAResponse& rsp = obs->response(); // Make sure that the observation holds a CTA event list. If this // is not the case then allocate and attach a CTA event list now. @@ -686,13 +704,24 @@ void ctobssim::simulate_source(GCTAObservation* obs, const GModels& models, GEnergy emin = events->ebounds().emin(ie); GEnergy emax = events->ebounds().emax(ie); + // Set true photon energy limits for simulation. If observation + // has energy dispersion then add margin + GEnergy e_true_min = emin; + GEnergy e_true_max = emax; + if (rsp.use_edisp()) { + e_true_min = rsp.ebounds_src(e_true_min).emin(); + e_true_max = rsp.ebounds_src(e_true_max).emax(); + } + // Dump energy range if (logNormal()) { if (events->ebounds().size() > 1) { indent++; wrklog->indent(indent); } - *wrklog << gammalib::parformat("Energy range", indent); + *wrklog << gammalib::parformat("Photon energy range", indent); + *wrklog << e_true_min << " - " << e_true_max << std::endl; + *wrklog << gammalib::parformat("Event energy range", indent); *wrklog << emin << " - " << emax << std::endl; } @@ -757,7 +786,7 @@ void ctobssim::simulate_source(GCTAObservation* obs, const GModels& models, // Get photons GPhotons photons = model->mc(m_area, dir, rad, - emin, emax, + e_true_min, e_true_max, tstart, tstop, ran); // Dump number of simulated photons @@ -785,7 +814,12 @@ void ctobssim::simulate_source(GCTAObservation* obs, const GModels& models, if (event != NULL) { // Use event only if it falls within ROI - if (events->roi().contains(*event)) { + // energy boundary and time slice + if (events->roi().contains(*event) && + event->energy() >= emin && + event->energy() <= emax && + event->time() >= tstart && + event->time() <= tstop) { event->event_id(m_event_id); events->append(*event); m_event_id++; @@ -998,7 +1032,7 @@ void ctobssim::init_members(void) m_tmax = 0.0; m_emin = 0.0; m_emax = 0.0; - m_apply_edisp = true; + m_apply_edisp = false; m_deadc = 1.0; // Initialise protected members diff --git a/src/ctobssim/ctobssim.par b/src/ctobssim/ctobssim.par index 7c5e750..c246ccf 100644 --- a/src/ctobssim/ctobssim.par +++ b/src/ctobssim/ctobssim.par @@ -17,13 +17,13 @@ # # File information #================= -infile, f, a, "test/data/crab.xml",,, "Model" -outfile,f, a, "events.fits",,, "Output event data file or observation definition file" -prefix, s, h, "sim_events_",,, "Prefix for event lists in observation definition file" +infile, f, a, "$CTOOLS/share/models/crab.xml",,, "Model" +outfile,f, a, "events.fits",,, "Output event data file or observation definition file" +prefix, s, h, "sim_events_",,, "Prefix for event lists in observation definition file" # # Observation parameters #======================= -caldb, s, a, "",,, "Calibration database" +caldb, s, a, "$CTOOLS/share/caldb/data/cta/dummy",,, "Calibration database" irf, s, a, "cta_dummy_irf",,, "Instrument response function" # # Simulation parameters