Lili Yang / ctools

Commit fb5380a82d3a9471c5f2fc092f3f305e0a7c9404

Authored by Jürgen Knödlseder
1 parent e0624ca2

Correct simulations for energy dispersion. 

The true energy range was used for event selection instead of the reconstructed energy range. This was a bug introduced in the recent refactoring. Now fixed.
Inline Side-by-side
Showing 2 changed files with 48 additions and 33 deletions
src/ctobssim/ctobssim.cpp
  View file @fb5380a
... ... @@ -39,8 +39,8 @@
39 39 "GModels&, GRan&, GLog*)"
40 40 #define G_SIMULATE_INTERVAL "ctobssim::simulate_interval(GCTAObservation*"\
41 41 ", GCTAEventList*, GCTAResponseIrf*, GModels&, GTime&"\
42   - ", GTime&, GEnergy&, GEnergy&, GSkyDir&, double&"\
43   - ", double&, GRan&, GLog*, int&)"
  42 + ", GTime&, GEnergy&, GEnergy&, GEnergy&, GEnergy&"\
  43 + ", GSkyDir&, double&, double&, GRan&, GLog*, int&)"
44 44 #define G_GET_AREA "ctobssim::get_area(GCTAObservation* obs, GEnergy&, "\
45 45 "GEnergy&)"
46 46  
... ... @@ -711,7 +711,8 @@ void ctobssim::simulate_source(GCTAObservation* obs,
711 711 GCTAEventList* events = static_cast<GCTAEventList*>(obs->events());
712 712  
713 713 // Get CTA response
714   - const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(obs->response());
  714 + const GCTAResponseIrf* rsp =
  715 + dynamic_cast<const GCTAResponseIrf*>(obs->response());
715 716 if (rsp == NULL) {
716 717 std::string cls = std::string(typeid(obs->response()).name());
717 718 std::string msg = "Response of type \""+cls+"\" is not a CTA "
... ... @@ -767,28 +768,28 @@ void ctobssim::simulate_source(GCTAObservation* obs,
767 768 // Loop over all energy boundaries
768 769 for (int ie = 0; ie < ebounds.size(); ++ie) {
769 770  
770   - // Set energy interval
771   - GEnergy emin = ebounds.emin(ie);
772   - GEnergy emax = ebounds.emax(ie);
  771 + // Set reconstructed energy interval
  772 + GEnergy ereco_min = ebounds.emin(ie);
  773 + GEnergy ereco_max = ebounds.emax(ie);
773 774  
774 775 // Set true photon energy limits for simulation. If the
775 776 // observation has energy dispersion then add a margin.
776   - GEnergy e_true_min = emin;
777   - GEnergy e_true_max = emax;
  777 + GEnergy etrue_min = ereco_min;
  778 + GEnergy etrue_max = ereco_max;
778 779 if (rsp->use_edisp()) {
779   - e_true_min = rsp->ebounds(e_true_min).emin();
780   - e_true_max = rsp->ebounds(e_true_max).emax();
  780 + etrue_min = rsp->ebounds(etrue_min).emin();
  781 + etrue_max = rsp->ebounds(etrue_max).emax();
781 782 }
782 783  
783 784 // Determine simulation area
784   - double area = get_area(obs, emin, emax);
  785 + double area = get_area(obs, ereco_min, ereco_max);
785 786  
786 787 // Log energy range and simulation area
787 788 if (logNormal()) {
788 789 *wrklog << gammalib::parformat("Photon energy range", indent);
789   - *wrklog << e_true_min << " - " << e_true_max << std::endl;
  790 + *wrklog << etrue_min << " - " << etrue_max << std::endl;
790 791 *wrklog << gammalib::parformat("Event energy range", indent);
791   - *wrklog << emin << " - " << emax << std::endl;
  792 + *wrklog << ereco_min << " - " << ereco_max << std::endl;
792 793 }
793 794  
794 795 // Increment indentation if there are several energy
... ... @@ -811,7 +812,8 @@ void ctobssim::simulate_source(GCTAObservation* obs,
811 812  
812 813 // Simulate events for this time and energy interval
813 814 simulate_interval(obs, rsp, events, models, tmin, tmax,
814   - e_true_min, e_true_max, dir, rad, area,
  815 + etrue_min, etrue_max, ereco_min, ereco_max,
  816 + dir, rad, area,
815 817 ran, wrklog, indent, nphotons, nevents);
816 818  
817 819 // Log simulation results
... ... @@ -887,8 +889,10 @@ void ctobssim::simulate_source(GCTAObservation* obs,
887 889 * @param[in] models Model list.
888 890 * @param[in] tmin Start time.
889 891 * @param[in] tmax Stop time.
890   - * @param[in] emin Minimum energy.
891   - * @param[in] emax Maximum energy.
  892 + * @param[in] etrue_min Minimum true energy.
  893 + * @param[in] etrue_max Maximum true energy.
  894 + * @param[in] ereco_min Minimum reconstructed energy.
  895 + * @param[in] ereco_max Maximum reconstructed energy.
892 896 * @param[in] dir Simulation cone centre.
893 897 * @param[in] rad Simulation cone radius (degrees).
894 898 * @param[in] area Simulation area (cm^2).
... ... @@ -906,8 +910,10 @@ void ctobssim::simulate_interval(GCTAObservation* obs,
906 910 const GModels& models,
907 911 const GTime& tmin,
908 912 const GTime& tmax,
909   - const GEnergy& emin,
910   - const GEnergy& emax,
  913 + const GEnergy& etrue_min,
  914 + const GEnergy& etrue_max,
  915 + const GEnergy& ereco_min,
  916 + const GEnergy& ereco_max,
911 917 const GSkyDir& dir,
912 918 const double& rad,
913 919 const double& area,
... ... @@ -938,7 +944,7 @@ void ctobssim::simulate_interval(GCTAObservation* obs,
938 944 // simulated photons to m_max_photons. The photon rate is estimated
939 945 // from the model flux and used to set the duration of the time
940 946 // slice.
941   - double flux = get_model_flux(model, emin, emax, dir, rad,
  947 + double flux = get_model_flux(model, etrue_min, etrue_max, dir, rad,
942 948 indent, wrklog);
943 949 double rate = flux * area;
944 950 double duration = 1800.0; // default: 1800 sec
... ... @@ -1017,7 +1023,8 @@ void ctobssim::simulate_interval(GCTAObservation* obs,
1017 1023  
1018 1024 // Simulate time slice
1019 1025 simulate_time_slice(obs, rsp, events, model, tstart, tstop,
1020   - emin, emax, dir, rad, area,
  1026 + etrue_min, etrue_max, ereco_min, ereco_max,
  1027 + dir, rad, area,
1021 1028 ran, wrklog, indent,
1022 1029 nphotons[i], nevents[i]);
1023 1030  
... ... @@ -1067,8 +1074,10 @@ void ctobssim::simulate_interval(GCTAObservation* obs,
1067 1074 * @param[in] models Model list.
1068 1075 * @param[in] tmin Start time.
1069 1076 * @param[in] tmax Stop time.
1070   - * @param[in] emin Minimum energy.
1071   - * @param[in] emax Maximum energy.
  1077 + * @param[in] etrue_min Minimum true energy.
  1078 + * @param[in] etrue_max Maximum true energy.
  1079 + * @param[in] ereco_min Minimum reconstructed energy.
  1080 + * @param[in] ereco_max Maximum reconstructed energy.
1072 1081 * @param[in] dir Simulation cone centre.
1073 1082 * @param[in] rad Simulation cone radius (degrees).
1074 1083 * @param[in] area Simulation area (cm^2).
... ... @@ -1086,8 +1095,10 @@ void ctobssim::simulate_time_slice(GCTAObservation* obs,
1086 1095 const GModelSky* model,
1087 1096 const GTime& tstart,
1088 1097 const GTime& tstop,
1089   - const GEnergy& emin,
1090   - const GEnergy& emax,
  1098 + const GEnergy& etrue_min,
  1099 + const GEnergy& etrue_max,
  1100 + const GEnergy& ereco_min,
  1101 + const GEnergy& ereco_max,
1091 1102 const GSkyDir& dir,
1092 1103 const double& rad,
1093 1104 const double& area,
... ... @@ -1098,8 +1109,8 @@ void ctobssim::simulate_time_slice(GCTAObservation* obs,
1098 1109 int& nevents)
1099 1110 {
1100 1111 // Get photons
1101   - GPhotons photons = model->mc(area, dir, rad, emin, emax, tstart, tstop,
1102   - ran);
  1112 + GPhotons photons = model->mc(area, dir, rad, etrue_min, etrue_max,
  1113 + tstart, tstop, ran);
1103 1114  
1104 1115 // Dump number of simulated photons
1105 1116 if (logExplicit()) {
... ... @@ -1122,11 +1133,11 @@ void ctobssim::simulate_time_slice(GCTAObservation* obs,
1122 1133 GCTAEventAtom* event = rsp->mc(area, photons[i], *obs, ran);
1123 1134  
1124 1135 // Use event only if it exists and if it falls within ROI, the
1125   - // energy interval and the time slice
  1136 + // reconstructed energy interval and the time slice
1126 1137 if (event != NULL) {
1127 1138 if (events->roi().contains(*event) &&
1128   - event->energy() >= emin &&
1129   - event->energy() <= emax &&
  1139 + event->energy() >= ereco_min &&
  1140 + event->energy() <= ereco_max &&
1130 1141 event->time() >= tstart &&
1131 1142 event->time() <= tstop) {
1132 1143 event->event_id(m_event_id);
... ...
src/ctobssim/ctobssim.hpp
  View file @fb5380a
... ... @@ -85,8 +85,10 @@ protected:
85 85 const GModels& models,
86 86 const GTime& tmin,
87 87 const GTime& tmax,
88   - const GEnergy& emin,
89   - const GEnergy& emax,
  88 + const GEnergy& etrue_min,
  89 + const GEnergy& etrue_max,
  90 + const GEnergy& ereco_min,
  91 + const GEnergy& ereco_max,
90 92 const GSkyDir& dir,
91 93 const double& rad,
92 94 const double& area,
... ... @@ -101,8 +103,10 @@ protected:
101 103 const GModelSky* model,
102 104 const GTime& tstart,
103 105 const GTime& tstop,
104   - const GEnergy& emin,
105   - const GEnergy& emax,
  106 + const GEnergy& etrue_min,
  107 + const GEnergy& etrue_max,
  108 + const GEnergy& ereco_min,
  109 + const GEnergy& ereco_max,
106 110 const GSkyDir& dir,
107 111 const double& rad,
108 112 const double& area,
... ...