GammaLib

#include <response/EffectiveAreaLookupIR.hh>

#include <utilities/StringTools.hh>

#include <muon/ArrayMuonParameters.hh>

#include <sash/HESSArray.hh>

#include <TH1.h>

#include <TCanvas.h>

#include <TLine.h>

#include <TGraph.h>

#include <TF1.h>

#include <TLegend.h>

#include <iostream>

#include <sstream>

int test_effArea() {

  // standard tests for fixed optical configuration

  Response::EffectiveAreaLookupIR lookup( "std", false );

  int muonEntry        = 100;

  int telescopePattern = 30;

  double azimuth       = 197.467;

  double zenith        = 29.9336;

  double offset        = 0.450788;

  double energy        = 0.404803;

  // lookup.PrintLookupArrays();

  double effArea;

  bool success = lookup.InterpolateEffectiveArea(energy,zenith,offset,muonEntry,telescopePattern,azimuth,effArea);

  std::cout << "EffArea: " << effArea << " success: " << success << std::endl;

  double thr,emax;

  lookup.GetRange(zenith,offset,muonEntry,telescopePattern,azimuth,thr,emax);

  std::cout << "Range: " << thr << " .. " << emax << std::endl;

  return 0;

}

int test_effArea_opticalConfiguration(){

  // test interpolation based on optical efficiencies

  Response::EffectiveAreaLookupIR lookup( "std", false );

  int telescopePattern = 30;

  double azimuth       = 359.704;

  double zenith        = 45.2264;

  double offset        = 0.446313;

  double energy        = 0.541777;

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

  ap.GetMcCorrectionFactorForSecondMuonEntry()  = 1.172;

  ap.GetDataCorrectionFactorForFirstMuonEntry() = 1.20635;

  ap.GetFirstMuonEntry()   = 101;

  ap.GetSecondMuonEntry()  = 102;

  ap.GetTelescopePattern() = telescopePattern;

  double spectralIndex = 2.5;

  double E1 = 1.0;

  double E2 = 100.0;

  double livetime = 1.0;

  double effArea;

  bool success = lookup.InterpolateEffectiveArea(&ap,energy,zenith,offset,azimuth,effArea);

  std::cout << "EffArea: " << effArea << " success: " << success << std::endl;

  double thr,emax;

  lookup.GetRange(&ap,zenith,offset,azimuth,thr,emax);

  std::cout << "Range: " << thr << " .. " << emax << std::endl;

  std::cout << "Energy was: " << energy << std::endl;

  int nbinsfailed;

  double xi = lookup.CalculateExpectedCounts(&ap,spectralIndex,E1,E2,zenith,offset,azimuth,livetime,nbinsfailed);

  std::cout << "Expected counts: " << xi << " (" << nbinsfailed << " bins failed)" << std::endl;

  return 0;

}

//

// Plot effective area as a function of the parameter given as first argument.

//

// example: plot_effArea(3,"muonRun=101,telPattern=28,azimuth=0.,offset=0.2,energy=40.0")

//          plot_effArea(1,"muonRun=101,azimuth=0.,zenith=8.3,offset=0.2,energy=40.0")

//          plot_effArea(5,"muonRun=101,azimuth=0.,zenith=8.3,offset=0.2,telPattern=28") (also plots threshold!)

//

