GammaLib

/*

*/

#include <TH1D.h>

#include <TString.h>

#include <TCanvas.h>

#include <TVector2.h>

#include <iostream>

#include <stdexcept>

#include <sash/HESSArray.hh>

#include <sash/MonitorBase.hh>

#include <sash/EnvelopeEntry.hh>

#include <tools/EventStats.hh>

#include "BgStats.hh"

#define DEBUG 0 // change to 1 to enable debugging here. Then, use the

                // DEBUG_OUT macro defined in "debugging.hh" instead

                // of "std::cout" to write out debug info.

#include <utilities/debugging.hh>

using std::cout;

using std::endl;

ENVELOPE_ENTRY(Sash::HESSArray,Background::BgStats);

Background::BgStats::BgStats(Sash::HESSArray& h, const char * name) 

  : MonitorBase("Handle<Background::BgStats>",name,Sash::MonitorBase::NonSplit),

    fRateStats(),

    fOnRegion(),

    fOffRegion(),

    fOnExclusionRegion(),

    fOffExclusionRegion(),

    fIdNumber(0),

    fNumberOfTelescopes(0),

    fTelescopePattern(0),

    fSkippedFlag(true),

    fBadFlag(false),

    fOnFlag(false),

    fOffFlag(false),

    fNumOn(0),

    fNumOff(0),

    fNumSkippedOn(0),

    fNumSkippedOff(0),

    fNumBadOn(0),

    fNumBadOff(0),

    fZenithDistOn(),

    fZenithDistOff(),

    fAzimuthDistOn(),

    fAzimuthDistOff(),

    fOffsetDistOn(),

    fOffsetDistOff(),

    fZenithOffsetDistOn(),

    fZenithOffsetDistOff(),

    fIsCumulative(false),

    fObservationPosition(0,0,Stash::System::GetDummySystem()),

    fMuonEntry(0),

    fArrayMuonParameters(h),

    fConfigName("")

{

  SetHistNames(name);

}

Background::BgStats::BgStats(Sash::HESSArray& h, const BgStats &stats, const char * name) 

  : MonitorBase("Handle<Background::BgStats>",name,Sash::MonitorBase::NonSplit),

    fObservationPosition(0,0,Stash::System::GetDummySystem()),

    fArrayMuonParameters(h)

{

  ReplaceStatsFrom(stats,false);

}

Background::BgStats::BgStats(const char* handleFunc, Sash::HESSArray& h, const char * name) 

  : MonitorBase(handleFunc,name,Sash::MonitorBase::NonSplit),

    fRateStats(),

    fOnRegion(),

    fOffRegion(),

    fOnExclusionRegion(),

    fOffExclusionRegion(),

    fIdNumber(0),

    fNumberOfTelescopes(0),

    fTelescopePattern(0),

    fSkippedFlag(true),

    fBadFlag(false),

    fOnFlag(false),

    fOffFlag(false),

    fNumOn(0),

    fNumOff(0),

    fNumSkippedOn(0),

    fNumSkippedOff(0),

    fNumBadOn(0),

    fNumBadOff(0),

    fZenithDistOn(),

    fZenithDistOff(),

    fAzimuthDistOn(),

    fAzimuthDistOff(),

    fOffsetDistOn(),

    fOffsetDistOff(),

    fZenithOffsetDistOn(),

    fZenithOffsetDistOff(),

    fIsCumulative(false),

    fObservationPosition(0,0,Stash::System::GetDummySystem()),

    fMuonEntry(0),

    fArrayMuonParameters(h),

    fConfigName("")

{

  SetHistNames(name);

}

Background::BgStats::~BgStats() {

  DEBUG_OUT << "GOODBYE" << endl;

}

void Background::BgStats::SetHistNames(const char* name)

