Commit f8b2c3603a9d508174ce7cd148c3c165f585c7e4

Authored by Jürgen Knödlseder
2 parents 1dbc3811 36d7e2d5

Merge branch 'nodes-mc' into integration

Showing 1 changed file with 40 additions and 8 deletions
src/ctobssim/ctobssim.cpp
1 /*************************************************************************** 1 /***************************************************************************
2 * ctobssim - Observation simulator tool * 2 * ctobssim - Observation simulator tool *
3 * ----------------------------------------------------------------------- * 3 * ----------------------------------------------------------------------- *
4 - * copyright (C) 2011-2016 by Juergen Knoedlseder * 4 + * copyright (C) 2011-2017 by Juergen Knoedlseder *
5 * ----------------------------------------------------------------------- * 5 * ----------------------------------------------------------------------- *
6 * * 6 * *
7 * This program is free software: you can redistribute it and/or modify * 7 * This program is free software: you can redistribute it and/or modify *
@@ -50,6 +50,8 @@ @@ -50,6 +50,8 @@
50 const double g_roi_margin = 0.5; //!< Simulation radius margin (degrees) 50 const double g_roi_margin = 0.5; //!< Simulation radius margin (degrees)
51 51
52 /* __ Debug definitions __________________________________________________ */ 52 /* __ Debug definitions __________________________________________________ */
  53 +//#define G_SOURCE_DEBUG
  54 +//#define G_BACKGROUND_DEBUG
