GammaLib

#include <iostream>

#include <TH1.h>

#include <TH2.h>

#include <TF1.h>

#include <TCanvas.h>

#include <TStyle.h>

#include <TFile.h>

#include <TLegend.h>

#include <sash/DataSet.hh>

#include <sash/HESSArray.hh>

#include <crash/EquatorSystem.hh>

#include <plotters/SkyHist.hh>

#include <plotters/FITSTools.hh>

#include <response/EffectiveAreaLookupIR.hh>

#include <background/BgStats.hh>

#include <background/BgMaps.hh>

#include <background/Utils.hh>

#include <background/AcceptanceLookupIR.hh>

#include <flux/PSFLookupIR.hh>

#include <flux/PSFInfo.hh>

//

// reads a PSF file created by the ExtractPSF tool and

// fits a parameterization to the PSF

//

int testfit( const char* filename = "psf.root" ){

  Flux::PSFInfo* psfinfo = Flux::PSFInfo::ReadFromFile(filename);

  if( !psfinfo ) return 1;

  TH1F* dist = psfinfo->CreateNormalizedDistribution();

  TF1* fit   = psfinfo->CreatePSFFunction("TRIPLEEXP");

  double radius = psfinfo->GetContainmentRadius();

  std::cout << "Containment radius: " << radius << std::endl;

  std::cout << "Integral: " << fit->Integral(0.0,0.45) << std::endl;

  TCanvas* psfCanvas = new TCanvas( "psfCanvas", "PSF (read)", 1000, 700 );

  psfCanvas->Divide(2,1);

  std::cout << *psfinfo << std::endl;

  psfCanvas->cd(1);

  dist->Draw();

  fit->Draw("L SAME");

  gPad->SetLogy();

  psfCanvas->cd(2);

  psfinfo->GetWeights().Draw();

  return 0;

}

//

// shows simple Fill functions of PSF class

//

int test_PSFInfo() {

  int muonRun    = 101;

  int telPattern = 30;

  double azimuth = 180.;

  double zenith  = 13.;

  double offset  = 0.5;

  double energy  = 1.0;

  double index = -2.0;

  // it's important to use the fullEnclosure config for effective areas! we don't want any theta cuts

  // since we also don't do them when filling the PSF lookup!

  Flux::PSFLookupIR* psfLookup = new Flux::PSFLookupIR( "std" );

  Response::EffectiveAreaLookupIR* effAreaLookup = 

    new Response::EffectiveAreaLookupIR( "std_fullEnclosure", true, true );

  Flux::PSFInfo* psfinfo_fixedE = new Flux::PSFInfo( Sash::HESSArray::GetHESSArray(), "test_fixedE", index );

  std::cout << *psfinfo_fixedE << std::endl;

  // fill PSF object for a single fixed energy

  psfinfo_fixedE->FillYourself(energy,muonRun,telPattern,azimuth,zenith,offset,1.0,psfLookup);

  // fill PSF object for an assumed E^(index) spectrum

  Flux::PSFInfo* psfinfo_spectrum = new Flux::PSFInfo( Sash::HESSArray::GetHESSArray(), "test_spectrum", index );

  std::cout << psfinfo_spectrum << std::endl;

  psfinfo_spectrum->FillYourself(muonRun,telPattern,azimuth,zenith,offset,1.0,true,psfLookup,effAreaLookup);

  TH1F& histo_fixedE = psfinfo_fixedE->GetPsf();

  TH1F& histo_spectrum = psfinfo_spectrum->GetPsf();

  TCanvas* testCanvas = new TCanvas("testCanvas", "test", 1100, 800 );

  testCanvas->Divide(2,2);

  testCanvas->cd(1);

  histo_fixedE.Draw();

  testCanvas->cd(3);

  psfinfo_fixedE->GetWeights().Draw();

  testCanvas->cd(2);

  histo_spectrum.Draw();

  testCanvas->cd(4);

  psfinfo_spectrum->GetWeights().Draw();

  // test fit

  TCanvas* psfCanvas = new TCanvas( "psfCanvas", "PSF fit", 1000, 700 );

  psfCanvas->cd();

  TH1F* dist = psfinfo_spectrum->CreateNormalizedDistribution();

  TF1* fit = psfinfo_spectrum->CreatePSFFunction("TRIPLEEXP");

  dist->Draw();

  fit->Draw("L SAME");

  // compare PSF for different optical efficiency

  Flux::PSFInfo* psfinfo_fixedE_mu1 = 

    new Flux::PSFInfo( Sash::HESSArray::GetHESSArray(), "test_fixedE_mu1", index );

  psfinfo_fixedE_mu1->FillYourself(energy,100,telPattern,azimuth,zenith,offset,1.0,psfLookup);

  std::cout << *psfinfo_fixedE_mu1 << std::endl;

  std::cout << psfinfo_fixedE_mu1->GetContainmentRadius() << std::endl;

  Flux::PSFInfo* psfinfo_fixedE_mu2 = 

    new Flux::PSFInfo( Sash::HESSArray::GetHESSArray(), "test_fixedE_mu2", index );

  psfinfo_fixedE_mu2->FillYourself(energy,102,telPattern,azimuth,zenith,offset,1.0,psfLookup);

  std::cout << *psfinfo_fixedE_mu2 << std::endl;

  std::cout << psfinfo_fixedE_mu2->GetContainmentRadius() << std::endl;

  TH1F* dist_mu1 = psfinfo_fixedE_mu1->CreateNormalizedDistribution();

  dist_mu1->SetLineColor(4);

  TH1F* dist_mu2 = psfinfo_fixedE_mu2->CreateNormalizedDistribution();

  dist_mu2->SetLineColor(2);

  dist_mu1->Scale(1.0/dist_mu1->GetMaximum());

  dist_mu2->Scale(1.0/dist_mu2->GetMaximum());

  TCanvas* optCanv = new TCanvas( "optCanv", "compare opt eff", 1000, 700 );

  optCanv->cd();

  dist_mu1->Draw();

  dist_mu2->Draw("same");

  gPad->SetLogy();

  return 0;

}