{

  std::string Name = name;

  GetOnRegion().SetName((Name+"_OnRegion").c_str());

  GetOffRegion().SetName((Name+"_OffRegion").c_str());

  GetOnExclusionRegion().SetName((Name+"_OnExclusionRegion").c_str());

  GetOffExclusionRegion().SetName((Name+"_OffExclusionRegion").c_str());

  TH1::AddDirectory(kFALSE);

  GetZenithDistOn().ResetBit(kCanDelete);

  GetZenithDistOff().ResetBit(kCanDelete);

  GetAzimuthDistOn().ResetBit(kCanDelete);

  GetAzimuthDistOff().ResetBit(kCanDelete);

  GetOffsetDistOn().ResetBit(kCanDelete);

  GetOffsetDistOff().ResetBit(kCanDelete);

  GetZenithOffsetDistOn().ResetBit(kCanDelete);

  GetZenithOffsetDistOff().ResetBit(kCanDelete);

  GetZenithDistOn().ResetBit(kMustCleanup);

  GetZenithDistOff().ResetBit(kMustCleanup);

  GetAzimuthDistOn().ResetBit(kMustCleanup);

  GetAzimuthDistOff().ResetBit(kMustCleanup);

  GetOffsetDistOn().ResetBit(kMustCleanup);

  GetOffsetDistOff().ResetBit(kMustCleanup);

  GetZenithOffsetDistOn().ResetBit(kMustCleanup);

  GetZenithOffsetDistOff().ResetBit(kMustCleanup);

  GetZenithDistOn().SetName((Name+"_ZenithDistOn").c_str());

  GetZenithDistOn().SetTitle("ON Zenith Distribution;zenith angle / deg");

  GetZenithDistOff().SetName((Name+ "_ZenithDistOff").c_str());

  GetZenithDistOff().SetTitle("OFF Zenith Distribution;zenith angle / deg");

  GetAzimuthDistOn().SetName((Name+"_AzimuthDistOn").c_str());

  GetAzimuthDistOn().SetTitle("ON Azimuth Distribution;azimuth angle / deg");

  GetAzimuthDistOff().SetName((Name+ "_AzimuthDistOff").c_str());

  GetAzimuthDistOff().SetTitle("OFF Azimuth Distribution;azimuth angle / deg");

  GetOffsetDistOn().SetName((Name+ "_OffsetDistOn").c_str());

  GetOffsetDistOn().SetTitle("ON Offset Distribution;offset / deg");

  GetOffsetDistOff().SetName((Name+ "_OffsetDistOff").c_str());

  GetOffsetDistOff().SetTitle("OFF Offset Distribution;offset / deg");

  GetZenithOffsetDistOn().SetName((Name+ "_ZenithOffsetDistOn").c_str());

  GetZenithOffsetDistOn().SetTitle("ON ZenithOffset Distribution;zenith angle / deg;offset / deg");

  GetZenithOffsetDistOff().SetName((Name+ "_ZenithOffsetDistOff").c_str());

  GetZenithOffsetDistOff().SetTitle("OFF ZenithOffset Distribution;zenith angle / deg;offset / deg");

  TH1::AddDirectory(kTRUE);

  return;

}

/**

 * clear all values in the BgStats

 */

void Background::BgStats::Clear(Option_t *option)

{

  DEBUG_OUT << "Clearing " << GetName() << " with option " << std::string(option) << std::endl;

  if(std::string(option) == std::string("online")) {

  } else {      

    GetOnRegion().Reset();

    GetOffRegion().Reset();

    GetOnExclusionRegion().Reset();

    GetOffExclusionRegion().Reset();

    GetIdNumber() = 0;

    GetNumberOfTelescopes() = 0;

    GetTelescopePattern() = 0;      

    GetOnFlag() = false;

    GetOffFlag() = false;

    GetMuonEntry() = 0;

    GetArrayMuonParameters().Clear();

    GetConfigName() = "";

  }

  GetRateStats().Clear(option);

  GetSkippedFlag() = false;

  GetBadFlag() = false;

  GetNumOn() = 0;

  GetNumOff() = 0;

  GetNumSkippedOn() = 0;

  GetNumSkippedOff() = 0;

  GetNumBadOn() = 0;

  GetNumBadOff() = 0;

  GetZenithDistOn().Reset();

  GetZenithDistOff().Reset();

  GetAzimuthDistOn().Reset();

  GetAzimuthDistOff().Reset();

  GetOffsetDistOn().Reset();

  GetOffsetDistOff().Reset();

  GetZenithOffsetDistOn().Reset();

  GetZenithOffsetDistOff().Reset();

  GetIsCumulative() = false;

  DEBUG_OUT << "ZenDistOff: " << GetZenithDistOff().GetName()

            << " entries: "<< GetZenithDistOff().GetEntries() << endl;

  return;

}

