Leonardo Di Venere / gammalib

Commit 455ed48861e8076febd7cdceae29c9160bf0d20d

Authored by knodlseder
1 parent aec4bcdd

Change GPointing argument to GObservation argument in GCTAResponse::mc method. T… 

…his allows taking into account the deadtime correction in the GCTAResponse::mc method.
Inline Side-by-side
Showing 3 changed files with 53 additions and 46 deletions
inst/cta/include/GCTAResponse.hpp
  View file @455ed48
... ... @@ -134,7 +134,7 @@ public:
134 134  
135 135 // Other Methods
136 136 GCTAEventAtom* mc(const double& area, const GPhoton& photon,
137   - const GPointing& pnt, GRan& ran) const;
  137 + const GObservation& obs, GRan& ran) const;
138 138 void caldb(const std::string& caldb);
139 139 std::string caldb(void) const { return m_caldb; }
140 140 void load(const std::string& rspname);
... ...
inst/cta/pyext/GCTAResponse.i
  View file @455ed48
... ... @@ -82,7 +82,7 @@ public:
82 82  
83 83 // Other Methods
84 84 GCTAEventAtom* mc(const double& area, const GPhoton& photon,
85   - const GPointing& pnt, GRan& ran) const;
  85 + const GObservation& obs, GRan& ran) const;
86 86 void caldb(const std::string& caldb);
87 87 std::string caldb(void) const;
88 88 void load(const std::string& rspname);
... ...
inst/cta/src/GCTAResponse.cpp
  View file @455ed48
... ... @@ -452,77 +452,84 @@ double GCTAResponse::npred(const GSkyDir& srcDir,
452 452 *
453 453 * @param[in] area Simulation surface area.
454 454 * @param[in] photon Photon.
455   - * @param[in] pnt Pointing.
  455 + * @param[in] obs Observation.
456 456 * @param[in] ran Random number generator.
457 457 *
458   - * @exception GCTAException::bad_pointing_type
459   - * Specified pointing is not a CTA pointing.
  458 + * @exception GCTAException::no_pointing
  459 + * No CTA pointing found in observation.
460 460 *
461 461 * Simulates a CTA event using the response function from an incident photon.
462 462 * If the event is not detected a NULL pointer is returned.
463 463 *
464   - * Note that this method does not take into account any deadtime correction.
465   - * This has to be applied externally (hence by adjusting the number of
466   - * photons to the livetime instead of the ontime).
  464 + * The method also applies a deadtime correction using a Monte Carlo process,
  465 + * taking into account temporal deadtime variations. For this purpose, the
  466 + * method makes use of the time dependent GObservation::deadc method.
467 467 *
468 468 * @todo Implement Phi dependence in CTA IRF
469 469 * @todo Implement energy dispersion
470 470 * @todo Implement the method for a 3 Gaussian PSF
471 471 ***************************************************************************/
472 472 GCTAEventAtom* GCTAResponse::mc(const double& area, const GPhoton& photon,
473   - const GPointing& pnt, GRan& ran) const
  473 + const GObservation& obs, GRan& ran) const
474 474 {
475 475 // Initialise event
476 476 GCTAEventAtom* event = NULL;
477 477  
478   - // Get pointer on CTA pointing
479   - const GCTAPointing *ctapnt = dynamic_cast<const GCTAPointing*>(&pnt);
480   - if (ctapnt == NULL) {
481   - throw GCTAException::bad_pointing_type(G_MC);
482   - }
  478 + // Continue only if the photon arrives at the moment when the detector
  479 + // is not dead. This implements the deadtime correction.
  480 + double deadc = obs.deadc(photon.time());
  481 + if (deadc >= 1.0 || ran.uniform() <= deadc) {
483 482  
484   - // Get pointing direction zenith angle and azimuth [radians]
485   - double zenith = ctapnt->zenith();
486   - double azimuth = ctapnt->azimuth();
  483 + // Get pointer on CTA pointing
  484 + GCTAPointing* pnt = dynamic_cast<GCTAPointing*>(obs.pointing());
  485 + if (pnt == NULL) {
  486 + throw GCTAException::no_pointing(G_MC);
  487 + }
487 488  
488   - // Get radial offset and polar angles of true photon in camera [radians]
489   - double theta_p = ctapnt->dir().dist(photon.dir());
490   - double phi_p = 0.0; //TODO Implement Phi dependence
  489 + // Get pointing direction zenith angle and azimuth [radians]
  490 + double zenith = pnt->zenith();
  491 + double azimuth = pnt->azimuth();
  492 +
  493 + // Get radial offset and polar angles of true photon in camera [radians]
  494 + double theta_p = pnt->dir().dist(photon.dir());
  495 + double phi_p = 0.0; //TODO Implement Phi dependence
  496 +
  497 + // Compute effective area for photon
  498 + double srcLogEng = photon.energy().log10TeV();
  499 + double effective_area = aeff(theta_p, phi_p, zenith, azimuth, srcLogEng);
491 500  
492   - // Compute effective area for photon
493   - double srcLogEng = photon.energy().log10TeV();
494   - double effective_area = aeff(theta_p, phi_p, zenith, azimuth, srcLogEng);
  501 + // Compute limiting value
  502 + double ulimite = effective_area / area;
495 503  
496   - // Compute limiting value
497   - double ulimite = effective_area / area;
  504 + // Continue only if event is detected
  505 + if (ran.uniform() <= ulimite) {
498 506  
499   - // Continue only if event is detected
500   - if (ran.uniform() <= ulimite) {
  507 + // Simulate offset from photon arrival direction
  508 + //TODO: Make a proper implementation depending on the response
  509 + // version. For now, the first Gaussian is used.
  510 + GCTAPsfPars pars = psf_dummy_sigma(srcLogEng, theta_p);
  511 + double theta = pars[1] * ran.chisq2() * rad2deg;
  512 + double phi = 360.0 * ran.uniform();
501 513  
502   - // Simulate offset from photon arrival direction
503   - //TODO: Make a proper implementation depending on the response
504   - // version. For now, the first Gaussian is used.
505   - GCTAPsfPars pars = psf_dummy_sigma(srcLogEng, theta_p);
506   - double theta = pars[1] * ran.chisq2() * rad2deg;
507   - double phi = 360.0 * ran.uniform();
  514 + // Rotate sky direction by offset
  515 + GSkyDir sky_dir = photon.dir();
  516 + sky_dir.rotate_deg(phi, theta);
508 517  
509   - // Rotate sky direction by offset
510   - GSkyDir sky_dir = photon.dir();
511   - sky_dir.rotate_deg(phi, theta);
  518 + // Set measured photon arrival direction
  519 + GCTAInstDir inst_dir;
  520 + inst_dir.skydir(sky_dir);
512 521  
513   - // Set measured photon arrival direction
514   - GCTAInstDir inst_dir;
515   - inst_dir.skydir(sky_dir);
  522 + // Allocate event
  523 + event = new GCTAEventAtom;
516 524  
517   - // Allocate event
518   - event = new GCTAEventAtom;
  525 + // Set event attributes
  526 + event->dir(inst_dir);
  527 + event->energy(photon.energy());
  528 + event->time(photon.time());
519 529  
520   - // Set event attributes
521   - event->dir(inst_dir);
522   - event->energy(photon.energy());
523   - event->time(photon.time());
  530 + } // endif: event was detected
524 531  
525   - } // endif: event was detected
  532 + } // endif: detector was alive
526 533  
527 534 // Return event
528 535 return event;
... ...