Commit 00c9c2013beb3bedb8d9d8954b4cd33d55b87555
1 parent
6517c486
Change GCTAEdisp::mc() interface (#2704)
Now use GEnergy instead of the log10 of the energy.
Showing
14 changed files
with
88 additions
and
87 deletions
NEWS
| ... | ... | @@ -142,7 +142,8 @@ The following methods have been renamed: |
| 142 | 142 | The arguments for the following methods have been changed: |
| 143 | 143 | - GCTAOnOffObservation(GCTAObservation&, GModelSpatial&, GEbounds&, GEbounds&, GSkyRegions&, GSkyRegions&, bool&) |
| 144 | 144 | - GCTAModelRadial::mc(GRan&) |
| 145 | -- GCTAEdisp::operator(GEnergy&, GEnergy&, double&, double&, double&, double&) | |
| 145 | +- GCTAEdisp::operator(GEnergy&, GEnergy&, double&, double&, double&, double&) | |
| 146 | +- GCTAEdisp::mc(GRan&, GEnergy&, double&, double&, double&, double&) | |
| 146 | 147 | - GCTAResponseIrf::edisp(GEnergy&, GEnergy&, double&, double&, double&, double&) |
| 147 | 148 | - GLATAeff(GFits&, std::string&) |
| 148 | 149 | - GLATEdisp(GFits&, std::string&) | ... | ... |
inst/cta/include/GCTAEdisp.hpp
| ... | ... | @@ -71,12 +71,12 @@ public: |
| 71 | 71 | virtual std::string classname(void) const = 0; |
| 72 | 72 | virtual void load(const GFilename& filename) = 0; |
| 73 | 73 | virtual GFilename filename(void) const = 0; |
| 74 | - virtual GEnergy mc(GRan& ran, | |
| 75 | - const double& logE, | |
| 76 | - const double& theta = 0.0, | |
| 77 | - const double& phi = 0.0, | |
| 78 | - const double& zenith = 0.0, | |
| 79 | - const double& azimuth = 0.0) const = 0; | |
| 74 | + virtual GEnergy mc(GRan& ran, | |
| 75 | + const GEnergy& etrue, | |
| 76 | + const double& theta = 0.0, | |
| 77 | + const double& phi = 0.0, | |
| 78 | + const double& zenith = 0.0, | |
| 79 | + const double& azimuth = 0.0) const = 0; | |
| 80 | 80 | virtual GEbounds ebounds_obs(const double& logEsrc, |
| 81 | 81 | const double& theta = 0.0, |
| 82 | 82 | const double& phi = 0.0, | ... | ... |
inst/cta/include/GCTAEdisp2D.hpp
| ... | ... | @@ -119,12 +119,12 @@ public: |
| 119 | 119 | std::string classname(void) const; |
| 120 | 120 | void load(const GFilename& filename); |
| 121 | 121 | GFilename filename(void) const; |
| 122 | - GEnergy mc(GRan& ran, | |
| 123 | - const double& logE, | |
| 124 | - const double& theta = 0.0, | |
| 125 | - const double& phi = 0.0, | |
| 126 | - const double& zenith = 0.0, | |
| 127 | - const double& azimuth = 0.0) const; | |
| 122 | + GEnergy mc(GRan& ran, | |
| 123 | + const GEnergy& etrue, | |
| 124 | + const double& theta = 0.0, | |
| 125 | + const double& phi = 0.0, | |
| 126 | + const double& zenith = 0.0, | |
| 127 | + const double& azimuth = 0.0) const; | |
| 128 | 128 | GEbounds ebounds_obs(const double& logEsrc, |
| 129 | 129 | const double& theta = 0.0, |
| 130 | 130 | const double& phi = 0.0, | ... | ... |
inst/cta/include/GCTAEdispPerfTable.hpp
| ... | ... | @@ -88,12 +88,12 @@ public: |
| 88 | 88 | std::string classname(void) const; |
| 89 | 89 | void load(const GFilename& filename); |
| 90 | 90 | GFilename filename(void) const; |
| 91 | - GEnergy mc(GRan& ran, | |
| 92 | - const double& logE, | |
| 93 | - const double& theta = 0.0, | |
| 94 | - const double& phi = 0.0, | |
| 95 | - const double& zenith = 0.0, | |
| 96 | - const double& azimuth = 0.0) const; | |
| 91 | + GEnergy mc(GRan& ran, | |
| 92 | + const GEnergy& etrue, | |
| 93 | + const double& theta = 0.0, | |
| 94 | + const double& phi = 0.0, | |
| 95 | + const double& zenith = 0.0, | |
| 96 | + const double& azimuth = 0.0) const; | |
| 97 | 97 | GEbounds ebounds_obs(const double& logEsrc, |
| 98 | 98 | const double& theta = 0.0, |
| 99 | 99 | const double& phi = 0.0, | ... | ... |
inst/cta/include/GCTAEdispRmf.hpp
| ... | ... | @@ -84,11 +84,11 @@ public: |
| 84 | 84 | void load(const GFilename& filename); |
| 85 | 85 | GFilename filename(void) const; |
| 86 | 86 | GEnergy mc(GRan& ran, |
| 87 | - const double& logEsrc, | |
| 88 | - const double& theta = 0.0, | |
| 89 | - const double& phi = 0.0, | |
| 90 | - const double& zenith = 0.0, | |
| 91 | - const double& azimuth = 0.0) const; | |
| 87 | + const GEnergy& etrue, | |
| 88 | + const double& theta = 0.0, | |
| 89 | + const double& phi = 0.0, | |
| 90 | + const double& zenith = 0.0, | |
| 91 | + const double& azimuth = 0.0) const; | |
| 92 | 92 | GEbounds ebounds_obs(const double& logEsrc, |
| 93 | 93 | const double& theta = 0.0, |
| 94 | 94 | const double& phi = 0.0, | ... | ... |
inst/cta/pyext/GCTAEdisp.i
| ... | ... | @@ -56,12 +56,12 @@ public: |
| 56 | 56 | virtual std::string classname(void) const = 0; |
| 57 | 57 | virtual void load(const GFilename& filename) = 0; |
| 58 | 58 | virtual GFilename filename(void) const = 0; |
| 59 | - virtual GEnergy mc(GRan& ran, | |
| 60 | - const double& logE, | |
| 61 | - const double& theta = 0.0, | |
| 62 | - const double& phi = 0.0, | |
| 63 | - const double& zenith = 0.0, | |
| 64 | - const double& azimuth = 0.0) const = 0; | |
| 59 | + virtual GEnergy mc(GRan& ran, | |
| 60 | + const GEnergy& etrue, | |
| 61 | + const double& theta = 0.0, | |
| 62 | + const double& phi = 0.0, | |
| 63 | + const double& zenith = 0.0, | |
| 64 | + const double& azimuth = 0.0) const = 0; | |
| 65 | 65 | virtual GEbounds ebounds_obs(const double& logEsrc, |
| 66 | 66 | const double& theta = 0.0, |
| 67 | 67 | const double& phi = 0.0, | ... | ... |
inst/cta/pyext/GCTAEdisp2D.i
| ... | ... | @@ -57,12 +57,12 @@ public: |
| 57 | 57 | std::string classname(void) const; |
| 58 | 58 | void load(const GFilename& filename); |
| 59 | 59 | GFilename filename(void) const; |
| 60 | - GEnergy mc(GRan& ran, | |
| 61 | - const double& logE, | |
| 62 | - const double& theta = 0.0, | |
| 63 | - const double& phi = 0.0, | |
| 64 | - const double& zenith = 0.0, | |
| 65 | - const double& azimuth = 0.0) const; | |
| 60 | + GEnergy mc(GRan& ran, | |
| 61 | + const GEnergy& etrue, | |
| 62 | + const double& theta = 0.0, | |
| 63 | + const double& phi = 0.0, | |
| 64 | + const double& zenith = 0.0, | |
| 65 | + const double& azimuth = 0.0) const; | |
| 66 | 66 | GEbounds ebounds_obs(const double& logEsrc, |
| 67 | 67 | const double& theta = 0.0, |
| 68 | 68 | const double& phi = 0.0, | ... | ... |
inst/cta/pyext/GCTAEdispPerfTable.i
| ... | ... | @@ -57,12 +57,12 @@ public: |
| 57 | 57 | std::string classname(void) const; |
| 58 | 58 | void load(const GFilename& filename); |
| 59 | 59 | GFilename filename(void) const; |
| 60 | - GEnergy mc(GRan& ran, | |
| 61 | - const double& logE, | |
| 62 | - const double& theta = 0.0, | |
| 63 | - const double& phi = 0.0, | |
| 64 | - const double& zenith = 0.0, | |
| 65 | - const double& azimuth = 0.0) const; | |
| 60 | + GEnergy mc(GRan& ran, | |
| 61 | + const GEnergy& etrue, | |
| 62 | + const double& theta = 0.0, | |
| 63 | + const double& phi = 0.0, | |
| 64 | + const double& zenith = 0.0, | |
| 65 | + const double& azimuth = 0.0) const; | |
| 66 | 66 | GEbounds ebounds_obs(const double& logEsrc, |
| 67 | 67 | const double& theta = 0.0, |
| 68 | 68 | const double& phi = 0.0, | ... | ... |
inst/cta/pyext/GCTAEdispRmf.i
| ... | ... | @@ -57,12 +57,12 @@ public: |
| 57 | 57 | std::string classname(void) const; |
| 58 | 58 | void load(const GFilename& filename); |
| 59 | 59 | GFilename filename(void) const; |
| 60 | - GEnergy mc(GRan& ran, | |
| 61 | - const double& logE, | |
| 62 | - const double& theta = 0.0, | |
| 63 | - const double& phi = 0.0, | |
| 64 | - const double& zenith = 0.0, | |
| 65 | - const double& azimuth = 0.0) const; | |
| 60 | + GEnergy mc(GRan& ran, | |
| 61 | + const GEnergy& etrue, | |
| 62 | + const double& theta = 0.0, | |
| 63 | + const double& phi = 0.0, | |
| 64 | + const double& zenith = 0.0, | |
| 65 | + const double& azimuth = 0.0) const; | |
| 66 | 66 | GEbounds ebounds_obs(const double& logEsrc, |
| 67 | 67 | const double& theta = 0.0, |
| 68 | 68 | const double& phi = 0.0, | ... | ... |
inst/cta/src/GCTAEdisp2D.cpp
| ... | ... | @@ -467,7 +467,7 @@ void GCTAEdisp2D::save(const GFilename& filename, const bool& clobber) const |
| 467 | 467 | * @brief Simulate energy dispersion |
| 468 | 468 | * |
| 469 | 469 | * @param[in] ran Random number generator. |
| 470 | - * @param[in] logEsrc Log10 of the true photon energy (\f$\log_{10}\f$ TeV). | |
| 470 | + * @param[in] etrue True photon energy. | |
| 471 | 471 | * @param[in] theta Offset angle in camera system (radians). |
| 472 | 472 | * @param[in] phi Azimuth angle in camera system (radians). Not used. |
| 473 | 473 | * @param[in] zenith Zenith angle in Earth system (radians). Not used. |
| ... | ... | @@ -482,23 +482,21 @@ void GCTAEdisp2D::save(const GFilename& filename, const bool& clobber) const |
| 482 | 482 | * angle @p theta. If no energy dispersion information is available the |
| 483 | 483 | * method will return the true photon energy. |
| 484 | 484 | ***************************************************************************/ |
| 485 | -GEnergy GCTAEdisp2D::mc(GRan& ran, | |
| 486 | - const double& logEsrc, | |
| 487 | - const double& theta, | |
| 488 | - const double& phi, | |
| 489 | - const double& zenith, | |
| 490 | - const double& azimuth) const | |
| 485 | +GEnergy GCTAEdisp2D::mc(GRan& ran, | |
| 486 | + const GEnergy& etrue, | |
| 487 | + const double& theta, | |
| 488 | + const double& phi, | |
| 489 | + const double& zenith, | |
| 490 | + const double& azimuth) const | |
| 491 | 491 | { |
| 492 | 492 | // Make sure that energy dispersion is online |
| 493 | 493 | fetch(); |
| 494 | 494 | |
| 495 | - // Initialise energies | |
| 496 | - GEnergy etrue; | |
| 495 | + // Initialise reconstructed photon energy | |
| 497 | 496 | GEnergy ereco; |
| 498 | - etrue.log10TeV(logEsrc); | |
| 499 | 497 | |
| 500 | 498 | // Get boundaries for observed energy |
| 501 | - GEbounds ebounds = ebounds_obs(logEsrc, theta, phi, zenith, azimuth); | |
| 499 | + GEbounds ebounds = ebounds_obs(etrue.log10TeV(), theta, phi, zenith, azimuth); | |
| 502 | 500 | double emin = ebounds.emin().log10TeV(); |
| 503 | 501 | double emax = ebounds.emax().log10TeV(); |
| 504 | 502 | |
| ... | ... | @@ -506,9 +504,8 @@ GEnergy GCTAEdisp2D::mc(GRan& ran, |
| 506 | 504 | // energy (this means that no energy dispersion information is available |
| 507 | 505 | // for that energy) |
| 508 | 506 | if (emin >= emax) { |
| 509 | - double eng = std::pow(10.0, logEsrc); | |
| 510 | 507 | std::string msg = "No valid energy dispersion information available " |
| 511 | - "for true photon energy "+gammalib::str(eng)+" TeV," | |
| 508 | + "for true photon energy "+etrue.print()+"," | |
| 512 | 509 | "offset angle "+ |
| 513 | 510 | gammalib::str(theta*gammalib::rad2deg)+" deg " |
| 514 | 511 | "and azimuth angle "+ |
| ... | ... | @@ -527,7 +524,7 @@ GEnergy GCTAEdisp2D::mc(GRan& ran, |
| 527 | 524 | |
| 528 | 525 | // Initialise rejection method |
| 529 | 526 | double ewidth = emax - emin; |
| 530 | - double logEobs = logEsrc; | |
| 527 | + double logEobs = etrue.log10TeV(); | |
| 531 | 528 | double f = 0.0; |
| 532 | 529 | double ftest = 1.0; |
| 533 | 530 | int zeros = 0; |
| ... | ... | @@ -549,10 +546,9 @@ GEnergy GCTAEdisp2D::mc(GRan& ran, |
| 549 | 546 | if (f == 0.0) { |
| 550 | 547 | zeros++; |
| 551 | 548 | if (zeros > max_subsequent_zeros) { |
| 552 | - double eng = std::pow(10.0, logEsrc); | |
| 553 | 549 | std::string msg = "No valid energy dispersion information " |
| 554 | 550 | "available for true photon energy "+ |
| 555 | - gammalib::str(eng)+" TeV, offset angle "+ | |
| 551 | + etrue.print()+", offset angle "+ | |
| 556 | 552 | gammalib::str(theta*gammalib::rad2deg)+ |
| 557 | 553 | " deg and azimuth angle "+ |
| 558 | 554 | gammalib::str(phi*gammalib::rad2deg)+ | ... | ... |
inst/cta/src/GCTAEdispPerfTable.cpp
| ... | ... | @@ -322,7 +322,7 @@ void GCTAEdispPerfTable::load(const GFilename& filename) |
| 322 | 322 | * @brief Simulate energy dispersion |
| 323 | 323 | * |
| 324 | 324 | * @param[in] ran Random number generator. |
| 325 | - * @param[in] logEsrc Log10 of the true photon energy (\f$\log_{10}\f$ TeV). | |
| 325 | + * @param[in] etrue True photon energy. | |
| 326 | 326 | * @param[in] theta Offset angle in camera system (radians). Not used. |
| 327 | 327 | * @param[in] phi Azimuth angle in camera system (radians). Not used. |
| 328 | 328 | * @param[in] zenith Zenith angle in Earth system (radians). Not used. |
| ... | ... | @@ -332,13 +332,16 @@ void GCTAEdispPerfTable::load(const GFilename& filename) |
| 332 | 332 | * Draws observed energy value from a normal distribution of width |
| 333 | 333 | * m_par_sigma around @p logE. |
| 334 | 334 | ***************************************************************************/ |
| 335 | -GEnergy GCTAEdispPerfTable::mc(GRan& ran, | |
| 336 | - const double& logEsrc, | |
| 337 | - const double& theta, | |
| 338 | - const double& phi, | |
| 339 | - const double& zenith, | |
| 340 | - const double& azimuth) const | |
| 335 | +GEnergy GCTAEdispPerfTable::mc(GRan& ran, | |
| 336 | + const GEnergy& etrue, | |
| 337 | + const double& theta, | |
| 338 | + const double& phi, | |
| 339 | + const double& zenith, | |
| 340 | + const double& azimuth) const | |
| 341 | 341 | { |
| 342 | + // Get log10 TeV of true photon energy | |
| 343 | + double logEsrc = etrue.log10TeV(); | |
| 344 | + | |
| 342 | 345 | // Update the parameter cache |
| 343 | 346 | update(logEsrc); |
| 344 | 347 | ... | ... |
inst/cta/src/GCTAEdispRmf.cpp
| ... | ... | @@ -287,7 +287,7 @@ void GCTAEdispRmf::load(const GFilename& filename) |
| 287 | 287 | * @brief Simulate energy dispersion |
| 288 | 288 | * |
| 289 | 289 | * @param[in] ran Random number generator. |
| 290 | - * @param[in] logEsrc Log10 of the true photon energy (\f$\log_{10}\f$ TeV). | |
| 290 | + * @param[in] etrue True photon energy. | |
| 291 | 291 | * @param[in] theta Offset angle in camera system (radians). Not used. |
| 292 | 292 | * @param[in] phi Azimuth angle in camera system (radians). Not used. |
| 293 | 293 | * @param[in] zenith Zenith angle in Earth system (radians). Not used. |
| ... | ... | @@ -296,20 +296,18 @@ void GCTAEdispRmf::load(const GFilename& filename) |
| 296 | 296 | * |
| 297 | 297 | * Draws observed energy value from RMF matrix. |
| 298 | 298 | ***************************************************************************/ |
| 299 | -GEnergy GCTAEdispRmf::mc(GRan& ran, | |
| 300 | - const double& logEsrc, | |
| 301 | - const double& theta, | |
| 302 | - const double& phi, | |
| 303 | - const double& zenith, | |
| 304 | - const double& azimuth) const | |
| 299 | +GEnergy GCTAEdispRmf::mc(GRan& ran, | |
| 300 | + const GEnergy& etrue, | |
| 301 | + const double& theta, | |
| 302 | + const double& phi, | |
| 303 | + const double& zenith, | |
| 304 | + const double& azimuth) const | |
| 305 | 305 | { |
| 306 | 306 | // Initialise energies |
| 307 | - GEnergy etrue; | |
| 308 | 307 | GEnergy ereco; |
| 309 | - etrue.log10TeV(logEsrc); | |
| 310 | 308 | |
| 311 | 309 | // Get boundaries for rejection method |
| 312 | - GEbounds ebounds = ebounds_obs(logEsrc, theta, phi, zenith, azimuth); | |
| 310 | + GEbounds ebounds = ebounds_obs(etrue.log10TeV(), theta, phi, zenith, azimuth); | |
| 313 | 311 | double emin = ebounds.emin().log10TeV(); |
| 314 | 312 | double emax = ebounds.emax().log10TeV(); |
| 315 | 313 | ... | ... |
inst/cta/src/GCTAResponseIrf.cpp
| ... | ... | @@ -578,8 +578,10 @@ GEbounds GCTAResponseIrf::ebounds(const GEnergy& obsEnergy) const |
| 578 | 578 | * @todo Set polar angle phi of photon in camera system |
| 579 | 579 | * @todo Implement energy dispersion |
| 580 | 580 | ***************************************************************************/ |
| 581 | -GCTAEventAtom* GCTAResponseIrf::mc(const double& area, const GPhoton& photon, | |
| 582 | - const GObservation& obs, GRan& ran) const | |
| 581 | +GCTAEventAtom* GCTAResponseIrf::mc(const double& area, | |
| 582 | + const GPhoton& photon, | |
| 583 | + const GObservation& obs, | |
| 584 | + GRan& ran) const | |
| 583 | 585 | { |
| 584 | 586 | // Initialise event |
| 585 | 587 | GCTAEventAtom* event = NULL; |
| ... | ... | @@ -639,7 +641,8 @@ GCTAEventAtom* GCTAResponseIrf::mc(const double& area, const GPhoton& photon, |
| 639 | 641 | // Set measured photon energy |
| 640 | 642 | GEnergy energy = photon.energy(); |
| 641 | 643 | if (use_edisp()) { |
| 642 | - energy = edisp()->mc(ran, srcLogEng, theta, phi, zenith, azimuth); | |
| 644 | + energy = edisp()->mc(ran, photon.energy(), theta, phi, | |
| 645 | + zenith, azimuth); | |
| 643 | 646 | } |
| 644 | 647 | |
| 645 | 648 | // Allocate event | ... | ... |
inst/cta/test/test_CTA.cpp
| ... | ... | @@ -917,7 +917,7 @@ void TestGCTAResponse::test_response_edisp_PerfTable(void) |
| 917 | 917 | |
| 918 | 918 | // Test mc method |
| 919 | 919 | GRan ran; |
| 920 | - test_value(edisp2.mc(ran, 0.0).TeV(), 0.88667418590027); | |
| 920 | + test_value(edisp2.mc(ran, etrue).TeV(), 0.88667418590027); | |
| 921 | 921 | |
| 922 | 922 | // Test normalisation |
| 923 | 923 | test_edisp_integration(edisp2); |
| ... | ... | @@ -973,7 +973,7 @@ void TestGCTAResponse::test_response_edisp_RMF(void) |
| 973 | 973 | |
| 974 | 974 | // Test mc method |
| 975 | 975 | GRan ran; |
| 976 | - test_value(edisp2.mc(ran, 0.0).TeV(), 0.930126476703161); | |
| 976 | + test_value(edisp2.mc(ran, etrue).TeV(), 0.930126476703161); | |
| 977 | 977 | |
| 978 | 978 | // Test normalisation |
| 979 | 979 | test_edisp_integration(edisp2, 0.5, 10.0); |
| ... | ... | @@ -1028,7 +1028,7 @@ void TestGCTAResponse::test_response_edisp_2D(void) |
| 1028 | 1028 | |
| 1029 | 1029 | // Test mc method |
| 1030 | 1030 | GRan ran; |
| 1031 | - test_value(edisp2.mc(ran, 0.0).TeV(), 0.872387491302767); | |
| 1031 | + test_value(edisp2.mc(ran, etrue).TeV(), 0.872387491302767); | |
| 1032 | 1032 | |
| 1033 | 1033 | // Test normalisation |
| 1034 | 1034 | test_edisp_integration(edisp2); | ... | ... |