TH1F* plot_effArea( unsigned int parameterToVary, const char* parameters = "" ){

  Response::EffectiveAreaLookupIR lookup( "std", true );

  std::string parstr(parameters);

  bool success = true;

  int muonRun      = ValueFromString(parstr,"muonRun",success);

  int telPattern   = ValueFromString(parstr,"telPattern",success);

  double azimuth   = ValueFromString(parstr,"azimuth",success);

  double zenith    = ValueFromString(parstr,"zenith",success);

  double offset    = ValueFromString(parstr,"offset",success);

  double energy    = ValueFromString(parstr,"energy",success);

  int nbins = 100;

  double xmin = 0.;

  double xmax = 0.;

  bool logx = false;

  switch( parameterToVary ){

  case 0: // optical efficiency

    nbins = 3;

    xmin = 100.;

    xmax = 103.;

    break;

  case 1: // telPattern

    nbins = 21;

    xmin = 11.5;

    xmax = 32.5;

    break;

  case 2: // azimuth

    nbins = 90;

    xmin = 0.;

    xmax = 360.;

    break;

  case 3: // zenith

    nbins = 70;

    xmin = 0.;

    xmax = 70.;

    break;

  case 4: // offset

    nbins = 50;

    xmin = 0.;

    xmax = 3.0;

    break;

  case 5: // energy

    nbins = 100;

    xmin = 0.1;

    xmax = 20.0;

    logx = true;

    break;

  default:

    return 0;

    break;

  }

  TH1F* histo = lookup.CalculateInterpolatedHistogram( energy, zenith, offset, muonRun, telPattern, azimuth, parameterToVary,

                                                       nbins, xmin, xmax, logx );

  histo->SetStats(false);

  TCanvas* c = new TCanvas( "c", "effArea", 900, 700 );

  c->cd();

  histo->Draw();

  if(parameterToVary==5){

    double threshold, energymax;

    lookup.GetRange( zenith, offset, 

                     muonRun,

                     telPattern, azimuth,

                     threshold, energymax );

    TLine* thrLine = new TLine(threshold,histo->GetMinimum(),

                               threshold,histo->GetMaximum()*1.05);

    thrLine->SetLineColor(2);

    thrLine->SetLineStyle(2);

    thrLine->Draw();

  }

  if( logx )

    gPad->SetLogx();

  return histo;

}

int compare_thresholds() {

  Response::EffectiveAreaLookupIR lookup( "std", true );

  int telPattern   = 30;

  double azimuth   = 180.;

  double zenith    = 20.;

  double offset    = 0.5;

  double dummy = 0.;

  int nbins = 400;

  double xmin = 0.1;

  double xmax = 20.0;

  bool logx = true;

  double ymin = 1.e4;

  double ymax = 3.e5;

  int muonRun = 0;

  double threshold, energymax;

  TH1F* histo100 = lookup.CalculateInterpolatedHistogram( dummy, zenith, offset,

                                                          muonRun = 100,

                                                          telPattern, azimuth, 5,

                                                          nbins, xmin, xmax, logx );

  histo100->SetLineColor(1);

  lookup.GetRange( zenith, offset, 

                   muonRun,

                   telPattern, azimuth,

                   threshold, energymax );

  TLine* thrLine100 = new TLine(threshold,ymin,

                                threshold,ymax*1.05);

  thrLine100->SetLineColor(1);

  thrLine100->SetLineWidth(2);

  std::cout << "Thr 100: " << threshold << std::endl;

  TH1F* histo101 = lookup.CalculateInterpolatedHistogram( dummy, zenith, offset,

                                                          muonRun = 101,

                                                          telPattern, azimuth, 5,

                                                          nbins, xmin, xmax, logx );

  histo101->SetLineColor(2);

  lookup.GetRange( zenith, offset, 

                   muonRun,

                   telPattern, azimuth,

                   threshold, energymax );

  TLine* thrLine101 = new TLine(threshold,ymin,

                                threshold,ymax*1.05);

  thrLine101->SetLineColor(2);

  thrLine101->SetLineWidth(2);

  std::cout << "Thr 101: " << threshold << std::endl;

  double corrFactor101 = 1.454;

  TLine* thrLine101Comparison = new TLine(corrFactor101*thrLine100->GetX1(),ymin,

                                          corrFactor101*thrLine100->GetX2(),ymax*1.05);

  thrLine101Comparison->SetLineColor(2);

  thrLine101Comparison->SetLineStyle(2);

  TH1F* histo102 = lookup.CalculateInterpolatedHistogram( dummy, zenith, offset,

                                                          muonRun = 102,

                                                          telPattern, azimuth, 5,

                                                          nbins, xmin, xmax, logx );

  histo102->SetLineColor(3);

  lookup.GetRange( zenith, offset, 

                   muonRun,

                   telPattern, azimuth,

                   threshold, energymax );

  TLine* thrLine102 = new TLine(threshold,ymin,

                                threshold,ymax*1.05);

  thrLine102->SetLineColor(3);

  thrLine102->SetLineWidth(2);

  std::cout << "Thr 102: " << threshold << std::endl;

  double corrFactor102 = 1.17;

  TLine* thrLine102Comparison = new TLine(corrFactor102*thrLine101->GetX1(),ymin,

                                          corrFactor102*thrLine101->GetX2(),ymax*1.05);

  thrLine102Comparison->SetLineColor(3);

  thrLine102Comparison->SetLineStyle(2);

  TCanvas* thrCanvas = new TCanvas( "thrCanvas", "thresholds", 1100, 800 );

  thrCanvas->cd();

  histo100->Draw("HIST");

  histo100->SetMinimum(ymin);

  histo100->SetMaximum(ymax);

  thrLine100->Draw();

  histo101->Draw("HIST SAME");

  thrLine101->Draw();

  histo102->Draw("HIST SAME");

  thrLine102->Draw();

  thrLine101Comparison->Draw();

  thrLine102Comparison->Draw();

  gPad->SetLogx();

  gPad->SetLogy();

  return 0;

}

