Commit 298cd84ab098d34a7e516005c1aab12166e2c97d
1 parent
f924c3ca
Use dispersion matrix instead of eres/ebias
Showing
4 changed files
with
624 additions
and
14 deletions
inst/cta/include/GCTAEdisp2D.hpp
| ... | ... | @@ -61,7 +61,7 @@ public: |
| 61 | 61 | |
| 62 | 62 | // Implemented methods |
| 63 | 63 | void clear(void); |
| 64 | - GCTAEdisp* clone(void) const; | |
| 64 | + GCTAEdisp2D* clone(void) const; | |
| 65 | 65 | std::string classname(void) const; |
| 66 | 66 | void load(const std::string& filename); |
| 67 | 67 | std::string filename(void) const; |
| ... | ... | @@ -107,6 +107,8 @@ private: |
| 107 | 107 | mutable double m_par_theta; |
| 108 | 108 | mutable double m_norm; |
| 109 | 109 | mutable double m_EobsOnEsrc; //!< linear ratio Eobs/Esrc |
| 110 | + mutable double m_eres; | |
| 111 | + mutable double m_ebias; | |
| 110 | 112 | }; |
| 111 | 113 | |
| 112 | 114 | /***********************************************************************//** |
| ... | ... | @@ -133,4 +135,4 @@ std::string GCTAEdisp2D::filename(void) const |
| 133 | 135 | return m_filename; |
| 134 | 136 | } |
| 135 | 137 | |
| 136 | -#endif /* GCTAEDISP2D_HPP */ | |
| 137 | 138 | \ No newline at end of file |
| 139 | +#endif /* GCTAEDISP2D_HPP */ | ... | ... |
inst/cta/src/GCTAEdisp2D.cpp
0 → 100644
| 1 | +/*************************************************************************** | |
| 2 | + * GCTAEdisp2D.cpp - CTA 2D energy dispersion class * | |
| 3 | + * ----------------------------------------------------------------------- * | |
| 4 | + * copyright (C) 2012-2015 by Juergen Knoedlseder * | |
| 5 | + * ----------------------------------------------------------------------- * | |
| 6 | + * * | |
| 7 | + * This program is free software: you can redistribute it and/or modify * | |
| 8 | + * it under the terms of the GNU General Public License as published by * | |
| 9 | + * the Free Software Foundation, either version 3 of the License, or * | |
| 10 | + * (at your option) any later version. * | |
| 11 | + * * | |
| 12 | + * This program is distributed in the hope that it will be useful, * | |
| 13 | + * but WITHOUT ANY WARRANTY; without even the implied warranty of * | |
| 14 | + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * | |
| 15 | + * GNU General Public License for more details. * | |
| 16 | + * * | |
| 17 | + * You should have received a copy of the GNU General Public License * | |
| 18 | + * along with this program. If not, see <http://www.gnu.org/licenses/>. * | |
| 19 | + * * | |
| 20 | + ***************************************************************************/ | |
| 21 | +/** | |
| 22 | + * @file GCTAEdisp2D.cpp | |
| 23 | + * @brief CTA 2D energy dispersion class implementation | |
| 24 | + * @author Florent Forest | |
| 25 | + */ | |
| 26 | + | |
| 27 | +/* __ Includes ___________________________________________________________ */ | |
| 28 | +#ifdef HAVE_CONFIG_H | |
| 29 | +#include <config.h> | |
| 30 | +#endif | |
| 31 | +#include <cmath> | |
| 32 | +#include "GTools.hpp" | |
| 33 | +#include "GFits.hpp" | |
| 34 | +#include "GFitsBinTable.hpp" | |
| 35 | +#include "GFitsTable.hpp" | |
| 36 | +#include "GCTAEdisp2D.hpp" | |
| 37 | +#include "GMath.hpp" | |
| 38 | + | |
| 39 | +/* __ Method name definitions ____________________________________________ */ | |
| 40 | + | |
| 41 | +/* __ Macros _____________________________________________________________ */ | |
| 42 | + | |
| 43 | +/* __ Coding definitions _________________________________________________ */ | |
| 44 | + | |
| 45 | +/* __ Debug definitions __________________________________________________ */ | |
| 46 | + | |
| 47 | +/* __ Constants __________________________________________________________ */ | |
| 48 | + | |
| 49 | + | |
| 50 | +/*========================================================================== | |
| 51 | + = = | |
| 52 | + = Constructors/destructors = | |
| 53 | + = = | |
| 54 | + ==========================================================================*/ | |
| 55 | + | |
| 56 | +/***********************************************************************//** | |
| 57 | + * @brief Void constructor | |
| 58 | + ***************************************************************************/ | |
| 59 | +GCTAEdisp2D::GCTAEdisp2D(void) : GCTAEdisp() | |
| 60 | +{ | |
| 61 | + // Initialise class members | |
| 62 | + init_members(); | |
| 63 | + | |
| 64 | + // Return | |
| 65 | + return; | |
| 66 | +} | |
| 67 | + | |
| 68 | + | |
| 69 | +/***********************************************************************//** | |
| 70 | + * @brief File constructor | |
| 71 | + * | |
| 72 | + * @param[in] filename FITS file name. | |
| 73 | + * | |
| 74 | + * Construct instance by loading the energy dispersion information from a FITS | |
| 75 | + * file. | |
| 76 | + ***************************************************************************/ | |
| 77 | +GCTAEdisp2D::GCTAEdisp2D(const std::string& filename) : GCTAEdisp() | |
| 78 | +{ | |
| 79 | + // Initialise class members | |
| 80 | + init_members(); | |
| 81 | + | |
| 82 | + // Load energy dispersion from file | |
| 83 | + load(filename); | |
| 84 | + | |
| 85 | + // Return | |
| 86 | + return; | |
| 87 | +} | |
| 88 | + | |
| 89 | + | |
| 90 | +/***********************************************************************//** | |
| 91 | + * @brief Copy constructor | |
| 92 | + * | |
| 93 | + * @param[in] edisp Energy dispersion. | |
| 94 | + ***************************************************************************/ | |
| 95 | +GCTAEdisp2D::GCTAEdisp2D(const GCTAEdisp2D& edisp) : GCTAEdisp(edisp) | |
| 96 | +{ | |
| 97 | + // Initialise class members | |
| 98 | + init_members(); | |
| 99 | + | |
| 100 | + // Copy members | |
| 101 | + copy_members(edisp); | |
| 102 | + | |
| 103 | + // Return | |
| 104 | + return; | |
| 105 | +} | |
| 106 | + | |
| 107 | + | |
| 108 | +/***********************************************************************//** | |
| 109 | + * @brief Destructor | |
| 110 | + ***************************************************************************/ | |
| 111 | +GCTAEdisp2D::~GCTAEdisp2D(void) | |
| 112 | +{ | |
| 113 | + // Free members | |
| 114 | + free_members(); | |
| 115 | + | |
| 116 | + // Return | |
| 117 | + return; | |
| 118 | +} | |
| 119 | + | |
| 120 | + | |
| 121 | +/*========================================================================== | |
| 122 | + = = | |
| 123 | + = Operators = | |
| 124 | + = = | |
| 125 | + ==========================================================================*/ | |
| 126 | + | |
| 127 | +/***********************************************************************//** | |
| 128 | + * @brief Assignment operator | |
| 129 | + * | |
| 130 | + * @param[in] edisp Energy dispersion. | |
| 131 | + * @return Energy dispersion. | |
| 132 | + ***************************************************************************/ | |
| 133 | +GCTAEdisp2D& GCTAEdisp2D::operator= (const GCTAEdisp2D& edisp) | |
| 134 | +{ | |
| 135 | + // Execute only if object is not identical | |
| 136 | + if (this != &edisp) { | |
| 137 | + | |
| 138 | + // Copy base class members | |
| 139 | + this->GCTAEdisp::operator=(edisp); | |
| 140 | + | |
| 141 | + // Free members | |
| 142 | + free_members(); | |
| 143 | + | |
| 144 | + // Initialise private members | |
| 145 | + init_members(); | |
| 146 | + | |
| 147 | + // Copy members | |
| 148 | + copy_members(edisp); | |
| 149 | + | |
| 150 | + } // endif: object was not identical | |
| 151 | + | |
| 152 | + // Return this object | |
| 153 | + return *this; | |
| 154 | +} | |
| 155 | + | |
| 156 | +/***********************************************************************//** | |
| 157 | + * @brief Return energy dispersion in units of s^-1 MeV^-1 | |
| 158 | + * | |
| 159 | + * @param[in] logEobs Log10 of the measured energy (TeV). | |
| 160 | + * @param[in] logEsrc Log10 of the true photon energy (TeV). | |
| 161 | + * @param[in] theta Offset angle in camera system (rad). Defaults to 0.0. | |
| 162 | + * @param[in] phi Azimuth angle in camera system (rad). Not used in this method. | |
| 163 | + * @param[in] zenith Zenith angle in Earth system (rad). Not used in this method. | |
| 164 | + * @param[in] azimuth Azimuth angle in Earth system (rad). Not used in this method. | |
| 165 | + * | |
| 166 | + * Returns the energy dispersion in units of s^-1 MeV^-1 for a given observed energy, | |
| 167 | + * true photon energy and offset angle. | |
| 168 | + ***************************************************************************/ | |
| 169 | +double GCTAEdisp2D::operator()(const double& logEobs, | |
| 170 | + const double& logEsrc, | |
| 171 | + const double& theta, | |
| 172 | + const double& phi, | |
| 173 | + const double& zenith, | |
| 174 | + const double& azimuth) const | |
| 175 | +{ | |
| 176 | + // Initalize edisp | |
| 177 | + double edisp = 0.0; | |
| 178 | + /* | |
| 179 | + // Update cache | |
| 180 | + update(logEobs, logEsrc, theta); | |
| 181 | + | |
| 182 | + // Continue if normalization is positive | |
| 183 | + if(m_norm > 0.0) | |
| 184 | + { | |
| 185 | + // Compute edisp | |
| 186 | + edisp = std::exp(-(m_EobsOnEsrc - m_ebias)*(m_EobsOnEsrc - m_ebias) / (m_eres*m_eres)) * m_norm; | |
| 187 | + } | |
| 188 | + */ | |
| 189 | + | |
| 190 | + edisp = m_edisp(0, logEsrc, logEobs - logEsrc, theta); | |
| 191 | + | |
| 192 | + // Return | |
| 193 | + return edisp; | |
| 194 | +} | |
| 195 | + | |
| 196 | +/*========================================================================== | |
| 197 | + = = | |
| 198 | + = Public methods = | |
| 199 | + = = | |
| 200 | + ==========================================================================*/ | |
| 201 | + | |
| 202 | +/***********************************************************************//** | |
| 203 | + * @brief Clear instance | |
| 204 | + * | |
| 205 | + * This method properly resets the object to an initial state. | |
| 206 | + ***************************************************************************/ | |
| 207 | +void GCTAEdisp2D::clear(void) | |
| 208 | +{ | |
| 209 | + // Free class members (base and derived classes, derived class first) | |
| 210 | + free_members(); | |
| 211 | + this->GCTAEdisp::free_members(); | |
| 212 | + | |
| 213 | + // Initialise members | |
| 214 | + this->GCTAEdisp::init_members(); | |
| 215 | + init_members(); | |
| 216 | + | |
| 217 | + // Return | |
| 218 | + return; | |
| 219 | +} | |
| 220 | + | |
| 221 | + | |
| 222 | +/***********************************************************************//** | |
| 223 | + * @brief Clone instance | |
| 224 | + * | |
| 225 | + * @return Deep copy of energy dispersion instance. | |
| 226 | + ***************************************************************************/ | |
| 227 | +GCTAEdisp2D* GCTAEdisp2D::clone(void) const | |
| 228 | +{ | |
| 229 | + return new GCTAEdisp2D(*this); | |
| 230 | +} | |
| 231 | + | |
| 232 | + | |
| 233 | +/***********************************************************************//** | |
| 234 | + * @brief Load energy dispersion from FITS file | |
| 235 | + * | |
| 236 | + * @param[in] filename FITS file. | |
| 237 | + * | |
| 238 | + * This method loads the energy dispersion information from a FITS file. | |
| 239 | + ***************************************************************************/ | |
| 240 | +void GCTAEdisp2D::load(const std::string& filename) | |
| 241 | +{ | |
| 242 | + // Open FITS file | |
| 243 | + GFits fits(filename); | |
| 244 | + | |
| 245 | + // Read energy dispersion from file | |
| 246 | + read(fits); | |
| 247 | + | |
| 248 | + // Close FITS file | |
| 249 | + fits.close(); | |
| 250 | + | |
| 251 | + // Store filename | |
| 252 | + m_filename = filename; | |
| 253 | + | |
| 254 | + // Return | |
| 255 | + return; | |
| 256 | +} | |
| 257 | + | |
| 258 | +/***********************************************************************//** | |
| 259 | + * @brief Write CTA energy dispersion table into FITS binary table object. | |
| 260 | + * | |
| 261 | + * @param[in] hdu FITS binary table. | |
| 262 | + * | |
| 263 | + * @todo Add necessary keywords. | |
| 264 | + ***************************************************************************/ | |
| 265 | +void GCTAEdisp2D::write(GFitsBinTable& hdu) const | |
| 266 | +{ | |
| 267 | + // Write background table | |
| 268 | + m_edisp.write(hdu); | |
| 269 | + | |
| 270 | + // Return | |
| 271 | + return; | |
| 272 | +} | |
| 273 | + | |
| 274 | +/***********************************************************************//** | |
| 275 | + * @brief Save energy dispersion table into FITS file | |
| 276 | + * | |
| 277 | + * @param[in] filename Energy dispersion table FITS file name. | |
| 278 | + * @param[in] clobber Overwrite existing file? (true=yes) | |
| 279 | + * | |
| 280 | + * Save the energy dispersion table into a FITS file. | |
| 281 | + ***************************************************************************/ | |
| 282 | +void GCTAEdisp2D::save(const std::string& filename, const bool& clobber) const | |
| 283 | +{ | |
| 284 | + // Create binary table | |
| 285 | + GFitsBinTable table; | |
| 286 | + table.extname("ENERGY DISPERSION"); | |
| 287 | + | |
| 288 | + // Write the Energy dispersion table | |
| 289 | + write(table); | |
| 290 | + | |
| 291 | + // Create FITS file, append table, and write into the file | |
| 292 | + GFits fits; | |
| 293 | + fits.append(table); | |
| 294 | + fits.saveto(filename, clobber); | |
| 295 | + | |
| 296 | + // Return | |
| 297 | + return; | |
| 298 | +} | |
| 299 | + | |
| 300 | + | |
| 301 | +/***********************************************************************//** | |
| 302 | + * @brief Read energy dispersion from FITS file | |
| 303 | + * | |
| 304 | + * @param[in] fits FITS file pointer. | |
| 305 | + * | |
| 306 | + * Reads the energy dispersion form the FITS file extension "ENERGY DISPERSION". | |
| 307 | + * The data are stored in m_edisp which is of type GCTAResponseTable. The | |
| 308 | + * energy axis will be set to log10, the offset angle axis to radians. | |
| 309 | + ***************************************************************************/ | |
| 310 | +void GCTAEdisp2D::read(const GFits& fits) | |
| 311 | +{ | |
| 312 | + // Clear response table | |
| 313 | + m_edisp.clear(); | |
| 314 | + | |
| 315 | + // Get energy dispersion table | |
| 316 | + const GFitsTable& table = *fits.table("ENERGY DISPERSION"); | |
| 317 | + | |
| 318 | + // Read energy dispersion table | |
| 319 | + m_edisp.read(table); | |
| 320 | + | |
| 321 | + int i_etrue = 0; | |
| 322 | + int i_migra = 0; | |
| 323 | + int i_theta = 0; | |
| 324 | + double sum = 0.0; | |
| 325 | + double deltaE = 0.0; | |
| 326 | + | |
| 327 | + double etrue_size = e_disp.axis(0); | |
| 328 | + double migra_size = e_disp.axis(1); | |
| 329 | + double theta_size = e_disp.axis(2); | |
| 330 | + | |
| 331 | + // Normalize vectors of migration matrix for all ETRUE and THETA | |
| 332 | + for(i_etrue = 0; i_etrue < etrue_size; ++i_etrue) | |
| 333 | + { | |
| 334 | + deltaE = e_disp.axis_hi(0, i_etrue) - e_disp.axis_lo(0, i_etrue); | |
| 335 | + for(i_theta = 0; i_theta < theta_size; ++i_theta) | |
| 336 | + { | |
| 337 | + // Compute sum | |
| 338 | + sum = 0.0; | |
| 339 | + for(i_migra = 0; i_migra < migra_size; ++i_migra) | |
| 340 | + { | |
| 341 | + sum += m_edisp(0, i_etrue + (i_migra + i_theta*migra_size)*etrue_size) * deltaE; | |
| 342 | + } | |
| 343 | + if(sum > 0.0) | |
| 344 | + { | |
| 345 | + for(i_migra = 0; i_migra < migra_size; ++i_migra) | |
| 346 | + { | |
| 347 | + m_edisp(0, i_etrue + (i_migra + i_theta*migra_size)*etrue_size) *= 1.0 / sum; | |
| 348 | + } | |
| 349 | + } | |
| 350 | + } | |
| 351 | + } | |
| 352 | + | |
| 353 | + // Set true energy axis to logarithmic scale | |
| 354 | + m_edisp.axis_log10(0); | |
| 355 | + | |
| 356 | + // Set migration energy ratio axis to logarithmic scale | |
| 357 | + m_edisp.axis_log10(1); | |
| 358 | + | |
| 359 | + // Set offset angle axis to radians | |
| 360 | + m_edisp.axis_radians(2); | |
| 361 | + | |
| 362 | + // Return | |
| 363 | + return; | |
| 364 | +} | |
| 365 | + | |
| 366 | +/***********************************************************************//** | |
| 367 | + * @brief Simulate energy dispersion | |
| 368 | + * | |
| 369 | + * @param[in] ran Random number generator. | |
| 370 | + * @param[in] logE Log10 of the true photon energy (TeV). | |
| 371 | + * @param[in] theta Offset angle in camera system (rad). Not used. | |
| 372 | + * @param[in] phi Azimuth angle in camera system (rad). Not used. | |
| 373 | + * @param[in] zenith Zenith angle in Earth system (rad). Not used. | |
| 374 | + * @param[in] azimuth Azimuth angle in Earth system (rad). Not used. | |
| 375 | + * | |
| 376 | + * Draws observed energy value from a normal distribution of width | |
| 377 | + * m_par_sigma around @p logE. | |
| 378 | + ***************************************************************************/ | |
| 379 | +GEnergy GCTAEdisp2D::mc(GRan& ran, | |
| 380 | + const double& logEsrc, | |
| 381 | + const double& theta, | |
| 382 | + const double& phi, | |
| 383 | + const double& zenith, | |
| 384 | + const double& azimuth) const | |
| 385 | +{ | |
| 386 | + // Update the parameter cache | |
| 387 | + update(m_par_logEobs, logEsrc, m_par_theta); | |
| 388 | + | |
| 389 | + // Draw log observed energy in TeV | |
| 390 | + double logEobs = m_eres * ran.normal() + logEsrc; | |
| 391 | + | |
| 392 | + // Set energy | |
| 393 | + GEnergy energy; | |
| 394 | + energy.log10TeV(logEobs); | |
| 395 | + | |
| 396 | + // Return energy | |
| 397 | + return energy; | |
| 398 | +} | |
| 399 | + | |
| 400 | + | |
| 401 | +/***********************************************************************//** | |
| 402 | + * @brief Return observed energy interval that contains the energy dispersion. | |
| 403 | + * | |
| 404 | + * @param[in] logEsrc Log10 of the true photon energy (TeV). | |
| 405 | + * @param[in] theta Offset angle in camera system (rad). Not used. | |
| 406 | + * @param[in] phi Azimuth angle in camera system (rad). Not used. | |
| 407 | + * @param[in] zenith Zenith angle in Earth system (rad). Not used. | |
| 408 | + * @param[in] azimuth Azimuth angle in Earth system (rad). Not used. | |
| 409 | + * | |
| 410 | + * Returns the band of observed energies outside of which the energy | |
| 411 | + * dispersion becomes negligible for a given true energy @p logEsrc. This | |
| 412 | + * band is set to \f$\pm 5 \times \sigma\f$, where \f$\sigma\f$ is the | |
| 413 | + * Gaussian width of the energy dispersion. | |
| 414 | + ***************************************************************************/ | |
| 415 | +GEbounds GCTAEdisp2D::ebounds_obs(const double& logEsrc, | |
| 416 | + const double& theta, | |
| 417 | + const double& phi, | |
| 418 | + const double& zenith, | |
| 419 | + const double& azimuth) const | |
| 420 | +{ | |
| 421 | + // Set energy band constant | |
| 422 | + const double number_of_sigmas = 5.0; | |
| 423 | + | |
| 424 | + // Update the parameter cache | |
| 425 | + update(m_par_logEobs, logEsrc, theta); | |
| 426 | + | |
| 427 | + // Compute energy boundaries | |
| 428 | + GEnergy emin; | |
| 429 | + GEnergy emax; | |
| 430 | + emin.log10TeV(logEsrc - number_of_sigmas * m_eres); | |
| 431 | + emax.log10TeV(logEsrc + number_of_sigmas * m_eres); | |
| 432 | + | |
| 433 | + // Return energy boundaries | |
| 434 | + return (GEbounds(emin, emax)); | |
| 435 | +} | |
| 436 | + | |
| 437 | + | |
| 438 | +/***********************************************************************//** | |
| 439 | + * @brief Return true energy interval that contains the energy dispersion. | |
| 440 | + * | |
| 441 | + * @param[in] logEobs Log10 of the observed event energy (TeV). | |
| 442 | + * @param[in] theta Offset angle in camera system (rad). Not used. | |
| 443 | + * @param[in] phi Azimuth angle in camera system (rad). Not used. | |
| 444 | + * @param[in] zenith Zenith angle in Earth system (rad). Not used. | |
| 445 | + * @param[in] azimuth Azimuth angle in Earth system (rad). Not used. | |
| 446 | + * | |
| 447 | + * Returns the band of true photon energies outside of which the energy | |
| 448 | + * dispersion becomes negligible for a given observed energy @p logEobs. This | |
| 449 | + * band is set to \f$\pm 5 \times \sigma\f$, where \f$\sigma\f$ is the | |
| 450 | + * Gaussian width of the energy dispersion. | |
| 451 | + ***************************************************************************/ | |
| 452 | +GEbounds GCTAEdisp2D::ebounds_src(const double& logEobs, | |
| 453 | + const double& theta, | |
| 454 | + const double& phi, | |
| 455 | + const double& zenith, | |
| 456 | + const double& azimuth) const | |
| 457 | +{ | |
| 458 | + // Set energy band constant | |
| 459 | + const double number_of_sigmas = 5.0; | |
| 460 | + | |
| 461 | + // Update the parameter cache | |
| 462 | + update(logEobs, m_par_logEsrc, theta); | |
| 463 | + | |
| 464 | + // Compute energy boundaries | |
| 465 | + GEnergy emin; | |
| 466 | + GEnergy emax; | |
| 467 | + emin.log10TeV(logEobs - number_of_sigmas * m_eres); | |
| 468 | + emax.log10TeV(logEobs + number_of_sigmas * m_eres); | |
| 469 | + | |
| 470 | + // Return energy boundaries | |
| 471 | + return (GEbounds(emin, emax)); | |
| 472 | +} | |
| 473 | + | |
| 474 | +/***********************************************************************//** | |
| 475 | + * @brief Print energy dispersion information | |
| 476 | + * | |
| 477 | + * @param[in] chatter Chattiness (defaults to NORMAL). | |
| 478 | + * @return String containing energy dispersion information. | |
| 479 | + ***************************************************************************/ | |
| 480 | +std::string GCTAEdisp2D::print(const GChatter& chatter) const | |
| 481 | +{ | |
| 482 | + // Initialise result string | |
| 483 | + std::string result; | |
| 484 | + | |
| 485 | + // Continue only if chatter is not silent | |
| 486 | + if (chatter != SILENT) { | |
| 487 | + | |
| 488 | + // Compute energy boundaries in TeV | |
| 489 | + double emin = m_edisp.axis_lo(0,0); | |
| 490 | + double emax = m_edisp.axis_hi(0,m_edisp.axis(0)-1); | |
| 491 | + | |
| 492 | + // Compute offset angle boundaries in deg | |
| 493 | + double omin = m_edisp.axis_lo(1,0); | |
| 494 | + double omax = m_edisp.axis_hi(1,m_edisp.axis(1)-1); | |
| 495 | + | |
| 496 | + // Append header | |
| 497 | + result.append("=== GCTAEdisp2D ==="); | |
| 498 | + | |
| 499 | + // Append information | |
| 500 | + result.append("\n"+gammalib::parformat("Filename")+m_filename); | |
| 501 | + result.append("\n"+gammalib::parformat("Number of energy bins") + | |
| 502 | + gammalib::str(m_edisp.axis(0))); | |
| 503 | + result.append("\n"+gammalib::parformat("Number of offset bins") + | |
| 504 | + gammalib::str(m_edisp.axis(1))); | |
| 505 | + result.append("\n"+gammalib::parformat("Log10(Energy) range")); | |
| 506 | + result.append(gammalib::str(emin)+" - "+gammalib::str(emax)+" TeV"); | |
| 507 | + result.append("\n"+gammalib::parformat("Offset angle range")); | |
| 508 | + result.append(gammalib::str(omin)+" - "+gammalib::str(omax)+" deg"); | |
| 509 | + | |
| 510 | + } // endif: chatter was not silent | |
| 511 | + | |
| 512 | + // Return result | |
| 513 | + return result; | |
| 514 | +} | |
| 515 | + | |
| 516 | + | |
| 517 | +/*========================================================================== | |
| 518 | + = = | |
| 519 | + = Private methods = | |
| 520 | + = = | |
| 521 | + ==========================================================================*/ | |
| 522 | + | |
| 523 | +/***********************************************************************//** | |
| 524 | + * @brief Initialise class members | |
| 525 | + ***************************************************************************/ | |
| 526 | +void GCTAEdisp2D::init_members(void) | |
| 527 | +{ | |
| 528 | + // Initialise members | |
| 529 | + m_filename.clear(); | |
| 530 | + m_edisp.clear(); | |
| 531 | + m_par_logEobs = -1.0e30; | |
| 532 | + m_par_logEsrc = -1.0e30; | |
| 533 | + m_par_theta = -1.0; | |
| 534 | + m_norm = 1.0; | |
| 535 | + m_EobsOnEsrc = 0.0; | |
| 536 | + m_eres = 0.0; | |
| 537 | + m_ebias = 0.0; | |
| 538 | + | |
| 539 | + // Return | |
| 540 | + return; | |
| 541 | +} | |
| 542 | + | |
| 543 | + | |
| 544 | +/***********************************************************************//** | |
| 545 | + * @brief Copy class members | |
| 546 | + * | |
| 547 | + * @param[in] edisp Energy dispersion. | |
| 548 | + ***************************************************************************/ | |
| 549 | +void GCTAEdisp2D::copy_members(const GCTAEdisp2D& edisp) | |
| 550 | +{ | |
| 551 | + // Copy members | |
| 552 | + m_filename = edisp.m_filename; | |
| 553 | + m_edisp = edisp.m_edisp; | |
| 554 | + | |
| 555 | + // Return | |
| 556 | + return; | |
| 557 | +} | |
| 558 | + | |
| 559 | + | |
| 560 | +/***********************************************************************//** | |
| 561 | + * @brief Delete class members | |
| 562 | + ***************************************************************************/ | |
| 563 | +void GCTAEdisp2D::free_members(void) | |
| 564 | +{ | |
| 565 | + // Return | |
| 566 | + return; | |
| 567 | +} | |
| 568 | + | |
| 569 | +/***********************************************************************//** | |
| 570 | + * @brief Update Energy dispersion parameter cache | |
| 571 | + * | |
| 572 | + * @param[in] logEobs Log10 of the observed photon energy (TeV). | |
| 573 | + * @param[in] logEsrc Log10 of the true photon energy (TeV). | |
| 574 | + * @param[in] theta Offset angle in camera system (rad). | |
| 575 | + * | |
| 576 | + * This method updates the Energy dispersion parameter cache. | |
| 577 | + ***************************************************************************/ | |
| 578 | +void GCTAEdisp2D::update(const double& logEobs, const double& logEsrc, const double& theta) const | |
| 579 | +{ | |
| 580 | + /* | |
| 581 | + // Only compute Edisp parameters if arguments have changed | |
| 582 | + if (logEobs != m_par_logEobs || logEsrc != m_par_logEsrc || theta != m_par_theta) { | |
| 583 | + | |
| 584 | + // Save parameters | |
| 585 | + m_par_logEobs = logEobs; | |
| 586 | + m_par_logEsrc = logEsrc; | |
| 587 | + m_par_theta = theta; | |
| 588 | + | |
| 589 | + // Compute linear Eobs/Esrc ratio | |
| 590 | + m_EobsOnEsrc = std::exp( (logEobs - logEsrc) * std::log(10)); | |
| 591 | + | |
| 592 | + | |
| 593 | + | |
| 594 | + | |
| 595 | + // Get ERES and EBIAS values from response table | |
| 596 | + m_eres = m_edisp(0, logEsrc, theta); | |
| 597 | + m_ebias = m_edisp(1, logEsrc, theta); | |
| 598 | + | |
| 599 | + // Compute global normalization parameter | |
| 600 | + m_norm = (m_eres > 0.0) ? gammalib::inv_sqrt2pi / m_eres : 0.0; | |
| 601 | + | |
| 602 | + } | |
| 603 | + */ | |
| 604 | + // Return | |
| 605 | + return; | |
| 606 | +} | ... | ... |
inst/cta/test/test_CTA.cpp
| ... | ... | @@ -36,6 +36,7 @@ |
| 36 | 36 | #include "GTools.hpp" |
| 37 | 37 | #include "GNodeArray.hpp" |
| 38 | 38 | #include "test_CTA.hpp" |
| 39 | +#include "GCTAEdisp2D.hpp" | |
| 39 | 40 | |
| 40 | 41 | /* __ Namespaces _________________________________________________________ */ |
| 41 | 42 | |
| ... | ... | @@ -59,6 +60,7 @@ const std::string cta_irf_king = "irf_file.fits"; |
| 59 | 60 | const std::string cta_edisp_rmf = PACKAGE_SOURCE"/inst/cta/test/caldb/dc1/rmf.fits"; |
| 60 | 61 | const std::string cta_modbck_fit = datadir+"/bg_test.fits"; |
| 61 | 62 | const std::string cta_point_table = datadir+"/crab_pointing.fits.gz"; |
| 63 | +const std::string cta_irf_matrix = PACKAGE_SOURCE"inst/cta/caldb/data/cta/e/bcf/IFAE20120510_50h/irf_file_matrix.fits"; | |
| 62 | 64 | |
| 63 | 65 | |
| 64 | 66 | /***********************************************************************//** |
| ... | ... | @@ -576,30 +578,29 @@ void TestGCTAResponse::test_response_edisp2D(void) |
| 576 | 578 | test_try_failure(e); |
| 577 | 579 | } |
| 578 | 580 | |
| 579 | - // Test GCTAEdisp2D copy constructor | |
| 580 | - test_try("GCTAEdisp2D copy constructor"); | |
| 581 | - try { | |
| 582 | - GCTAEdisp2D edisp2D_2(edisp2D); | |
| 583 | - test_try_success(); | |
| 584 | - } | |
| 585 | - catch (std::exception &e) { | |
| 586 | - test_try_failure(e); | |
| 587 | - } | |
| 588 | - | |
| 589 | 581 | // Test GCTAEdisp2D file constructor |
| 590 | 582 | test_try("GCTAEdisp2D file constructor"); |
| 591 | 583 | try { |
| 592 | - GCTAEdisp2D edisp2D_3(cta_irf_king); | |
| 584 | + GCTAEdisp2D edisp2D_2(cta_irf_matrix); | |
| 593 | 585 | test_try_success(); |
| 594 | 586 | } |
| 595 | 587 | catch (std::exception &e) { |
| 596 | 588 | test_try_failure(e); |
| 597 | 589 | } |
| 598 | - | |
| 590 | + | |
| 599 | 591 | // Return |
| 600 | 592 | return; |
| 601 | 593 | } |
| 602 | 594 | |
| 595 | +/***********************************************************************//** | |
| 596 | + * @brief Test CTA Edisp Performance Table computation | |
| 597 | + ***************************************************************************/ | |
| 598 | +void test_response_edispPerfTable(void) | |
| 599 | +{ | |
| 600 | + | |
| 601 | + return; | |
| 602 | +} | |
| 603 | + | |
| 603 | 604 | |
| 604 | 605 | /***********************************************************************//** |
| 605 | 606 | * @brief Test CTA IRF computation for diffuse source model | ... | ... |
inst/cta/test/test_CTA.hpp
| ... | ... | @@ -56,6 +56,7 @@ public: |
| 56 | 56 | void test_response_edisp(void); |
| 57 | 57 | void test_response_edispRMF(void); |
| 58 | 58 | void test_response_edisp2D(void); |
| 59 | + void test_response_edispPerfTable(void); | |
| 59 | 60 | void test_response_irf_diffuse(void); |
| 60 | 61 | void test_response_npred_diffuse(void); |
| 61 | 62 | void test_response_expcube(void); | ... | ... |