Lili Yang / ctools

Commit f1ce5f252cfe20e6b3ab286362f4784b98a43831

Authored by Jürgen Knödlseder
1 parent 58ac71c8

Introduce energy slices in ctobssim to accelerate computations, in particular at low energies.

Inline Side-by-side
Showing 5 changed files with 517 additions and 247 deletions
ChangeLog
  View file @f1ce5f2
1   -2016-02-11 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
  1 +2016-02-17 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
2 2  
3 3 * Version 1.1.0 released
4 4 ========================
5 5  
  6 + Subdivide simulation energy range in ctobssim (#1673)
6 7 Delete existing events before simulation in ctobssim (#1663)
7 8 Add csiactcopy script (#1661)
8 9 ctselect allows usage of non-default events and GTI extensions (#1598)
... ...
... ... @@ -2,7 +2,7 @@ New Features and Important Changes in ctools 1.1.0
2 2  
3 3 Author(s) : Juergen Knoedlseder
4 4  
5   -11 February 2016
  5 +17 February 2016
6 6  
7 7  
8 8 1. Introduction
... ... @@ -28,6 +28,16 @@ None
28 28  
29 29 5. ctobssim - CTA observation simulation
30 30 ----------------------------------------
  31 +The simulation energy range is subdivided internally to avoid having a too
  32 +large difference between simulation area and effective area. Before, there
  33 +was a simulation area of 19634954 m^2 that was applied to all energies, which
  34 +led to the simulation of a huge number of low energy photons that were thrown
  35 +away later due to the small effective area. Now, the energy range is subdivided
  36 +(controlled using the "eslices" parameter), and for each energy range, the
  37 +simulation area is computed from the maximum effective area in that energy
  38 +range multiplied by a security factor of 2. This leads to a much more efficient
  39 +simulation of low-energy events (#1673).
  40 +
31 41 Any existing events in an event list are now removed before starting the
32 42 simulation to enable providing an existing event file on input from which
33 43 only the observation attributes are extracted (for example GTI, ROI, energy
... ...
src/ctobssim/ctobssim.cpp
  View file @f1ce5f2
... ... @@ -29,6 +29,7 @@
29 29 #include <config.h>
30 30 #endif
31 31 #include <cstdio>
  32 +#include <typeinfo>
32 33 #include "ctobssim.hpp"
33 34 #include "GTools.hpp"
34 35  
... ... @@ -36,6 +37,12 @@
36 37 #define G_GET_PARAMETERS "ctobssim::get_parameters()"
37 38 #define G_SIMULATE_SOURCE "ctobssim::simulate_source(GCTAObservation*, "\
38 39 "GModels&, GRan&, GLog*)"
  40 +#define G_SIMULATE_INTERVAL "ctobssim::simulate_interval(GCTAObservation*"\
  41 + ", GCTAEventList*, GCTAResponseIrf*, GModels&, GTime&"\
  42 + ", GTime&, GEnergy&, GEnergy&, GSkyDir&, double&"\
  43 + ", double&, GRan&, GLog*, int&)"
  44 +#define G_GET_AREA "ctobssim::get_area(GCTAObservation* obs, GEnergy&, "\
  45 + "GEnergy&)"
39 46  
40 47 /* __ Constants __________________________________________________________ */
41 48 const double g_roi_margin = 0.5; //!< Simulation radius margin (degrees)
... ... @@ -382,7 +389,7 @@ void ctobssim::run(void)
382 389 simulate_background(&obs_clone, models, m_rans[i], &wrklog);
383 390  
384 391 // Dump simulation results
385   - if (logNormal()) {
  392 + if (logTerse()) {
386 393 wrklog << gammalib::parformat("MC events");
387 394 wrklog << obs_clone.events()->size() - events_before;
388 395 wrklog << " (all models)";
... ... @@ -512,6 +519,7 @@ void ctobssim::init_members(void)
512 519 m_outevents.clear();
513 520 m_prefix.clear();
514 521 m_seed = 1;
  522 + m_eslices = 10;
515 523 m_apply_edisp = false;
516 524 m_max_rate = 1.0e6;
517 525  
... ... @@ -521,8 +529,7 @@ void ctobssim::init_members(void)
521 529 m_save_and_dispose = false;
522 530  
523 531 // Set fixed parameters
524   - m_area = 19634954.0 * 1.0e4; //!< pi*(2500^2) m^2
525   - m_max_photons = 1000000; //!< Maximum number of photons / time slice
  532 + m_max_photons = 1000000; //!< Maximum number of photons / time slice
526 533  
527 534 // Initialise first event identifier
528 535 m_event_id = 1;
... ... @@ -543,11 +550,11 @@ void ctobssim::copy_members(const ctobssim&amp; app)
543 550 m_outevents = app.m_outevents;
544 551 m_prefix = app.m_prefix;
545 552 m_seed = app.m_seed;
  553 + m_eslices = app.m_eslices;
546 554 m_apply_edisp = app.m_apply_edisp;
547 555 m_max_rate = app.m_max_rate;
548 556  
549 557 // Copy protected members
550   - m_area = app.m_area;
551 558 m_max_photons = app.m_max_photons;
552 559 m_rans = app.m_rans;
553 560 m_obs = app.m_obs;
... ... @@ -608,6 +615,7 @@ void ctobssim::get_parameters(void)
608 615  
609 616 // Get other parameters
610 617 m_seed = (*this)["seed"].integer();
  618 + m_eslices = (*this)["eslices"].integer();
611 619 m_apply_edisp = (*this)["edisp"].boolean();
612 620 m_max_rate = (*this)["maxrate"].real();
613 621  
... ... @@ -673,17 +681,29 @@ void ctobssim::get_parameters(void)
673 681 *
674 682 * @param[in] obs Pointer on CTA observation.
675 683 * @param[in] models Model list.
676   - * @param[in] ran Random number generator.
  684 + * @param[in,out] ran Random number generator.
677 685 * @param[in] wrklog Pointer to logger.
678 686 *
679   - * Simulate source events from a photon list for a given CTA observation.
680   - * The events are stored in as event list in the observation.
  687 + * Simulate source events from a photon list for a given CTA observation and
  688 + * all source models. The events are stored in form of an event list in the
  689 + * observation.
681 690 *
682   - * This method does nothing if the observation pointer is NULL. It also
683   - * verifies if the observation has a valid pointing and response.
  691 + * The method will loop over all Good Time Intervals (GTI) to perform the
  692 + * simulations. Within a given GTI the requested energy range is split into
  693 + * a number of energy slices so that the simulation area can be adapted to
  694 + * the effective area of the instrument (this avoids simulating a large
  695 + * number of sources photons at low energies while accepting only few of
  696 + * them do to the small effective area). In case that the flux of a source
  697 + * is large, the simulation may be done within time slices so that the
  698 + * memory requirements won't get too large.
  699 + *
  700 + * This method does nothing if the observation pointer is NULL. It verifies
  701 + * if the observation has a CTA IRF response.
684 702 ***************************************************************************/
685   -void ctobssim::simulate_source(GCTAObservation* obs, const GModels& models,
686   - GRan& ran, GLog* wrklog)
  703 +void ctobssim::simulate_source(GCTAObservation* obs,
  704 + const GModels& models,
  705 + GRan& ran,
  706 + GLog* wrklog)
687 707 {
688 708 // Continue only if observation pointer is valid
689 709 if (obs != NULL) {
... ... @@ -693,36 +713,42 @@ void ctobssim::simulate_source(GCTAObservation* obs, const GModels&amp; models,
693 713 wrklog = &log;
694 714 }
695 715  
  716 + // Get pointer on event list
  717 + GCTAEventList* events = static_cast<GCTAEventList*>(obs->events());
  718 +
696 719 // Get CTA response
697   - const GCTAResponseIrf* rsp =
698   - dynamic_cast<const GCTAResponseIrf*>(obs->response());
  720 + const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(obs->response());
699 721 if (rsp == NULL) {
700   - std::string msg = "Response is not an IRF response.\n" +
701   - obs->response()->print();
  722 + std::string cls = std::string(typeid(obs->response()).name());
  723 + std::string msg = "Response of type \""+cls+"\" is not a CTA "
  724 + "IRF response. Please make sure that a CTA "
  725 + "IRF response is contained in the CTA "
  726 + "observation.";
702 727 throw GException::invalid_value(G_SIMULATE_SOURCE, msg);
703 728 }
704 729  
705   - // Get pointer on event list (circumvent const correctness)
706   - GCTAEventList* events =
707   - static_cast<GCTAEventList*>(const_cast<GEvents*>(obs->events()));
708   -
709   - // Extract simulation region.
  730 + // Extract simulation region from event list ROI
710 731 GSkyDir dir = events->roi().centre().dir();
711 732 double rad = events->roi().radius() + g_roi_margin;
712 733  
713   - // Dump simulation cone information
714   - if (logNormal()) {
715   - *wrklog << gammalib::parformat("Simulation area");
716   - *wrklog << m_area << " cm2" << std::endl;
  734 + // Determine energy boundaries for simulation
  735 + GEbounds ebounds = get_ebounds(events->ebounds());
  736 +
  737 + // Log simulation cone information
  738 + if (logTerse()) {
717 739 *wrklog << gammalib::parformat("Simulation cone");
718 740 *wrklog << "RA=" << dir.ra_deg() << " deg";
719 741 *wrklog << ", Dec=" << dir.dec_deg() << " deg";
720   - *wrklog << ", r=" << rad << " deg" << std::endl;
  742 + *wrklog << ", radius=" << rad << " deg" << std::endl;
721 743 }
722 744  
723 745 // Initialise indentation for logging
724 746 int indent = 0;
725 747  
  748 + // Initialise photon and event counters
  749 + std::vector<int> nphotons(models.size(),0);
  750 + std::vector<int> nevents(models.size(),0);
  751 +
726 752 // Loop over all Good Time Intervals
727 753 for (int it = 0; it < events->gti().size(); ++it) {
728 754  
... ... @@ -730,7 +756,8 @@ void ctobssim::simulate_source(GCTAObservation* obs, const GModels&amp; models,
730 756 GTime tmin = events->gti().tstart(it);
731 757 GTime tmax = events->gti().tstop(it);
732 758  
733   - // Dump time interval
  759 + // Log time interval. Increment indentation if there are
  760 + // several Good Time Intervals.
734 761 if (logNormal()) {
735 762 if (events->gti().size() > 1) {
736 763 indent++;
... ... @@ -744,14 +771,14 @@ void ctobssim::simulate_source(GCTAObservation* obs, const GModels&amp; models,
744 771 }
745 772  
746 773 // Loop over all energy boundaries
747   - for (int ie = 0; ie < events->ebounds().size(); ++ie) {
  774 + for (int ie = 0; ie < ebounds.size(); ++ie) {
748 775  
749   - // Extract energy boundaries
750   - GEnergy emin = events->ebounds().emin(ie);
751   - GEnergy emax = events->ebounds().emax(ie);
  776 + // Set energy interval
  777 + GEnergy emin = ebounds.emin(ie);
  778 + GEnergy emax = ebounds.emax(ie);
752 779  
753   - // Set true photon energy limits for simulation. If observation
754   - // has energy dispersion then add margin
  780 + // Set true photon energy limits for simulation. If the
  781 + // observation has energy dispersion then add a margin.
755 782 GEnergy e_true_min = emin;
756 783 GEnergy e_true_max = emax;
757 784 if (rsp->use_edisp()) {
... ... @@ -759,210 +786,51 @@ void ctobssim::simulate_source(GCTAObservation* obs, const GModels&amp; models,
759 786 e_true_max = rsp->ebounds(e_true_max).emax();
760 787 }
761 788  
762   - // Dump energy range
  789 + // Determine simulation area
  790 + double area = get_area(obs, emin, emax);
  791 +
  792 + // Log energy range and simulation area
763 793 if (logNormal()) {
764   - if (events->ebounds().size() > 1) {
765   - indent++;
766   - wrklog->indent(indent);
767   - }
768 794 *wrklog << gammalib::parformat("Photon energy range", indent);
769 795 *wrklog << e_true_min << " - " << e_true_max << std::endl;
770 796 *wrklog << gammalib::parformat("Event energy range", indent);
771 797 *wrklog << emin << " - " << emax << std::endl;
772 798 }
773 799  
774   - // Loop over all sky models
775   - for (int i = 0; i < models.size(); ++i) {
776   -
777   - // Get sky model (NULL if not a sky model)
778   - const GModelSky* model =
779   - dynamic_cast<const GModelSky*>(models[i]);
780   -
781   - // If we have a sky model that applies to the present
782   - // observation then simulate events
783   - if (model != NULL &&
784   - model->is_valid(obs->instrument(), obs->id())) {
785   -
786   - // Determine duration of a time slice by limiting the
787   - // number of simulated photons to m_max_photons.
788   - // The photon rate is estimated from the model flux
789   - // and used to set the duration of the time slice.
790   - double flux = get_model_flux(model, emin, emax,
791   - dir, rad,
792   - indent, wrklog);
793   - double rate = flux * m_area;
794   - double duration = 1800.0; // default: 1800 sec
795   - if (rate > 0.0) {
796   - duration = m_max_photons / rate;
797   - if (duration < 1.0) { // not <1 sec
798   - duration = 1.0;
799   - }
800   - else if (duration > 180000.0) { // not >50 hr
801   - duration = 180000.0;
802   - }
803   - }
804   - GTime tslice(duration, "sec");
805   -
806   - // Skip model if photon rate is 0
807   - if (rate <= 0.0) {
808   - continue;
809   - }
810   -
811   - // Dump photon rate
812   - if (logNormal()) {
813   - *wrklog << gammalib::parformat("Photon rate", indent);
814   - *wrklog << rate << " photons/s";
815   - if (model->name().length() > 0) {
816   - *wrklog << " [" << model->name() << "]";
817   - }
818   - *wrklog << std::endl;
819   - }
  800 + // Increment indentation if there are several energy
  801 + // boundaries
  802 + if (logNormal()) {
  803 + if (ebounds.size() > 1) {
  804 + indent++;
  805 + wrklog->indent(indent);
  806 + }
  807 + }
820 808  
821   - // If photon rate exceeds the maximum photon rate
822   - // then throw an exception
823   - if (rate > m_max_rate) {
824   - std::string modnam = (model->name().length() > 0) ?
825   - model->name() : "Unknown";
826   - std::string msg = "Photon rate "+
827   - gammalib::str(rate)+
828   - " photons/s for model \""+
829   - modnam+"\" exceeds maximum"
830   - " allowed photon rate of "+
831   - gammalib::str(m_max_rate)+
832   - " photons/s. Please check"
833   - " the model parameters for"
834   - " model \""+modnam+"\" or"
835   - " increase the value of the"
836   - " hidden \"maxrate\""
837   - " parameter.";
838   - throw GException::invalid_value(G_SIMULATE_SOURCE, msg);
839   - }
  809 + // Log simulation area
  810 + if (logNormal()) {
  811 + *wrklog << gammalib::parformat("Simulation area", indent);
  812 + *wrklog << area << " cm2" << std::endl;
  813 + }
840 814  
841   - // To reduce memory requirements we split long time
842   - // intervals into several slices.
843   - GTime tstart = tmin;
844   - GTime tstop = tstart + tslice;
845   -
846   - // Initialise cumulative photon counters
847   - int nphotons = 0;
848   -
849   - // Loop over time slices
850   - while (tstart < tmax) {
851   -
852   - // Make sure that tstop <= tmax
853   - if (tstop > tmax) {
854   - tstop = tmax;
855   - }
856   -
857   - // Dump time slice
858   - if (logExplicit()) {
859   - if (tmax - tmin > tslice) {
860   - indent++;
861   - wrklog->indent(indent);
862   - }
863   - *wrklog << gammalib::parformat("Time slice", indent);
864   - *wrklog << tstart.convert(m_cta_ref);
865   - *wrklog << " - ";
866   - *wrklog << tstop.convert(m_cta_ref);
867   - *wrklog << " s";
868   - if (model->name().length() > 0) {
869   - *wrklog << " [" << model->name() << "]";
870   - }
871   - *wrklog << std::endl;
872   - }
873   -
874   - // Get photons
875   - GPhotons photons = model->mc(m_area, dir, rad,
876   - e_true_min, e_true_max,
877   - tstart, tstop, ran);
878   -
879   - // Dump number of simulated photons
880   - if (logExplicit()) {
881   - *wrklog << gammalib::parformat("MC source photons/slice", indent);
882   - *wrklog << photons.size();
883   - if (model->name().length() > 0) {
884   - *wrklog << " [" << model->name() << "]";
885   - }
886   - *wrklog << std::endl;
887   - }
888   -
889   - // Simulate events from photons
890   - for (int i = 0; i < photons.size(); ++i) {
891   -
892   - // Increment photon counter
893   - nphotons++;
894   -
895   - // Simulate event. Note that this method
896   - // includes the deadtime correction.
897   - GCTAEventAtom* event = rsp->mc(m_area,
898   - photons[i],
899   - *obs,
900   - ran);
901   -
902   - if (event != NULL) {
903   -
904   - // Use event only if it falls within ROI
905   - // energy boundary and time slice
906   - if (events->roi().contains(*event) &&
907   - event->energy() >= emin &&
908   - event->energy() <= emax &&
909   - event->time() >= tstart &&
910   - event->time() <= tstop) {
911   - event->event_id(m_event_id);
912   - events->append(*event);
913   - m_event_id++;
914   - }
915   - delete event;
916   - }
917   -
918   - } // endfor: looped over events
919   -
920   - // Go to next time slice
921   - tstart = tstop;
922   - tstop = tstart + tslice;
923   -
924   - // Reset indentation
925   - if (logExplicit()) {
926   - if (tmax - tmin > tslice) {
927   - indent--;
928   - wrklog->indent(indent);
929   - }
930   - }
931   -
932   - } // endwhile: looped over time slices
933   -
934   - // Dump simulation results
935   - if (logNormal()) {
936   - *wrklog << gammalib::parformat("MC source photons", indent);
937   - *wrklog << nphotons;
938   - if (model->name().length() > 0) {
939   - *wrklog << " [" << model->name() << "]";
940   - }
941   - *wrklog << std::endl;
942   - *wrklog << gammalib::parformat("MC source events", indent);
943   - *wrklog << events->size();
944   - if (model->name().length() > 0) {
945   - *wrklog << " [" << model->name() << "]";
946   - }
947   - *wrklog << std::endl;
  815 + // Save state of event counter before doing the simulation
  816 + int nevents_before = events->size();
948 817  
949   - }
950   -
951   - } // endif: model was a sky model
  818 + // Simulate events for this time and energy interval
  819 + simulate_interval(obs, rsp, events, models, tmin, tmax,
  820 + e_true_min, e_true_max, dir, rad, area,
  821 + ran, wrklog, indent, nphotons, nevents);
952 822  
953   - } // endfor: looped over models
954   -
955   - // Dump simulation results
  823 + // Log simulation results
956 824 if (logNormal()) {
957 825 *wrklog << gammalib::parformat("MC source events", indent);
958   - *wrklog << events->size();
  826 + *wrklog << events->size() - nevents_before;
959 827 *wrklog << " (all source models)";
960 828 *wrklog << std::endl;
961 829 }
962 830  
963   - // Reset indentation
  831 + // Reset indentation if there were several energy boundaries
964 832 if (logNormal()) {
965   - if (events->ebounds().size() > 1) {
  833 + if (ebounds.size() > 1) {
966 834 indent--;
967 835 wrklog->indent(indent);
968 836  
... ... @@ -971,7 +839,7 @@ void ctobssim::simulate_source(GCTAObservation* obs, const GModels&amp; models,
971 839  
972 840 } // endfor: looped over all energy boundaries
973 841  
974   - // Reset indentation
  842 + // Reset indentation if there were several Good Time Intervals
975 843 if (logNormal()) {
976 844 if (events->gti().size() > 1) {
977 845 indent--;
... ... @@ -984,6 +852,31 @@ void ctobssim::simulate_source(GCTAObservation* obs, const GModels&amp; models,
984 852 // Reset indentation
985 853 wrklog->indent(0);
986 854  
  855 + // Log simulation summary
  856 + if (logTerse()) {
  857 + for (int i = 0; i < models.size(); ++i) {
  858 + const GModelSky* model = dynamic_cast<const GModelSky*>(models[i]);
  859 + if (model == NULL) {
  860 + continue;
  861 + }
  862 + if (!model->is_valid(obs->instrument(), obs->id())) {
  863 + continue;
  864 + }
  865 + *wrklog << gammalib::parformat("MC source photons");
  866 + *wrklog << nphotons[i];
  867 + if (model->name().length() > 0) {
  868 + *wrklog << " [" << model->name() << "]";
  869 + }
  870 + *wrklog << std::endl;
  871 + *wrklog << gammalib::parformat("MC source events");
  872 + *wrklog << nevents[i];
  873 + if (model->name().length() > 0) {
  874 + *wrklog << " [" << model->name() << "]";
  875 + }
  876 + *wrklog << std::endl;
  877 + }
  878 + }
  879 +
987 880 } // endif: observation pointer was valid
988 881  
989 882 // Return
... ... @@ -992,6 +885,335 @@ void ctobssim::simulate_source(GCTAObservation* obs, const GModels&amp; models,
992 885  
993 886  
994 887 /***********************************************************************//**
  888 + * @brief Simulate source events for a time and energy interval
  889 + *
  890 + * @param[in] obs Pointer on CTA observation.
  891 + * @param[in] rsp Pointer on CTA IRF response.
  892 + * @param[in,out] events Pointer on CTA event list.
  893 + * @param[in] models Model list.
  894 + * @param[in] tmin Start time.
  895 + * @param[in] tmax Stop time.
  896 + * @param[in] emin Minimum energy.
  897 + * @param[in] emax Maximum energy.
  898 + * @param[in] dir Simulation cone centre.
  899 + * @param[in] rad Simulation cone radius (degrees).
  900 + * @param[in] area Simulation area (cm^2).
  901 + * @param[in,out] ran Random number generator.
  902 + * @param[in] wrklog Pointer to logger.
  903 + * @param[in,out] indent Logger indent.
  904 + * @param[in,out] nphotons Number of photons for all models.
  905 + * @param[in,out] nevents Number of events for all models.
  906 + *
  907 + * Simulate source events for a time and energy interval.
  908 + ***************************************************************************/
  909 +void ctobssim::simulate_interval(GCTAObservation* obs,
  910 + const GCTAResponseIrf* rsp,
  911 + GCTAEventList* events,
  912 + const GModels& models,
  913 + const GTime& tmin,
  914 + const GTime& tmax,
  915 + const GEnergy& emin,
  916 + const GEnergy& emax,
  917 + const GSkyDir& dir,
  918 + const double& rad,
  919 + const double& area,
  920 + GRan& ran,
  921 + GLog* wrklog,
  922 + int& indent,
  923 + std::vector<int>& nphotons,
  924 + std::vector<int>& nevents)
  925 +{
  926 + // Loop over all models
  927 + for (int i = 0; i < models.size(); ++i) {
  928 +
  929 + // Get sky model (NULL if not a sky model)
  930 + const GModelSky* model = dynamic_cast<const GModelSky*>(models[i]);
  931 +
  932 + // If the model is not a sky model then skip the model
  933 + if (model == NULL) {
  934 + continue;
  935 + }
  936 +
  937 + // If the model does not apply to the instrument and observation
  938 + // identifier then skip the model
  939 + if (!model->is_valid(obs->instrument(), obs->id())) {
  940 + continue;
  941 + }
  942 +
  943 + // Determine duration of a time slice by limiting the number of
  944 + // simulated photons to m_max_photons. The photon rate is estimated
  945 + // from the model flux and used to set the duration of the time
  946 + // slice.
  947 + double flux = get_model_flux(model, emin, emax, dir, rad,
  948 + indent, wrklog);
  949 + double rate = flux * area;
  950 + double duration = 1800.0; // default: 1800 sec
  951 + if (rate > 0.0) {
  952 + duration = m_max_photons / rate;
  953 + if (duration < 1.0) { // not <1 sec
  954 + duration = 1.0;
  955 + }
  956 + else if (duration > 180000.0) { // not >50 hr
  957 + duration = 180000.0;
  958 + }
  959 + }
  960 + GTime tslice(duration, "sec");
  961 +
  962 + // Skip model if photon rate is 0
  963 + if (rate <= 0.0) {
  964 + continue;
  965 + }
  966 +
  967 + // Log photon rate
  968 + if (logNormal()) {
  969 + *wrklog << gammalib::parformat("Photon rate", indent);
  970 + *wrklog << rate << " photons/s";
  971 + if (model->name().length() > 0) {
  972 + *wrklog << " [" << model->name() << "]";
  973 + }
  974 + *wrklog << std::endl;
  975 + }
  976 +
  977 + // If photon rate exceeds the maximum photon rate that is allowed
  978 + // then throw an exception
  979 + if (rate > m_max_rate) {
  980 + std::string mod = (model->name().length() > 0) ?
  981 + model->name() : "Unknown";
  982 + std::string msg = "Photon rate "+gammalib::str(rate)+
  983 + " photons/s for model \""+mod+"\" exceeds "
  984 + "maximum allowed photon rate of "+
  985 + gammalib::str(m_max_rate)+" photons/s. "
  986 + "Please check the model parameters for "
  987 + "model \""+mod+"\" or increase the value "
  988 + "of the \"maxrate\" parameter.";
  989 + throw GException::invalid_value(G_SIMULATE_INTERVAL, msg);
  990 + }
  991 +
  992 + // To reduce memory requirements we split long time intervals into
  993 + // several time slices
  994 + GTime tstart = tmin;
  995 + GTime tstop = tstart + tslice;
  996 +
  997 + // Save state of photon and event counters before doing the
  998 + // simulation
  999 + int nphotons_before = nphotons[i];
  1000 + int nevents_before = nevents[i];
  1001 +
  1002 + // Loop over time slices
  1003 + while (tstart < tmax) {
  1004 +
  1005 + // Make sure that tstop <= tmax
  1006 + if (tstop > tmax) {
  1007 + tstop = tmax;
  1008 + }
  1009 +
  1010 + // Log time slice
  1011 + if (logExplicit()) {
  1012 + if (tmax - tmin > tslice) {
  1013 + indent++;
  1014 + wrklog->indent(indent);
  1015 + }
  1016 + *wrklog << gammalib::parformat("Time slice", indent);
  1017 + *wrklog << tstart.convert(m_cta_ref) << " - ";
  1018 + *wrklog << tstop.convert(m_cta_ref) << " s";
  1019 + if (model->name().length() > 0) {
  1020 + *wrklog << " [" << model->name() << "]";
  1021 + }
  1022 + *wrklog << std::endl;
  1023 + }
  1024 +
  1025 + // Simulate time slice
  1026 + simulate_time_slice(obs, rsp, events, model, tstart, tstop,
  1027 + emin, emax, dir, rad, area,
  1028 + ran, wrklog, indent,
  1029 + nphotons[i], nevents[i]);
  1030 +
  1031 + // Go to next time slice
  1032 + tstart = tstop;
  1033 + tstop = tstart + tslice;
  1034 +
  1035 + // Reset indentation
  1036 + if (logExplicit()) {
  1037 + if (tmax - tmin > tslice) {
  1038 + indent--;
  1039 + wrklog->indent(indent);
  1040 + }
  1041 + }
  1042 +
  1043 + } // endwhile: looped over time slices
  1044 +
  1045 + // Log simulation results
  1046 + if (logNormal()) {
  1047 + *wrklog << gammalib::parformat("MC source photons", indent);
  1048 + *wrklog << nphotons[i] - nphotons_before;
  1049 + if (model->name().length() > 0) {
  1050 + *wrklog << " [" << model->name() << "]";
  1051 + }
  1052 + *wrklog << std::endl;
  1053 + *wrklog << gammalib::parformat("MC source events", indent);
  1054 + *wrklog << nevents[i] - nevents_before;
  1055 + if (model->name().length() > 0) {
  1056 + *wrklog << " [" << model->name() << "]";
  1057 + }
  1058 + *wrklog << std::endl;
  1059 +
  1060 + }
  1061 +
  1062 + } // endfor: looped over models
  1063 +
  1064 + // Return
  1065 + return;
  1066 +}
  1067 +
  1068 +
  1069 +/***********************************************************************//**
  1070 + * @brief Simulate source events for a time slice
  1071 + *
  1072 + * @param[in] obs Pointer on CTA observation.
  1073 + * @param[in,out] events Pointer on CTA event list.
  1074 + * @param[in] models Model list.
  1075 + * @param[in] tmin Start time.
  1076 + * @param[in] tmax Stop time.
  1077 + * @param[in] emin Minimum energy.
  1078 + * @param[in] emax Maximum energy.
  1079 + * @param[in] dir Simulation cone centre.
  1080 + * @param[in] rad Simulation cone radius (degrees).
  1081 + * @param[in] area Simulation area (cm^2).
  1082 + * @param[in,out] ran Random number generator.
  1083 + * @param[in] wrklog Pointer to logger.
  1084 + * @param[in,out] indent Logger indent.
  1085 + * @param[in,out] nphotons Number of photons.
  1086 + * @param[in,out] nevents Number of events.
  1087 + *
  1088 + * Simulate source events for a time slice.
  1089 + ***************************************************************************/
  1090 +void ctobssim::simulate_time_slice(GCTAObservation* obs,
  1091 + const GCTAResponseIrf* rsp,
  1092 + GCTAEventList* events,
  1093 + const GModelSky* model,
  1094 + const GTime& tstart,
  1095 + const GTime& tstop,
  1096 + const GEnergy& emin,
  1097 + const GEnergy& emax,
  1098 + const GSkyDir& dir,
  1099 + const double& rad,
  1100 + const double& area,
  1101 + GRan& ran,
  1102 + GLog* wrklog,
  1103 + int& indent,
  1104 + int& nphotons,
  1105 + int& nevents)
  1106 +{
  1107 + // Get photons
  1108 + GPhotons photons = model->mc(area, dir, rad, emin, emax, tstart, tstop,
  1109 + ran);
  1110 +
  1111 + // Dump number of simulated photons
  1112 + if (logExplicit()) {
  1113 + *wrklog << gammalib::parformat("MC source photons/slice", indent);
  1114 + *wrklog << photons.size();
  1115 + if (model->name().length() > 0) {
  1116 + *wrklog << " [" << model->name() << "]";
  1117 + }
  1118 + *wrklog << std::endl;
  1119 + }
  1120 +
  1121 + // Simulate events from photons
  1122 + for (int i = 0; i < photons.size(); ++i) {
  1123 +
  1124 + // Increment photon counter
  1125 + nphotons++;
  1126 +
  1127 + // Simulate event. Note that this method includes the deadtime
  1128 + // correction.
  1129 + GCTAEventAtom* event = rsp->mc(area, photons[i], *obs, ran);
  1130 +
  1131 + // Use event only if it exists and if it falls within ROI, the
  1132 + // energy interval and the time slice
  1133 + if (event != NULL) {
  1134 + if (events->roi().contains(*event) &&
  1135 + event->energy() >= emin &&
  1136 + event->energy() <= emax &&
  1137 + event->time() >= tstart &&
  1138 + event->time() <= tstop) {
  1139 + event->event_id(m_event_id);
  1140 + events->append(*event);
  1141 + m_event_id++;
  1142 + nevents++;
  1143 + }
  1144 + delete event;
  1145 + }
  1146 +
  1147 + } // endfor: looped over events
  1148 +
  1149 + // Return
  1150 + return;
  1151 +}
  1152 +
  1153 +
  1154 +/***********************************************************************//**
  1155 + * @brief Get energy boundaries
  1156 + *
  1157 + * @param[in] ebounds Energy boundaries of events.
  1158 + * @return Energy boundaries for simulation.
  1159 + *
  1160 + * Get the energy boundaries for the source simulation.
  1161 + *
  1162 + * @todo Need to take care of the (rare) case that there are distinct
  1163 + * energy boundaries in the event list. In that case we need to
  1164 + * create mixed energy boundaries.
  1165 + ***************************************************************************/
  1166 +GEbounds ctobssim::get_ebounds(const GEbounds& ebounds) const
  1167 +{
  1168 + // Set energy bins
  1169 + GEbounds ebins(m_eslices, ebounds.emin(), ebounds.emax());
  1170 +
  1171 + // Return energy bins
  1172 + return (ebins);
  1173 +}
  1174 +
  1175 +
  1176 +/***********************************************************************//**
  1177 + * @brief Get simulation area (cm^2)
  1178 + *
  1179 + * @param[in] obs Pointer on CTA observation.
  1180 + * @param[in] emin Minimum energy.
  1181 + * @param[in] emax Maximum energy.
  1182 + * @return Simulation area (cm^2).
  1183 + *
  1184 + * Get the simulation area for an energy interval in units of cm^2.
  1185 + ***************************************************************************/
  1186 +double ctobssim::get_area(GCTAObservation* obs,
  1187 + const GEnergy& emin,
  1188 + const GEnergy& emax) const
  1189 +{
  1190 + // Get CTA response
  1191 + const GCTAResponseIrf* rsp =
  1192 + dynamic_cast<const GCTAResponseIrf*>(obs->response());
  1193 + if (rsp == NULL) {
  1194 + std::string cls = std::string(typeid(rsp).name());
  1195 + std::string msg = "Response of type \""+cls+"\" is not a CTA IRF "
  1196 + "response. Please make sure that a CTA IRF "
  1197 + "response is contained in the CTA observation.";
  1198 + throw GException::invalid_value(G_GET_AREA, msg);
  1199 + }
  1200 +
  1201 + // Compute effective area at minimum and maximum energy
  1202 + double area_emin = rsp->aeff()->max(emin.log10TeV(), 0.0, 0.0);
  1203 + double area_emax = rsp->aeff()->max(emax.log10TeV(), 0.0, 0.0);
  1204 +
  1205 + // Use maximum of both
  1206 + double area = (area_emin > area_emax) ? area_emin : area_emax;
  1207 +
  1208 + // Multiply by security factor
  1209 + area *= 2.0;
  1210 +
  1211 + // Return simulation area
  1212 + return area;
  1213 +}
  1214 +
  1215 +
  1216 +/***********************************************************************//**
995 1217 * @brief Determine sky model flux
996 1218 *
997 1219 * @param[in] model Sky model.
... ... @@ -1000,7 +1222,7 @@ void ctobssim::simulate_source(GCTAObservation* obs, const GModels&amp; models,
1000 1222 * @param[in] centre Centre of region for photon rate determination.
1001 1223 * @param[in] radius Radius of region for photon rate determination (degrees).
1002 1224 * @param[in] indent Indent for logging.
1003   - * @param[in,out] wrklog Pointer to logger (default: NULL).
  1225 + * @param[in,out] wrklog Pointer to logger.
1004 1226 * @return Model flux (photons/cm2/sec).
1005 1227 ***************************************************************************/
1006 1228 double ctobssim::get_model_flux(const GModelSky* model,
... ...
src/ctobssim/ctobssim.hpp
  View file @f1ce5f2
1 1 /***************************************************************************
2 2 * ctobssim - Observation simulator tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2016 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -34,7 +34,7 @@
34 34  
35 35 /* __Definitions _________________________________________________________ */
36 36 #define CTOBSSIM_NAME "ctobssim"
37   -#define CTOBSSIM_VERSION "1.0.0"
  37 +#define CTOBSSIM_VERSION "1.1.0"
38 38  
39 39  
40 40 /***********************************************************************//**
... ... @@ -71,38 +71,74 @@ public:
71 71  
72 72 protected:
73 73 // Protected methods
74   - void init_members(void);
75   - void copy_members(const ctobssim& app);
76   - void free_members(void);
77   - void get_parameters(void);
78   - void simulate_source(GCTAObservation* obs,
79   - const GModels& models,
80   - GRan& ran,
81   - GLog* wrklog = NULL);
82   - double get_model_flux(const GModelSky* model,
83   - const GEnergy& emin,
84   - const GEnergy& emax,
85   - const GSkyDir& centre,
86   - const double& radius,
87   - const int& indent,
88   - GLog* wrklog = NULL);
89   - void simulate_background(GCTAObservation* obs,
90   - const GModels& models,
91   - GRan& ran,
92   - GLog* wrklog = NULL);
93   - void save_fits(void);
94   - void save_xml(void);
  74 + void init_members(void);
  75 + void copy_members(const ctobssim& app);
  76 + void free_members(void);
  77 + void get_parameters(void);
  78 + void simulate_source(GCTAObservation* obs,
  79 + const GModels& models,
  80 + GRan& ran,
  81 + GLog* wrklog = NULL);
  82 + void simulate_interval(GCTAObservation* obs,
  83 + const GCTAResponseIrf* rsp,
  84 + GCTAEventList* events,
  85 + const GModels& models,
  86 + const GTime& tmin,
  87 + const GTime& tmax,
  88 + const GEnergy& emin,
  89 + const GEnergy& emax,
  90 + const GSkyDir& dir,
  91 + const double& rad,
  92 + const double& area,
  93 + GRan& ran,
  94 + GLog* wrklog,
  95 + int& indent,
  96 + std::vector<int>& nphotons,
  97 + std::vector<int>& nevents);
  98 + void simulate_time_slice(GCTAObservation* obs,
  99 + const GCTAResponseIrf* rsp,
  100 + GCTAEventList* events,
  101 + const GModelSky* model,
  102 + const GTime& tstart,
  103 + const GTime& tstop,
  104 + const GEnergy& emin,
  105 + const GEnergy& emax,
  106 + const GSkyDir& dir,
  107 + const double& rad,
  108 + const double& area,
  109 + GRan& ran,
  110 + GLog* wrklog,
  111 + int& indent,
  112 + int& nphotons,
  113 + int& nevents);
  114 + GEbounds get_ebounds(const GEbounds& ebounds) const;
  115 + double get_area(GCTAObservation* obs,
  116 + const GEnergy& emin,
  117 + const GEnergy& emax) const;
  118 + double get_model_flux(const GModelSky* model,
  119 + const GEnergy& emin,
  120 + const GEnergy& emax,
  121 + const GSkyDir& centre,
  122 + const double& radius,
  123 + const int& indent,
  124 + GLog* wrklog);
  125 + void simulate_background(GCTAObservation* obs,
  126 + const GModels& models,
  127 + GRan& ran,
  128 + GLog* wrklog = NULL);
  129 + void save_fits(void);
  130 + void save_xml(void);
95 131  
96 132 // User parameters
97 133 std::string m_outevents; //!< Output events file
98 134 std::string m_prefix; //!< Prefix for multiple event lists
99 135 int m_seed; //!< Random number generator seed
  136 + int m_eslices; //!< Number of energy slices
100 137 bool m_apply_edisp; //!< Apply energy dispersion?
101 138  
102 139 // Protected members
103 140 mutable GObservations m_obs; //!< Observation container
104 141 mutable bool m_save_and_dispose; //!< Save and dispose immediately
105   - double m_area; //!< Surface area for simulation (cm2)
106 142 int m_max_photons; //!< Maximum number of photons/slice
107 143 double m_max_rate; //!< Maximum photon rate
108 144 std::vector<GRan> m_rans; //!< Random number generators
... ...
src/ctobssim/ctobssim.par
  View file @f1ce5f2
... ... @@ -5,7 +5,7 @@
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # File: ctobssim.par #
8   -# Version: 1.0.0 #
  8 +# Version: 1.1.0 #
9 9 # #
10 10 # Author: Jurgen Knodlseder #
11 11 # Jurgen.Knodlseder@irap.omp.eu #
... ... @@ -38,6 +38,7 @@ emin, r, a, 0.1,,, &quot;Lower energy limit (TeV)&quot;
38 38 emax, r, a, 100.0,,, "Upper energy limit (TeV)"
39 39 deadc, r, h, 0.95,0,1, "Average deadtime correction factor"
40 40 maxrate, r, h, 1.0e6,,, "Maximum photon rate"
  41 +eslices, i, h, 10,1,100,"Number of energy slices"
41 42  
42 43 #
43 44 # Standard parameters
... ...