//

// open a hap result file and compare the PSF to data

//

// needs a hap result file and the psf file created by the ExtractPSF task

//

// arguments:

//    hap_result_file:    hap result file containing a BgMaps object

//    psf_file:           PSF as produced by the ExtractPSF task

//    run:                only use specified run (0=all runs) from result file

//    thetaSqrBins:       number of bins for the final theta^2 comparison plot

//    integrationRange:   theta^2 value up to which data and PSF model are integrated

//                        in the calculation of the overall normalization

//    cor_radius:         correlation radius to be used in the calculation of the correlated excess map

//    minPhi, maxPhi:     phi range used for the extraction of the theta^2 plot, in case of nearby sources

//

int compare_to_data( const char* hap_result_file, const char* psf_file, int run = 0,

                     int thetaSqrBins = 100, double integrationRange = 0.05,

                     double cor_radius = 0.08, double minPhi = 0., double maxPhi = 360. ){

  double thetaSqrMax = 0.3;

  Flux::PSFInfo* psfinfo = Flux::PSFInfo::ReadFromFile(psf_file);

  if( !psfinfo ) return 10;

  TFile* file = new TFile( hap_result_file );

  if( !file->IsOpen() ) return 1;

  Sash::HESSArray* hess = &(Sash::HESSArray::GetHESSArray());

  Background::BgStats *stats = 0;

  Sash::DataSet* ds = 0;

  if( run == 0 ) ds = Background::FindTotalDataSet();

  else ds = Background::FindRunDataSet();

  if( !ds ) return 2;

  if( run == 0 ){

    ds->GetEntry(0);

    std::vector<const Background::BgStats*> statsvec = 

      hess->GetAll<Background::BgStats>();

    for( unsigned int j=0 ; j<statsvec.size() ; ++j )

      std::cout << "BgStats " << j << ": " << statsvec[j]->GetName() << std::endl;

    stats = statsvec.size() > 0 ? const_cast<Background::BgStats*>(statsvec.front()) : 0;

  }

  else{

    for( int i=0 ; i<ds->GetEntries() ; ++i ){

      ds->GetEntry(i);

      std::vector<const Background::BgStats*> statsvec = 

        hess->GetAll<Background::BgStats>();

      for( unsigned int j=0 ; j<statsvec.size() ; ++j )

        std::cout << "BgStats " << j << ": " << statsvec[j]->GetName() << std::endl;

      stats = statsvec.size() > 0 ? const_cast<Background::BgStats*>(statsvec.front()) : 0;

      if( stats->GetIdNumber() == run ) break;

      else stats = 0;

    }

  }

  if( !stats ) return 3;

  std::cout << *stats << std::endl;

  std::vector<const Background::BgMaps*> mapsvec = 

    hess->GetAll<Background::BgMaps>();

  for( unsigned int j=0 ; j<mapsvec.size() ; ++j ){

    std::cout << "BgMaps  " << j << ": " << mapsvec[j]->GetName() << std::endl;

  }

  Background::BgMaps *maps = mapsvec.size() > 0 ? const_cast<Background::BgMaps*>(mapsvec.front()) : 0;

  if( !maps ) return 4;

  TH1F* psfdist = psfinfo->CreateNormalizedDistribution();

  TF1* psffit   = psfinfo->CreatePSFFunction("TRIPLEEXP");

  double radius = psfinfo->GetContainmentRadius();

  std::cout << "Containment radius: " << radius << std::endl;

  TCanvas* psfCanvas = new TCanvas( "psfCanvas", "PSF", 800, 600 );

  psfCanvas->Divide(2,1);

  std::cout << *psfinfo << std::endl;

  psfCanvas->cd(1);

  psfdist->Draw();

  psffit->Draw("L SAME");

  gPad->SetLogy();

  psfCanvas->cd(2);

  psfinfo->GetWeights().Draw();

  TCanvas* excessCanvas = new TCanvas( "excessCanvas", "excess map", 800, 800 );

  excessCanvas->cd();

  std::cout << "Making the Excess map..." << std::endl;

  Plotters::SkyHist *excess = maps->MakeExcessMap();

  std::cout << "Smoothing the Excess map..." << std::endl;

  Plotters::SkyHist *smooth = excess->GetCorrelatedSkyMap( cor_radius );

  smooth->SetVisFov(1.0); //2.0

  smooth->Draw("colz");

  smooth->DrawPSFFunction(psffit,cor_radius,false);

  // add a wedge showing the range of the theta^2 plot if that differs from the full circle

  if( minPhi > 0.0 || maxPhi < 360. ) {

    Plotters::SkyHist* sh = smooth;

    Tools::CircleSegment* wedge = new Tools::CircleSegment( sh->GetPosition(),

                                                            0., sqrt(thetaSqrMax),

                                                            minPhi, maxPhi, "wedge" );

    wedge->SetLineWidth(1);

    wedge->SetLineColor(kWhite);

    sh->DrawRegion(wedge);

  }

  gPad->Modified();

  TCanvas* azmCanvas = new TCanvas( "azmCanvas", "azimuth-related plots", 800, 600 );

  azmCanvas->cd();

  std::cout << "Making azimuth profile..." << std::endl;

  TH1F* azmProfile = excess->MakeAzimuthProfile(excess->GetPosition(), 0., sqrt(thetaSqrMax), 12, false );

  azmProfile->SetStats(0);

  azmProfile->Draw();

  if( azmProfile->GetMinimum() > 0.0 )

    azmProfile->SetMinimum(0);

  gPad->SetGridy();

  TCanvas* thetaCanvas = new TCanvas( "thetaCanvas", "theta plots", 800, 600 );

  thetaCanvas->cd();

  std::cout << "Making th2 for ON..." << std::endl;

  TH1F* thon = maps->GetOnMap().MakeThetaSq(thetaSqrMax,thetaSqrBins,minPhi,maxPhi,true);

  thon->Draw("hist");

  Plotters::SkyHist *alphaoff = (Plotters::SkyHist*) maps->GetOffMap().Clone("AlphaOff");

  std::cout << "Making alpha map..." << std::endl;

  Plotters::SkyHist *alpha = maps->MakeAlphaMap( );

  alphaoff->Multiply( alpha );

  std::cout << "Making th2 for alpha*OFF..." << std::endl;

  TH1F *thoff = alphaoff->MakeThetaSq(thetaSqrMax,thetaSqrBins,minPhi,maxPhi,true);

  thoff->SetLineColor(kRed);

  thoff->Draw("hist same");

  gPad->Modified();

  TF1* bgLevel = new TF1( "bgLevel", "pol0", 0., thetaSqrMax );

  thoff->Fit( bgLevel, "RW0" );

  bgLevel->Draw("L SAME");

  TCanvas* compCanvas = new TCanvas( "compCanvas", "comparison of data and PSF", 800, 800 );

  compCanvas->cd();

  TH1F* excessThetaSq = excess->MakeThetaSq(thetaSqrMax,thetaSqrBins,minPhi,maxPhi,true);

  excessThetaSq->SetTitle("comparison of PSF to data");

  excessThetaSq->Draw();

  excessThetaSq->GetYaxis()->SetTickLength(-0.01);

  TH1F* psf_rescaled = (TH1F*)psfinfo->GetPsf().Clone("psf_rescaled");

  psf_rescaled->SetLineColor(2);

  psf_rescaled->SetMarkerColor(2);

  // calculate scale factor from integral over first few bins

  // start with coarse histogram

  int dataBin = excessThetaSq->FindBin(integrationRange);

  double exactIntegrationRange = excessThetaSq->GetXaxis()->GetBinUpEdge(dataBin);

  double dataIntegral = excessThetaSq->Integral(1,dataBin,"width");

  int psfBin = psf_rescaled->FindBin(exactIntegrationRange-0.0001);

  double psfIntegral = psf_rescaled->Integral(1,psfBin,"width");

  double scaleFactor = dataIntegral/psfIntegral;

  std::cout << "dataBin: " << dataBin << " psfBin: " << psfBin 

            << " exact range: " << exactIntegrationRange << std::endl;

  psf_rescaled->Scale(scaleFactor);

  psf_rescaled->Draw("same");

  TLegend* leg = new TLegend( 0.6,0.6,0.89,0.89 );

  leg->AddEntry(excessThetaSq, "data", "L");

  leg->AddEntry(psf_rescaled, "PSF model", "L" );

  leg->SetFillColor(kWhite);

  leg->SetLineColor(kWhite);

  leg->Draw();

  // Crab fit (similar fit as in Crab paper, for comparison)

  TF1* func = new TF1( "psf_function",

                       "[0]*(exp(-((x)/(2*[1]*[1])))+[2]*exp(-((x)/(2*[3]*[3]))))",

                        0.0005, 0.2 );

  func->SetParameter(0,excessThetaSq->GetMaximum());

  func->SetParameter(1,0.046);

  func->SetParameter(2,0.15);

  func->SetParameter(3,0.12);

  excessThetaSq->Fit(func,"RI0");

  func->SetLineColor(kGreen);

  func->SetLineWidth(1);

  // func->Draw("L SAME");

  std::cout << "Simple fit: Chi2 / ndf = " << func->GetChisquare() << " / " << func->GetNDF() 

            << " prob: " << func->GetProb() << std::endl;

  excessThetaSq->GetXaxis()->SetRangeUser(0.,0.1);

  // create histograms with identical binning for Kolmogorov test

  // we assume that the PSF model has a finer binning than the data histogram

  // with the number of bins in a given range being an integer multiple of the data bins

  double maxK = 0.1;

  int nbinsK = excessThetaSq->FindBin(maxK);

  double xmaxK = excessThetaSq->GetXaxis()->GetBinUpEdge(nbinsK);

  double ratiod = psf_rescaled->FindBin(xmaxK-0.0001)/nbinsK;

  int ratio = int(round(ratiod));

  std::cout << "Ratio: " << ratiod << " -> " << ratio << std::endl;

  TH1F* dataHistoK = new TH1F( "dataHistoK", "comparison of data and PSF model;#theta^{2} / deg^{2}", 

                               nbinsK, 0., xmaxK );

  TH1F* modelHistoK = new TH1F( "modelHistoK", "modelHistoK", nbinsK, 0., xmaxK );

  for( int i=1 ; i<=nbinsK ; ++i ){

    dataHistoK->SetBinContent(i,excessThetaSq->GetBinContent(i));

    dataHistoK->SetBinError(i,excessThetaSq->GetBinError(i));

    double contentSum = 0.0;

    double varSum = 0.0;

    for( int j=1+(i-1)*ratio ; j<=i*ratio ; ++j ){

      double binc = psf_rescaled->GetBinContent(j);

      double bine = psf_rescaled->GetBinError(j);

      contentSum += binc;

      varSum += bine*bine;

    }

    modelHistoK->SetBinContent(i,contentSum/ratiod);

    modelHistoK->SetBinError(i,sqrt(varSum)/ratiod);

  }

  TCanvas* kolmoCanv = new TCanvas("kolmoCanv","Kolmogorov test", 900, 900 );

  kolmoCanv->cd(1);

  dataHistoK->Draw();

  dataHistoK->GetYaxis()->SetTickLength(-0.01);

  modelHistoK->SetLineColor(psf_rescaled->GetLineColor());

  modelHistoK->Draw("SAME");

  dataHistoK->SetMaximum(1.05*max(dataHistoK->GetMaximum(),modelHistoK->GetMaximum()));

//   quickPsf->Draw("hist same");

  dataHistoK->GetXaxis()->SetRangeUser(0.,0.036);

  TLegend* leg2 = new TLegend( 0.6,0.6,0.89,0.89 );

  leg2->AddEntry(dataHistoK, "data", "L");

  leg2->AddEntry(modelHistoK, "PSF model", "L" );

  leg2->SetFillColor(kWhite);

  leg2->SetLineColor(kWhite);

  leg2->Draw();

  double prob = dataHistoK->KolmogorovTest(modelHistoK);

  std::cout << " Kolmogorov prob: " << prob << std::endl;

  return 0;

}

