#ifndef DETRESPONSE_INCLUDED_HH
#define DETRESPONSE_INCLUDED_HH

#include "rootutil.hh"
#include "ana/search/EffectiveArea.hh"
#include "TRandom.h"


struct DetResponse
{
  string name;

  EffectiveArea Aeff;

  DetResponse() {}

  // for background

  virtual double bg_dP_dlogErec_dOmega( double logE_rec, double cos_theta_rec ) {return 0;}

  // the following are for signal

  virtual double dP_dlogErec_dOmega_zenith( double E_rec, double cos_theta_rec, double E_true, double cos_theta_true) {return 0;};

  virtual double dP_dlogErec_dOmega_reco_error( double E_rec, double reco_error, double E_true, double cos_theta_true) {return 0;};

  virtual double dP_dlogErec_dOmega_zenith_bins(int bin_E_rec, int bin_cos_theta_rec, int bin_E_true, int bin_cos_theta_true) 
  {return 0;}

  virtual ~DetResponse() {}




  void dus();

  ClassDef( DetResponse, 1);
};


/*! read object from root file... reasonably */
template<typename T>
void read( TFile& f , T& dst , string name )
{
  TObject* p = f.Get(name.c_str());
  if (!p) fatal("failed to find object called", name, "in TFile", f.GetName());
  dst = *((T*) p);
}

  

struct DefaultDetResponse : public DetResponse
{
  TH3D h_COSZENREC; // x = etrue, y=costhetatrue, z =costheta_rec
  TH3D h_EERR;      // x = etrue, y=costheattrue, z =E_rec
  TH3D h_ANGERR;    // x = etrue, y=costhetatrue, z =log(angle-muon-rec-track)
  

  TH2D h_ATMNU;        // x = etrue, y=costhetatrue
  TH2D h_ATMNUREC;     // x = erec   y=costhetarec
  double n_background; // number of background events from h_ATMMU

  TH2D A; // Aeff

  DefaultDetResponse() {}

  DefaultDetResponse( string ingredients_file, string flavour = "numuCC")
  {
    TFile f( ingredients_file.c_str() );

    read( f, A  , "h_EFFAREA" + flavour );
    Aeff.Haeff = &A;
    read( f, h_COSZENREC ,"h_COSZENREC"+flavour );
    read( f, h_ANGERR    ,"h_ANGERR"+flavour    );
    read( f, h_EERR      ,"h_EERR"+flavour      );
    read( f, h_ATMNU     ,"h_ATMNU"+flavour     );

    // these are pdfs of the quantity on the z-axis. In other words
    // they are dP/ ( dlogE dcostheta).
    // So they must be normalized to that sum ( content * binsize ) = 1

    normalize_z( h_COSZENREC ,1.0, true );
    normalize_z( h_ANGERR    ,1.0, true );
    normalize_z( h_EERR      ,1.0, true );

    check( h_COSZENREC );
    check( h_ANGERR    );
    check( h_EERR      );
    check( h_ATMNU     );

    h_ATMNUREC = smeared_bin( h_ATMNU, [this](int a, int b, int c, int d) {
      double x = this->dP_dlogErec_dOmega_zenith_bins(a,b,c,d);
      //print( a,b,c,d, x );
      return x;
    });

    n_background = h_ATMNU.Integral();
    h_ATMNUREC.SetName("h_ATMNUREC");
    normalize( h_ATMNUREC );
  }

  /*!
  the probability P( Erec, cos(theta)_rec | Etrue , cos(theta)_true ) is factorised as
   P( cos(theta)_rec | Etrue , cos(theta)_true ) *  P( Erec | Etrue , cos(theta)_true )
  */
  virtual double dP_dlogErec_dOmega_zenith(double E_rec, double cos_theta_rec, double E_true, double cos_theta_true)
  {
    double a = interpolate( &h_EERR     , E_true, cos_theta_true, E_rec ) ;
    double b = interpolate( &h_COSZENREC, E_true, cos_theta_true, cos_theta_rec );
    return  a* b / (2*pi);      
  }

  virtual double dP_dlogErec_dOmega_zenith_bins(int bin_E_rec, int bin_cos_theta_rec, int bin_E_true, int bin_cos_theta_true)
  {

    double bb = h_EERR.GetBin(bin_E_true, bin_cos_theta_true, bin_E_rec );
    double a = h_EERR     .GetBinContent(bb ) ;

    if (!isfinite(a)) fatal("h_errr failed: ", &h_EERR, bin_E_true, bin_cos_theta_true, bin_E_rec , a );

    double b = h_COSZENREC.GetBinContent( bin_E_true, bin_cos_theta_true, bin_cos_theta_rec );

    if (!isfinite(b)) fatal("h_cos failed: ", bin_E_true, bin_cos_theta_true, bin_cos_theta_rec);

    return a* b / (2*pi);      
  }


  virtual double dP_dlogErec_dOmega_reco_error( double E_rec, double reco_error_deg, double E_true, double cos_theta_true)
  {
    const double log10_reco_error_deg = log10( reco_error_deg );

    double dP_dlog10alpha = interpolate( &h_ANGERR, E_true, cos_theta_true, log10_reco_error_deg );
    double dP_dalpha = dP_dlog10alpha / ( reco_error_deg * log(10) );
    double dP_dOmega = dP_dalpha / (2*pi * sin( reco_error_deg*pi/180) );

    return interpolate( &h_EERR, E_true, cos_theta_true, E_rec ) *
           dP_dOmega;

  }

  virtual double dP_dlogErec_dOmega_reco_error_log( double E_rec, double log_reco_error_deg, double E_true, double cos_theta_true)
  {
    return dP_dlogErec_dOmega_reco_error( E_rec, pow(10, log_reco_error_deg), E_true, cos_theta_true );
  }

  // dP/dOmega = dP/(dcostheta * 2*pi)
  virtual double bg_dP_dlogErec_dOmega( double logE_rec, double cos_theta_rec )
  {
    int bin = h_ATMNUREC.FindBin( logE_rec, cos_theta_rec );
    return h_ATMNUREC.GetBinContent( bin ) / (2*pi);
  }

  virtual ~DefaultDetResponse() {}

  ClassDef( DefaultDetResponse, 1);
};
#endif