Commit 36d7e2d589866ea286dfe52fd41785218fa12988
1 parent
9d4325e6
Add debug code to ctobssim
Showing
1 changed file
with
40 additions
and
8 deletions
src/ctobssim/ctobssim.cpp
| 1 | 1 | /*************************************************************************** |
| 2 | 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 | 7 | * This program is free software: you can redistribute it and/or modify * |
| ... | ... | @@ -50,6 +50,8 @@ |
| 50 | 50 | const double g_roi_margin = 0.5; //!< Simulation radius margin (degrees) |
| 51 | 51 | |
| 52 | 52 | /* __ Debug definitions __________________________________________________ */ |
| 53 | +//#define G_SOURCE_DEBUG | |
| 54 | +//#define G_BACKGROUND_DEBUG | |
| 53 | 55 | |
| 54 | 56 | /* __ Coding definitions _________________________________________________ */ |
| 55 | 57 | |
| ... | ... | @@ -609,7 +611,7 @@ void ctobssim::get_parameters(void) |
| 609 | 611 | |
| 610 | 612 | // Get observation container |
| 611 | 613 | m_obs = get_observations(); |
| 612 | - | |
| 614 | + | |
| 613 | 615 | } |
| 614 | 616 | |
| 615 | 617 | // ... otherwise add response information and energy boundaries in case |
| ... | ... | @@ -1079,7 +1081,7 @@ void ctobssim::simulate_interval(GCTAObservation* obs, |
| 1079 | 1081 | *wrklog << std::endl; |
| 1080 | 1082 | |
| 1081 | 1083 | } |
| 1082 | - | |
| 1084 | + | |
| 1083 | 1085 | } // endfor: looped over models |
| 1084 | 1086 | |
| 1085 | 1087 | // Return |
| ... | ... | @@ -1130,10 +1132,20 @@ void ctobssim::simulate_time_slice(GCTAObservation* obs, |
| 1130 | 1132 | int& nphotons, |
| 1131 | 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 | 1140 | // Get photons |
| 1134 | 1141 | GPhotons photons = model->mc(area, dir, rad, etrue_min, etrue_max, |
| 1135 | 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 | 1149 | // Dump number of simulated photons |
| 1138 | 1150 | if (logExplicit()) { |
| 1139 | 1151 | *wrklog << gammalib::parformat("MC source photons/slice", indent); |
| ... | ... | @@ -1150,10 +1162,20 @@ void ctobssim::simulate_time_slice(GCTAObservation* obs, |
| 1150 | 1162 | // Increment photon counter |
| 1151 | 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 | 1170 | // Simulate event. Note that this method includes the deadtime |
| 1154 | 1171 | // correction. |
| 1155 | 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 | 1179 | // Use event only if it exists and if it falls within ROI, the |
| 1158 | 1180 | // reconstructed energy interval and the time slice |
| 1159 | 1181 | if (event != NULL) { |
| ... | ... | @@ -1193,7 +1215,7 @@ GEbounds ctobssim::get_ebounds(const GEbounds& ebounds) const |
| 1193 | 1215 | { |
| 1194 | 1216 | // Set energy bins |
| 1195 | 1217 | GEbounds ebins(m_eslices, ebounds.emin(), ebounds.emax()); |
| 1196 | - | |
| 1218 | + | |
| 1197 | 1219 | // Return energy bins |
| 1198 | 1220 | return (ebins); |
| 1199 | 1221 | } |
| ... | ... | @@ -1322,7 +1344,7 @@ double ctobssim::get_model_flux(const GModelSky* model, |
| 1322 | 1344 | spectral = nodes; |
| 1323 | 1345 | |
| 1324 | 1346 | } // endif: spatial model was a diffuse cube |
| 1325 | - | |
| 1347 | + | |
| 1326 | 1348 | // Compute flux within [emin, emax] in model from spectral |
| 1327 | 1349 | // component (units: ph/cm2/s) |
| 1328 | 1350 | double flux0 = spectral->flux(emin, emax); |
| ... | ... | @@ -1396,11 +1418,21 @@ void ctobssim::simulate_background(GCTAObservation* obs, |
| 1396 | 1418 | if (model != NULL && |
| 1397 | 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 | 1426 | // Get simulated CTA event list. Note that this method |
| 1400 | 1427 | // includes the deadtime correction. |
| 1401 | 1428 | GCTAEventList* list = |
| 1402 | 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 | 1436 | // Continue only if we got a CTA event list |
| 1405 | 1437 | if (list != NULL) { |
| 1406 | 1438 | |
| ... | ... | @@ -1410,7 +1442,7 @@ void ctobssim::simulate_background(GCTAObservation* obs, |
| 1410 | 1442 | // Initialise statistics |
| 1411 | 1443 | int n_appended = 0; |
| 1412 | 1444 | int n_outside_roi = 0; |
| 1413 | - | |
| 1445 | + | |
| 1414 | 1446 | // Append events |
| 1415 | 1447 | for (int k = 0; k < list->size(); k++) { |
| 1416 | 1448 | |
| ... | ... | @@ -1589,10 +1621,10 @@ std::string ctobssim::outfile(const int& index) |
| 1589 | 1621 | |
| 1590 | 1622 | // ... otherwise use the outfile parameter |
| 1591 | 1623 | else { |
| 1592 | - | |
| 1624 | + | |
| 1593 | 1625 | // Get output event list file name |
| 1594 | 1626 | m_outevents = (*this)["outevents"].filename(); |
| 1595 | - | |
| 1627 | + | |
| 1596 | 1628 | // Set output filename |
| 1597 | 1629 | outfile = std::string(m_outevents); |
| 1598 | 1630 | } | ... | ... |