int psf_to_fits( const char* psffile, const char* fitsfile ){

  Flux::PSFInfo* psfinfo = Flux::PSFInfo::ReadFromFile(psffile);

  if( !psfinfo ) return 1;

  TF1* psffunc = psfinfo->CreatePSFFunction( "TRIPLEEXP" );

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

                         Stash::Beta( "0d0'0.0" ),

                         Crash::GetRADecJ2000System() );

  Stash::Beta fov(0.5,Stash::Angle::Degrees);

  Stash::SmallAngle binsize(0.002,Stash::Angle::Degrees);

  Plotters::SkyHist* sh = new Plotters::SkyHist("sh","PSF",pos,fov,binsize);

  Background::AcceptanceLookupIR::FillAcceptanceMap( *sh, pos, 1.0, psffunc );

  std::cout << "Writing map \"" << sh->GetName() << "\" to file " << fitsfile << std::endl;

  Plotters::SkyHistToFITS(sh, fitsfile, sh->GetName() );

  return 0;

}

//

//

//

//   THE FOLLOWING FUNCTIONS ARE FOR EXPERT TESTS ONLY!!!

//

//

//

// test the --fast option of ExtractPSF

//

int test_fast( const char* file1 = "psf.root", const char* file2 = "psf_fast.root" ) {

  Flux::PSFInfo* psfinfo1 = Flux::PSFInfo::ReadFromFile(file1);

  if( !psfinfo1 ) return 1;

  Flux::PSFInfo* psfinfo2 = Flux::PSFInfo::ReadFromFile(file2);

  if( !psfinfo2 ) return 2;

  std::cout << "1:\n" << *psfinfo1 << std::endl;

  std::cout << "2:\n" << *psfinfo2 << std::endl;

  TH1F* dist1 = psfinfo1->CreateNormalizedDistribution();

  TH1F* dist2 = psfinfo2->CreateNormalizedDistribution();

  dist2->SetLineColor(2);

  TCanvas* psfCanvas = new TCanvas( "psfFastCanvas", "PSF test", 1000, 700 );

  psfCanvas->cd();

  dist1->Draw();

  dist2->Draw("SAME");

  gPad->SetLogy();

  return 0;

}