void Background::BgStats::print(std::ostream& os) const

{

  os << "[" << GetName() << "]" << std::endl;

  if (GetIsCumulative() == false) {

    os << "\tId            = " << GetIdNumber() << std::endl;

    os << "\tNtel          = " << GetNumberOfTelescopes() << std::endl;

    os << "\tTelPattern    = " << GetTelescopePattern() << std::endl;

    os << "\tMuonEntry     = " << GetMuonEntry() << std::endl;

    os << "\tArrayMuonPars = " << GetArrayMuonParameters() << std::endl;

    os << "\tIsOnRun       = " << GetOnFlag() << std::endl;

    os << "\tIsOffRun      = " << GetOffFlag() << std::endl;

  }  

  else {

    os << "\tNumOn         = " << GetNumOn() << " runs"<<std::endl

       << "\tNumSkippedOn  = " << GetNumSkippedOn() << " runs" << std::endl

       << "\tNumBadOn      = " << GetNumBadOn() << " runs"<< std::endl;

    os << "\tNumOff        = " << GetNumOff() << " runs" << std::endl

       << "\tNumSkippedOff = " << GetNumSkippedOff() << " runs" << std::endl

       << "\tNumBadOff     = " << GetNumBadOff() << " runs" << std::endl;

  }

  if (GetSkippedFlag())

    os << "\tSkipped       = true"<< std::endl;

  else 

    os << "\tSkipped       = false"<< std::endl;

  if (GetBadFlag())

    os << "\tBad           = true"<< std::endl;

  else 

    os << "\tBad           = false"<< std::endl;

  if (GetSkippedFlag() == false){

    os << "\tConfigName    = " << GetConfigName() << std::endl;

    os << "\tObservationPos= " << GetObservationPosition() << std::endl;

    os << "\tMeanZenithOn  = " << GetMeanZenithOn() << " deg" << std::endl;

    os << "\tMeanZenithOff = " << GetMeanZenithOff() << " deg" << std::endl;

    os << "\tMeanOffsetOn  = " << GetMeanOffsetOn() << " deg" << std::endl;

    os << "\tMeanOffsetOff = " << GetMeanOffsetOff() << " deg" << std::endl;

    os << "\tMeanAzimuthOn = " << GetMeanAzimuthOn() << " deg" << std::endl;

    os << "\tMeanAzimuthOff= " << GetMeanAzimuthOff() << " deg" << std::endl;

    os << GetRateStats() << std::endl;

  }

}

/**

 * Sum up stats from another BgStats object, making this one a

 * cumulative BgStats.

 *

 * \todo: check if ON region is the same and throw exception if it's

 * not. Also, "or" together the EXCLUSION regions.

 *

 * \todo: sum OFF regions too

 */

void

Background::BgStats::

