GammaLib

#ifndef ROOT_BgMaps

#define ROOT_BgMaps

#include <sash/EnvelopeEntry.hh>

#include <sash/HESSArray.hh>

#include <plotters/SkyHist.hh>

namespace Background {

  // BgMaps contains fEnergyBound, which is of the type EnergyBound.

  // The first parameter is the lower energy bound of an energy range,

  // the second is the upper bound.

  typedef std::pair<Double_t,Double_t> EnergyBound;

  /**

   * \brief Data class containing all 2D background output.

   *

   * This, along with a list of RunStats objects,

   * is the primary output of a BgMaker.

   *

   * BgMaps provides:

   *

   *  - storage of "fine-binned" raw maps generated by BgMakers

   *

   *  - methods for combining these maps to produce e.g. a 

   *    significance or excess map 

   *

   *  - all the logic needed to accumulate

   * 

   *  - a series of BgMaps for each run. BgMaps may be added which do

   *  not have the same geometry (System, field-of-view, center, or

   *  binsize ) as the target, and in that case the various maps are

   *  resampled and transformed. Note that this may introduce some

   *  binning errors or coordinate distortion on small scales.

   * 

   *

   * For example, you can combine a list of BgMaps into a final

   * "summed" BgMap via the following:

   *

   * \code

   *    Background::BgMaps totals;

   *    totals.Init( center, binsize, fov_x, fov_y );

   *

   *

   *    Background::BgMaps *maps;

   *        maps = hess->Get<Background::BgMaps>("RingBgMaker_CurrentMaps"),

   *

   *    for (int ii=0; ii<RunResult->GetEntries(); ii++) {

   *       RunResult->GetEntry(ii);

   *       totals.AddMapsFrom( *maps );

   *    }

   *

   *    totals.MakeSignificanceMap( sqrt(0.05) )->Draw();

   * \endcode

   *

   * \note Exposure is the *background exposure* and is defined as:

   * Expected number of events at a particular position, given the

   * acceptance at that point, such that the integral of the exposure

   * map excluding exclusions == total number of ON events.  I.e., the

   * exposure map is equivalent to the FOV background map. The alpha

   * value (used in calculating significance and excess) is defined as

   * the OnExposureMap/OffExposureMap.  ExposureON/OFF should be

   * filled always by the BgMaker base class, with an option to use

   * real events to make the exposure maps. TemplateBgMaker will have

   * to fill these in a differnt way, since these exposure maps assume

   * gamma-ray-like events.

   *

   *  \ingroup background

   */

  class BgMaps : public Sash::MonitorBase

  {

    // used by FluxMapFillFunction

    enum FluxMapType { FLUX_POINTLIKE, UPPER_LIMIT_POINTLIKE,

                       SENSITIVITY_POINTLIKE,

                       SURFACE_BRIGHTNESS, FLUX_EXTENDED,

                       UPPER_LIMIT_EXTENDED, SENSITIVITY_EXTENDED,

                       UPPER_LIMIT_EXTENDED_AVG };

  public:

    BgMaps(Sash::HESSArray& hessarray

#if defined(__CINT__) || defined(CINTOBJECT)

            = hessdummy

#endif

            ,const char * name="");

    //this constructors allows creation of a fully initilized copy with a different name

    BgMaps(Sash::HESSArray& hessarray, const BgMaps &maps, const char * name);

    // this constructor allows to inherited classes

    BgMaps(const std::string &handle,

           Sash::HESSArray& hessarray

#if defined(__CINT__) || defined(CINTOBJECT)

            = hessdummy

#endif

            ,const char * name="");

    void Clear(Option_t *option="");

    // \todo make it virtual and const

    void print(std::ostream& os) const;

    // sub-classes may add new features to these methods

    virtual void Init(Stash::Coordinate center, double binsize,

                      double fov_x, double fov_y, bool strictlyKeepFov);

    virtual void AddMapsFrom( BgMaps &maps );

    virtual void ReplaceMapsFrom( const BgMaps &maps );

    virtual void Rebin( int factor );

    virtual void CutAwayBorders( double radius, Stash::Coordinate &centerpos );    

    static void InitSkyHist(Plotters::SkyHist& hist, const char* name,

                            Stash::Coordinate centre, double binsize, 

                            double fov_x, double fov_y, bool strictlyKeepFov);

    bool InitFluxMap();

    GETMEMBER( OnMap,          Plotters::SkyHist);

    GETMEMBER( OffMap,         Plotters::SkyHist);      

    GETMEMBER( OnExposureMap,  Plotters::SkyHist);

    GETMEMBER( OffExposureMap, Plotters::SkyHist);

    GETMEMBER( ExclusionMap,   Plotters::SkyHist);

    GETMEMBER( OffMapIsCorrelated, Bool_t);

    GETMEMBER( IsFluxMap,    Bool_t);