//

// test uncertainty handling

//

int test_uncertainty( const char* filename = "psf.root" ) {

  Flux::PSFInfo* psfinfo = Flux::PSFInfo::ReadFromFile(filename);

  if( !psfinfo ) return 1;

  std::cout << *psfinfo << std::endl;

  std::cout << "Fitting psf with DOUBLEEXP" << std::endl;

  TF1* fit21 = psfinfo->CreatePSFFunction("DOUBLEEXP");

  std::cout << "Fitting psf with TRIPLEEXP" << std::endl;

  TF1* fit22 = psfinfo->CreatePSFFunction("TRIPLEEXP");

  TH1F*  dist  = psfinfo->CreateNormalizedDistributionOldStyle();

  TH1F*  dist2 = psfinfo->CreateNormalizedDistribution();

  std::cout << "name: " << dist2->GetName() << std::endl;

  TCanvas* testCanvas = new TCanvas("uncertaintyCanvas", "test uncertainties", 1100, 800 );

  testCanvas->Divide(2,2);

  testCanvas->cd(1);

  dist->SetTitle("old-style");

  dist->Draw();

  gPad->SetLogy();

  testCanvas->cd(2);

  dist2->SetTitle("new style, double exp");

  dist2->Draw();

  fit21->Draw("L SAME");

  gPad->SetLogy();

  testCanvas->cd(4);

  TH1D* dist2c = (TH1D*)dist2->Clone("dist2c");

  dist2c->SetTitle("new style, triple exp");

  dist2c->Draw();

  fit22->Draw("L SAME");

  gPad->SetLogy();

  return(0);

}