AddStatsFrom(const BgStats &stats, bool online){

  DEBUG_OUT << "adding stats from: "<< stats.GetName() 

            << endl; 

  GetIsCumulative() = true;

  GetSkippedFlag() = false;

  // Don't accumulate from skipped or bad runs:

  if (stats.GetSkippedFlag() == true) {

    if (stats.GetOnFlag() == true )   GetNumSkippedOn()++;

    if (stats.GetOffFlag() == true )   GetNumSkippedOff()++;

    DEBUG_OUT << "Run "<< stats.GetIdNumber() << " was skipped" << endl;

    return;

  }

  if (stats.GetBadFlag() == true) {

    if (stats.GetOnFlag() == true )   GetNumBadOn()++;

    if (stats.GetOffFlag() == true )  GetNumBadOff()++;

    DEBUG_OUT << "Run "<< stats.GetIdNumber() << " was skipped as BAD" << endl;

    return;

  }

  // count number of accumulated on and off runs

  if (stats.GetOnFlag() == true )   GetNumOn()++;

  if (stats.GetOffFlag() == true )  GetNumOff()++;

  // update all rate stats:

  GetRateStats() += stats.GetRateStats();

  // other info

  GetIdNumber() = -1; // special ID for cumulative BgStats

  // number of telescopes

  int nTel = GetNumberOfTelescopes();

  if( nTel==0 ){

    GetNumberOfTelescopes() = stats.GetNumberOfTelescopes();

    GetTelescopePattern()   = stats.GetTelescopePattern();

  }

  else if( nTel == -1 ){

    // do nothing

  }

  else{

    if( nTel != stats.GetNumberOfTelescopes() ){

      GetNumberOfTelescopes() = -1; // from now on, mixed number of telescopes

      GetTelescopePattern()   = -1;

    }

  }

  // on/off/excl regions:

  DEBUG_OUT << "Add regions" << std::endl;

  // accumulate ON region (want to keep just the first one, so we

  // don't get hundreds of copies of the same one)

  if(!online) {

    if (GetOnRegion().GetNumRegions() == 0) {

      DEBUG_OUT << "Accumulating ON region..." << endl;

      GetOnRegion().AddRegion( &(stats.GetOnRegion()));

    }

    DEBUG_OUT << "Accumulating OFF regions..." << endl;

    GetOffRegion().AddRegion( &(stats.GetOffRegion() ) );

    DEBUG_OUT << "Accumulating EXCLUSION regions..." << endl;

    GetOnExclusionRegion().

      AddRegion( &(stats.GetOnExclusionRegion() ) );

    GetOffExclusionRegion().

      AddRegion( &(stats.GetOffExclusionRegion() ) );

  }

  // set the Observation Position:

  GetObservationPosition() = stats.GetObservationPosition();

  // Zenith and Offset Distributions:

  DEBUG_OUT << "Adding zenith and offset distributions..." 

            << endl;

  if (GetZenithDistOn().GetNbinsX() == stats.GetZenithDistOn().GetNbinsX()) 

    GetZenithDistOn().Add( &stats.GetZenithDistOn() );

  else if(GetZenithDistOn().GetNbinsX() <= 1)

    GetZenithDistOn() = stats.GetZenithDistOn();

  else

    throw std::range_error("Can't add ZenithDistOn with incompatible size");

  if (GetZenithDistOff().GetNbinsX() == stats.GetZenithDistOff().GetNbinsX()) 

    GetZenithDistOff().Add( &stats.GetZenithDistOff() );

  else if(GetZenithDistOff().GetNbinsX() <= 1)

    GetZenithDistOff() = stats.GetZenithDistOff();

  else 

    throw std::range_error("Can't add ZenithDistOff with incompatible size");

  if (GetAzimuthDistOn().GetNbinsX() == stats.GetAzimuthDistOn().GetNbinsX()) 

    GetAzimuthDistOn().Add( &stats.GetAzimuthDistOn() );

  else if(GetAzimuthDistOn().GetNbinsX() <= 1)

    GetAzimuthDistOn() = stats.GetAzimuthDistOn();

  else

    throw std::range_error("Can't add AzimuthDistOn with incompatible size");

  if (GetAzimuthDistOff().GetNbinsX() == stats.GetAzimuthDistOff().GetNbinsX()) 

    GetAzimuthDistOff().Add( &stats.GetAzimuthDistOff() );

  else if(GetAzimuthDistOff().GetNbinsX() <= 1)

    GetAzimuthDistOff() = stats.GetAzimuthDistOff();

  else 

    throw std::range_error("Can't add AzimuthDistOff with incompatible size");

  if (GetOffsetDistOn().GetNbinsX() == stats.GetOffsetDistOn().GetNbinsX()) 

    GetOffsetDistOn().Add( &stats.GetOffsetDistOn() );

  else if(GetOffsetDistOn().GetNbinsX() <= 1)

    GetOffsetDistOn() = stats.GetOffsetDistOn();

  else 

    throw std::range_error("Can't add OffsetDistOn with incompatible size");

  if (GetOffsetDistOff().GetNbinsX() == stats.GetOffsetDistOff().GetNbinsX()) 

    GetOffsetDistOff().Add( &stats.GetOffsetDistOff() );

  else if(GetOffsetDistOff().GetNbinsX() <= 1)

    GetOffsetDistOff() = stats.GetOffsetDistOff();

  else 

    throw std::range_error("Can't add OffsetDistOff with incompatible size");

  if (GetZenithOffsetDistOn().GetNbinsX() == stats.GetZenithOffsetDistOn().GetNbinsX() &&

      GetZenithOffsetDistOn().GetNbinsY() == stats.GetZenithOffsetDistOn().GetNbinsY()) 

    GetZenithOffsetDistOn().Add( &stats.GetZenithOffsetDistOn() );

  else if((GetZenithOffsetDistOn().GetNbinsX() <= 1) &&

          (GetZenithOffsetDistOn().GetNbinsY() <= 1))

    GetZenithOffsetDistOn() = stats.GetZenithOffsetDistOn();

  else 

    throw std::range_error("ZenithOffsetDistOn out of range");

  if (GetZenithOffsetDistOff().GetNbinsX() == stats.GetZenithOffsetDistOff().GetNbinsX() &&

      GetZenithOffsetDistOff().GetNbinsY() == stats.GetZenithOffsetDistOff().GetNbinsY()) 

    GetZenithOffsetDistOff().Add( &stats.GetZenithOffsetDistOff() );

  else if((GetZenithOffsetDistOff().GetNbinsX() <= 1) &&

          (GetZenithOffsetDistOff().GetNbinsY() <= 1))

    GetZenithOffsetDistOff() = stats.GetZenithOffsetDistOff();

  else 

    throw std::range_error("ZenithOffsetDistOff out of range");

  // check and set the ConfigName:

  if (GetConfigName() =="") {

    GetConfigName() = stats.GetConfigName();

  }

  else {

    if (GetConfigName().find(stats.GetConfigName())==std::string::npos) {

      WARN_OUT << "You have added statistics using " 

               << stats.GetConfigName() <<  " to a run which was analyzed "

               << "with "<<GetConfigName()<<"! Are you sure that's wise?"

               << endl;

      GetConfigName() += "+"+stats.GetConfigName();

    }

  }

  DEBUG_OUT << "Done adding histograms" << endl;

}