int comparePatterns(){

  TH1F* histo = plot_effArea(1,"muonRun=103,azimuth=0.,zenith=25.,offset=0.5,energy=1.0");

//   TH1F* histo = plot_effArea(1,"muonRun=103,azimuth=0.,zenith=45.,offset=0.5,energy=0.8");

  TH1F* dhisto = (TH1F*)histo->Clone("dhisto");

  dhisto->Reset();

  int patterns[5] = { 14, 22, 26, 28, 30 };

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

    int bin = histo->FindBin(patterns[i]);

    dhisto->SetBinContent(bin,histo->GetBinContent(bin));

    dhisto->SetBinError(bin,histo->GetBinError(bin));

  }

  dhisto->Draw();

  return 0;

}

int compareExpectedCounts(){

  Response::EffectiveAreaLookupIR lookup( "std", false );

  TH1F* histo = new TH1F( "histo", "expected counts;pattern;expected counts / a.u.", 25, 10.5, 35.5 );

  int muonEntry        = 103;

  double azimuth       = 0.;

  double zenith        = 55.;

  double offset        = 0.5;

  double spectralIndex = 2.5;

  double E2            = 90.;

  int telPatterns[5] = { 14, 22, 26, 28, 30 };

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

    int pattern = telPatterns[i];

    double thr,emax;

    int nBinsFailed;

    lookup.GetRange(zenith,offset,muonEntry,pattern,azimuth,thr,emax);

    std::cout << "Pattern: " << pattern << " Range: " << thr << " .. " << emax << std::endl;

    double expCounts = lookup.CalculateExpectedCounts( spectralIndex, thr, E2,

                                                       zenith, offset, azimuth,

                                                       muonEntry, pattern,

                                                       1., nBinsFailed );

    std::cout << "Exp counts: " << expCounts << std::endl;

    histo->Fill(pattern,expCounts);

  }

  TCanvas* c = new TCanvas( "expCountsCanvas", "expected counts", 1100, 800 );

  c->cd();

  histo->SetStats(0);

  histo->Draw();

  return 0;

}

int calcExpCounts(){

  Response::EffectiveAreaLookupIR lookup( "hard", true );

  int nb;

  double expCounts = lookup.CalculateExpectedCounts(2.5,1.,50.,25.,0.5,335.,103,30,3600.,nb);

  std::cout << "EXP: " << expCounts << std::endl;

  return 0;

}

