From 455ed48861e8076febd7cdceae29c9160bf0d20d Mon Sep 17 00:00:00 2001 From: knodlseder <> Date: Fri, 9 Mar 2012 10:50:02 +0000 Subject: [PATCH] Change GPointing argument to GObservation argument in GCTAResponse::mc method. This allows taking into account the deadtime correction in the GCTAResponse::mc method. --- inst/cta/include/GCTAResponse.hpp | 2 +- inst/cta/pyext/GCTAResponse.i | 2 +- inst/cta/src/GCTAResponse.cpp | 95 +++++++++++++++++++++++++++++++++++++++++++++++++++-------------------------------------------- 3 files changed, 53 insertions(+), 46 deletions(-) diff --git a/inst/cta/include/GCTAResponse.hpp b/inst/cta/include/GCTAResponse.hpp index c763fca..53a3bed 100644 --- a/inst/cta/include/GCTAResponse.hpp +++ b/inst/cta/include/GCTAResponse.hpp @@ -134,7 +134,7 @@ public: // Other Methods GCTAEventAtom* mc(const double& area, const GPhoton& photon, - const GPointing& pnt, GRan& ran) const; + const GObservation& obs, GRan& ran) const; void caldb(const std::string& caldb); std::string caldb(void) const { return m_caldb; } void load(const std::string& rspname); diff --git a/inst/cta/pyext/GCTAResponse.i b/inst/cta/pyext/GCTAResponse.i index 173463b..10071ad 100644 --- a/inst/cta/pyext/GCTAResponse.i +++ b/inst/cta/pyext/GCTAResponse.i @@ -82,7 +82,7 @@ public: // Other Methods GCTAEventAtom* mc(const double& area, const GPhoton& photon, - const GPointing& pnt, GRan& ran) const; + const GObservation& obs, GRan& ran) const; void caldb(const std::string& caldb); std::string caldb(void) const; void load(const std::string& rspname); diff --git a/inst/cta/src/GCTAResponse.cpp b/inst/cta/src/GCTAResponse.cpp index 0d19a63..a31bc81 100644 --- a/inst/cta/src/GCTAResponse.cpp +++ b/inst/cta/src/GCTAResponse.cpp @@ -452,77 +452,84 @@ double GCTAResponse::npred(const GSkyDir& srcDir, * * @param[in] area Simulation surface area. * @param[in] photon Photon. - * @param[in] pnt Pointing. + * @param[in] obs Observation. * @param[in] ran Random number generator. * - * @exception GCTAException::bad_pointing_type - * Specified pointing is not a CTA pointing. + * @exception GCTAException::no_pointing + * No CTA pointing found in observation. * * Simulates a CTA event using the response function from an incident photon. * If the event is not detected a NULL pointer is returned. * - * Note that this method does not take into account any deadtime correction. - * This has to be applied externally (hence by adjusting the number of - * photons to the livetime instead of the ontime). + * The method also applies a deadtime correction using a Monte Carlo process, + * taking into account temporal deadtime variations. For this purpose, the + * method makes use of the time dependent GObservation::deadc method. * * @todo Implement Phi dependence in CTA IRF * @todo Implement energy dispersion * @todo Implement the method for a 3 Gaussian PSF ***************************************************************************/ GCTAEventAtom* GCTAResponse::mc(const double& area, const GPhoton& photon, - const GPointing& pnt, GRan& ran) const + const GObservation& obs, GRan& ran) const { // Initialise event GCTAEventAtom* event = NULL; - // Get pointer on CTA pointing - const GCTAPointing *ctapnt = dynamic_cast(&pnt); - if (ctapnt == NULL) { - throw GCTAException::bad_pointing_type(G_MC); - } + // Continue only if the photon arrives at the moment when the detector + // is not dead. This implements the deadtime correction. + double deadc = obs.deadc(photon.time()); + if (deadc >= 1.0 || ran.uniform() <= deadc) { - // Get pointing direction zenith angle and azimuth [radians] - double zenith = ctapnt->zenith(); - double azimuth = ctapnt->azimuth(); + // Get pointer on CTA pointing + GCTAPointing* pnt = dynamic_cast(obs.pointing()); + if (pnt == NULL) { + throw GCTAException::no_pointing(G_MC); + } - // Get radial offset and polar angles of true photon in camera [radians] - double theta_p = ctapnt->dir().dist(photon.dir()); - double phi_p = 0.0; //TODO Implement Phi dependence + // Get pointing direction zenith angle and azimuth [radians] + double zenith = pnt->zenith(); + double azimuth = pnt->azimuth(); + + // Get radial offset and polar angles of true photon in camera [radians] + double theta_p = pnt->dir().dist(photon.dir()); + double phi_p = 0.0; //TODO Implement Phi dependence + + // Compute effective area for photon + double srcLogEng = photon.energy().log10TeV(); + double effective_area = aeff(theta_p, phi_p, zenith, azimuth, srcLogEng); - // Compute effective area for photon - double srcLogEng = photon.energy().log10TeV(); - double effective_area = aeff(theta_p, phi_p, zenith, azimuth, srcLogEng); + // Compute limiting value + double ulimite = effective_area / area; - // Compute limiting value - double ulimite = effective_area / area; + // Continue only if event is detected + if (ran.uniform() <= ulimite) { - // Continue only if event is detected - if (ran.uniform() <= ulimite) { + // Simulate offset from photon arrival direction + //TODO: Make a proper implementation depending on the response + // version. For now, the first Gaussian is used. + GCTAPsfPars pars = psf_dummy_sigma(srcLogEng, theta_p); + double theta = pars[1] * ran.chisq2() * rad2deg; + double phi = 360.0 * ran.uniform(); - // Simulate offset from photon arrival direction - //TODO: Make a proper implementation depending on the response - // version. For now, the first Gaussian is used. - GCTAPsfPars pars = psf_dummy_sigma(srcLogEng, theta_p); - double theta = pars[1] * ran.chisq2() * rad2deg; - double phi = 360.0 * ran.uniform(); + // Rotate sky direction by offset + GSkyDir sky_dir = photon.dir(); + sky_dir.rotate_deg(phi, theta); - // Rotate sky direction by offset - GSkyDir sky_dir = photon.dir(); - sky_dir.rotate_deg(phi, theta); + // Set measured photon arrival direction + GCTAInstDir inst_dir; + inst_dir.skydir(sky_dir); - // Set measured photon arrival direction - GCTAInstDir inst_dir; - inst_dir.skydir(sky_dir); + // Allocate event + event = new GCTAEventAtom; - // Allocate event - event = new GCTAEventAtom; + // Set event attributes + event->dir(inst_dir); + event->energy(photon.energy()); + event->time(photon.time()); - // Set event attributes - event->dir(inst_dir); - event->energy(photon.energy()); - event->time(photon.time()); + } // endif: event was detected - } // endif: event was detected + } // endif: detector was alive // Return event return event; -- libgit2 0.21.4