diff --git a/src/ctbin/ctbin.cpp b/src/ctbin/ctbin.cpp index a7eda26..b25cbbf 100644 --- a/src/ctbin/ctbin.cpp +++ b/src/ctbin/ctbin.cpp @@ -227,19 +227,21 @@ void ctbin::bin(void) log << num_outside_ebds << std::endl; } - // Log observation + // Log map if (logTerse()) { log << std::endl; log.header1("Counts map"); log << map << std::endl; } - - // Save counts map - GFits file; - map.write(&file); - obs.ebounds().write(&file); - obs.gti().write(&file); - file.saveto(m_outfile, clobber()); + + // Create events cube from sky map + GCTAEventCube cube(map); + + // Replace event list by event cube in observation + obs.events(&cube); + + // Save observation + obs.save(m_outfile, clobber()); // Return return; diff --git a/src/ctlike/ctlike.cpp b/src/ctlike/ctlike.cpp index 4282361..104d934 100644 --- a/src/ctlike/ctlike.cpp +++ b/src/ctlike/ctlike.cpp @@ -104,6 +104,7 @@ ctlike::~ctlike(void) ***************************************************************************/ void ctlike::run(void) { + log.cout(true); // Get parameters get_parameters(); @@ -119,6 +120,13 @@ void ctlike::run(void) else binned(m_cntmap); + // Write observation into logger + if (logTerse()) { + log << std::endl; + log.header1("Observation"); + log << m_obs << std::endl; + } + // Setup models for optimizing. m_models = GModels(m_srcmdl); m_obs.models(m_models); @@ -265,7 +273,7 @@ void ctlike::unbinned(const std::string& evfile) // for analysis run.load_unbinned(evfile); run.response(m_irf, m_caldb); - run.roi(roi); + run.roi(&roi); run.ebounds(ebds); m_obs.append(run); diff --git a/src/ctobssim/ctobssim.cpp b/src/ctobssim/ctobssim.cpp index 4ec91d1..b7573a0 100644 --- a/src/ctobssim/ctobssim.cpp +++ b/src/ctobssim/ctobssim.cpp @@ -108,6 +108,9 @@ void ctobssim::run(void) if (logTerse()) { log_parameters(); log << std::endl; + log << m_pnt << std::endl; + log << m_obs << std::endl; + log << m_models << std::endl; } // Write header @@ -174,6 +177,36 @@ void ctobssim::get_parameters(void) // Initialise random number generator m_ran.seed(m_seed); + // Set CTA pointing + GSkyDir dir; + dir.radec_deg(m_ra, m_dec); + m_pnt.dir(dir); + + // Set GTI + GGti gti; + GTime tstart; + GTime tstop; + tstart.met(m_tmin); + tstop.met(m_tmax); + gti.append(tstart, tstop); + + // Set energy boundaries + GEbounds ebds; + GEnergy emin; + GEnergy emax; + emin.TeV(m_emin); + emax.TeV(m_emax); + ebds.append(emin, emax); + + // Set observation + m_obs.ebounds(ebds); + m_obs.gti(gti); + m_obs.response(m_irf, m_caldb); + m_obs.pointing(m_pnt); + + // Load models + m_models.load(m_infile); + // Return return; } @@ -198,17 +231,14 @@ GPhotons ctobssim::simulate_photons(void) tmin.met(m_tmin); tmax.met(m_tmax); - // Load models - GModels models(m_infile); - // Allocate photons GPhotons photons; // Simulate photons for all sky models - for (int i = 0; i < models.size(); ++i) { + for (int i = 0; i < m_models.size(); ++i) { // Get sky model (NULL if not a sky model) - GModelSky* model = dynamic_cast(models(i)); + GModelSky* model = dynamic_cast(m_models(i)); // If we have a sky model then simulate photons if (model != NULL) { @@ -260,21 +290,15 @@ GPhotons ctobssim::simulate_photons(void) ***************************************************************************/ GCTAEventList* ctobssim::simulate_events(const GPhotons& photons) { - // Setup CTA response - GCTAResponse rsp(m_irf, m_caldb); - - // Setup CTA pointing - GCTAPointing pnt; - // Allocate CTA event list GCTAEventList* events = new GCTAEventList; // Reserves space for events events->reserve(photons.size()); - // Simulate events + // Simulate events from photons for (int i = 0; i < photons.size(); ++i) { - GCTAEventAtom* event = rsp.mc(m_area, photons[i], pnt, m_ran); + GCTAEventAtom* event = m_obs.response()->mc(m_area, photons[i], m_pnt, m_ran); if (event != NULL) events->append(*event); delete event; @@ -282,9 +306,43 @@ GCTAEventList* ctobssim::simulate_events(const GPhotons& photons) // Dump simulation results if (logNormal()) { - log << parformat("Simulated events") << str(events->size()) << std::endl; + log << parformat("MC source events"); + log << str(events->size()) << std::endl; } + // Simulate background events + for (int i = 0; i < m_models.size(); ++i) { + + // Get model (NULL if not a radial acceptance model) + GCTAModelRadialAcceptance* model = + dynamic_cast(m_models(i)); + + // If we have a radial acceptance model then simulate events + if (model != NULL) { + + // Get simulated event list + GCTAEventList* list = model->mc(m_obs, m_ran); + + // Reserves space for events + events->reserve(list->size()+events->size()); + + // Append events + for (int k = 0; k < list->size(); k++) + events->append(*(list->pointer(k))); + + // Dump simulation results + if (logNormal()) { + log << parformat("MC background events"); + log << str(list->size()) << std::endl; + } + + // Free event list + delete list; + + } // endif: model was valid + + } // endfor: looped over all models + // Return event list return events; } @@ -353,9 +411,13 @@ void ctobssim::set_events_keywords(GFits* file) hdu->card("ONTIME", m_tmax-m_tmin, "[s] Total good time including deadtime"); hdu->card("LIVETIME", m_tmax-m_tmin, "[s] Total livetime"); + // Set target information + //hdu->card("RA_OBJ", m_ra, "[deg] Target Right Ascension"); + //hdu->card("DEC_OBJ", m_dec, "[deg] Target Declination"); + // Set pointing information - hdu->card("RA_OBJ", m_ra, "[deg] Target Right Ascension"); - hdu->card("DEC_OBJ", m_dec, "[deg] Target Declination"); + hdu->card("RA_PNT", m_ra, "[deg] Pointing Right Ascension"); + hdu->card("DEC_PNT", m_dec, "[deg] Pointing Declination"); // Return return; @@ -378,6 +440,8 @@ void ctobssim::init_members(void) m_outfile.clear(); m_caldb.clear(); m_irf.clear(); + m_obs.clear(); + m_pnt.clear(); m_seed = 1; m_ra = 0.0; m_dec = 0.0; diff --git a/src/ctobssim/ctobssim.hpp b/src/ctobssim/ctobssim.hpp index a70de6c..37f0d93 100644 --- a/src/ctobssim/ctobssim.hpp +++ b/src/ctobssim/ctobssim.hpp @@ -54,20 +54,23 @@ protected: void free_members(void); // User parameters - std::string m_infile; //!< Input model - std::string m_outfile; //!< Output events file - std::string m_caldb; //!< Calibration database repository - std::string m_irf; //!< Instrument response function - int m_seed; //!< Random number generator seed - double m_ra; //!< RA of ROI centre - double m_dec; //!< DEC of ROI centre - double m_rad; //!< ROI radius - double m_tmin; //!< Start time (MET) - double m_tmax; //!< Stop time (MET) - double m_emin; //!< Lower energy (TeV) - double m_emax; //!< Upper energy (TeV) - double m_area; //!< Surface area for simulation (cm2) - GRan m_ran; //!< Random number generator + std::string m_infile; //!< Input model + std::string m_outfile; //!< Output events file + std::string m_caldb; //!< Calibration database repository + std::string m_irf; //!< Instrument response function + int m_seed; //!< Random number generator seed + double m_ra; //!< RA of pointing direction + double m_dec; //!< DEC of pointing direction + double m_rad; //!< FOV radius + double m_tmin; //!< Start time (MET) + double m_tmax; //!< Stop time (MET) + double m_emin; //!< Lower energy (TeV) + double m_emax; //!< Upper energy (TeV) + double m_area; //!< Surface area for simulation (cm2) + GRan m_ran; //!< Random number generator + GCTAPointing m_pnt; //!< CTA pointing + GCTAObservation m_obs; //!< CTA observation + GModels m_models; //!< Models }; #endif /* CTOBSSIM_HPP */