Commit 44f157c1b0dcbeaf668a34094b0b44d921e16b08
1 parent
437ae542
support simulation of multiple CTA observations; update code tu support standard…
…ized Python interface
Showing
1 changed file
with
421 additions
and
160 deletions
src/ctobssim/ctobssim.cpp
| ... | ... | @@ -24,6 +24,9 @@ |
| 24 | 24 | #include "ctobssim.hpp" |
| 25 | 25 | #include "GTools.hpp" |
| 26 | 26 | |
| 27 | +/* __ Method name definitions ____________________________________________ */ | |
| 28 | +#define G_SIMULATE_SOURCE "ctobssim::simulate_source(GCTAObservation*, GPhotons&)" | |
| 29 | + | |
| 27 | 30 | /* __ Debug definitions __________________________________________________ */ |
| 28 | 31 | |
| 29 | 32 | /* __ Coding definitions _________________________________________________ */ |
| ... | ... | @@ -52,6 +55,28 @@ ctobssim::ctobssim(void) : GApplication(CTOBSSIM_NAME, CTOBSSIM_VERSION) |
| 52 | 55 | |
| 53 | 56 | |
| 54 | 57 | /***********************************************************************//** |
| 58 | + * @brief Observations constructor | |
| 59 | + * | |
| 60 | + * This method creates an instance of the class by copying an existing | |
| 61 | + * observations container. | |
| 62 | + ***************************************************************************/ | |
| 63 | +ctobssim::ctobssim(GObservations obs) : GApplication(CTOBSSIM_NAME, CTOBSSIM_VERSION) | |
| 64 | +{ | |
| 65 | + // Initialise members | |
| 66 | + init_members(); | |
| 67 | + | |
| 68 | + // Set observations | |
| 69 | + m_obs = obs; | |
| 70 | + | |
| 71 | + // Write header into logger | |
| 72 | + log_header(); | |
| 73 | + | |
| 74 | + // Return | |
| 75 | + return; | |
| 76 | +} | |
| 77 | + | |
| 78 | + | |
| 79 | +/***********************************************************************//** | |
| 55 | 80 | * @brief Command line constructor |
| 56 | 81 | * |
| 57 | 82 | * @param[in] argc Number of arguments in command line. |
| ... | ... | @@ -97,6 +122,24 @@ ctobssim::~ctobssim(void) |
| 97 | 122 | ==========================================================================*/ |
| 98 | 123 | |
| 99 | 124 | /***********************************************************************//** |
| 125 | + * @brief Clear instance | |
| 126 | + ***************************************************************************/ | |
| 127 | +void ctobssim::clear(void) | |
| 128 | +{ | |
| 129 | + // Free members | |
| 130 | + free_members(); | |
| 131 | + this->GApplication::free_members(); | |
| 132 | + | |
| 133 | + // Initialise members | |
| 134 | + this->GApplication::init_members(); | |
| 135 | + init_members(); | |
| 136 | + | |
| 137 | + // Return | |
| 138 | + return; | |
| 139 | +} | |
| 140 | + | |
| 141 | + | |
| 142 | +/***********************************************************************//** | |
| 100 | 143 | * @brief Simulate event data |
| 101 | 144 | ***************************************************************************/ |
| 102 | 145 | void ctobssim::run(void) |
| ... | ... | @@ -108,48 +151,116 @@ void ctobssim::run(void) |
| 108 | 151 | // Get parameters |
| 109 | 152 | get_parameters(); |
| 110 | 153 | |
| 111 | - // Write parameters into logger | |
| 154 | + // Write input parameters into logger | |
| 112 | 155 | if (logTerse()) { |
| 113 | 156 | log_parameters(); |
| 114 | 157 | log << std::endl; |
| 115 | - log << m_pnt << std::endl; | |
| 116 | - log << m_obs << std::endl; | |
| 117 | - log << m_models << std::endl; | |
| 118 | 158 | } |
| 119 | 159 | |
| 120 | - // Write header | |
| 160 | + // Write observation(s) into logger | |
| 121 | 161 | if (logTerse()) { |
| 122 | 162 | log << std::endl; |
| 123 | - log.header1("Simulate photons"); | |
| 163 | + if (m_obs.size() > 1) | |
| 164 | + log.header1("Observations"); | |
| 165 | + else | |
| 166 | + log.header1("Observation"); | |
| 167 | + log << m_obs << std::endl; | |
| 124 | 168 | } |
| 125 | 169 | |
| 126 | - // Simulate photons | |
| 127 | - GPhotons photons = simulate_photons(); | |
| 128 | - | |
| 129 | 170 | // Write header |
| 130 | 171 | if (logTerse()) { |
| 131 | 172 | log << std::endl; |
| 132 | - log.header1("Simulate events"); | |
| 173 | + if (m_obs.size() > 1) | |
| 174 | + log.header1("Simulate observations"); | |
| 175 | + else | |
| 176 | + log.header1("Simulate observation"); | |
| 133 | 177 | } |
| 134 | 178 | |
| 135 | - // Simulate events | |
| 136 | - GCTAEventList* events = simulate_events(photons); | |
| 179 | + // Loop over all observation in the container | |
| 180 | + for (int i = 0; i < m_obs.size(); ++i) { | |
| 137 | 181 | |
| 182 | + // Get CTA observation | |
| 183 | + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs(i)); | |
| 184 | + | |
| 185 | + // Continue only if observation is a CTA observation | |
| 186 | + if (obs != NULL) { | |
| 187 | + | |
| 188 | + // Write header for observation | |
| 189 | + if (logTerse()) | |
| 190 | + log.header3("Observation "+obs->obsname()); | |
| 191 | + | |
| 192 | + // Allocate CTA event list and attach it to the observation. | |
| 193 | + // Note that the events() method removes any existing event | |
| 194 | + // list before attaching a new event list. | |
| 195 | + GCTAEventList* events = new GCTAEventList; | |
| 196 | + obs->events(events); | |
| 197 | + | |
| 198 | + // Write header for photon simulation | |
| 199 | + //if (logTerse()) | |
| 200 | + // log.header3("Simulate photons"); | |
| 201 | + | |
| 202 | + // Simulate photons | |
| 203 | + GPhotons photons = simulate_photons(obs, m_obs.models()); | |
| 204 | + | |
| 205 | + // Write header for sky event simulation | |
| 206 | + //if (logTerse()) | |
| 207 | + // log.header3("Simulate sky events"); | |
| 208 | + | |
| 209 | + // Simulate source events | |
| 210 | + simulate_source(obs, photons); | |
| 211 | + | |
| 212 | + // Write header for bacgkround event simulation | |
| 213 | + //if (logTerse()) | |
| 214 | + // log.header3("Simulate background events"); | |
| 215 | + | |
| 216 | + // Simulate source events | |
| 217 | + simulate_background(obs, m_obs.models()); | |
| 218 | + | |
| 219 | + } // endif: CTA observation found | |
| 220 | + | |
| 221 | + } // endfor: looped over observations | |
| 222 | + | |
| 223 | + // Return | |
| 224 | + return; | |
| 225 | +} | |
| 226 | + | |
| 227 | + | |
| 228 | +/***********************************************************************//** | |
| 229 | + * @brief Save simulated observation | |
| 230 | + ***************************************************************************/ | |
| 231 | +void ctobssim::save(void) | |
| 232 | +{ | |
| 138 | 233 | // Write header |
| 139 | 234 | if (logTerse()) { |
| 140 | 235 | log << std::endl; |
| 141 | - log.header1("Save events in FITS file"); | |
| 236 | + if (m_obs.size() > 1) | |
| 237 | + log.header1("Save observations"); | |
| 238 | + else | |
| 239 | + log.header1("Save observation"); | |
| 142 | 240 | } |
| 143 | 241 | |
| 144 | - // Save events | |
| 145 | - save_events(events); | |
| 242 | + // Loop over all observation in the container | |
| 243 | + int file_num = 0; | |
| 244 | + for (int i = 0; i < m_obs.size(); ++i) { | |
| 146 | 245 | |
| 147 | - // Release events | |
| 148 | - delete events; | |
| 246 | + // Get CTA observation | |
| 247 | + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs(i)); | |
| 149 | 248 | |
| 150 | - // Write separator into logger | |
| 151 | - if (logTerse()) | |
| 152 | - log << std::endl; | |
| 249 | + // Save only if observation is a CTA observation | |
| 250 | + if (obs != NULL) { | |
| 251 | + | |
| 252 | + // Set filename. If more than one file will be created an | |
| 253 | + // index "_xxx" will be appended. | |
| 254 | + std::string filename = m_outfile; | |
| 255 | + if (file_num > 0) | |
| 256 | + filename += "_"+str(file_num); | |
| 257 | + | |
| 258 | + // Save file | |
| 259 | + obs->save(filename, clobber()); | |
| 260 | + | |
| 261 | + } // endif: observation was a CTA observation | |
| 262 | + | |
| 263 | + } // endfor: looped over files | |
| 153 | 264 | |
| 154 | 265 | // Return |
| 155 | 266 | return; |
| ... | ... | @@ -160,195 +271,340 @@ void ctobssim::run(void) |
| 160 | 271 | * @brief Get application parameters |
| 161 | 272 | * |
| 162 | 273 | * Get all task parameters from parameter file or (if required) by querying |
| 163 | - * the user. | |
| 274 | + * the user. Most parameters are only required if no observation exists so | |
| 275 | + * far in the observation container. In this case, a single CTA observation | |
| 276 | + * will be added to the container, using the definition provided in the | |
| 277 | + * parameter file. | |
| 164 | 278 | ***************************************************************************/ |
| 165 | 279 | void ctobssim::get_parameters(void) |
| 166 | 280 | { |
| 167 | - // Get parameters | |
| 168 | - m_infile = par("infile")->value(); | |
| 281 | + // If there are no observations in container then add a single CTA | |
| 282 | + // observation using the parameters from the parameter file | |
| 283 | + if (m_obs.size() == 0) { | |
| 284 | + | |
| 285 | + // Get CTA observation parameters | |
| 286 | + m_infile = par("infile")->value(); | |
| 287 | + m_caldb = par("caldb")->value(); | |
| 288 | + m_irf = par("irf")->value(); | |
| 289 | + m_ra = par("ra")->real(); | |
| 290 | + m_dec = par("dec")->real(); | |
| 291 | + m_rad = par("rad")->real(); | |
| 292 | + m_tmin = par("tmin")->real(); | |
| 293 | + m_tmax = par("tmax")->real(); | |
| 294 | + m_emin = par("emin")->real(); | |
| 295 | + m_emax = par("emax")->real(); | |
| 296 | + | |
| 297 | + // Allocate CTA observation | |
| 298 | + GCTAObservation obs; | |
| 299 | + | |
| 300 | + // Set pointing direction | |
| 301 | + GCTAPointing pnt; | |
| 302 | + GSkyDir skydir; | |
| 303 | + skydir.radec_deg(m_ra, m_dec); | |
| 304 | + pnt.dir(skydir); | |
| 305 | + | |
| 306 | + // Set ROI | |
| 307 | + GCTARoi roi; | |
| 308 | + GCTAInstDir instdir; | |
| 309 | + instdir.radec_deg(m_ra, m_dec); | |
| 310 | + roi.centre(instdir); | |
| 311 | + roi.radius(m_rad); | |
| 312 | + | |
| 313 | + // Set GTI | |
| 314 | + GGti gti; | |
| 315 | + GTime tstart; | |
| 316 | + GTime tstop; | |
| 317 | + tstart.met(m_tmin); | |
| 318 | + tstop.met(m_tmax); | |
| 319 | + gti.append(tstart, tstop); | |
| 320 | + | |
| 321 | + // Set energy boundaries | |
| 322 | + GEbounds ebounds; | |
| 323 | + GEnergy emin; | |
| 324 | + GEnergy emax; | |
| 325 | + emin.TeV(m_emin); | |
| 326 | + emax.TeV(m_emax); | |
| 327 | + ebounds.append(emin, emax); | |
| 328 | + | |
| 329 | + // Set CTA observation attributes | |
| 330 | + obs.pointing(pnt); | |
| 331 | + obs.response(m_irf, m_caldb); | |
| 332 | + obs.roi(&roi); | |
| 333 | + obs.gti(gti); | |
| 334 | + obs.ebounds(ebounds); | |
| 335 | + | |
| 336 | + // Append CTA observation to container | |
| 337 | + m_obs.append(obs); | |
| 338 | + | |
| 339 | + // Load models into container | |
| 340 | + m_obs.models(m_infile); | |
| 341 | + | |
| 342 | + } // endif: there was no observation in the container | |
| 343 | + | |
| 344 | + // Get other parameters | |
| 169 | 345 | m_outfile = par("outfile")->value(); |
| 170 | - m_caldb = par("caldb")->value(); | |
| 171 | - m_irf = par("irf")->value(); | |
| 172 | 346 | m_seed = par("seed")->integer(); |
| 173 | - m_ra = par("ra")->real(); | |
| 174 | - m_dec = par("dec")->real(); | |
| 175 | - m_rad = par("rad")->real(); | |
| 176 | - m_tmin = par("tmin")->real(); | |
| 177 | - m_tmax = par("tmax")->real(); | |
| 178 | - m_emin = par("emin")->real(); | |
| 179 | - m_emax = par("emax")->real(); | |
| 180 | 347 | |
| 181 | 348 | // Initialise random number generator |
| 182 | 349 | m_ran.seed(m_seed); |
| 183 | 350 | |
| 184 | - // Set CTA pointing | |
| 185 | - GSkyDir dir; | |
| 186 | - dir.radec_deg(m_ra, m_dec); | |
| 187 | - m_pnt.dir(dir); | |
| 188 | - | |
| 189 | - // Set GTI | |
| 190 | - GGti gti; | |
| 191 | - GTime tstart; | |
| 192 | - GTime tstop; | |
| 193 | - tstart.met(m_tmin); | |
| 194 | - tstop.met(m_tmax); | |
| 195 | - gti.append(tstart, tstop); | |
| 196 | - | |
| 197 | - // Set energy boundaries | |
| 198 | - GEbounds ebds; | |
| 199 | - GEnergy emin; | |
| 200 | - GEnergy emax; | |
| 201 | - emin.TeV(m_emin); | |
| 202 | - emax.TeV(m_emax); | |
| 203 | - ebds.append(emin, emax); | |
| 204 | - | |
| 205 | - // Set observation | |
| 206 | - m_obs.ebounds(ebds); | |
| 207 | - m_obs.gti(gti); | |
| 208 | - m_obs.response(m_irf, m_caldb); | |
| 209 | - m_obs.pointing(m_pnt); | |
| 210 | - | |
| 211 | - // Load models | |
| 212 | - m_models.load(m_infile); | |
| 213 | - | |
| 214 | 351 | // Return |
| 215 | 352 | return; |
| 216 | 353 | } |
| 217 | 354 | |
| 218 | 355 | |
| 219 | 356 | /***********************************************************************//** |
| 220 | - * @brief Simulate photons | |
| 357 | + * @brief Simulate photons for a given CTA observation | |
| 358 | + * | |
| 359 | + * @param[in] obs Pointer on CTA observation. | |
| 360 | + * @param[in] models Models. | |
| 361 | + * | |
| 362 | + * This method simulates photons for a given CTA observation. It will loop | |
| 363 | + * over all energy boundaries, all Good Time Intervals and all sky models | |
| 364 | + * to generate a list of photons. | |
| 365 | + * | |
| 366 | + * @todo Add margins to the data selections so that we generate also photons | |
| 367 | + * that fall slightly outside the ROI and the energy interval so that | |
| 368 | + * they can be scatter within the ROI and energy interval by the | |
| 369 | + * PSF and energy dispersion. | |
| 221 | 370 | ***************************************************************************/ |
| 222 | -GPhotons ctobssim::simulate_photons(void) | |
| 371 | +GPhotons ctobssim::simulate_photons(const GCTAObservation* obs, | |
| 372 | + const GModels& models) | |
| 223 | 373 | { |
| 224 | - // Allocate MC parameters | |
| 225 | - GSkyDir dir; | |
| 226 | - GEnergy emin; | |
| 227 | - GEnergy emax; | |
| 228 | - GTime tmin; | |
| 229 | - GTime tmax; | |
| 230 | - | |
| 231 | - // Define MC parameters | |
| 232 | - dir.radec_deg(m_ra, m_dec); | |
| 233 | - emin.TeV(m_emin); | |
| 234 | - emax.TeV(m_emax); | |
| 235 | - tmin.met(m_tmin); | |
| 236 | - tmax.met(m_tmax); | |
| 374 | + // Initialise indentation for logging | |
| 375 | + int indent = 0; | |
| 376 | + | |
| 377 | + // Extract MC parameters from observation | |
| 378 | + const GCTARoi* roi = static_cast<const GCTARoi*>(obs->roi()); | |
| 379 | + GSkyDir dir = roi->centre().skydir(); | |
| 380 | + double rad = roi->radius(); | |
| 381 | + | |
| 382 | + // Dump simulation cone information | |
| 383 | + if (logNormal()) { | |
| 384 | + log << parformat("Random number seed"); | |
| 385 | + log << str(m_seed) << std::endl; | |
| 386 | + log << parformat("Simulation area"); | |
| 387 | + log << str(m_area) << " cm2" << std::endl; | |
| 388 | + log << parformat("Simulation cone"); | |
| 389 | + log << "RA=" << dir.ra_deg() << " deg"; | |
| 390 | + log << ", Dec=" << dir.dec_deg() << "deg"; | |
| 391 | + log << ", r=" << rad << " deg" << std::endl; | |
| 392 | + } | |
| 237 | 393 | |
| 238 | 394 | // Allocate photons |
| 239 | 395 | GPhotons photons; |
| 240 | 396 | |
| 241 | - // Simulate photons for all sky models | |
| 242 | - for (int i = 0; i < m_models.size(); ++i) { | |
| 397 | + // Loop over all Good Time Intervals | |
| 398 | + for (int it = 0; it < obs->gti().size(); ++it) { | |
| 399 | + | |
| 400 | + // Extract time interval | |
| 401 | + GTime tmin = obs->gti().tstart(it); | |
| 402 | + GTime tmax = obs->gti().tstop(it); | |
| 243 | 403 | |
| 244 | - // Get sky model (NULL if not a sky model) | |
| 245 | - GModelSky* model = dynamic_cast<GModelSky*>(m_models(i)); | |
| 404 | + // Dump time interval | |
| 405 | + if (logNormal()) { | |
| 406 | + if (obs->gti().size() > 1) { | |
| 407 | + indent++; | |
| 408 | + log.indent(indent); | |
| 409 | + } | |
| 410 | + log << parformat("Time interval"); | |
| 411 | + log << tmin << " - " << tmax << std::endl; | |
| 412 | + } | |
| 246 | 413 | |
| 247 | - // If we have a sky model then simulate photons | |
| 248 | - if (model != NULL) { | |
| 414 | + // Loop over all energy boundaries | |
| 415 | + for (int ie = 0; ie < obs->ebounds().size(); ++ie) { | |
| 249 | 416 | |
| 250 | - // Get photons | |
| 251 | - GPhotons p = model->mc(m_area, dir, m_rad, emin, emax, tmin, tmax, m_ran); | |
| 417 | + // Extract energy boundaries | |
| 418 | + GEnergy emin = obs->ebounds().emin(ie); | |
| 419 | + GEnergy emax = obs->ebounds().emax(ie); | |
| 252 | 420 | |
| 253 | - // Reserve new space for photons | |
| 254 | - photons.reserve(photons.size() + p.size()); | |
| 255 | - | |
| 256 | - // Add photons to list | |
| 257 | - for (int k = 0; k < p.size(); ++k) { | |
| 258 | - p[k].mcid(i); | |
| 259 | - photons.push_back(p[k]); | |
| 421 | + // Dump energy range | |
| 422 | + if (logNormal()) { | |
| 423 | + if (obs->ebounds().size() > 1) { | |
| 424 | + indent++; | |
| 425 | + log.indent(indent); | |
| 426 | + } | |
| 427 | + log << parformat("Energy range"); | |
| 428 | + log << emin << " - " << emax << std::endl; | |
| 260 | 429 | } |
| 261 | 430 | |
| 262 | - } // endif: model was a sky model | |
| 431 | + // Loop over all sky models | |
| 432 | + for (int i = 0; i < models.size(); ++i) { | |
| 263 | 433 | |
| 264 | - } // endfor: looped over models | |
| 434 | + // Get sky model (NULL if not a sky model) | |
| 435 | + const GModelSky* model = dynamic_cast<const GModelSky*>(models(i)); | |
| 436 | + | |
| 437 | + // If we have a sky model then simulate photons | |
| 438 | + if (model != NULL) { | |
| 439 | + | |
| 440 | + // Get photons | |
| 441 | + GPhotons p = model->mc(m_area, dir, rad, emin, emax, | |
| 442 | + tmin, tmax, m_ran); | |
| 443 | + | |
| 444 | + // Reserve new space for photons | |
| 445 | + photons.reserve(photons.size() + p.size()); | |
| 446 | + | |
| 447 | + // Add photons to list | |
| 448 | + for (int k = 0; k < p.size(); ++k) { | |
| 449 | + p[k].mcid(i); | |
| 450 | + photons.push_back(p[k]); | |
| 451 | + } | |
| 452 | + | |
| 453 | + // Dump simulation result for model | |
| 454 | + if (logNormal()) { | |
| 455 | + log << parformat("MC source photons"); | |
| 456 | + log << str(p.size()); | |
| 457 | + if (model->name().length() > 0) | |
| 458 | + log << " [" << model->name() << "]"; | |
| 459 | + log << std::endl; | |
| 460 | + } | |
| 461 | + | |
| 462 | + } // endif: model was a sky model | |
| 463 | + | |
| 464 | + } // endfor: looped over models | |
| 465 | + | |
| 466 | + } // endfor: looped over all energy boundaries | |
| 467 | + | |
| 468 | + } // endfor: looped over all time intervals | |
| 469 | + | |
| 470 | + // Reset indentation | |
| 471 | + log.indent(0); | |
| 265 | 472 | |
| 266 | - // Dump simulation results | |
| 267 | - if (logNormal()) { | |
| 268 | - log << parformat("Simulated photons") << str(photons.size()) << std::endl; | |
| 269 | - log << parformat("Simulation area"); | |
| 270 | - log << str(m_area) << " cm2" << std::endl; | |
| 271 | - log << parformat("Simulation cone"); | |
| 272 | - log << "RA=" << dir.ra_deg() << " deg"; | |
| 273 | - log << ", Dec=" << dir.dec_deg() << "deg"; | |
| 274 | - log << ", r=" << m_rad << " deg" << std::endl; | |
| 275 | - log << parformat("Time interval"); | |
| 276 | - log << tmin << " - " << tmax << std::endl; | |
| 277 | - log << parformat("Energy range"); | |
| 278 | - log << emin << " - " << emax << std::endl; | |
| 279 | - log << parformat("Random number seed") << str(m_seed) << std::endl; | |
| 280 | - } | |
| 281 | - | |
| 282 | 473 | // Return photons |
| 283 | 474 | return photons; |
| 284 | 475 | } |
| 285 | 476 | |
| 286 | 477 | |
| 287 | 478 | /***********************************************************************//** |
| 288 | - * @brief Simulate events | |
| 479 | + * @brief Simulate source events from photon list | |
| 289 | 480 | * |
| 481 | + * @param[in] obs Pointer on CTA observation. | |
| 290 | 482 | * @param[in] photons Photon list. |
| 291 | 483 | * |
| 292 | - * Simulate events from a photon list using a given response function and | |
| 293 | - * a given pointing. | |
| 484 | + * @exception GCTAException::no_pointing | |
| 485 | + * No valid pointing found in CTA observation | |
| 486 | + * @exception GException::no_response | |
| 487 | + * No valid response found in CTA observation | |
| 488 | + * | |
| 489 | + * Simulate source events from a photon list for a given CTA observation. | |
| 490 | + * The events are stored in as event list in the observation. | |
| 491 | + * | |
| 492 | + * This method does nothing if the observation pointer is NULL. It also | |
| 493 | + * verifies if the observation has a valid pointing and response. | |
| 294 | 494 | ***************************************************************************/ |
| 295 | -GCTAEventList* ctobssim::simulate_events(const GPhotons& photons) | |
| 495 | +void ctobssim::simulate_source(GCTAObservation* obs, const GPhotons& photons) | |
| 296 | 496 | { |
| 297 | - // Allocate CTA event list | |
| 298 | - GCTAEventList* events = new GCTAEventList; | |
| 299 | - | |
| 300 | - // Reserves space for events | |
| 301 | - events->reserve(photons.size()); | |
| 302 | - | |
| 303 | - // Simulate events from photons | |
| 304 | - for (int i = 0; i < photons.size(); ++i) { | |
| 305 | - GCTAEventAtom* event = m_obs.response()->mc(m_area, photons[i], m_pnt, m_ran); | |
| 306 | - if (event != NULL) | |
| 307 | - events->append(*event); | |
| 308 | - delete event; | |
| 309 | - } | |
| 497 | + // Continue only if observation pointer is valid | |
| 498 | + if (obs != NULL) { | |
| 499 | + | |
| 500 | + // Get pointer on CTA pointing. As we have no time dependent | |
| 501 | + // pointings implement we just determine the pointing from the | |
| 502 | + // beginning of the observation. Throw an exception if the pointing | |
| 503 | + // is not defined. | |
| 504 | + GCTAPointing* pnt = obs->pointing(obs->tstart()); | |
| 505 | + if (pnt == NULL) | |
| 506 | + throw GCTAException::no_pointing(G_SIMULATE_SOURCE); | |
| 507 | + | |
| 508 | + // Get pointer on CTA response. Throw an exception if the response | |
| 509 | + // is not defined. | |
| 510 | + GCTAResponse* rsp = obs->response(); | |
| 511 | + if (rsp == NULL) | |
| 512 | + throw GException::no_response(G_SIMULATE_SOURCE); | |
| 513 | + | |
| 514 | + // Make sure that the observation holds a CTA event list. If this | |
| 515 | + // is not the case then allocate and attach a CTA event list now. | |
| 516 | + if (dynamic_cast<const GCTAEventList*>(obs->events()) == NULL) | |
| 517 | + obs->events(new GCTAEventList); | |
| 518 | + | |
| 519 | + // Get pointer on event list (circumvent const correctness) | |
| 520 | + GCTAEventList* events = (GCTAEventList*)(obs->events()); | |
| 521 | + | |
| 522 | + // Reserves space for source events in event list | |
| 523 | + events->reserve(photons.size()); | |
| 524 | + | |
| 525 | + // Simulate events from photons | |
| 526 | + for (int i = 0; i < photons.size(); ++i) { | |
| 527 | + GCTAEventAtom* event = rsp->mc(m_area, photons[i], *pnt, m_ran); | |
| 528 | + if (event != NULL) { | |
| 529 | + events->append(*event); | |
| 530 | + delete event; | |
| 531 | + } | |
| 532 | + } | |
| 310 | 533 | |
| 311 | - // Dump simulation results | |
| 312 | - if (logNormal()) { | |
| 313 | - log << parformat("MC source events"); | |
| 314 | - log << str(events->size()) << std::endl; | |
| 315 | - } | |
| 534 | + // Dump simulation results | |
| 535 | + if (logNormal()) { | |
| 536 | + log << parformat("MC source events"); | |
| 537 | + log << str(events->size()) << std::endl; | |
| 538 | + } | |
| 539 | + | |
| 540 | + } // endif: observation pointer was valid | |
| 541 | + | |
| 542 | + // Return | |
| 543 | + return; | |
| 544 | +} | |
| 545 | + | |
| 546 | + | |
| 547 | +/***********************************************************************//** | |
| 548 | + * @brief Simulate background events from model | |
| 549 | + * | |
| 550 | + * @param[in] obs Pointer on CTA observation. | |
| 551 | + * @param[in] models Models. | |
| 552 | + * | |
| 553 | + * Simulate background events from models. The events are stored as event | |
| 554 | + * list in the observation. | |
| 555 | + * | |
| 556 | + * This method does nothing if the observation pointer is NULL. | |
| 557 | + ***************************************************************************/ | |
| 558 | +void ctobssim::simulate_background(GCTAObservation* obs, const GModels& models) | |
| 559 | +{ | |
| 560 | + // Continue only if observation pointer is valid | |
| 561 | + if (obs != NULL) { | |
| 316 | 562 | |
| 317 | - // Simulate background events | |
| 318 | - for (int i = 0; i < m_models.size(); ++i) { | |
| 563 | + // Make sure that the observation holds a CTA event list. If this | |
| 564 | + // is not the case then allocate and attach a CTA event list now. | |
| 565 | + if (dynamic_cast<const GCTAEventList*>(obs->events()) == NULL) | |
| 566 | + obs->events(new GCTAEventList); | |
| 567 | + | |
| 568 | + // Get pointer on event list (circumvent const correctness) | |
| 569 | + GCTAEventList* events = (GCTAEventList*)(obs->events()); | |
| 319 | 570 | |
| 320 | - // Get model (NULL if not a radial acceptance model) | |
| 321 | - GCTAModelRadialAcceptance* model = | |
| 322 | - dynamic_cast<GCTAModelRadialAcceptance*>(m_models(i)); | |
| 571 | + // Loop over all models | |
| 572 | + for (int i = 0; i < models.size(); ++i) { | |
| 323 | 573 | |
| 324 | - // If we have a radial acceptance model then simulate events | |
| 325 | - if (model != NULL) { | |
| 574 | + // Get model (NULL if not a radial acceptance model) | |
| 575 | + const GCTAModelRadialAcceptance* model = | |
| 576 | + dynamic_cast<const GCTAModelRadialAcceptance*>(models(i)); | |
| 326 | 577 | |
| 327 | - // Get simulated event list | |
| 328 | - GCTAEventList* list = model->mc(m_obs, m_ran); | |
| 578 | + // If we have a radial acceptance model then simulate events | |
| 579 | + if (model != NULL) { | |
| 329 | 580 | |
| 330 | - // Reserves space for events | |
| 331 | - events->reserve(list->size()+events->size()); | |
| 581 | + // Get simulated event list | |
| 582 | + GCTAEventList* list = model->mc(*obs, m_ran); | |
| 332 | 583 | |
| 333 | - // Append events | |
| 334 | - for (int k = 0; k < list->size(); k++) | |
| 335 | - events->append(*(list->pointer(k))); | |
| 584 | + // Reserves space for events | |
| 585 | + events->reserve(list->size()+events->size()); | |
| 336 | 586 | |
| 337 | - // Dump simulation results | |
| 338 | - if (logNormal()) { | |
| 339 | - log << parformat("MC background events"); | |
| 340 | - log << str(list->size()) << std::endl; | |
| 341 | - } | |
| 587 | + // Append events | |
| 588 | + for (int k = 0; k < list->size(); k++) | |
| 589 | + events->append(*(list->pointer(k))); | |
| 590 | + | |
| 591 | + // Dump simulation results | |
| 592 | + if (logNormal()) { | |
| 593 | + log << parformat("MC background events"); | |
| 594 | + log << str(list->size()) << std::endl; | |
| 595 | + } | |
| 596 | + | |
| 597 | + // Free event list | |
| 598 | + delete list; | |
| 342 | 599 | |
| 343 | - // Free event list | |
| 344 | - delete list; | |
| 600 | + } // endif: model was valid | |
| 345 | 601 | |
| 346 | - } // endif: model was valid | |
| 602 | + } // endfor: looped over all models | |
| 347 | 603 | |
| 348 | - } // endfor: looped over all models | |
| 604 | + } // endif: observation pointer was valid | |
| 349 | 605 | |
| 350 | - // Return event list | |
| 351 | - return events; | |
| 606 | + // Return | |
| 607 | + return; | |
| 352 | 608 | } |
| 353 | 609 | |
| 354 | 610 | |
| ... | ... | @@ -357,6 +613,7 @@ GCTAEventList* ctobssim::simulate_events(const GPhotons& photons) |
| 357 | 613 | * |
| 358 | 614 | * @todo Set header keywords |
| 359 | 615 | ***************************************************************************/ |
| 616 | +/* | |
| 360 | 617 | void ctobssim::save_events(const GCTAEventList* events) |
| 361 | 618 | { |
| 362 | 619 | // Allocate events file |
| ... | ... | @@ -396,13 +653,14 @@ void ctobssim::save_events(const GCTAEventList* events) |
| 396 | 653 | // Return |
| 397 | 654 | return; |
| 398 | 655 | } |
| 399 | - | |
| 656 | +*/ | |
| 400 | 657 | |
| 401 | 658 | /***********************************************************************//** |
| 402 | 659 | * @brief Set keywords of EVENTS header |
| 403 | 660 | * |
| 404 | 661 | * @param[in] file FITS file. |
| 405 | 662 | ***************************************************************************/ |
| 663 | +/* | |
| 406 | 664 | void ctobssim::set_events_keywords(GFits* file) |
| 407 | 665 | { |
| 408 | 666 | // Get EVENTS HDU |
| ... | ... | @@ -426,7 +684,7 @@ void ctobssim::set_events_keywords(GFits* file) |
| 426 | 684 | // Return |
| 427 | 685 | return; |
| 428 | 686 | } |
| 429 | - | |
| 687 | +*/ | |
| 430 | 688 | |
| 431 | 689 | /*========================================================================== |
| 432 | 690 | = = |
| ... | ... | @@ -445,7 +703,6 @@ void ctobssim::init_members(void) |
| 445 | 703 | m_caldb.clear(); |
| 446 | 704 | m_irf.clear(); |
| 447 | 705 | m_obs.clear(); |
| 448 | - m_pnt.clear(); | |
| 449 | 706 | m_seed = 1; |
| 450 | 707 | m_ra = 0.0; |
| 451 | 708 | m_dec = 0.0; |
| ... | ... | @@ -471,6 +728,10 @@ void ctobssim::init_members(void) |
| 471 | 728 | ***************************************************************************/ |
| 472 | 729 | void ctobssim::free_members(void) |
| 473 | 730 | { |
| 731 | + // Write separator into logger | |
| 732 | + if (logTerse()) | |
| 733 | + log << std::endl; | |
| 734 | + | |
| 474 | 735 | // Return |
| 475 | 736 | return; |
| 476 | 737 | } | ... | ... |