int expCountsVsAzm(){

  Response::EffectiveAreaLookupIR lookup( "hard", true );

  int nb;

  TH1D* histo = new TH1D( "histo", "exp counts vs azm", 360, 0., 360. );

  for( double azm = 0.0 ; azm < 360. ; azm += 5.0 ){

    double expCounts = lookup.CalculateExpectedCounts(2.5,1.,50.,25.,0.5,azm,103,30,3600.,nb);

    histo->Fill( azm, expCounts );

  }

  histo->SetStats(0);

  new TCanvas;

  histo->Draw();

  return 0;

}

int indexDependenceOfIntegralFlux( bool useRawEnergy = true ){

  gROOT->SetStyle("Plain");

  Response::EffectiveAreaLookupIR lookup( "std", true, useRawEnergy );

  double referenceIndex = 2.5; //3.0;

  double indexRange = 0.6;

  double indexStep = 0.05;

  double thresholdEnergies[2] = { 0.2, 1.0 };

  double zeniths[3] = { 20., 40., 60. };

  double offsets[3] = { 0.5, 0.75, 1.0};

  double azimuth       = 180.;

  int muonEntry        = 101;

  int telescopePattern = 30;

  double Emax = 100.;

  TGraph* graphs[2][3][3];

  for( int ene = 0 ; ene < 2 ; ++ene ){

    double E1 = thresholdEnergies[ene];

    for( int zen = 0 ; zen < 3 ; ++zen ){

      double zenith = zeniths[zen];

      for( int off = 0 ; off < 3 ; ++off ){

        double offset = offsets[off];

        int nBinsFailed;

        double expCountsRef = lookup.CalculateExpectedCounts( referenceIndex, E1, Emax,

                                                              zenith, offset, azimuth,

                                                              muonEntry, telescopePattern,

                                                              1.0, nBinsFailed );

        TGraph*& graph = graphs[ene][zen][off];

        graph = new TGraph;

        double index = referenceIndex - indexRange;

        int i=0;

        while( index <= referenceIndex + indexRange ){

          double expCounts = lookup.CalculateExpectedCounts( index, E1, Emax,

                                                             zenith, offset, azimuth,

                                                             muonEntry, telescopePattern,

                                                             1.0, nBinsFailed );

          double fluxRatio = expCountsRef / expCounts;

          double deviation = 100.0 * (fluxRatio - 1.0);

          graph->SetPoint(i,index,deviation);

          index += indexStep;

          ++i;

        }

      }

    }

  }

  TCanvas* canv = new TCanvas( "indexDependenceOfIntegralFluxCanvas", "index dependence of integral flux", 900, 600 );

  canv->Divide(2,1);

  TLegend* legend = 0;

  for( int ene = 0 ; ene < 2 ; ++ene ){

    canv->cd(ene+1);

    legend = new TLegend(0.2,0.6,0.4,0.85);

    graphs[ene][0][0]->SetLineColor(2);

    graphs[ene][0][0]->Draw("AL");

    std::stringstream title;

    title << "E > " << thresholdEnergies[ene] << " TeV";

    graphs[ene][0][0]->GetHistogram()->SetTitle(title.str().c_str());

    graphs[ene][0][0]->GetHistogram()->SetMaximum(200.);

    graphs[ene][0][0]->GetHistogram()->SetMinimum(-200.);

    graphs[ene][0][0]->GetHistogram()->GetXaxis()->SetTitle("index");

    graphs[ene][0][0]->GetHistogram()->GetYaxis()->SetTitle("deviation (%)");

    for( int zen = 0 ; zen < 3 ; ++zen ){

      for( int off = 0 ; off < 3 ; ++off ){

        if( !zen && !off ) continue;

        graphs[ene][zen][off]->SetLineColor(zen+2);

        graphs[ene][zen][off]->SetLineStyle(off+1);

        graphs[ene][zen][off]->Draw("L");

      }

      std::stringstream text;

      text << zeniths[zen] << " deg zenith";

      legend->AddEntry(graphs[ene][zen][0],text.str().c_str(),"L");

    }

    if( thresholdEnergies[ene] == 1.0 ){

      TF1* expectation = new TF1( "expectation", "100*((x-1)/([0]-1)-1)", referenceIndex-indexRange,  referenceIndex+indexRange);

      expectation->SetParameter(0,referenceIndex);

      expectation->SetLineStyle(2);

      expectation->SetLineColor(1);

      expectation->Draw("L SAME");

      cout << expectation->Eval(2.) << " " << expectation->Eval(3.) << " " << expectation->Eval(4.) << endl;

      cout << "Drawing expectation." << endl;

      legend->AddEntry(expectation,"expectation","L");

    }

    legend->SetFillColor(10);

    legend->SetLineColor(10);

    legend->Draw();

  }

  return 0;

}

