#ifndef EFFECTIVEAREA_INCLUDED_HH
#define EFFECTIVEAREA_INCLUDED_HH

#include "io_util.hh"
#include "TH2D.h"
#include "TMath.h"
#include "Flux.hh"
#include "Astro.hh"
#include "definitions.hh"

/*! Effective area vs energy */

struct EffectiveArea1D
{
  TH1D H;
  flavor flav;
  channel chan;


  EffectiveArea1D() {}


  EffectiveArea1D( TH1D h , flavor flav, channel chan ) : H(h), flav(flav), chan(chan)  {}


  /*! Return the amount of GeV's in a logarithmic bin */
  double GetWidthE (int i) 
  { 
    double Emin = pow(10.,H.GetXaxis()->GetBinLowEdge(i));
    double Emax = pow(10.,H.GetXaxis()->GetBinUpEdge(i));
    return Emax-Emin;
  }

  void prnt()
  {
    print ( " H has ", H.GetNbinsX(), H.Integral() );
  }


  /*! Integrate a function(logE) (e.g. a flux) over the effective area. Func must have units GeV-1 */

  template< typename F > 
  double integrate ( const F& func )
  {
    double c = 0;

    for (int bx = 1; bx <= H.GetNbinsX() ; bx++ )
    {
      double aeff  = H.GetBinContent(bx);
      double loge  = H.GetBinCenter(bx);
      double width = GetWidthE( bx );
      double f     = func( loge );

      //print ( aeff, loge, width, f );

      c += aeff * width * f;
    }    
    return c;
  }



  /*! Integrate a Aeff* flux to get the total event rate. 
      The flux should provide the method dNdE(flavour,log_E) */

  double GetTotalRate( Flux& F, double livetime = 365*24*3600 )
  {
    double r = integrate( [&](double loge) { return F.dNdE( flav, loge ); } ) * livetime ;
    return r;
  }


  /*! Return a histogratm of the Rate per energy bin. */

  TH1D GetRate( Flux& flux , double livetime = 365*24*3600 )
  {
    TH1D r = H;
    r.SetName("r");

    for (int bx = 1; bx <= r.GetNbinsX() ; bx++ )
    {
      double aeff     = H.GetBinContent(bx);
      double loge     = r.GetXaxis()->GetBinCenter(bx);
      double f        = flux.dNdE( flav, loge );

      //print( "xxx", f, aeff, GetWidthE(bx) , livetime);

      r.SetBinContent(bx, aeff * f * GetWidthE(bx) * livetime );
    }
    return r;
  }


  /*! Version of GetRate that returns 4 vectors: bin-center Energy, Flux, Aeff, Rate */

  vector< vector< double > > GetRateInfo( Flux& flux, int nseconds = 365*24*3600 )
  {
    vector< vector< double >> r;
    r.resize(4);

    TH1D h = GetRate( flux, nseconds );

    for (int bx = 1; bx <= h.GetNbinsX() ; bx++ ) {
      r[0].push_back( h.GetBinCenter(bx) );
      r[1].push_back( flux.dNdE( flav,  h.GetBinCenter(bx) ));
      r[2].push_back( H.GetBinContent(bx) ); // big H = Aeff 
      r[3].push_back( h.GetBinContent(bx) ); // rate
    }

    return r;
  }

  ~EffectiveArea1D() {}
};


/*! Effective Area vs E_true and costheta_true */

struct EffectiveArea
{
  TH2D Haeff_costheta;
  TH2D Haeff_sindecl;
  bool Haeff_sindecl_init = false;

  vector< EffectiveArea1D > Haeff_slices_sindecl; // the 1-d slices  of Haeff_sindecl
  bool GetAreaAtDecl_init = false ;

  Det det; // ugly but....  todo: make one global detector 

  flavor flav;
  channel chan;
  

  EffectiveArea () {
    //print("warning default Aeff ctor" );
    det.latitude = -100;
  }

  // EffectiveArea(TH2D H_costheta, TH2D H_sindecl) : 
  // Haeff_costheta(H_costheta),Haeff_sindecl(H_sindecl) {}

  EffectiveArea(TH2D H_costheta, Det& det): Haeff_costheta(H_costheta), det(det) 
  {

    if (det.latitude < 0.4 || det.latitude > 0.7 ) fatal ("wrong lattitude");
  }


  
  /* compute Aeff vs sin(declination) */
  void ensure_Haeff_sindecl(); // in cc
 

  //functions to extract information from Haeff and Haeff_sindec
  double GetPhaseSpace_costheta (int i) { 
    return GetWidthE_costheta(i) * 2*pi * Haeff_costheta.GetYaxis()->GetBinWidth(1); 
  }

  double GetPhaseSpace_sindecl  (int i) { 
    ensure_Haeff_sindecl();
    return GetWidthE_sindecl(i) * 2*pi * Haeff_sindecl.GetYaxis()->GetBinWidth(1); 
  }

  double GetWidthE_costheta (int i) { 
    double Emin = TMath::Power(10.,Haeff_costheta.GetXaxis()->GetBinLowEdge(i));
    double Emax = TMath::Power(10.,Haeff_costheta.GetXaxis()->GetBinUpEdge(i));
    return Emax-Emin;
  }