void Background::BgStats::ReplaceStatsFrom(const BgStats &stats, bool online)

{

  std::string Name = GetName();

  GetRateStats() = stats.GetRateStats();

  if(!online) {

    GetOnRegion() = *(stats.GetOnRegion().Clone((Name+"_OnRegion").c_str()));

    GetOffRegion() = *(stats.GetOffRegion().Clone((Name+"_OffRegion").c_str()));

    GetOnExclusionRegion() = *(stats.GetOnExclusionRegion().Clone((Name+"_OnExclusionRegion").c_str()));

    GetOffExclusionRegion() = *(stats.GetOffExclusionRegion().Clone((Name+"_OffExclusionRegion").c_str()));

  }

  GetIdNumber() = stats.GetIdNumber();

  GetNumberOfTelescopes() = stats.GetNumberOfTelescopes();

  GetTelescopePattern() = stats.GetTelescopePattern();

  GetSkippedFlag() = stats.GetSkippedFlag();

  GetBadFlag() = stats.GetBadFlag();

  GetOnFlag() = stats.GetOnFlag();

  GetOffFlag() = stats.GetOffFlag();

  GetNumOn() = stats.GetNumOn();

  GetNumOff() = stats.GetNumOff();

  GetNumSkippedOn() = stats.GetNumSkippedOn();

  GetNumSkippedOff() = stats.GetNumSkippedOff();

  GetNumBadOn() = stats.GetNumBadOn();

  GetNumBadOff() = stats.GetNumBadOff();

  GetZenithDistOn() = stats.GetZenithDistOn();

  GetZenithDistOn().SetName((Name+"_ZenithDistOn").c_str());

  GetZenithDistOff() = stats.GetZenithDistOff();

  GetZenithDistOff().SetName((Name+"_ZenithDistOff").c_str());

  GetAzimuthDistOn() = stats.GetAzimuthDistOn();

  GetAzimuthDistOn().SetName((Name+"_AzimuthDistOn").c_str());

  GetAzimuthDistOff() = stats.GetAzimuthDistOff();

  GetAzimuthDistOff().SetName((Name+"_AzimuthDistOff").c_str());

  GetOffsetDistOn() = stats.GetOffsetDistOn();

  GetOffsetDistOn().SetName((Name+"_OffsetDistOn").c_str());

  GetOffsetDistOff() = stats.GetOffsetDistOff();

  GetOffsetDistOff().SetName((Name+"_OffsetDistOff").c_str());

  GetZenithOffsetDistOn() = stats.GetZenithOffsetDistOn();

  GetZenithOffsetDistOn().SetName((Name+"_ZenithOffsetDistOn").c_str());

  GetZenithOffsetDistOff() = stats.GetZenithOffsetDistOff();

  GetZenithOffsetDistOff().SetName((Name+"_ZenithOffsetDistOff").c_str());

  GetIsCumulative() = stats.GetIsCumulative();

  GetObservationPosition() = stats.GetObservationPosition();

  GetMuonEntry() = stats.GetMuonEntry();

  GetArrayMuonParameters() = stats.GetArrayMuonParameters();

  GetConfigName() = stats.GetConfigName();

}