//

// access the PSF lookup file directly

//

int draw_fixed_PSF() {

  gStyle->SetOptFit(1111);

  TFile* f = new TFile( "$HESSCONFIG/std/PSF.root" , "READ" );

  TH2F* psf = (TH2F*)f->Get("ThetaSq_muon101_30telp_180azm_20deg_0.5off");

  double logE = log10(1.0);

  int bin = psf->GetXaxis()->FindBin(logE);

  TH1D* proj = psf->ProjectionY("_py",bin,bin);

  proj->SetStats(0);

  TCanvas* canvas = new TCanvas("canvas", "PSF", 1100, 800 );

  canvas->cd();

  proj->Draw("hist");

  gPad->SetLogy();

  Stash::Coordinate pos( Stash::Lambda( "18h55m52.63s" ),

                         Stash::Beta( "-10d0'52.5" ),

                         Crash::GetRADecJ2000System() );

  Stash::Beta fov(0.3,Stash::Angle::Degrees);

  Stash::SmallAngle binsize(0.002,Stash::Angle::Degrees);

  Plotters::SkyHist* sh = new Plotters::SkyHist("sh","PSF",pos,fov,binsize);

  Plotters::SkyHist* shf = new Plotters::SkyHist("shf","PSF func",pos,fov,binsize);

  Background::AcceptanceLookupIR::FillAcceptanceMap( *sh, pos, 1.0, proj );

  sh->SetContour(50);

  TCanvas* psfCanvas = new TCanvas("psfCanvas", "PSF", 1100, 800 );

  psfCanvas->Divide(2,2);

  psfCanvas->cd(1);

  psf->Draw("colz");

  psf->SetStats(0);

  psf->SetMaximum(100.);

  psfCanvas->cd(2);

  proj->Draw();

  gPad->SetLogy();

  TF1* func = new TF1( "psf_function",

                       "[0]*exp(-((x)/(2*[1]*[1])))+[2]*exp(-((x)/(2*[3]*[3])))",

                       0.0005, 0.2 );

  double s = proj->GetMaximum();

  func->SetParameters( s*0.3, 0.2, s, 0.05 );

  func->SetParLimits(0,0.001,1000000.);

  func->SetParLimits(1,0.05,0.8);

  func->SetParLimits(2,0.001,1000000.);

  func->SetParLimits(3,0.01,0.2);

  proj->Fit( func, "", "", 0.000, 0.04 );

  proj->SetStats(1);

  Background::AcceptanceLookupIR::FillAcceptanceMap( *shf, pos, 1.0, func );

  shf->SetContour(50);

  psfCanvas->cd(3);

  sh->Draw("colz");

  gPad->Modified();

  psfCanvas->cd(4);

  shf->Draw("colz");

  gPad->Modified();

  return 0;

}