  double GetWidthE_sindecl (int i) { 
    ensure_Haeff_sindecl();
    double Emin = TMath::Power(10.,Haeff_sindecl.GetXaxis()->GetBinLowEdge(i));
    double Emax = TMath::Power(10.,Haeff_sindecl.GetXaxis()->GetBinUpEdge(i));
    return Emax-Emin;
  }


  /*! Integrate over energy at the specified cos(zenith). Func must have units GeV-1 */
  template< typename F> 
  double integrate ( double cos_zenith, F func )
  {
    double c = 0;

    int by = Haeff_costheta.GetYaxis()->FindBin( cos_zenith );

    for (int bx = 1; bx <= Haeff_costheta.GetNbinsX() ; bx++ )
    {
      double aeff  = Haeff_costheta.GetBinContent(bx,by);
      double loge  = Haeff_costheta.GetXaxis()->GetBinCenter(bx);
      double width = GetWidthE_costheta( bx );
      c += aeff * width * func( loge );
    }
    
    return c;
  }


  /*! return the total rate. The function will call the flux object with the correct
      flavor */

  double rate( double cos_zenith, Flux* PF )
  {
    return integrate( cos_zenith, [&] (double logEtrue ) { return PF->dNdE( flav , logEtrue); } );
  }

  /*! return the rate (th2d) for the given flux and the given time-span */
  TH2D GetRate( DiffuseFlux* F , int nseconds = 365*24*3600 )
  {

    TH2D r = Haeff_costheta;
    r.SetName("rate");

    for (int bx = 1; bx <= Haeff_costheta.GetNbinsX(); bx++ )
    {
      double loge = Haeff_costheta.GetXaxis()->GetBinCenter(bx);
      for (int by = 1; by <= Haeff_costheta.GetNbinsY(); by++ )
          {
            double costheta = Haeff_costheta.GetYaxis()->GetBinCenter(by);
            double flux     = F->dNdEdOmega(flav, loge, costheta );
            double aeff     = Haeff_costheta.GetBinContent(bx,by);
            r.SetBinContent( bx, by, aeff * flux * GetPhaseSpace_costheta(bx) );
          }
    }    
    r.Scale(nseconds);
    return r;
  }


  /*! return the total rate for the given flux and the given time-span */
  double GetTotalRate( DiffuseFlux* F , int nseconds = 365*24*3600 )
  {
    double rate = 0;
    for (int bx = 1; bx <= Haeff_costheta.GetNbinsX(); bx++ )
    {
      double loge = Haeff_costheta.GetXaxis()->GetBinCenter(bx);
      for (int by = 1; by <= Haeff_costheta.GetNbinsY(); by++ )
          {
            double costheta = Haeff_costheta.GetYaxis()->GetBinCenter(by);
            double aeff     = Haeff_costheta.GetBinContent(bx,by);
            double flux     = F->dNdEdOmega(flav, loge, costheta ); 
            rate += aeff * flux * GetPhaseSpace_costheta(bx);
          }
    }    
    return rate * nseconds;
  }


  EffectiveArea1D GetAreaAtDecl( double decl )
  {
    if (!GetAreaAtDecl_init) {

      ensure_Haeff_sindecl();
      
      int n = Haeff_sindecl.GetYaxis()->GetNbins();
      Haeff_slices_sindecl.resize( n + 1 );

      print ("initializing GetAreaAtDecl for",n,"slices", flav, chan );

      for (int b = 1 ; b <= n ; b++ ) {
          TH1D* h = Haeff_sindecl.ProjectionX( ("hsindeclslice"+str(b)+str(chan)+str(flav)).c_str(),b,b);
          h->SetDirectory(0);
          Haeff_slices_sindecl[b] = EffectiveArea1D( *h, flav, chan );
      }

      GetAreaAtDecl_init = true;
    }
    

    int b = Haeff_sindecl.GetYaxis()->FindBin( sin(decl) );


    return Haeff_slices_sindecl[b];
  }
 

  EffectiveArea1D GetAverage()
  {
    print ("EffectiveArea1D GetAverage() -- warning: this leaks "); // don't remove print before fixing the leak, please

    EffectiveArea1D r( *Haeff_costheta.ProjectionX("aeff_average"), flav, chan );
    r.H.Scale( 1.0 / Haeff_costheta.GetNbinsY() );
    return r;
  }


  //! return total rate of point source at declination= decl
  double GetRateDecl( double decl, Flux* F, int nseconds = 365*24*3600 )
  {
    ensure_Haeff_sindecl();

    double sd = sin(decl);
    double r =0;

    int by = Haeff_sindecl.GetYaxis()->FindBin( sd );

    for (int bx = 1; bx <= Haeff_sindecl.GetNbinsX() ; bx++ )
    {
      double loge     = Haeff_sindecl.GetXaxis()->GetBinCenter(bx);
      double flux     = F->dNdE( flav , loge ); 
      double aeff     = Haeff_sindecl.GetBinContent(bx,by);

      r += flux * aeff * GetWidthE_sindecl( bx );
    }    

    return r * nseconds;
  }

  virtual EffectiveArea *clone() const
  {
    return new EffectiveArea(*this);
  }

  virtual ~EffectiveArea() {}
  
  ClassDef( EffectiveArea, 1)
};

#endif