/**

 * Fills zenith and offset distributions, assuming the ON and OFF

 * regions have been set already. If you only want to fill ON or

 * OFF distributions use the flags 'fillon' and 'filloff'. This

 * might be necessary, if ON and OFF data are taken from different

 * runs.

 *

 * \param obspos - observation position

 * \param starttime - starting time of the run

 * \param endtime - ending time of the run

 * \param psicut - fill offset distribut only with region fullfilling

 *                 the psi cut; ignore the psi, if parameter < 0

 * \param fillon - fill ON distributions (default: true)

 * \param filloff - fill OFF distributions (default: true)

 *

 * \note this assumes the correct ON and OFF regions have already been

 * set up.

 */

void

Background::BgStats::

FillDistributions( const Stash::Coordinate &obspos,

                   const Sash::Time starttime, 

                   const Sash::Time endtime,

                   double psicut, bool fillon, bool filloff) {

  const double MAX_OFFSET = 4.0;

  const double MAX_ZENITH = 90.0;

  const int NBINS_OFFSET = 80;

  const int NBINS_ZENITH = 90;

  const int NBINS_AZIMUTH = 90;

  const double MAX_AZIMUTH = 360.0;

  DEBUG_OUT << "Filling Zenith and Offset distributions " << endl;

  if (!fillon && filloff) {

    DEBUG_OUT << "only fill OFF histograms" << endl;

  }

  if (fillon && !filloff) {

    DEBUG_OUT << "only fill ON histograms" << endl;

  }

  if (!fillon && !filloff) {

    WARN_OUT << "don't fill either ON nor OFF histograms" << endl;

    return;

  }

  if (GetOnRegion().GetArea() < 1.0e-9) {

    WARN_OUT << GetName()

             << ": There is no ON region: ON zenith distribution will "

             << "not be calculated correctly!" << endl;

  }

  if (GetOffRegion().GetArea() < 1.0e-9 ) {

    if( TString(GetName()).Contains("AcceptanceLookupBgMaker") )

      filloff = false;

    else

      WARN_OUT << GetName()

               << ": There is no OFF region: OFF zenith distribution will "

               << "not be calculated correctly!" << endl;

  }

  if (fillon) {

    GetZenithDistOn().SetBins( NBINS_ZENITH, 0.0, MAX_ZENITH );

    GetOffsetDistOn().SetBins( NBINS_OFFSET, 0.0, MAX_OFFSET );

    GetZenithOffsetDistOn().SetBins( NBINS_ZENITH,0.0,MAX_ZENITH,

                                     NBINS_OFFSET,0.0,MAX_OFFSET );

    GetAzimuthDistOn().SetBins( NBINS_AZIMUTH, 0.0, MAX_AZIMUTH );

    GetOnRegion().FillZenithOrAzimuthDistribution( GetZenithDistOn(),starttime,endtime,true);

    GetOnRegion().FillZenithOrAzimuthDistribution( GetAzimuthDistOn(),starttime,endtime,false);

    GetOnRegion().FillOffsetDistribution( obspos, GetOffsetDistOn(), psicut );

    GetOnRegion().FillZenithOffsetDistribution( obspos, 

                                                GetZenithOffsetDistOn(),

                                                starttime,endtime,

                                                psicut );

  }

  if (filloff) {

    GetZenithDistOff().SetBins( NBINS_ZENITH, 0.0, MAX_ZENITH );

    GetOffsetDistOff().SetBins( NBINS_OFFSET, 0.0, MAX_OFFSET );

    GetZenithOffsetDistOff().SetBins( NBINS_ZENITH,0.0,MAX_ZENITH,

                                      NBINS_OFFSET,0.0,MAX_OFFSET );

    GetAzimuthDistOff().SetBins( NBINS_AZIMUTH, 0.0, MAX_AZIMUTH );

    GetOffRegion().FillZenithOrAzimuthDistribution(GetZenithDistOff(),starttime,endtime,true);

    GetOffRegion().FillZenithOrAzimuthDistribution(GetAzimuthDistOff(),starttime,endtime,false);

    GetOffRegion().FillOffsetDistribution( obspos, GetOffsetDistOff(), psicut );

    GetOffRegion().FillZenithOffsetDistribution( obspos, 

                                                 GetZenithOffsetDistOff(),

                                                 starttime, endtime, psicut );

  }

}