//

// compare azimuth dependence of PSF

//

int azimuth_dependence() {

  double spectralIndex = 3.5;

  Flux::PSFLookupIR* psfLookup = new Flux::PSFLookupIR( "std" );

  Response::EffectiveAreaLookupIR* effAreaLookup = 

    new Response::EffectiveAreaLookupIR( "std_fullEnclosure", true, true );

  Flux::PSFInfo* psfinfo_north = 

    new Flux::PSFInfo( Sash::HESSArray::GetHESSArray(), "psf_north", spectralIndex );

  std::cout << *psfinfo_north << std::endl;

  Flux::PSFInfo* psfinfo_south = 

    new Flux::PSFInfo( Sash::HESSArray::GetHESSArray(), "psf_south", spectralIndex );

  std::cout << *psfinfo_south << std::endl;

  Flux::PSFInfo* psfinfo_ipol = 

    new Flux::PSFInfo( Sash::HESSArray::GetHESSArray(), "psf_ipol", spectralIndex );

  std::cout << *psfinfo_ipol << std::endl;

  psfinfo_north->FillYourself(101,30,0.,30.,0.5,1.0,true,psfLookup,effAreaLookup);

  psfinfo_south->FillYourself(101,30,180.,30.,0.5,1.0,true,psfLookup,effAreaLookup);

  psfinfo_ipol->FillYourself( 101,30,90.,30.,0.5,1.0,true,psfLookup,effAreaLookup);

  TH1F* psfdist_north = psfinfo_north->CreateNormalizedDistribution();

  TH1F* psfdist_south = psfinfo_south->CreateNormalizedDistribution();

  TH1F* psfdist_ipol  = psfinfo_ipol->CreateNormalizedDistribution();

  psfdist_north->Rebin(6);

  psfdist_north->Scale(1./6.);

  psfdist_south->Rebin(6);

  psfdist_south->Scale(1./6.);

  psfdist_ipol->Rebin(6);

  psfdist_ipol->Scale(1./6.);

  TCanvas* azmCanv = new TCanvas("azmCanv","azimuth dependence of PSF", 900, 900 );

  azmCanv->cd();

  psfdist_south->SetStats(0);

  psfdist_south->SetLineColor(4);

  psfdist_south->Draw();

  psfdist_south->GetYaxis()->SetTickLength(-0.01);

  psfdist_south->GetXaxis()->SetRangeUser(0.,0.1);

  psfdist_north->SetLineColor(1);

  psfdist_north->Draw("same");

  psfdist_ipol->SetLineColor(2);

  psfdist_ipol->Draw("same");

  return 0;

}

