GammaLib

#ifndef FLUX_SPECTRUM_HH

#define FLUX_SPECTRUM_HH

#include <TH1D.h>

#include <TH2D.h>

#include <TString.h>

#include <sash/MonitorBase.hh>

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

#include <sash/EnvelopeEntry.hh>

#include <sash/HESSArray.hh>

#endif

#include "SpectrumConstants.hh"

#include "SpectrumStats.hh"

#include "EnergyAxis.hh"

namespace Flux {

  typedef std::pair<int, int> BinRange;

  class FitFunction;

  /*! \class Flux::Spectrum

    \brief Stores count and exposure information in energy bins.

    \ingroup Flux

    \author Heidelberg team

    \todo Bkg -> Off in documentation

    Filled by SpectrumMaker, read by SpectrumFitter.

    use this to loop over bins:

    \code

    for(int bin=1; bin<=GetNbins(); bin++){ 

    ... 

    }

    \endcode

  */

  class Spectrum : public Sash::MonitorBase {

  public:

    //

    // Constructors and destructors

    //

    Spectrum(Sash::HESSArray&

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

             = hessdummy

#endif

             , const char * name ="");

    Spectrum(const char *handleFunc,

             Sash::HESSArray&

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

             = hessdummy

#endif

             , const char * name="");

    Spectrum(const Spectrum& other);

    Spectrum(const Spectrum& other, const char* newname);

    /** Clone this object.

        \todo Why doesn't the Sash::MonitorBase:Clone method actually clone an object?

    */

    virtual TObject* clone(const char* /*newname*/) const

    { return new Spectrum(*this/*,newname*/); }

    inline virtual ~Spectrum() {}

    virtual void Clear(Option_t *option="R");

    void Init();

    void InitBinning(double Emin = DEFAULT_EMIN,

                     double Emax = DEFAULT_EMAX, 

                     int nbinspdec = DEFAULT_BINS_PER_DEC, 

                     bool reinitialise = true);

    //

    // Operators

    //

    bool operator==(const Spectrum& rhs) const;

    bool operator!=(const Spectrum& rhs) const;

    void Add(const Spectrum& rhs);

    void operator+=(const Spectrum& rhs);

    Spectrum operator+(const Spectrum& rhs) const;

    //

    // GETMEMBERs

    //

    GETMEMBER(Moniker, TString);

    GETMEMBER(EnergyAxis, Flux::EnergyAxis);

    GETMEMBER(SafeBinRange, BinRange);

    GETMEMBER(Cells, std::vector<Flux::SpectrumStats>);

    GETMEMBER(EnergyResolution, TH2D);

    GETMEMBER(EnergyResolution2, TH2D);

    GETMEMBER(AreaTimeEnergyResolution, TH2D);

    GETMEMBER(AreaTimeEnergyResolution2, TH2D);

    //

    // Moniker ("title": Using Title conflicts with MonitorBase!)

    //

    void SetMoniker(TString title){ GetMoniker() = title; }

    void SetMoniker(const char* title){ GetMoniker() = TString(title); }

    //! Number of bins of axis.    

    int GetNbins() const {return GetEnergyAxis().GetNbins(); } 

    //

    // bin finders & helpers

    //

    void FindFineBinRange(const int bin, 

                          int& fine_bmin, int& fine_bmax);

    //

    // data finders & helpers

    //

    bool HasContentLess(const std::string field, const int bin, 

                        const double value = 0) const;

    bool AllHasContentLess(const std::string field, 

                           const double value = 0) const;

    void GetContentLessBinRange(const std::string field, const double value,

                                int &minbin, int &maxbin) const;

    void GetContentLessEnergyRange(const std::string field, const double value, 

                                   double &emin, double &emax) const;

    bool HasNoCounts(const int bin) const;

    bool AllHasNoCounts() const;

    void GetCountsBinRange(int &minbin, int &maxbin) const;

    void GetCountsEnergyRange(double &emin, double &emax) const;

    void GetMinLivetimeBinRange(double minLivetimeFraction, 

                                int &minbin, int &maxbin) const;

    void GetMinLivetimeEnergyRange(double minLivetimeFraction, 

                                   double &emin, double &emax) const;

    bool HasNoData(const int bin) const;

    bool AllHasNoData() const;

    void GetDataBinRange(int &minbin, int &maxbin) const;

    void GetDataEnergyRange(double &emin, double &emax) const;

    bool IsSafe(const double energy) const;

    void GetSafeEnergyRange(double &emin, double &emax) const;

    void SetSafeEnergyRange(const double emin, const double emax);

    // Clear methods

    void ClearBinsOutsideBinRange(const int bmin, const int bmax);

    void ClearBinsOutsideEnergyRange(const double emin, const double emax);

