GammaLib

#include <background/AcceptanceLookupIR.hh>

#include <background/BgMaps.hh>

#include <plotters/SkyHist.hh>

#include <crash/EquatorSystem.hh>

#include <response/Smoothing.hh>

#include <TCanvas.h>

#include <TLegend.h>

#include <sstream>

int test_radialAcceptance() {

  Stash::Coordinate pos( Stash::Lambda( "18h0m0s" ),

                         Stash::Beta( "0d" ),

                         Crash::GetRADecJ2000System() );

  double binning = 0.02;

  double fov = 3.0;

  double zenith = 25.;

  Background::BgMaps* maps = new Background::BgMaps(Sash::HESSArray::GetHESSArray());

  maps->Init( pos, binning, fov, fov, false );

  Background::AcceptanceLookupIR lookup("std");

  Plotters::SkyHist* exposure = &(maps->GetOnExposureMap());

  lookup.FillAcceptanceMap( *exposure, pos, zenith );

  TCanvas* canv = new TCanvas("canv","canvas",800,800);

  canv->cd();

  exposure->Draw();

  gPad->Modified();

  return 0;

}

//

// test of energy-dependent acceptance lookup

//

int test_edep( const char* config = "std", double Emin = 2., double Emax = 10.,

               double zenith = 25., double psiSqrMax = 12.25,

               int degree = 6, int nbins = 150, int stepsize = 1) {

  Background::AcceptanceLookupIR lookup( config, "RadialAcceptance.root" );

//   Background::AcceptanceLookupIR lookup( config, "RadialAcceptanceOff.root" );

  TCanvas* canv = new TCanvas("edepcanv","edep",1000,800);

  canv->cd();

  TH1F* ac = lookup.CreateAcceptanceCurve(zenith, Emin, Emax);

  if( !ac ) return 1;

  ac->Draw();

  TH1F* smooth = lookup.CreateSmoothedAcceptanceCurve(zenith, Emin, Emax, psiSqrMax,

                                                      degree, nbins, stepsize );

  if( !smooth ) return 1;

  smooth->SetLineColor(4);

  smooth->Draw("same");

  return 0;

}

//

// experiment with different paramters for the polynomial smoothing of acceptance curves

// a sliding-window technique is used for the smoothing

// standard parameters to try are

// degree=3, nbins=30 (fine smoothing)

// degree=6, nbins=150 (coarse smoothing)

//

// or try the default values of degree=0, nbins=0: different parameters are tried and

// the best smoothing is chosen based on chi^2

//

// the trick is to find parameters that give a good description of the acceptance curves

// without following the statistical fluctuations more than necessary.

//

int test_all_edep( const char* config = "std", double Emin = 2., double Emax = 10.,

                   int degree = 6, int nbins = 150, int stepsize = 1,

                   double psiSqrMax = 12.25 ) {

  Background::AcceptanceLookupIR lookup( config, "RadialAcceptance.root" );

  TCanvas* canv = new TCanvas("lookupCanvas","energy-dependent acceptances",1100,900);

  canv->Divide(3,2);

  TCanvas* rescanv = new TCanvas("residualsCanvas","smoothing residulas",1200,900);

  rescanv->Divide(1,6);

  double zeniths[6] = { 10., 25., 35., 42., 47., 60. };

  for( int i=0 ; i<6 ; ++i ){

    double zenith = zeniths[i];

    canv->cd(i+1);

    TH1F* ac = lookup.CreateAcceptanceCurve(zenith, Emin, Emax);

    if( !ac ) return 1;

    ac->Draw();

    TH1F* smooth = lookup.CreateSmoothedAcceptanceCurve(zenith, Emin, Emax, psiSqrMax,

                                                        degree, nbins, stepsize );

    if( !smooth ) return 1;

    smooth->SetLineColor(4);

    smooth->Draw("same");

    rescanv->cd(i+1);

    std::stringstream resname;

    resname << "residuals_" << i;

    TH1F* residuals = (TH1F*)smooth->Clone(resname.str().c_str());

    residuals->Reset();

    for( int bin=1 ; bin<=smooth->GetNbinsX() ; ++bin ){

      double pull = 0.0;

      int acbin = ac->FindBin(smooth->GetBinCenter(bin));

      if( ac->GetBinError(acbin) > 0.0 )

        pull = ( smooth->GetBinContent(bin)-ac->GetBinContent(acbin) ) / ac->GetBinError(acbin);

      residuals->SetBinContent(bin,pull);

      residuals->SetBinError(bin,1.);

    }

    residuals->Draw();

    residuals->SetMinimum(-3.);

    residuals->SetMaximum(3.);

    residuals->SetStats(0);

    residuals->GetYaxis()->SetTickLength(-0.005);

    gPad->SetGridy();

  }

  return 0;

}

int compare_edep() {

  const char* config = "std";

  double Emin_1 = 0.1;

  double Emin_2 = 1.5;

  double Emax = 50.;

  double zenith = 10.;

  double psiSqrMax = 9.;

  int degree = 6;

  int nbins = 150;

  int stepsize = 1;

  Background::AcceptanceLookupIR lookup( config, "RadialAcceptance.root" );

  TCanvas* canv = new TCanvas("edepcanv","edep",1000,800);

  canv->cd();

  TH1F* smooth1 = lookup.CreateSmoothedAcceptanceCurve(zenith, Emin_1, Emax, psiSqrMax,

                                                       degree, nbins, stepsize );

  if( !smooth1 ) return 1;

  smooth1->SetLineColor(2);

  TH1F* smooth2 = lookup.CreateSmoothedAcceptanceCurve(zenith, Emin_2, Emax, psiSqrMax,

                                                       degree, nbins, stepsize );

  if( !smooth2 ) return 2;

  smooth2->SetLineColor(4);

  smooth2->Scale(smooth1->GetMaximum()/smooth2->GetMaximum());

  smooth1->Draw();

  smooth2->Draw("same");

  std::stringstream le1, le2;

  le1 << "E > " << Emin_1 << " TeV";

  le2 << "E > " << Emin_2 << " TeV";

  TLegend* leg = new TLegend(0.8,0.8,0.99,0.99);

  leg->SetFillColor(kWhite);

  leg->AddEntry(smooth1,le1.str().c_str(),"L");

  leg->AddEntry(smooth2,le2.str().c_str(),"L");

  leg->Draw();

  return 0;

}