//

// compare zenith dependence of PSF

//

int zenith_dependence() {

  double spectralIndex = 3.5;

  Flux::PSFLookupIR* psfLookup = new Flux::PSFLookupIR( "std" );

  Response::EffectiveAreaLookupIR* effAreaLookup = 

    new Response::EffectiveAreaLookupIR( "std_fullEnclosure", true, true );

  Flux::PSFInfo* psfinfo_20 = 

    new Flux::PSFInfo( Sash::HESSArray::GetHESSArray(), "psf_20", spectralIndex );

  std::cout << *psfinfo_20 << std::endl;

  Flux::PSFInfo* psfinfo_30 = 

    new Flux::PSFInfo( Sash::HESSArray::GetHESSArray(), "psf_30", spectralIndex );

  std::cout << *psfinfo_30 << std::endl;

  Flux::PSFInfo* psfinfo_ipol = 

    new Flux::PSFInfo( Sash::HESSArray::GetHESSArray(), "psf_ipol", spectralIndex );

  std::cout << *psfinfo_ipol << std::endl;

  psfinfo_20->FillYourself(101,30,180.,20.,0.5,1.0,true,psfLookup,effAreaLookup);

  psfinfo_30->FillYourself(101,30,180.,30.,0.5,1.0,true,psfLookup,effAreaLookup);

  psfinfo_ipol->FillYourself(101,30,180.,25.,0.5,1.0,true,psfLookup,effAreaLookup);

  TH1F* psfdist_20 = psfinfo_20->CreateNormalizedDistribution();

  TH1F* psfdist_30 = psfinfo_30->CreateNormalizedDistribution();

  TH1F* psfdist_ipol  = psfinfo_ipol->CreateNormalizedDistribution();

  psfdist_20->Rebin(6);

  psfdist_20->Scale(1./6.);

  psfdist_30->Rebin(6);

  psfdist_30->Scale(1./6.);

  psfdist_ipol->Rebin(6);

  psfdist_ipol->Scale(1./6.);

  TCanvas* zenCanv = new TCanvas("zenCanv","zenith dependence of PSF", 900, 900 );

  zenCanv->cd();

  psfdist_30->SetStats(0);

  psfdist_30->SetLineColor(4);

  psfdist_30->Draw();

  psfdist_30->GetYaxis()->SetTickLength(-0.01);

  psfdist_30->GetXaxis()->SetRangeUser(0.,0.1);

  psfdist_20->SetLineColor(1);

  psfdist_20->Draw("same");

  psfdist_ipol->SetLineColor(2);

  psfdist_ipol->Draw("same");

  return 0;

}

//

// compare dependence of PSF on spectral index

//

int spectral_dependence() {

  Flux::PSFLookupIR* psfLookup = new Flux::PSFLookupIR( "std" );

  Response::EffectiveAreaLookupIR* effAreaLookup = 

    new Response::EffectiveAreaLookupIR( "std_fullEnclosure", true, true );

  Flux::PSFInfo* psfinfo_20 = new Flux::PSFInfo( Sash::HESSArray::GetHESSArray(), "psf_20", 2.0 );

  std::cout << *psfinfo_20 << std::endl;

  Flux::PSFInfo* psfinfo_25 = new Flux::PSFInfo( Sash::HESSArray::GetHESSArray(), "psf_25", 2.5 );

  std::cout << *psfinfo_25 << std::endl;

  Flux::PSFInfo* psfinfo_30 = new Flux::PSFInfo( Sash::HESSArray::GetHESSArray(), "psf_30", 3.0 );

  std::cout << *psfinfo_30 << std::endl;

  psfinfo_20->FillYourself(101,30,180.,20.,0.5,1.0,true,psfLookup,effAreaLookup);

  psfinfo_25->FillYourself(101,30,180.,20.,0.5,1.0,true,psfLookup,effAreaLookup);

  psfinfo_30->FillYourself(101,30,180.,20.,0.5,1.0,true,psfLookup,effAreaLookup);

  TH1F* psfdist_20 = psfinfo_20->CreateNormalizedDistribution();

  TH1F* psfdist_25 = psfinfo_25->CreateNormalizedDistribution();

  TH1F* psfdist_30 = psfinfo_30->CreateNormalizedDistribution();

  psfdist_20->Rebin(6);

  psfdist_20->Scale(1./6.);

  psfdist_30->Rebin(6);

  psfdist_30->Scale(1./6.);

  psfdist_25->Rebin(6);

  psfdist_25->Scale(1./6.);

  TCanvas* specCanv = new TCanvas("specCanv","dependence of PSF on spectral index", 900, 900 );

  specCanv->cd();

  psfdist_30->SetStats(0);

  psfdist_30->SetLineColor(4);

  psfdist_30->Draw();

  psfdist_30->GetYaxis()->SetTickLength(-0.01);

  psfdist_30->GetXaxis()->SetRangeUser(0.,0.1);

  psfdist_20->SetLineColor(1);

  psfdist_20->Draw("same");

  psfdist_25->SetLineColor(2);

  psfdist_25->Draw("same");

  return 0;

}