Stash::Lambda

Background::BgStats::

GetMaximumOffset() {

  if (GetOffsetDistOn().GetEntries() == 0 ) {

    WARN_OUT << "No entries in ON offset distribution!" << endl;

    return Stash::Lambda(180.0, Stash::Angle::Degrees);

  }

  // check overflow bin - /todo is this bin filled ?

  if (GetOffsetDistOn().GetBinContent( GetOffsetDistOn().GetNbinsX()+1 ) > 0) {

    // low edge of overflow bin

    Stash::Lambda ret = Stash::Lambda( GetOffsetDistOn().

                                       GetBinLowEdge( GetOffsetDistOn().GetNbinsX() ),

                                       Stash::Angle::Degrees );

    WARN_OUT << "overflow bin of OffsetDistOn is not empty, "

             << "the maximum offset might be larger than " << ret << endl;

    return ret;

  }

  for (int i=GetOffsetDistOn().GetNbinsX()-1; i>0; i--) {

    if ( GetOffsetDistOn().GetBinContent( i ) > 0 ) {

      return Stash::Lambda( GetOffsetDistOn().GetBinLowEdge(i+1),

                            Stash::Angle::Degrees );

    }

  }

  // not found, this code shouldn't be reached

  WARN_OUT << "Cannot find a valid entry (entries: " 

           << GetOffsetDistOn().GetEntries()

           << ")" << endl;

  return Stash::Lambda(180, Stash::Angle::Degrees);

}