int test_GaussianEnergyResolution() {

  // standard tests for fixed optical configuration

  Response::EffectiveAreaLookupIR lookup( "std", false );

  int muonEntry        = 100;

  int telescopePattern = 30;

  double azimuth       = 197.467;

  double zenith        = 29.9336;

  double offset        = 0.450788;

  double energy        = 0.404803;

//   lookup.PrintLookupArrays();

  double energyResolution, energyBias;

  bool success1 = lookup.EnergyResolution(energy,zenith,offset,muonEntry,telescopePattern,azimuth,energyResolution);

  bool success2 = lookup.EnergyBias(energy,zenith,offset,muonEntry,telescopePattern,azimuth,energyBias);

  std::cout << "EnergyResolution: " << energyResolution << " success: " << success1 << std::endl;

  std::cout << "EnergyBias:       " << energyBias << " success: " << success2 << std::endl;

  return 0;

}

int test_GaussianEnergyResolution_opticalConfiguration(){

  // test interpolation based on optical efficiencies

  Response::EffectiveAreaLookupIR lookup( "std", false );

  int telescopePattern = 30;

  double azimuth       = 0.0;

  double zenith        = 13.0;

  double offset        = 0.4;

  double energy        = 1.4;

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

  ap.GetMcCorrectionFactorForSecondMuonEntry()  = 1.4;

  ap.GetDataCorrectionFactorForFirstMuonEntry() = 1.1;

  ap.GetFirstMuonEntry()   = 101;

  ap.GetSecondMuonEntry()  = 102;

  ap.GetTelescopePattern() = telescopePattern;

  double energyResolution, energyBias;

  bool success1 = lookup.EnergyResolution(&ap,energy,zenith,offset,azimuth,energyResolution);

  bool success2 = lookup.EnergyBias(&ap,energy,zenith,offset,azimuth,energyBias);

  std::cout << "EnergyResolution: " << energyResolution << " success: " << success1 << std::endl;

  std::cout << "EnergyBias:       " << energyBias << " success: " << success2 << std::endl;

  return 0;

}

int plot_GaussianEnergyResolution( double zenith = 10., double azimuth = 42., double offset = 0.42,

                                   int muonEntry = 100, int telPattern = 30 ) {

  Response::EffectiveAreaLookupIR lookup( "std", false );

  TH1* eres = new TH1D("eres", "Energy Resolution", 100, -2, 2);

  eres->GetXaxis()->SetTitle("log_{10}(E)");

  eres->GetYaxis()->SetTitle("relative energy resolution");

  TH1* ebias = new TH1D("ebias", "Energy Resolution", 100, -2, 2);

  ebias->GetXaxis()->SetTitle("log_{10}(E)");

  ebias->GetYaxis()->SetTitle("relative energy bias");

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

    double energy = pow10(eres->GetBinCenter(bin));

    double energyResolution, energyBias;

    lookup.EnergyResolution(energy,zenith,offset,muonEntry,telPattern,azimuth,energyResolution);

    lookup.EnergyBias(energy,zenith,offset,muonEntry,telPattern,azimuth,energyBias);

    eres->SetBinContent(bin, energyResolution);

    ebias->SetBinContent(bin, energyBias);

  }

  TCanvas* eresCanvas = new TCanvas( "eresCanvas", "energy resolution", 1200, 800 );

  eresCanvas->Divide(2,1);

  eresCanvas->cd(1);

  ebias->Draw();

  eresCanvas->cd(2);

  eres->Draw();

  return 0;

}