//

// compare zenith interpolation of PSF

// set HESSCONFIG to right directory !!!

// ~/scratch/lookups/std_10deg_fullEnclosure/result

//

int fancy_zenith_dependence() {

  double spectralIndex = 3.5;

  Flux::PSFLookupIR* psfLookup = new Flux::PSFLookupIR( "std_fullEnclosure" );

  Response::EffectiveAreaLookupIR* effAreaLookup = 

    new Response::EffectiveAreaLookupIR( "std_fullEnclosure", true, true );

  Flux::PSFLookupIR* psfLookup2 = new Flux::PSFLookupIR( "std_10deg_fullEnclosure" );

  Response::EffectiveAreaLookupIR* effAreaLookup2 = 

    new Response::EffectiveAreaLookupIR( "std_10deg_fullEnclosure", true, true );

  Flux::PSFInfo* psfinfo_0 = 

    new Flux::PSFInfo( Sash::HESSArray::GetHESSArray(), "psf_0", spectralIndex );

  std::cout << *psfinfo_0 << std::endl;

  Flux::PSFInfo* psfinfo_10_corr = 

    new Flux::PSFInfo( Sash::HESSArray::GetHESSArray(), "psf_10_corr", spectralIndex );

  std::cout << *psfinfo_10_corr << std::endl;

  Flux::PSFInfo* psfinfo_10_ipol = 

    new Flux::PSFInfo( Sash::HESSArray::GetHESSArray(), "psf_10_ipol", spectralIndex );

  std::cout << *psfinfo_10_ipol << std::endl;

  Flux::PSFInfo* psfinfo_20 = 

    new Flux::PSFInfo( Sash::HESSArray::GetHESSArray(), "psf_20", spectralIndex );

  std::cout << *psfinfo_20 << std::endl;

  psfinfo_0->FillYourself(101,30,180.,0.,0.5,1.0,true,psfLookup,effAreaLookup);

  psfinfo_10_corr->FillYourself(101,30,180.,10.,0.5,1.0,true,psfLookup2,effAreaLookup2);

  psfinfo_10_ipol->FillYourself( 101,30,180.,10.,0.5,1.0,true,psfLookup,effAreaLookup);

  psfinfo_20->FillYourself(101,30,180.,20.,0.5,1.0,true,psfLookup,effAreaLookup);

  TH1F* psfdist_0 = psfinfo_0->CreateNormalizedDistribution();

  TH1F* psfdist_10_corr = psfinfo_10_corr->CreateNormalizedDistribution();

  TH1F* psfdist_10_ipol = psfinfo_10_ipol->CreateNormalizedDistribution();

  TH1F* psfdist_20 = psfinfo_20->CreateNormalizedDistribution();

  psfdist_0->Rebin(6);

  psfdist_0->Scale(1./6.);

  psfdist_10_corr->Rebin(6);

  psfdist_10_corr->Scale(1./6.);

  psfdist_10_ipol->Rebin(6);

  psfdist_10_ipol->Scale(1./6.);

  psfdist_20->Rebin(6);

  psfdist_20->Scale(1./6.);

  TCanvas* zenCanv = new TCanvas("zenCanv","zenith dependence of PSF", 900, 900 );

  zenCanv->cd();

  psfdist_0->SetStats(0);

  psfdist_0->SetLineColor(4);

  psfdist_0->Draw();

  psfdist_0->GetYaxis()->SetTickLength(-0.01);

  psfdist_0->GetXaxis()->SetRangeUser(0.,0.1);

  psfdist_20->SetLineColor(1);

  psfdist_20->Draw("same");

  psfdist_10_ipol->SetLineColor(2);

  psfdist_10_ipol->Draw("same");

  psfdist_10_corr->SetLineColor(3);

  psfdist_10_corr->Draw("same");

  return 0;

}