GammaLib

#include <iostream>

#include <TH1.h>

#include <TCanvas.h>

#include <TLine.h>

#include <sash/HESSArray.hh>

#include <muon/ArrayMuonParameters.hh>

#include <flux/EnergyResolutionLookupIR.hh>

#include <flux/EnergyResolution2LookupIR.hh>

int test_energyResolution() {

  // standard test for fixed optical configuration

  Flux::EnergyResolutionLookupIR lookup( "std" );

  Flux::EnergyResolution2LookupIR lookup2( "std" );

  int muonEntry        = 100;

  int telescopePattern = 30;

  double azimuth       = 197.467;

  double zenith        = 29.9336;

  double offset        = 0.450788;

  double energy        = 0.404803;

  double relativeError = 0.1;

  double reco_energy   = (relativeError + 1) * energy;

  double pdf, pdf_err;

  bool success = lookup.InterpolatePdf(energy,relativeError,zenith,offset,muonEntry,

                                       telescopePattern,azimuth,pdf,pdf_err);

  std::cout << "pdf: " << pdf << " +- " << pdf_err << " success: " << success << std::endl;

  double pdf2, pdf_err2;

  bool success2 = lookup2.InterpolatePdf(energy,reco_energy,zenith,offset,muonEntry,

                                         telescopePattern,azimuth,pdf2,pdf_err2);

  std::cout << "pdf2: " << pdf2 << " +- " << pdf_err2 << " success2: " << success2 << std::endl;

  return 0;

}

int test_energyResolution_opticalConfiguration(){

  // test interpolation based on optical efficiencies

  Flux::EnergyResolutionLookupIR lookup( "std" );

  Flux::EnergyResolution2LookupIR lookup2( "std" );

  int telescopePattern = 30;

  double azimuth       = 0.0;

  double zenith        = 13.0;

  double offset        = 0.4;

  double energy        = 1.4;

  double relativeError = 0.1;

  double reco_energy   = (relativeError + 1) * energy;

  Muon::ArrayMuonParameters ap(Sash::HESSArray::GetHESSArray());

  ap.GetMcCorrectionFactorForSecondMuonEntry()  = 1.4;

  ap.GetDataCorrectionFactorForFirstMuonEntry() = 1.1;

  ap.GetFirstMuonEntry()   = 101;

  ap.GetSecondMuonEntry()  = 102;

  ap.GetTelescopePattern() = telescopePattern;

  double pdf, pdf_err;

  bool success = lookup.InterpolatePdf(&ap,energy,relativeError,zenith,offset,azimuth,pdf,pdf_err);

  std::cout << "pdf: " << pdf << " +- " << pdf_err << " success: " << success << std::endl;

  double pdf2, pdf_err2;

  bool success2 = lookup2.InterpolatePdf(&ap,energy,reco_energy,zenith,offset,azimuth,pdf2,pdf_err2);

  std::cout << "pdf2: " << pdf2 << " +- " << pdf_err2 << " success2: " << success2 << std::endl;

  return 0;

}

int plot_slice( double energy = 1.0, double zenith = 25., double offset = 1.3, double azimuth = 223.,

                int muonEntry = 100, int telPattern = 30 ) {

  Flux::EnergyResolutionLookupIR lookup( "std" );

  TH1F* pdf = lookup.CalculatePdfHistogram( energy, zenith, offset, muonEntry, telPattern, azimuth );

  pdf->Draw();

  return 0;

}

int plot_slice_opticalConfiguration( double energy = 1.0, double zenith = 25., 

                                     double offset = 1.3, double azimuth = 223. ) {

  Flux::EnergyResolutionLookupIR lookup( "std" );

  Muon::ArrayMuonParameters ap(Sash::HESSArray::GetHESSArray());

  ap.GetMcCorrectionFactorForSecondMuonEntry()  = 1.4;

  ap.GetDataCorrectionFactorForFirstMuonEntry() = 1.1;

  ap.GetFirstMuonEntry()   = 101;

  ap.GetSecondMuonEntry()  = 102;

  ap.GetTelescopePattern() = 30;

  TH1F* pdf = lookup.CalculatePdfHistogram(&ap, energy, zenith, offset, azimuth );

  pdf->Draw();

  return 0;

}

/*

int plot_pdf( double zenith = 25., double offset = 1.3, double azimuth = 223.,

              int muonEntry = 100, int telPattern = 30 ) {

  Flux::EnergyResolutionLookupIR lookup( "std" );

  TH1F* pdf = lookup.CalculatePdfHistogram( energy, zenith, offset, muonEntry, telPattern, azimuth );

  TCanvas* canvas = new TCanvas("EnergyResolution");

  canvas->Divide(2,1);

  canvas->cd(1);

  pdf->Draw("colz");

  canvas->cd(2);

  return 0;

}

*/