Peter Deiml / gammalib

Commit 1457b1cb1da5a9b7a4c7710cf181bd0679914019

Authored by Jürgen Knödlseder
2 parents abbc9eb4 030a746c

Merge branch 'pull-request-1988-Change-Effective-Area' into integration 

# Conflicts:
#	ChangeLog
#	NEWS
Inline Side-by-side
Showing 5 changed files with 77 additions and 22 deletions
... ... @@ -4,6 +4,7 @@ Joshua Cardenzana (jcardenzana@irap.omp.eu)
4 4 Jean-Baptiste Cayrou (cayrou@etud.insa-toulouse.fr)
5 5 Christoph Deil (Christoph.Deil@mpi-hd.mpg.de)
6 6 Alex Donath (Axel.Donath@mpi-hd.mpg.de)
  7 +Stefan Eschbach (stefan.eschbach@fau.de)
7 8 Florent Forest (florent.forest@isae.fr)
8 9 Lucie Gerard (lucie.gerard@desy.de)
9 10 Tarek Hassan (thassan@ifae.es)
... ...
ChangeLog
  View file @1457b1c
... ... @@ -7,6 +7,7 @@
7 7 Add GCTAAeff::ebounds() method (#1989)
8 8 Add GEbounds::emin() and GEbounds::emax() setter methods (#1989)
9 9 Protect GEbounds::set_lin() and set_log() methods from invalid arguments (#1989)
  10 + Avoid hang-up of GCTAModelAeffBackground::mc() method (#1988)
10 11 Copy shared Python libraries instead of creating a symbolic link (#1987)
11 12 Added FITS image conversion operators (#1986)
12 13 If code originates from Git, always compile the Python bindings (#1985)
... ...
... ... @@ -192,6 +192,11 @@ An GCTAAeff::ebounds() method has been added to the abstract interface for
192 192 effective area responses to return the valid energy range of the effective
193 193 area (#1989).
194 194  
  195 +The GCTAModelAeffBackground::mc() could get stuck in the rejection loop if the
  196 +effective area returned is always zero. A counter has been added that leads to
  197 +the exit of the loop in case that 1000 subsequent zero values are returned
  198 +(#1988).
  199 +
195 200  
196 201 21. COMPTEL interface
197 202 ---------------------
... ...
inst/cta/src/GCTAModelAeffBackground.cpp
  View file @1457b1c
1 1 /***************************************************************************
2 2 * GCTAModelAeffBackground.cpp - CTA Aeff background model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2015-2016 by Michael Mayer *
  4 + * copyright (C) 2015-2017 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -509,13 +509,14 @@ GCTAEventList* GCTAModelAeffBackground::mc(const GObservation& obs,
509 509 // Retrieve CTA observation
510 510 const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>(&obs);
511 511 if (cta == NULL) {
512   - std::string msg = "Specified observation is not a CTA observation.\n" +
513   - obs.print();
  512 + std::string msg = "Specified observation is not a CTA "
  513 + "observation.\n" + obs.print();
514 514 throw GException::invalid_argument(G_MC, msg);
515 515 }
516 516  
517 517 // Get pointer on CTA IRF response
518   - const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(cta->response());
  518 + const GCTAResponseIrf* rsp =
  519 + dynamic_cast<const GCTAResponseIrf*>(cta->response());
519 520 if (rsp == NULL) {
520 521 std::string msg = "Specified observation does not contain"
521 522 " an IRF response.\n" + obs.print();
... ... @@ -534,7 +535,8 @@ GCTAEventList* GCTAModelAeffBackground::mc(const GObservation&amp; obs,
534 535 }
535 536  
536 537 // Retrieve event list to access the ROI, energy boundaries and GTIs
537   - const GCTAEventList* events = dynamic_cast<const GCTAEventList*>(obs.events());
  538 + const GCTAEventList* events =
  539 + dynamic_cast<const GCTAEventList*>(obs.events());
538 540 if (events == NULL) {
539 541 std::string msg = "No CTA event list found in observation.\n" +
540 542 obs.print();
... ... @@ -559,16 +561,21 @@ GCTAEventList* GCTAModelAeffBackground::mc(const GObservation&amp; obs,
559 561 // Set up spectral model to draw random energies from. Here we use
560 562 // a fixed energy sampling for an instance of GModelSpectralNodes.
561 563 // This is analogous to to the GCTAModelIrfBackground::mc method.
562   - GEbounds spectral_ebounds = GEbounds(m_n_mc_energies, ebounds.emin(), ebounds.emax(), true);
  564 + // We make sure that only non-negative nodes get appended.
  565 + GEbounds spectral_ebounds =
  566 + GEbounds(m_n_mc_energies, ebounds.emin(), ebounds.emax(), true);
563 567 GModelSpectralNodes spectral;
564 568 for (int i = 0; i < spectral_ebounds.size(); ++i) {
565   - GEnergy energy = spectral_ebounds.elogmean(i);
566   - double intensity = aeff_integral(obs, energy.log10TeV());
567   - double norm = m_spectral->eval(energy, events->tstart());
568   - spectral.append(energy, norm * intensity);
  569 + GEnergy energy = spectral_ebounds.elogmean(i);
  570 + double intensity = aeff_integral(obs, energy.log10TeV());
  571 + double norm = m_spectral->eval(energy, events->tstart());
  572 + double arg = norm * intensity;
  573 + if (arg > 0.0) {
  574 + spectral.append(energy, arg);
  575 + }
569 576 }
570 577  
571   - // Loop over all energy boundaries
  578 + // Loop over all energy bins
572 579 for (int ieng = 0; ieng < ebounds.size(); ++ieng) {
573 580  
574 581 // Compute the background rate in model within the energy
... ... @@ -583,6 +590,11 @@ GCTAEventList* GCTAModelAeffBackground::mc(const GObservation&amp; obs,
583 590 std::cout << "rate=" << rate << " cts/s)" << std::endl;
584 591 #endif
585 592  
  593 + // If the rate is not positive then skip this energy bins
  594 + if (rate <= 0.0) {
  595 + continue;
  596 + }
  597 +
586 598 // Loop over all good time intervals
587 599 for (int itime = 0; itime < gti.size(); ++itime) {
588 600  
... ... @@ -630,14 +642,22 @@ GCTAEventList* GCTAModelAeffBackground::mc(const GObservation&amp; obs,
630 642 ran);
631 643  
632 644 // Get maximum effective area for rejection method
633   - double max_aeff = aeff->max(energy.log10TeV(), pnt.zenith(), pnt.azimuth(), false);
  645 + double max_aeff = aeff->max(energy.log10TeV(), pnt.zenith(),
  646 + pnt.azimuth(), false);
  647 +
  648 + // Skip event if the maximum effective area is not positive
  649 + if (max_aeff <= 0.0) {
  650 + continue;
  651 + }
634 652  
635 653 // Initialise randomised coordinates
636 654 double offset = 0.0;
637 655 double phi = 0.0;
638 656  
639   - // Initialise acceptance fraction for rejection method
  657 + // Initialise acceptance fraction and counter of zeros for
  658 + // rejection method
640 659 double acceptance_fraction = 0.0;
  660 + int zeros = 0;
641 661  
642 662 // Start rejection method loop
643 663 do {
... ... @@ -649,18 +669,39 @@ GCTAEventList* GCTAModelAeffBackground::mc(const GObservation&amp; obs,
649 669 phi = ran.uniform() * gammalib::twopi;
650 670  
651 671 // Compute function value at this offset angle
652   - double value = (*aeff)(energy.log10TeV(), offset, phi, pnt.zenith(), pnt.azimuth(), false);
  672 + double value = (*aeff)(energy.log10TeV(), offset, phi,
  673 + pnt.zenith(), pnt.azimuth(),
  674 + false);
  675 +
  676 + // If the value is not positive then increment the
  677 + // zeros counter and fall through. The counter assures
  678 + // that this loop does not lock up.
  679 + if (value <= 0.0) {
  680 + zeros++;
  681 + continue;
  682 + }
  683 +
  684 + // Value is non-zero so reset the zeros counter
  685 + zeros = 0;
653 686  
654 687 // Compute acceptance fraction
655 688 acceptance_fraction = value / max_aeff;
656   -
657   - } while (ran.uniform() > acceptance_fraction);
  689 +
  690 + } while ((ran.uniform() > acceptance_fraction) &&
  691 + (zeros < 1000));
  692 +
  693 + // If the zeros counter is non-zero then the loop was
  694 + // exited due to exhaustion and the event is skipped
  695 + if (zeros > 0) {
  696 + continue;
  697 + }
658 698  
659 699 // Convert CTA pointing direction in instrument system
660 700 GCTAInstDir mc_dir(pnt.dir());
661 701  
662 702 // Rotate pointing direction by offset and azimuth angle
663   - mc_dir.dir().rotate_deg(phi * gammalib::rad2deg, offset * gammalib::rad2deg);
  703 + mc_dir.dir().rotate_deg(phi * gammalib::rad2deg,
  704 + offset * gammalib::rad2deg);
664 705  
665 706 // Compute DETX and DETY coordinates
666 707 double detx(0.0);
... ...
inst/cta/src/GCTAModelIrfBackground.cpp
  View file @1457b1c
1 1 /***************************************************************************
2 2 * GCTAModelIrfBackground.cpp - CTA IRF background model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2014-2017 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -551,13 +551,14 @@ GCTAEventList* GCTAModelIrfBackground::mc(const GObservation&amp; obs, GRan&amp; ran) co
551 551 // Retrieve CTA observation
552 552 const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>(&obs);
553 553 if (cta == NULL) {
554   - std::string msg = "Specified observation is not a CTA observation.\n" +
555   - obs.print();
  554 + std::string msg = "Specified observation is not a CTA "
  555 + "observation.\n" + obs.print();
556 556 throw GException::invalid_argument(G_MC, msg);
557 557 }
558 558  
559 559 // Get pointer on CTA IRF response
560   - const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(cta->response());
  560 + const GCTAResponseIrf* rsp =
  561 + dynamic_cast<const GCTAResponseIrf*>(cta->response());
561 562 if (rsp == NULL) {
562 563 std::string msg = "Specified observation does not contain"
563 564 " an IRF response.\n" + obs.print();
... ... @@ -576,7 +577,8 @@ GCTAEventList* GCTAModelIrfBackground::mc(const GObservation&amp; obs, GRan&amp; ran) co
576 577 }
577 578  
578 579 // Retrieve event list to access the ROI, energy boundaries and GTIs
579   - const GCTAEventList* events = dynamic_cast<const GCTAEventList*>(obs.events());
  580 + const GCTAEventList* events =
  581 + dynamic_cast<const GCTAEventList*>(obs.events());
580 582 if (events == NULL) {
581 583 std::string msg = "No CTA event list found in observation.\n" +
582 584 obs.print();
... ... @@ -618,6 +620,11 @@ GCTAEventList* GCTAModelIrfBackground::mc(const GObservation&amp; obs, GRan&amp; ran) co
618 620 std::cout << "rate=" << rate << " cts/s)" << std::endl;
619 621 #endif
620 622  
  623 + // If the rate is not positive then skip this energy bins
  624 + if (rate <= 0.0) {
  625 + continue;
  626 + }
  627 +
621 628 // Loop over all good time intervals
622 629 for (int itime = 0; itime < gti.size(); ++itime) {
623 630  
... ...