43 #define G_GET_PARAMETERS "ctobssim::get_parameters()"
44 #define G_SIMULATE_SOURCE "ctobssim::simulate_source(GCTAObservation*, "\
45 "GModels&, GRan&, GLog*)"
46 #define G_SIMULATE_INTERVAL "ctobssim::simulate_interval(GCTAObservation*"\
47 ", GCTAEventList*, GCTAResponseIrf*, GModels&, GTime&"\
48 ", GTime&, GEnergy&, GEnergy&, GEnergy&, GEnergy&"\
49 ", GSkyDir&, double&, double&, GRan&, GLog*, int&)"
50 #define G_GET_AREA "ctobssim::get_area(GCTAObservation* obs, GEnergy&, "\
214 this->GApplication::clear();
252 int n_observations = 0;
253 for (
int i = 0; i <
m_obs.size(); ++i) {
254 GCTAObservation*
obs =
dynamic_cast<GCTAObservation*
>(
m_obs[i]);
262 if (n_observations > 1) {
267 log_header1(NORMAL,
"Execution mode");
269 ?
"Save and dispose (reduces memory needs)"
270 :
"Keep events in memory";
272 ?
"Write Observation Definition XML file"
273 :
"Write single event list FITS file";
274 log_value(NORMAL,
"Event list management", mode);
275 log_value(NORMAL,
"Output format", xml);
278 log_header1(NORMAL,
"Seed values");
279 for (
int i = 0; i <
m_rans.size(); ++i) {
280 log_value(NORMAL,
"Seed "+gammalib::str(i),
281 gammalib::str(
m_rans[i].seed()));
288 log_header1(TERSE, gammalib::number(
"Simulate observation",
m_obs.size()));
305 wrklog.date(log.date());
306 wrklog.name(log.name());
309 wrklog.buffer_size(10000000);
314 for (
int i = 0; i <
m_obs.size(); ++i) {
322 GCTAObservation*
obs =
dynamic_cast<GCTAObservation*
>(
m_obs[i]);
327 wrklog <<
" Skipping ";
328 wrklog <<
m_obs[i]->instrument();
329 wrklog <<
" observation" << std::endl;
335 if (obs->eventtype() ==
"CountsCube") {
337 wrklog <<
" Skipping binned ";
338 wrklog << obs->instrument();
339 wrklog <<
" observation" << std::endl;
345 if (obs->id().empty()) {
348 obs->id(std::string(buffer));
355 GCTAEventList* events =
static_cast<GCTAEventList*
>(obs->events());
356 events->remove(0, events->size());
361 GCTAObservation obs_clone = *
obs;
364 int events_before = obs_clone.events()->size();
377 wrklog << gammalib::parformat(
"MC events");
378 wrklog << obs_clone.events()->size() - events_before;
379 wrklog <<
" (all models)";
384 obs->events(*(obs_clone.events()));
393 obs->eventfile(outfile);
397 #pragma omp critical(ctobssim_run)
400 obs->save(outfile, clobber());
404 obs->dispose_events();
412 #pragma omp critical (log)
423 if ((*
this)[
"publish"].
boolean()) {
454 log_header1(TERSE, gammalib::number(
"Save observation",
m_obs.size()));
479 log_header1(TERSE, gammalib::number(
"Publish event list",
m_obs.size()));
482 for (
int i = 0; i <
m_obs.size(); ++i) {
485 GCTAObservation*
obs =
dynamic_cast<GCTAObservation*
>(
m_obs[i]);
491 if (obs->events()->size() != 0) {
494 std::string user_name(name);
495 if (user_name.empty()) {
501 user_name += gammalib::str(i);
505 log_value(NORMAL,
"Event list name", user_name);
512 fits.publish(gammalib::extname_cta_events, user_name);
618 if (
m_obs.size() == 0) {
645 if (
m_obs.models().size() == 0) {
646 std::string inmodel = (*this)[
"inmodel"].filename();
647 m_obs.models(inmodel);
651 m_seed = (*this)[
"seed"].integer();
652 m_eslices = (*this)[
"eslices"].integer();
659 (*this)[
"outevents"].query();
660 (*this)[
"prefix"].query();
676 std::vector<unsigned long long int> seeds;
679 for (
int i = 0; i <
m_obs.size(); ++i) {
682 unsigned long long int new_seed;
688 new_seed = (
unsigned long long int)(master.uniform() * 1.0e10) +
691 for (
int j = 0; j < seeds.size(); ++j) {
692 if (new_seed == seeds[j]) {
700 seeds.push_back(new_seed);
704 m_rans.push_back(GRan(new_seed));
710 int nthreads = (*this)[
"nthreads"].integer();
712 omp_set_num_threads(nthreads);
749 const GModels& models,
754 #if defined(G_SOURCE_DEBUG)
755 std::cout <<
"ctobssim::simulate_source: in" << std::endl;
762 if (wrklog == NULL) {
767 GCTAEventList* events =
static_cast<GCTAEventList*
>(obs->events());
770 const GCTAResponseIrf* rsp =
771 dynamic_cast<const GCTAResponseIrf*
>(obs->response());
773 std::string cls = std::string(
typeid(obs->response()).name());
774 std::string msg =
"Response of type \""+cls+
"\" is not a CTA "
775 "IRF response. Please make sure that a CTA "
776 "IRF response is contained in the CTA "
782 GSkyDir dir = events->roi().centre().dir();
790 *wrklog << gammalib::parformat(
"Simulation cone");
791 *wrklog <<
"RA=" << dir.ra_deg() <<
" deg";
792 *wrklog <<
", Dec=" << dir.dec_deg() <<
" deg";
793 *wrklog <<
", radius=" << rad <<
" deg" << std::endl;
800 std::vector<int> nphotons(models.size(),0);
801 std::vector<int> nevents(models.size(),0);
804 for (
int it = 0; it < events->gti().size(); ++it) {
807 GTime tmin = events->gti().tstart(it);
808 GTime tmax = events->gti().tstop(it);
809 GTimeReference tref = events->gti().reference();
814 if (events->gti().size() > 1) {
816 wrklog->indent(indent);
818 *wrklog << gammalib::parformat(
"Time interval", indent);
819 *wrklog << tmin.convert(tref);
821 *wrklog << tmax.convert(tref);
822 *wrklog <<
" s" << std::endl;
826 for (
int ie = 0; ie < ebounds.size(); ++ie) {
829 GEnergy ereco_min = ebounds.emin(ie);
830 GEnergy ereco_max = ebounds.emax(ie);
834 GEnergy etrue_min = ereco_min;
835 GEnergy etrue_max = ereco_max;
836 if (rsp->use_edisp()) {
837 etrue_min = rsp->ebounds(etrue_min).emin();
838 etrue_max = rsp->ebounds(etrue_max).emax();
842 double area =
get_area(obs, etrue_min, etrue_max);
846 *wrklog << gammalib::parformat(
"Photon energy range", indent);
847 *wrklog << etrue_min <<
" - " << etrue_max << std::endl;
848 *wrklog << gammalib::parformat(
"Event energy range", indent);
849 *wrklog << ereco_min <<
" - " << ereco_max << std::endl;
855 if (ebounds.size() > 1) {
857 wrklog->indent(indent);
863 *wrklog << gammalib::parformat(
"Simulation area", indent);
864 *wrklog << area <<
" cm2" << std::endl;
868 int nevents_before = events->size();
872 etrue_min, etrue_max, ereco_min, ereco_max,
874 ran, wrklog, indent, nphotons, nevents);
878 *wrklog << gammalib::parformat(
"MC source events", indent);
879 *wrklog << events->size() - nevents_before;
880 *wrklog <<
" (all source models)";
881 *wrklog << std::endl;
886 if (ebounds.size() > 1) {
888 wrklog->indent(indent);
897 if (events->gti().size() > 1) {
899 wrklog->indent(indent);
910 for (
int i = 0; i < models.size(); ++i) {
911 const GModelSky* model =
dynamic_cast<const GModelSky*
>(models[i]);
915 if (!model->is_valid(obs->instrument(), obs->id())) {
918 *wrklog << gammalib::parformat(
"MC source photons");
919 *wrklog << nphotons[i];
920 if (model->name().length() > 0) {
921 *wrklog <<
" [" << model->name() <<
"]";
923 *wrklog << std::endl;
924 *wrklog << gammalib::parformat(
"MC source events");
925 *wrklog << nevents[i];
926 if (model->name().length() > 0) {
927 *wrklog <<
" [" << model->name() <<
"]";
929 *wrklog << std::endl;
936 #if defined(G_SOURCE_DEBUG)
937 std::cout <<
"ctobssim::simulate_source: out" << std::endl;
970 const GCTAResponseIrf* rsp,
971 GCTAEventList* events,
972 const GModels& models,
975 const GEnergy& etrue_min,
976 const GEnergy& etrue_max,
977 const GEnergy& ereco_min,
978 const GEnergy& ereco_max,
985 std::vector<int>& nphotons,
986 std::vector<int>& nevents)
989 #if defined(G_SOURCE_DEBUG)
990 std::cout <<
"ctobssim::simulate_interval: in";
991 std::cout <<
" Etrue=[" << etrue_min <<
"," << etrue_max <<
"]";
992 std::cout <<
" Ereco=[" << ereco_min <<
"," << ereco_max <<
"]";
993 std::cout << std::endl;
997 for (
int i = 0; i < models.size(); ++i) {
1000 const GModelSky* model =
dynamic_cast<const GModelSky*
>(models[i]);
1003 if (model == NULL) {
1009 if (!model->is_valid(obs->instrument(), obs->id())) {
1017 double flux =
get_model_flux(model, etrue_min, etrue_max, dir, rad,
1019 double rate = flux * area;
1020 double duration = 1800.0;
1023 if (duration < 1.0) {
1026 else if (duration > 180000.0) {
1027 duration = 180000.0;
1038 *wrklog << gammalib::parformat(
"Photon rate", indent);
1039 *wrklog << rate <<
" photons/s";
1040 if (model->name().length() > 0) {
1041 *wrklog <<
" [" << model->name() <<
"]";
1043 *wrklog << std::endl;
1049 std::string mod = (model->name().length() > 0) ?
1050 model->name() :
"Unknown";
1051 std::string msg =
"Photon rate "+gammalib::str(rate)+
1052 " photons/s for model \""+mod+
"\" exceeds "
1053 "maximum allowed photon rate of "+
1055 "Please check the parameters of model "
1056 "\""+mod+
"\" or increase the value of the "
1057 "\"maxrate\" parameter.";
1063 GTime tstart = tmin;
1064 GTime tstop = tstart + duration;
1068 int nphotons_before = nphotons[i];
1069 int nevents_before = nevents[i];
1072 while (tstart < tmax) {
1080 if (logExplicit()) {
1081 if (tmax - tmin > duration) {
1083 wrklog->indent(indent);
1085 *wrklog << gammalib::parformat(
"Time slice", indent);
1086 *wrklog << tstart.convert(obs->gti().reference()) <<
" - ";
1087 *wrklog << tstop.convert(obs->gti().reference()) <<
" s";
1088 if (model->name().length() > 0) {
1089 *wrklog <<
" [" << model->name() <<
"]";
1091 *wrklog << std::endl;
1096 etrue_min, etrue_max, ereco_min, ereco_max,
1098 ran, wrklog, indent,
1099 nphotons[i], nevents[i]);
1103 tstop = tstart + duration;
1106 if (logExplicit()) {
1107 if (tmax - tmin > duration) {
1109 wrklog->indent(indent);
1117 *wrklog << gammalib::parformat(
"MC source photons", indent);
1118 *wrklog << nphotons[i] - nphotons_before;
1119 if (model->name().length() > 0) {
1120 *wrklog <<
" [" << model->name() <<
"]";
1122 *wrklog << std::endl;
1123 *wrklog << gammalib::parformat(
"MC source events", indent);
1124 *wrklog << nevents[i] - nevents_before;
1125 if (model->name().length() > 0) {
1126 *wrklog <<
" [" << model->name() <<
"]";
1128 *wrklog << std::endl;
1135 #if defined(G_SOURCE_DEBUG)
1136 std::cout <<
"ctobssim::simulate_interval: out" << std::endl;
1170 const GCTAResponseIrf* rsp,
1171 GCTAEventList* events,
1172 const GModelSky* model,
1174 const GTime& tstart,
1176 const GEnergy& etrue_min,
1177 const GEnergy& etrue_max,
1178 const GEnergy& ereco_min,
1179 const GEnergy& ereco_max,
1190 #if defined(G_SOURCE_DEBUG)
1191 std::cout <<
"ctobssim::simulate_time_slice: model->mc in" << std::endl;
1195 GPhotons photons = model->mc(area, dir, rad, etrue_min, etrue_max,
1196 tstart, tstop, ran);
1199 #if defined(G_SOURCE_DEBUG)
1200 std::cout <<
"ctobssim::simulate_time_slice: model->mc out" << std::endl;
1204 if (logExplicit()) {
1205 *wrklog << gammalib::parformat(
"MC source photons/slice", indent);
1206 *wrklog << photons.size();
1207 if (model->name().length() > 0) {
1208 *wrklog <<
" [" << model->name() <<
"]";
1210 *wrklog << std::endl;
1214 for (
int i = 0; i < photons.size(); ++i) {
1220 #if defined(G_SOURCE_DEBUG)
1221 std::cout <<
"ctobssim::simulate_time_slice: rsp->mc in" << std::endl;
1226 GCTAEventAtom*
event = rsp->mc(area, photons[i], *obs, ran);
1229 #if defined(G_SOURCE_DEBUG)
1230 std::cout <<
"ctobssim::simulate_time_slice: rsp->mc out" << std::endl;
1235 if (event != NULL) {
1236 if (events->roi().contains(*event) &&
1237 event->energy() >= ereco_min &&
1238 event->energy() <= ereco_max &&
1239 event->time() >= tstart &&
1240 event->time() <= tstop) {
1246 event->mc_id(mc_id);
1249 events->append(*event);
1261 events->has_mc_id(
true);
1283 GEbounds ebins(
m_eslices, ebounds.emin(), ebounds.emax());
1304 const GEnergy& emin,
1305 const GEnergy& emax)
const
1308 const GCTAResponseIrf* rsp =
1309 dynamic_cast<const GCTAResponseIrf*
>(obs->response());
1311 std::string cls = std::string(
typeid(rsp).name());
1312 std::string msg =
"Response of type \""+cls+
"\" is not a CTA IRF "
1313 "response. Please make sure that a CTA IRF "
1314 "response is contained in the CTA observation.";
1315 throw GException::invalid_value(
G_GET_AREA, msg);
1319 const int nbins = 10;
1320 double logE = emin.log10TeV();
1321 double logEbin = (emax.log10TeV() - logE)/
double(nbins-1);
1323 for (
int i = 0; i < nbins; ++i, logE += logEbin) {
1324 double aeff = rsp->aeff()->max(logE, 0.0, 0.0);
1351 const GEnergy& emin,
1352 const GEnergy& emax,
1353 const GSkyDir& centre,
1354 const double& radius,
1359 #if defined(G_SOURCE_DEBUG)
1360 std::cout <<
"ctobssim::get_model_flux: in";
1361 std::cout <<
" Etrue=[" << emin <<
"," << emax <<
"]";
1362 std::cout << std::endl;
1371 double norm = model->spatial()->mc_norm(centre, radius);
1372 bool use_model = (norm > 0.0) ?
true :
false;
1375 *wrklog << gammalib::parformat(
"Use model", indent);
1378 *wrklog << gammalib::parformat(
"Skip model", indent);
1380 if (model->name().length() > 0) {
1381 *wrklog << model->name();
1383 *wrklog << std::endl;
1384 *wrklog << gammalib::parformat(
"Normalization", indent);
1386 if (model->name().length() > 0) {
1387 *wrklog <<
" [" << model->name() <<
"]";
1389 *wrklog << std::endl;
1396 bool free_spectral =
false;
1399 const GModelSpectral* spectral = model->spectral();
1405 GModelSpatialDiffuseCube* cube =
1406 dynamic_cast<GModelSpatialDiffuseCube*
>(model->spatial());
1410 cube->mc_cone(GSkyRegionCircle(centre, radius));
1413 GModelSpectralNodes* nodes =
1414 new GModelSpectralNodes(cube->spectrum());
1415 for (
int i = 0; i < nodes->nodes(); ++i) {
1416 GEnergy energy = nodes->energy(i);
1417 double intensity = nodes->intensity(i);
1418 double value = spectral->eval(energy);
1419 nodes->intensity(i, value * intensity);
1423 free_spectral =
true;
1431 if (nodes->nodes() == 0) {
1439 double flux0 = (use_model) ? spectral->flux(emin, emax) : 0.0;
1440 flux = flux0 * norm;
1444 *wrklog << gammalib::parformat(
"Flux", indent);
1446 if (model->name().length() > 0) {
1447 *wrklog <<
" [" << model->name() <<
"]";
1449 *wrklog <<
" photons/cm2/s" << std::endl;
1450 *wrklog << gammalib::parformat(
"Normalized flux", indent);
1452 if (model->name().length() > 0) {
1453 *wrklog <<
" [" << model->name() <<
"]";
1455 *wrklog <<
" photons/cm2/s" << std::endl;
1459 if (free_spectral)
delete spectral;
1464 #if defined(G_SOURCE_DEBUG)
1465 std::cout <<
"ctobssim::get_model_flux: out" << std::endl;
1487 const GModels& models,
1495 if (wrklog == NULL) {
1500 GCTAEventList* events =
1501 static_cast<GCTAEventList*
>(
const_cast<GEvents*
>(obs->events()));
1504 for (
int i = 0; i < models.size(); ++i) {
1507 const GModelData* model =
1508 dynamic_cast<const GModelData*
>(models[i]);
1512 if (model != NULL &&
1513 model->is_valid(obs->instrument(), obs->id())) {
1516 #if defined(G_BACKGROUND_DEBUG)
1517 std::cout <<
"ctobssim::simulate_background: model->mc in"
1522 GCTAEventList* list =
1523 dynamic_cast<GCTAEventList*
>(model->mc(*obs, ran));
1526 #if defined(G_BACKGROUND_DEBUG)
1527 std::cout <<
"ctobssim::simulate_background: model->mc out"
1535 events->reserve(list->size()+events->size());
1539 int n_outside_roi = 0;
1542 for (
int k = 0; k < list->size(); k++) {
1545 GCTAEventAtom*
event = (*list)[k];
1548 if (events->roi().contains(*event)) {
1558 events->append(*event);
1574 *wrklog << gammalib::parformat(
"MC events outside ROI");
1575 *wrklog << n_outside_roi << std::endl;
1576 *wrklog << gammalib::parformat(
"MC background events");
1577 *wrklog << n_appended << std::endl;
1590 events->has_mc_id(
true);
1612 const GModels& models,
1619 if (wrklog == NULL) {
1624 GCTAEventList* events =
1625 static_cast<GCTAEventList*
>(
const_cast<GEvents*
>(obs->events()));
1628 std::vector<int> ids;
1629 std::vector<std::string> names;
1632 for (
int i = 0; i < events->size(); ++i) {
1635 int id = (*events)[i]->mc_id();
1641 names.push_back(models[
id-1]->name());
1647 bool id_found =
false;
1648 for (
int k = 0; k < ids.size(); ++k) {
1656 names.push_back(models[
id-1]->name());
1663 if ((ids.size() > 0) && (names.size() > 0)) {
1666 events->set_mc_id_names(ids, names);
1670 for (
int k = 0; k < ids.size(); ++k) {
1671 *wrklog << gammalib::parformat(
"MC identifier " +
1672 gammalib::str(ids[k]));
1673 *wrklog << names[k] << std::endl;
1702 GCTAObservation*
obs =
dynamic_cast<GCTAObservation*
>(
m_obs[0]);
1705 log_value(NORMAL,
"Event list file",
m_outevents);
1742 for (
int i = 0; i <
m_obs.size(); ++i) {
1745 GCTAObservation*
obs =
dynamic_cast<GCTAObservation*
>(
m_obs[i]);
1751 if (obs->events()->size() != 0) {
1757 obs->eventfile(outfile);
1760 log_value(NORMAL,
"Event list file", outfile);
1763 obs->save(outfile, clobber());
1803 m_prefix = (*this)[
"prefix"].string();
1808 std::sprintf(buffer,
"%s%6.6d.fits",
m_prefix.c_str(),
1812 outfile = std::string(buffer);
std::string m_outevents
Output events file.
void save(void)
Save the selected event list(s)
int m_startindex
Start index for multiple event lists.
bool m_apply_edisp
Apply energy dispersion?
virtual ~ctobssim(void)
Destructor.
void process(void)
Process the ctobssim tool.
#define G_SIMULATE_SOURCE
const GObservations & obs(void) const
Return observation container.
std::string m_prefix
Prefix for multiple event lists.
void clear(void)
Clear ctobssim tool.
std::vector< GRan > m_rans
Random number generators.
double get_model_flux(const GModelSky *model, const GEnergy &emin, const GEnergy &emax, const GSkyDir ¢re, const double &radius, const int &indent, GLog *wrklog)
Determine sky model flux.
void save_fits(void)
Save event list in FITS format.
double get_area(GCTAObservation *obs, const GEnergy &emin, const GEnergy &emax) const
Get simulation area (cm^2)
std::string outfile(const int &index)
Return output filename.
Observation simulator tool definition.
GModels m_models
Optionally provided models.
Observation simulator tool.
bool m_save_and_dispose
Save and dispose immediately.
void get_parameters(void)
Get application parameters.
int m_event_id
Event identifier.
Base class for observation tools.
void publish(const std::string &name="")
Publish event lists.
void simulate_source(GCTAObservation *obs, const GModels &models, GRan &ran, GLog *wrklog=NULL)
Simulate source events from photon list.
void free_members(void)
Delete class members.
void simulate_interval(GCTAObservation *obs, const GCTAResponseIrf *rsp, GCTAEventList *events, const GModels &models, const GTime &tmin, const GTime &tmax, const GEnergy &etrue_min, const GEnergy &etrue_max, const GEnergy &ereco_min, const GEnergy &ereco_max, const GSkyDir &dir, const double &rad, const double &area, GRan &ran, GLog *wrklog, int &indent, std::vector< int > &nphotons, std::vector< int > &nevents)
Simulate source events for a time and energy interval.
int m_seed
Random number generator seed.
void set_obs_bounds()
Set observation boundaries for CTA observations.
void init_members(void)
Initialise class members.
ctobssim & operator=(const ctobssim &app)
Assignment operator.
const double g_roi_margin
Simulation radius margin (degrees)
double m_max_rate
Maximum photon rate.
#define G_SIMULATE_INTERVAL
ctobservation & operator=(const ctobservation &app)
Assignment operator.
void copy_members(const ctobssim &app)
Copy class members.
ctobssim(void)
Void constructor.
void simulate_time_slice(GCTAObservation *obs, const GCTAResponseIrf *rsp, GCTAEventList *events, const GModelSky *model, const int &mc_id, const GTime &tstart, const GTime &tstop, const GEnergy &etrue_min, const GEnergy &etrue_max, const GEnergy &ereco_min, const GEnergy &ereco_max, const GSkyDir &dir, const double &rad, const double &area, GRan &ran, GLog *wrklog, int &indent, int &nphotons, int &nevents)
Simulate source events for a time slice.
void models(const GModels &models)
Set models for simulation.
int m_eslices
Number of energy slices.
void simulate_background(GCTAObservation *obs, const GModels &models, GRan &ran, GLog *wrklog=NULL)
Simulate background events from model.
GObservations m_obs
Observation container.
void save_xml(void)
Save event list(s) in XML format.
int m_max_photons
Maximum number of photons/slice.
void free_members(void)
Delete class members.
void set_mc_id_names(GCTAObservation *obs, const GModels &models, GLog *wrklog=NULL)
Set correspondance between Monte Carlo identifier and model names.
void init_members(void)
Initialise class members.