GammaLib

#include <iostream>

#include <cmath>

#include "TRolke.h"

#include <utilities/Statistics.hh>

#include <utilities/debugging.hh>

#include "SpectrumStats.hh"

using std::endl;

namespace Flux {

SpectrumStats::SpectrumStats() 

  : nOn(0),

    nOff(0),

    exposureOn(0),

    exposureOff(0),

    trueAreaTimeOn(0),

    trueAreaTimeOff(0),

    recoAreaTimeOn(0),

    recoAreaTimeOff(0),

    liveTime(0),

    offLiveTime(0),

    invRecoAreaOn(0),

    invRecoAreaOff(0),

    fPrintPrefix("\t")

{}

void SpectrumStats::Clear()

{

  nOn             = 0;

  nOff            = 0;

  exposureOn      = 0;

  exposureOff     = 0;

  trueAreaTimeOn  = 0;

  trueAreaTimeOff = 0;

  recoAreaTimeOn  = 0;

  recoAreaTimeOff = 0;

  liveTime        = 0;

  offLiveTime     = 0;

  invRecoAreaOn   = 0;

  invRecoAreaOff  = 0;

}

double SpectrumStats::GetField(std::string field) const

{

  if(field=="nOn") return nOn;

  else if(field=="exposureOn") return exposureOn;

  else if(field=="nOff") return nOff;

  else if(field=="exposureOff") return exposureOff;

  else if(field=="trueAreaTimeOn") return trueAreaTimeOn;

  else if(field=="trueAreaTimeOff") return trueAreaTimeOff;

  else if(field=="recoAreaTimeOn") return recoAreaTimeOn;

  else if(field=="recoAreaTimeOff") return recoAreaTimeOff;

  else if(field=="liveTime") return liveTime;

  else if(field=="offLiveTime") return offLiveTime;

  else if(field=="invRecoAreaOn") return invRecoAreaOn;

  else if(field=="invRecoAreaOff") return invRecoAreaOff;

  else if(field=="Alpha") return Alpha();

  else if(field=="Significance") return Significance();

  else if(field=="Excess") return Excess();

  else if(field=="ExcessError") return ExcessError();

  else if(field=="ExcessErrorUp") return ExcessErrorUp();

  else if(field=="ExcessErrorDown") return ExcessErrorDown();

  else if(field=="FracError") return FracError();

  else if(field=="FracErrorUp") return FracErrorUp();

  else if(field=="FracErrorDown") return FracErrorDown();

  else if(field=="Background") return Background();

  else{

    WARN_OUT << "Unknown field string = " << field << endl;

    return 0;

  }

}

void SpectrumStats::SetField(std::string field, double value)

{

  if(field=="nOn") nOn = value;

  else if(field=="exposureOn") exposureOn = value;

  else if(field=="nOff") nOff = value;

  else if(field=="exposureOff") exposureOff = value;

  else if(field=="trueAreaTimeOn") trueAreaTimeOn = value;

  else if(field=="trueAreaTimeOff") trueAreaTimeOff = value;

  else if(field=="recoAreaTimeOn") recoAreaTimeOn = value;

  else if(field=="recoAreaTimeOff") recoAreaTimeOff = value;

  else if(field=="liveTime") liveTime = value;

  else if(field=="offLiveTime") offLiveTime = value;

  else if(field=="invRecoAreaOn") invRecoAreaOn = value;

  else if(field=="invRecoAreaOff") invRecoAreaOff = value;

  else{

    WARN_OUT << "Unknown field string = " << field << endl;

    return;

  }

}

/** \return true only if all fields are equal */

bool SpectrumStats::Compare(const SpectrumStats& rhs) const

{

  if(nOn != rhs.nOn) return false;

  if(exposureOn != rhs.exposureOn) return false;

  if(nOff != rhs.nOff) return false;

  if(exposureOff != rhs.exposureOff) return false;

  if(trueAreaTimeOn != rhs.trueAreaTimeOn) return false;

  if(trueAreaTimeOff != rhs.trueAreaTimeOff) return false;

  if(recoAreaTimeOn != rhs.recoAreaTimeOn) return false;

  if(recoAreaTimeOff != rhs.recoAreaTimeOff) return false;

  if(offLiveTime != rhs.offLiveTime) return false;

  if(invRecoAreaOn != rhs.invRecoAreaOn) return false;

  if(invRecoAreaOff != rhs.invRecoAreaOff) return false;

  return true;

}

bool SpectrumStats::operator==(const SpectrumStats& rhs) const

{

  return Compare(rhs);

}

bool SpectrumStats::operator!=(const SpectrumStats& rhs) const

{

  return !Compare(rhs);

}

/** Performs the operation: this = this + c1*rhs.  */

void SpectrumStats::Add(const SpectrumStats& rhs)

{

  nOn             += rhs.nOn;

  nOff            += rhs.nOff;

  exposureOn      += rhs.exposureOn;

  exposureOff     += rhs.exposureOff;

  trueAreaTimeOn  += rhs.trueAreaTimeOn;

  trueAreaTimeOff += rhs.trueAreaTimeOff;

  recoAreaTimeOn  += rhs.recoAreaTimeOn;

  recoAreaTimeOff += rhs.recoAreaTimeOff;

  liveTime        += rhs.liveTime;

  offLiveTime     += rhs.offLiveTime;

  invRecoAreaOn   += rhs.invRecoAreaOn;

  invRecoAreaOff  += rhs.invRecoAreaOff;

}

/** Add the rhs object to this one, return a reference to this */

SpectrumStats& SpectrumStats::operator+=(const SpectrumStats& rhs)

{

  Add(rhs);

  return *this;

}

/** Return a new SpectrumStats object which is the sum of this one

 * and the rhs.

 */

SpectrumStats SpectrumStats::operator+(const SpectrumStats& rhs) const

{

  SpectrumStats newstats;

  newstats += *this;

  newstats += rhs;

  return newstats;

}

/**

 * Returns alpha, defined as the ratio of ON/OFF exposure.

 */

double SpectrumStats::Alpha() const

{

  return Utilities::Statistics::Alpha(exposureOn,exposureOff);

}

double SpectrumStats::AlphaUnchecked() const

{

  if (exposureOff>1e-10) 

    return exposureOn/exposureOff;

  else

    return 0;

}

/**

 * Returns the significance (Li and Ma formula) 

 */

double SpectrumStats::Significance() const

{

  double alpha = Alpha();

  return Utilities::Statistics::Significance(nOn, nOff, alpha);

}

/**

 * Same as Significance but throws no errors (for scripts)

 */

double SpectrumStats::SignificanceUnchecked() const {

  double signif=0;

  try {

    signif = Significance();

  }

  catch (Utilities::Statistics::InvalidResult &res) { 

    signif = 0;

  }

  return signif;

}

/**

 * Returns the excess counts

 */

double SpectrumStats::Excess() const

{

  double excess=0;

  try {

    excess = (double)nOn - Alpha()*(double)nOff;

  }

  catch( Utilities::Statistics::InvalidResult &res ){

    // do nothing - leave excess at 0 if alpha is not calculable

  }

  return excess;

}

/**

 * Arithmetic mean of Li & Ma up and down errors.

 *

 * Note that the nOff == 0 case is not treated separately

 * as in ExcessError(Down|Up)

 */

double SpectrumStats::ExcessError() const

{

  double error_up

    = Utilities::Statistics::LiMa_dExcess_Up(nOn,nOff,Alpha(),1.); 

  double error_down

    = Utilities::Statistics::LiMa_dExcess_Down(nOn,nOff,Alpha(),1.); 

  return (error_down + error_up) / 2;

}

/**

 * Returns the error on the excess, using Li & Ma

 */

double SpectrumStats::ExcessErrorUp() const

{

  double excesserrorup

    = Utilities::Statistics::LiMa_dExcess_Up(nOn,nOff,Alpha(),1.); 

  if (nOff < 1e-10)

    excesserrorup = Utilities::Statistics::ErrorUp(nOn);

  return fabs(excesserrorup);

}

/**

 * Returns the error on the excess, using Li & Ma

 */

double SpectrumStats::ExcessErrorDown() const

{

  double excesserrordown

    = Utilities::Statistics::LiMa_dExcess_Down(nOn,nOff,Alpha(),1.); 

  if (nOff < 1e-10)

    excesserrordown = Utilities::Statistics::ErrorDown(nOn);

  return fabs(excesserrordown);

}

/**

 * Returns the fractional error on the excess, using Li & Ma

 */

double SpectrumStats::FracErrorUp() const

{

  return fabs(ExcessErrorUp() / Excess());

}

/**

 * Returns the fractional error on the excess, using Li & Ma

 */

double SpectrumStats::FracErrorDown() const

{

  return fabs(ExcessErrorDown() / Excess());

}

/**

 * Returns the fractional error on the excess, using Li & Ma

 */

double SpectrumStats::FracError() const

{

  return fabs(ExcessError() / Excess());

}

/**

 * Returns the excess sensitivity for a desired significance level.

 */

double SpectrumStats::ExcessSensitivity(double significance) const {

  if( significance <= 0.0 ) return 0;

  double result = 0;

  try {

    result = Utilities::Statistics::LiMa_requiredExcess(significance, nOff,

                                                        AlphaUnchecked());

  } catch (Utilities::Statistics::InvalidResult error) {

    std::cerr << error.what() << std::endl;

    result = 0;

  }

  return result;

}

/**

 * Rolke confidence limits on the excess.

 * \param cl Confidence level in %

 */

void SpectrumStats::ExcessLimits(double& excess_down, double& excess_up, 

                                 double cl) const

{

  // \note Creating and deleting a new TRolke for every bin is inefficient,

  // but I don't know how else to implement it here.

  // Where speed is a problem you can just make a TRolke yourself instead

  // of calling this method for many bins.

  TRolke* rolke = new TRolke(cl / 100);

  rolke->SetBounding(true); //use method 2 of Rolke et al (2005)

  rolke->SetPoissonBkgKnownEff(nOn, nOff, 1 / Alpha(), 1);

  rolke->GetLimits(excess_down, excess_up);

  delete rolke;

}

/**

 * Rolke confidence limits on the excess.

 * \param cl Confidence level in %

 */

void SpectrumStats::ZeroExcessLimits(double& excess_down, double& excess_up, 

                                     double cl) const

{

  // \note Creating and deleting a new TRolke for every bin is inefficient,

  // but I don't know how else to implement it here.

  // Where speed is a problem you can just make a TRolke yourself instead

  // of calling this method for many bins.

  TRolke* rolke = new TRolke(cl / 100);

  rolke->SetBounding(true); //use method 2 of Rolke et al (2005)

  rolke->SetPoissonBkgKnownEff(Background(), nOff, 1 / Alpha(), 1);

  rolke->GetLimits(excess_down, excess_up);

  delete rolke;

}

double SpectrumStats::Background() const

{

  return AlphaUnchecked() * nOff;

}

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

{

  try {

    os << Form("%sN_on          = %d", 

               fPrintPrefix.c_str(), nOn) << endl;

    os << Form("%sN_off         = %d", 

               fPrintPrefix.c_str(), nOff) << endl;

    os << Form("%sExposure_on   = %.1f", 

               fPrintPrefix.c_str(), exposureOn) << endl;

    os << Form("%sExposure_off  = %.1f", 

               fPrintPrefix.c_str(), exposureOff) << endl;

    os << Form("%sAlpha         = %.2f", 

               fPrintPrefix.c_str(), Alpha()) << endl;

    os << Form("%sExcess        = %.2f", 

               fPrintPrefix.c_str(), Excess()) << endl;

    os << Form("%sSignificance  = %.2f", 

               fPrintPrefix.c_str(), Significance()) << endl;

    os << Form("%strueAT_on     = %.2e cm^2 s",

               fPrintPrefix.c_str(), trueAreaTimeOn) << endl;

    os << Form("%strueAT_off    = %.2e cm^2 s",

               fPrintPrefix.c_str(), trueAreaTimeOff) << endl;

    os << Form("%srecoAT_on     = %.2e cm^2 s", 

               fPrintPrefix.c_str(), recoAreaTimeOn) << endl;

    os << Form("%srecoAT_off    = %.2e cm^2 s", 

               fPrintPrefix.c_str(), recoAreaTimeOff) << endl;

    os << Form("%sLiveTime      = %.2f s", 

               fPrintPrefix.c_str(), liveTime) << endl;

    os << Form("%sLiveTime_off  = %.2f s", 

               fPrintPrefix.c_str(), offLiveTime) << endl;

    os << Form("%sinvArea_on    = %.2e cm^-2",

               fPrintPrefix.c_str(), invRecoAreaOn) << endl;

    os << Form("%sinvArea_off   = %.2e cm^-2",

               fPrintPrefix.c_str(), invRecoAreaOff) << endl;

  }

  catch( Utilities::Statistics::InvalidResult &res ){

    os << "# Couldn't calculate all statistics because: " << res.what() << endl;

  }

}

} // end namespace Flux

std::ostream& operator<<(std::ostream& stream,

                         const Flux::SpectrumStats& stats)

{

  stats.print(stream);

  return stream;

}

ClassImp(Flux::SpectrumStats)