    GETMEMBER( ExpGammaMap,     Plotters::SkyHist);

    GETMEMBER( EnergyBounds,    EnergyBound);

    GETMEMBER( ExpGammaMapInfoString, std::string);

    GETMEMBER( AssumedRefIndex,  Double_t );

    Plotters::SkyHist* MakeSignificanceMap( double correlation_radius,

                                            bool exclusions=false,

                                            std::string name="Significance");

    Plotters::SkyHist* MakeExtendedSignificanceMap( double correlation_radius); // test

    Plotters::SkyHist* MakeExcludedSignificanceMap(double correlation_radius,

                                                   std::string name="Significance_Excluded");

    Plotters::SkyHist* MakeExcessMap();

    Plotters::SkyHist* MakeExposureCorrectedExcessMap();

    Plotters::SkyHist* MakeAlphaMap( double correlation_radius=0.0, bool exclusions = false );

    Plotters::SkyHist* MakeBackgroundMapExcessStyle();

    Plotters::SkyHist* MakeBackgroundMapSignificanceStyle(double correlation_radius);

    Plotters::SkyHist* MakePointsourceFluxMap(double E1 = 1.0, double E2 = 0., double maxSensitivity = 0.0);

    Plotters::SkyHist* MakePointsourceSensitivityMap(double threshold = 5.0, double E1 = 1.0, double E2 = 0.);

    Plotters::SkyHist* MakePointsourceUpperLimitMap( double confidenceLevel = 0.95,

                                                     double E1 = 1.0, double E2 = 0., double maxSensitivity = 0.0);

    Plotters::SkyHist* MakeSurfaceBrightnessMap(double E1 = 1.0, double E2 = 0.,

                                                double maxSensitivity = 0.0, double correlationRadius = 0.1);

    Plotters::SkyHist* MakeExtendedFluxMap(double radius, double E1 = 1.0, double E2 = 0., double maxSensitivity = 0.0);

    Plotters::SkyHist* MakeExtendedSensitivityMap(double radius, double threshold = 5.0, double E1 = 1.0, double E2 = 0.);

    Plotters::SkyHist* MakeExtendedUpperLimitMap(double radius,  double confidenceLevel = 0.95,

                                                 double E1 = 1.0, double E2 = 0., double maxSensitivity = 0.0);

    Plotters::SkyHist* MakeExtendedUpperLimitAvgMap(double radius,  double confidenceLevel = 0.95,

                                                    double E1 = 1.0, double E2 = 0.);

    void SaveAsFITS( std::string basename );

    int GetResampleFactor() const { return fResampleFactor; }

    void SetResampleFactor(int factor=-1) { fResampleFactor=factor; }

    void FillExclusionMap( Tools::Compound &exclusionregion, bool init = true );

  private:

    Plotters::SkyHist* FluxMapFillFunction(double radius, double E1, double E2,

                                           FluxMapType maptype,

                                           std::string mapname = "dummy",

                                           std::string maptitle = "dummy",

                                           std::string zaxis = "dummy",

                                           double confidence_level = 0.95);

  protected:

    Plotters::SkyHist fOnMap;  //< Event map with independent bins

    Plotters::SkyHist fOffMap; //< Event map with independent bins 

    Plotters::SkyHist fOnExposureMap;  //< ON exposure

    Plotters::SkyHist fOffExposureMap; //< OFF exposure

    /// Exclusion regions map -

    /// definition: data from a set bin (eq. 1)

    /// has been used to calculate the Off data.

    /// In constrast, the information whether

    /// the Off map (or On map) contains data

    /// is stored in the Exposure map. The ring 

    /// background is a good example for this

    /// definition. The Off map includes data

    /// at the excluded region.

    Plotters::SkyHist fExclusionMap; 

    /// Optionally filled by ExposureMapMaker

    /// calculated such that you can get

    /// Norm = (Excess / ExpGamma)

    /// ExpGamma is the expected Excess from a

    /// source with the fixed reference flux fully

    /// contained in the On region used to

    /// calculate the Excess.

    Plotters::SkyHist fExpGammaMap;

    /// Energy cuts can be applied to

    /// Excess and ExpGamma to reduce

    /// systematic errors near the safe

    /// energy thresholds:

    EnergyBound fEnergyBounds;

    std::string fExpGammaMapInfoString; //< pair=value string containing info needed for flux map calculation

    /// Spectral index assumed for filling of ExpGammaMap

    Double_t fAssumedRefIndex;

    Bool_t fOffMapIsCorrelated; /// set if fOffMap is correlated (--> RingBg)

    Bool_t fIsFluxMap;    /// flag will be set when fExpGammaMap is initialised in InitFluxMap.

    int fResampleFactor;  //! subpixel resampling factor used when adding (default -1)

    ClassDef(Background::BgMaps,6);

  };

}

#endif