53 55
54 /* __ Coding definitions _________________________________________________ */ 56 /* __ Coding definitions _________________________________________________ */
55 57
@@ -609,7 +611,7 @@ void ctobssim::get_parameters(void) @@ -609,7 +611,7 @@ void ctobssim::get_parameters(void)
609 611
610 // Get observation container 612 // Get observation container
611 m_obs = get_observations(); 613 m_obs = get_observations();
612 - 614 +
613 } 615 }
614 616
615 // ... otherwise add response information and energy boundaries in case 617 // ... otherwise add response information and energy boundaries in case
@@ -1079,7 +1081,7 @@ void ctobssim::simulate_interval(GCTAObservation* obs, @@ -1079,7 +1081,7 @@ void ctobssim::simulate_interval(GCTAObservation* obs,
1079 *wrklog << std::endl; 1081 *wrklog << std::endl;
1080 1082
1081 } 1083 }
1082 - 1084 +
1083 } // endfor: looped over models 1085 } // endfor: looped over models
1084 1086
1085 // Return 1087 // Return
@@ -1130,10 +1132,20 @@ void ctobssim::simulate_time_slice(GCTAObservation* obs, @@ -1130,10 +1132,20 @@ void ctobssim::simulate_time_slice(GCTAObservation* obs,
1130 int& nphotons, 1132 int& nphotons,
1131 int& nevents) 1133 int& nevents)
1132 { 1134 {
  1135 + // Debug code: signal that we step into the model MC method
  1136 + #if defined(G_SOURCE_DEBUG)
  1137 + std::cout << "ctobssim::simulate_time_slice: model->mc in" << std::endl;
  1138 + #endif
  1139 +
1133 // Get photons 1140 // Get photons
1134 GPhotons photons = model->mc(area, dir, rad, etrue_min, etrue_max, 1141 GPhotons photons = model->mc(area, dir, rad, etrue_min, etrue_max,
1135 tstart, tstop, ran); 1142 tstart, tstop, ran);
1136 1143
  1144 + // Debug code: signal that we came back
  1145 + #if defined(G_SOURCE_DEBUG)
  1146 + std::cout << "ctobssim::simulate_time_slice: model->mc out" << std::endl;
  1147 + #endif
  1148 +
1137 // Dump number of simulated photons 1149 // Dump number of simulated photons
1138 if (logExplicit()) { 1150 if (logExplicit()) {
1139 *wrklog << gammalib::parformat("MC source photons/slice", indent); 1151 *wrklog << gammalib::parformat("MC source photons/slice", indent);
@@ -1150,10 +1162,20 @@ void ctobssim::simulate_time_slice(GCTAObservation* obs, @@ -1150,10 +1162,20 @@ void ctobssim::simulate_time_slice(GCTAObservation* obs,
1150 // Increment photon counter 1162 // Increment photon counter
1151 nphotons++; 1163 nphotons++;
1152 1164
  1165 + // Debug code: signal that we step into the response MC method
  1166 + #if defined(G_SOURCE_DEBUG)
  1167 + std::cout << "ctobssim::simulate_time_slice: rsp->mc in" << std::endl;
  1168 + #endif
  1169 +
1153 // Simulate event. Note that this method includes the deadtime 1170 // Simulate event. Note that this method includes the deadtime
1154 // correction. 1171 // correction.
1155 GCTAEventAtom* event = rsp->mc(area, photons[i], *obs, ran); 1172 GCTAEventAtom* event = rsp->mc(area, photons[i], *obs, ran);
1156 1173
  1174 + // Debug code: signal that we came back
  1175 + #if defined(G_SOURCE_DEBUG)
  1176 + std::cout << "ctobssim::simulate_time_slice: rsp->mc out" << std::endl;
  1177 + #endif
  1178 +
1157 // Use event only if it exists and if it falls within ROI, the 1179 // Use event only if it exists and if it falls within ROI, the
1158 // reconstructed energy interval and the time slice 1180 // reconstructed energy interval and the time slice
1159 if (event != NULL) { 1181 if (event != NULL) {
@@ -1193,7 +1215,7 @@ GEbounds ctobssim::get_ebounds(const GEbounds&amp; ebounds) const @@ -1193,7 +1215,7 @@ GEbounds ctobssim::get_ebounds(const GEbounds&amp; ebounds) const
1193 { 1215 {
1194 // Set energy bins 1216 // Set energy bins
1195 GEbounds ebins(m_eslices, ebounds.emin(), ebounds.emax()); 1217 GEbounds ebins(m_eslices, ebounds.emin(), ebounds.emax());
1196 - 1218 +
1197 // Return energy bins 1219 // Return energy bins
1198 return (ebins); 1220 return (ebins);
1199 } 1221 }
@@ -1322,7 +1344,7 @@ double ctobssim::get_model_flux(const GModelSky* model, @@ -1322,7 +1344,7 @@ double ctobssim::get_model_flux(const GModelSky* model,
1322 spectral = nodes; 1344 spectral = nodes;
1323 1345
1324 } // endif: spatial model was a diffuse cube 1346 } // endif: spatial model was a diffuse cube
1325 - 1347 +
1326 // Compute flux within [emin, emax] in model from spectral 1348 // Compute flux within [emin, emax] in model from spectral
1327 // component (units: ph/cm2/s) 1349 // component (units: ph/cm2/s)
1328 double flux0 = spectral->flux(emin, emax); 1350 double flux0 = spectral->flux(emin, emax);
@@ -1396,11 +1418,21 @@ void ctobssim::simulate_background(GCTAObservation* obs, @@ -1396,11 +1418,21 @@ void ctobssim::simulate_background(GCTAObservation* obs,
1396 if (model != NULL && 1418 if (model != NULL &&
1397 model->is_valid(obs->instrument(), obs->id())) { 1419 model->is_valid(obs->instrument(), obs->id())) {
1398 1420
  1421 + // Debug code: signal that we step into the response MC method
  1422 + #if defined(G_BACKGROUND_DEBUG)
  1423 + std::cout << "ctobssim::simulate_background: model->mc in" << std::endl;
  1424 + #endif
  1425 +
1399 // Get simulated CTA event list. Note that this method 1426 // Get simulated CTA event list. Note that this method
1400 // includes the deadtime correction. 1427 // includes the deadtime correction.
1401 GCTAEventList* list = 1428 GCTAEventList* list =
1402 dynamic_cast<GCTAEventList*>(model->mc(*obs, ran)); 1429 dynamic_cast<GCTAEventList*>(model->mc(*obs, ran));
1403 1430
  1431 + // Debug code: signal that we came back
  1432 + #if defined(G_BACKGROUND_DEBUG)
  1433 + std::cout << "ctobssim::simulate_background: model->mc out" << std::endl;
  1434 + #endif
  1435 +
1404 // Continue only if we got a CTA event list 1436 // Continue only if we got a CTA event list
1405 if (list != NULL) { 1437 if (list != NULL) {
1406 1438
@@ -1410,7 +1442,7 @@ void ctobssim::simulate_background(GCTAObservation* obs, @@ -1410,7 +1442,7 @@ void ctobssim::simulate_background(GCTAObservation* obs,
1410 // Initialise statistics 1442 // Initialise statistics
1411 int n_appended = 0; 1443 int n_appended = 0;
1412 int n_outside_roi = 0; 1444 int n_outside_roi = 0;
1413 - 1445 +
1414 // Append events 1446 // Append events
1415 for (int k = 0; k < list->size(); k++) { 1447 for (int k = 0; k < list->size(); k++) {
1416 1448
@@ -1589,10 +1621,10 @@ std::string ctobssim::outfile(const int&amp; index) @@ -1589,10 +1621,10 @@ std::string ctobssim::outfile(const int&amp; index)
1589 1621
1590 // ... otherwise use the outfile parameter 1622 // ... otherwise use the outfile parameter
1591 else { 1623 else {
1592 - 1624 +
1593 // Get output event list file name 1625 // Get output event list file name
1594 m_outevents = (*this)["outevents"].filename(); 1626 m_outevents = (*this)["outevents"].filename();
1595 - 1627 +
1596 // Set output filename 1628 // Set output filename
1597 outfile = std::string(m_outevents); 1629 outfile = std::string(m_outevents);
1598 } 1630 }