    void ClearBinsOutsideSafeRange();

    // Get range & clear methods

    void ClearAndGetDataEnergyRange(double &emin, double &emax);

    void ClearAndGetCountsEnergyRange(double &emin, double &emax);

    void ClearAndGetMinLivetimeEnergyRange(double minLivetimeFraction, 

                                           double &emin, double &emax);

    //

    // filling

    //

    void Fill(double energy, double recoarea, bool isOn = true);

    //

    // setting / getting

    //

    void SetBinContent(int bin, Flux::SpectrumStats stats);

    void SetBinContent(std::string field, int bin, double value);

    void SetBinContents(std::vector<Flux::SpectrumStats> points);

    void SetBinContents(std::string field, double value);

    void SetBinContents(std::string field, std::vector<double> values);

    void ClearBin(int bin);

    //! Reference to a bin that let's you get and set without making 

    //! local copies and two calls (get, set).

    Flux::SpectrumStats& GetCell(int bin) { return GetCells()[bin]; }

    std::vector<Flux::SpectrumStats> GetBinContents() const;

    Flux::SpectrumStats GetBinContent(int bin) const;

    std::vector<double> GetBinContents(const std::string field,

                                       const int bmin = 1, const int bmax = N_BINS) const;

    double GetBinContent(std::string field, int bin) const;

    int GetNDataBins() const;

    std::vector<double> GetEmptyVector() const;

    //

    // integration tools

    //

    Flux::SpectrumStats Integral() const;

    Flux::SpectrumStats Integral(int bmin, int bmax) const;

    //! Return integral (sum total) of all bin contents. 

    Flux::SpectrumStats GetTotalStats() const { return Integral(); }

    Flux::SpectrumStats GetTotalStats(int bmin, int bmax) const;

    Flux::SpectrumStats GetTotalStats(double emin, double emax) const;

    double Integral(std::string field) const;

    double Integral(std::string field, int bmax, int bmin) const;

    //! Return integral (sum total) of all bin contents of a field. 

    double GetStatTotal(std::string field) const { return Integral(field); }

    double GetStatTotal(std::string field, 

                        int bmin, int bmax) const;

    double GetStatTotal(std::string field, 

                        double emin, double emax) const;

    //

    // Expected Counts

    //

    double ExpectedCounts(int bin, 

                          FitFunction* model, bool integrate = DEFAULT_INTEGRATE) const;

    double ExpectedCounts(int bmin, int bmax,

                          FitFunction* model, bool integrate = DEFAULT_INTEGRATE) const;

    double ExpectedCounts(double emin, double emax, 

                          FitFunction* model, bool integrate = DEFAULT_INTEGRATE) const;

    double ExpectedCountsPL(int bin,

                            double gamma = DEFAULT_GAMMA, 

                            double norm = DEFAULT_NORM, 

                            double eref = DEFAULT_EREF, 

                            bool integrate = DEFAULT_INTEGRATE) const;

    double ExpectedCountsPL(int bmin, int bmax,

                            double gamma = DEFAULT_GAMMA, 

                            double norm = DEFAULT_NORM, 

                            double eref = DEFAULT_EREF, 

                            bool integrate = DEFAULT_INTEGRATE) const;

    double ExpectedCountsPL(double emin, double emax,

                            double gamma = DEFAULT_GAMMA, 

                            double norm = DEFAULT_NORM, 

                            double eref = DEFAULT_EREF, 

                            bool integrate = DEFAULT_INTEGRATE) const;

    //

    // Expected Hadrons

    //

    static double ExpectedHadron(int nOn, int nOff,

                                 double alpha, double expGamma);

    double ExpectedHadron(int bin, double expCounts=-1.) const;

    //

    // Flux, limits and sensitivity

    //

    void FluxAndLimits(double& flux, double& flux_down, double& flux_up,

                        double energy, const double emin, const double emax, FitFunction* model,

                        double cl = DEFAULT_CL, bool integrate = DEFAULT_INTEGRATE) const;

    void FluxAndLimitsPL(double& flux, double& flux_down, double& flux_up,

                        double energy, double emin, double emax, double gamma = DEFAULT_GAMMA,

                        double cl = DEFAULT_CL, bool integrate = DEFAULT_INTEGRATE) const;

    double Sensitivity(double energy, double emin, double emax,

                        FitFunction* model, double cl = DEFAULT_CL, bool integrate =

                                        DEFAULT_INTEGRATE) const;

    double SensitivityPL(double energy, double emin, double emax, double gamma =

                        DEFAULT_GAMMA, double cl = DEFAULT_CL,

                        bool integrate = DEFAULT_INTEGRATE) const;

    //

    // weighted quantities

    //