/**

 * \returns fraction of time the ON region spent above the given offset. 

 * This is a replacement for the old "PsiSqr" run cut, 

 * where runs further away than a given radius were thrown out. Instead, we check that the

 */

double

Background::BgStats::

GetFractionOfTimeAboveOffset( double offs_deg ) {

  // KPK: removed the offset_deg check below because it breaks the

  // skip-run logic. If we return 0 here, the run will never be

  // skipped, which isn't what is intended by this function... we

  // should fall back to an estimate instead.

  // if ( offs_deg > GetOffsetDistOn().GetBinLowEdge( GetOffsetDistOn().GetNbinsX() ) ) {

  //   WARN_OUT << "cannot test for offset angles of " << offs_deg << "deg; " 

  //              << "OffsetDistOn stops at " 

  //              << GetOffsetDistOn().GetBinLowEdge( GetOffsetDistOn().GetNbinsX() )

  //              << endl;

  //       return 0.0;

  // }

  Int_t bin = GetOffsetDistOn().FindBin( offs_deg );

  Int_t maxbin = GetOffsetDistOn().GetNbinsX();

  Double_t integral_above, integral;

  // note Integral() integrates the histogram while GetIntegral()

  // gives you the sum of all points filled including the overflow!

  integral_above = GetOffsetDistOn().Integral( bin, maxbin );

  integral = GetOffsetDistOn().Integral();

  double frac;

  if (integral>1.0e-10) 

    frac = (double)integral_above/(double)integral;

  else

    frac = 1.0; // all bad, since no data!

  return frac;

}

double 

Background::BgStats::

GetMeanMuonCorrection () {

  int count = 0;

  double avg = 0.;

  for(TelMap_t::iterator j = GetMuonCorrections().begin(); 

      j !=  GetMuonCorrections().end(); ++j) {        

    avg += (*j).second;

    count++;

  }

  if (count) avg /= count;

  return avg;

}

void 

Background::BgStats::

Draw(const Option_t* opt) {

  TCanvas *canv = new TCanvas("BgStatsCanvas","BgStatsCanvas",900,600);

  canv->Divide(3,2);

  canv->cd(1);

  GetZenithDistOn().Draw(opt);

  GetZenithDistOff().SetLineColor( kRed );

  GetZenithDistOff().Draw("same");

  canv->cd(2);

  GetAzimuthDistOn().Draw(opt);

  GetAzimuthDistOff().SetLineColor( kRed );

  GetAzimuthDistOff().Draw("same");

  canv->cd(3);

  GetOffsetDistOn().Draw(opt);

  GetOffsetDistOff().SetLineColor( kRed );

  GetOffsetDistOff().Draw("same");

  canv->cd(4);

  GetZenithOffsetDistOn().Draw("colz");

  canv->cd(5);

  GetZenithOffsetDistOff().Draw("colz");

}

/**

 * \brief Calculate mean azimuth angle taking the discontinuity at 0 deg = 360 deg into account.

 */

double Background::BgStats::MeanAzimuth(const TH1& azimuthDist) {

  // idea: sum up vectors on the unit circle for each histogram entry,

  // then calculate mean vector and transform back

  TVector2 weightedVector(0.,0.);

  double sumOfWeights = 0.0;

  for( int i=1 ; i<=azimuthDist.GetNbinsX() ; ++i ){

    double w = azimuthDist.GetBinContent(i);

    if( w == 0.0 ) continue;

    double d = azimuthDist.GetBinCenter(i);

    TVector2 v(TMath::Cos(d*TMath::DegToRad()),TMath::Sin(d*TMath::DegToRad()));

    weightedVector += (w*v);

    sumOfWeights   +=  w;

  }

  if( sumOfWeights <= 0.0 ){

    WARN_OUT << "Empty azimuth distribution??? Returning south." << std::endl;

    return 180.;

  }

  weightedVector /= sumOfWeights;

  double phi = weightedVector.Phi()*TMath::RadToDeg();

  return phi;

}

ClassImp(Background::BgStats)