    // \todo Make weighting more flexible, e.g. via string argument.

    double WeightedMean(const std::vector<double> values, const int bmin = 0,

                        const int bmax = N_BINS, FitFunction* model = 0,

                        const bool weighted = true) const;

    double WeightedMean(const std::string field, const int bmin = 0,

                        const int bmax = N_BINS, FitFunction* model = 0,

                        const bool weighted = true) const;

    double WeightedMeanEnergy(int bmin = 0, int bmax = N_BINS,

                              FitFunction* model = 0, bool weighted = true) const;

    //

    // rebinning 

    //

    static void ListRebinAlgos();

    Spectrum * Rebin(const char* algorithm = DEFAULT_REBIN_ALGO.c_str(),

                     double parameter = DEFAULT_REBIN_PARAM,

                     const char* newname = "");

    Spectrum * PairBinsInRange(int bmin, int bmax, const char* newname = "");

    Spectrum * RegroupBins(int Ngroup = 1, const char* newname = "");

    //

    // Fitting

    //

    int FitFlux(FitFunction* f1,  

                Option_t* algorithm, Option_t* option);

    //

    // histograms

    //

    TH1D * MakeEmptyHist(std::string name, std::string title = "") const;

    TH1D * MakeHist(std::string field, std::string name, 

                    std::string title = "") const;

    //

    // plotting

    //

    void DiagnosticsPlots();

    TCanvas * StatisticsPlots();

    TCanvas * EnergyResolutionPlots();

    void Draw(Option_t* opt = "");

    //

    // printing

    //

    void DoPrint(std::ostream& os, Option_t* opt = "") const;

    void Print(Option_t* opt = "") const; // *MENU*

    void print(std::ostream&) const;

    void print_TotalStatsForBins(std::ostream& os, int bmin = NONE,

                                 int bmax = NONE, bool verbose = true) const;

    void print_TotalStats(std::ostream& os, double Emin = DEFAULT_EMIN,

                          double Emax = DEFAULT_EMAX, bool verbose = true) const;

    void print_Table(std::ostream& os, Option_t* opt = "") const;

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

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

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

    //

    // Convenience static functions

    //

    static Spectrum* GetEmptySpectrum(const char* name = "");

    static Spectrum* GetSpectrumFromFile(const char* filename,

                                         bool filelist = false,

                                         std::string specname = "");

    static Spectrum* FindSpectrum(std::string specname = "");

    static Spectrum* GetSpectrumNoHAP(const char* filename = "spectrum.root",

                                      const char* specname = "Spectrum");

  private : 

    // \todo Why do we need a moniker? Can't we simply use GetName()?

    TString fMoniker;              //< A moniker ("title")

    Flux::EnergyAxis fEnergyAxis;  //< The energy axis of this spectrum object

    BinRange fSafeBinRange;        //< Bin range within the safe energy thresholds.

    std::vector<Flux::SpectrumStats> fCells;  //< vector of statistics (counts and exposure) for each energy bin

    // Q: Why four resolutions and not one?

    // A: Because we are investigating / debugging this at the moment.

    //    EnergyResolution is only meaningful for each run, it is only used for debugging

    //    AreaTimeEnergyResolution is meaningful for each run and total, it is used for spectral fitting

    //    Version 1 or 2 are different formats to do the same thing.

    //    Both have advantages and disadvantages.

    //    You should get consistent results if you use either one.

    TH2D fEnergyResolution;  //< Energy resolution format 1 with y = (E_reco - E) / E

    TH2D fEnergyResolution2;  //< Energy resolution format 2 with y = log(E_reco)

    TH2D fAreaTimeEnergyResolution;  //< Area x Time x Energy resolution format 1 with y = (E_reco - E) / E

    TH2D fAreaTimeEnergyResolution2;  //< Area x Time x Energy resolution format 2 with y = log(E_reco)

    // Helper functions

    static bool getFromFile(const char* filename, Spectrum **spec, std::string specname = "");

    static bool getFromFiles(const char* filename, Spectrum **spec);

    ClassDef(Flux::Spectrum, 9);

    // 2 -> 3: new flux module

    // 3 -> 4: add fSafeEnergyRange and remove

    // fUBEnergyOn, fUBEffAreaOn, fUBEnergyOff, fUBEffAreaOff

    // 4 -> 5: add TTree* fEvents

    // 5 -> 6: change fSafeEnergyRange to fSafeBinRange to avoid off by one errors

    // 6 -> 7: replace TTree with vector and make the code simpler.

    // 7 -> 8: major clean-up

    // 8 -> 9: finally, more energy resolution histos!

  };

}// end namespace flux

std::ostream& operator<<(std::ostream& os, const Flux::Spectrum& spectrum);

#endif