#ifndef DETRESPONSE_INCLUDED_HH
#define DETRESPONSE_INCLUDED_HH
#include "rootutil.hh"
#include "ana/search/EffectiveArea.hh"
#include "TRandom.h"
#include "PsfParam.hh"
#include "definitions.hh"
#include "histogram_conversion.hh"
#include "statutil.hh"

/*! 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("\n\nfailed to find object called", name, "\nin TFile", f.GetName(), "\n\n");
  }
  dst = *((T *)p);
  dst.SetDirectory(0);
}

struct DetResponse
{
  string name;
  channel chan;
  flavor flav;
  EffectiveArea Aeff;
  PsfParam psf; // smooth parametrisation.
  Det det;

  DetResponse() {}

  // 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);
};

struct DefaultDetResponse : public DetResponse
{
  string name;

  // detector resolution
  TH3D h_COSZENREC;    // x = etrue, y=costhetatrue, z =costheta_rec
  TH3D h_EREC;         // x = etrue, y=costheattrue, z =E_rec
  TH3D h_ANGERR;       // x = etrue, y=costhetatrue, z =log(angle-muon-rec-track)
  TH3D h_EREC_sindecl; // like h_EREC bug with sindecl on the y-axis

  // expected number of events
  map<string, TH3D> h_BKGREC_3D; // x = erec   y=costhetarec (Background contributions) z = azimuth
  TH3D h_BKGREC_azi_sindecl;

  map<string, TH2D> h_BKGREC; // x = erec   y=costhetarec (Background contributions)
  TH2D h_BKGREC_sindecl;

  TH3D h_EFFAREA_3D; // Aeff vs energy vs zenith vs azimuth
  TH2D h_EFFAREA;    // Aeff vs energy and zenith angle

  // maps to store histograms for NC and CC interactions ["CC/NC"]
  map<string, TH3D> h_COSZENREC_ALL;
  map<string, TH3D> h_EREC_ALL;
  map<string, TH3D> h_ANGERR_ALL;
  map<string, TH3D> h_BKGREC_3D_ALL;
  map<string, TH2D> h_BKGREC_ALL;
  map<string, TH3D> h_EFFAREA_3D_ALL;
  map<string, TH2D> h_EFFAREA_ALL;

  DefaultDetResponse()
  {
    init_det();
  }

  DefaultDetResponse(string ingredients_file,
                     flavor fl = numu,
                     vector<flavor_interaction> fl_int = {numuCC}, // you can do NC, CC or ALL here
                     string rec_type = "track",
                     double src_width = 0.,
                     int psf_mergebins = 2,
                     bool gaus = true,
                     double cos_zenith_cut = 999 )
  {
    print ("DefaultDetResponse ctor");

    init_det();

    TFile f(ingredients_file.c_str());

    // check if we have 2D or 3D histograms for EFFAREA and BKGREC
    TObject *tmp = f.Get(("h_EFFAREA_numuCC_" + rec_type).c_str());
    if (!tmp) fatal("failed to get ","h_EFFAREA_numuCC_" + rec_type )

    bool use_3D_hist = tmp->InheritsFrom("TH3");
    print ("DefaultDetResponse 1");  


    if (fl != muon)
    {

      name = "DefaultDetResponse for " + flavor_name(fl) + "/" + rec_type + ", using 3D histograms " + use_3D_hist;

      // extract histograms from file and store
      for (auto fl_int_tmp : fl_int)
      {
        cout << "before read_histograms " << flavor_name(fl) << " " << flavor_interaction_name(fl_int_tmp) << endl;

        read_histograms(f, flavor_name(fl), fl_int_tmp, rec_type, use_3D_hist);
      }

      // add CC and NC for ALL, or copy histograms if only 1 interaction type
      add_histograms(flavor_name(fl), fl_int, rec_type);

      f.Close();



      // zero out all bins in h_EFFARA wtih cos_zenith > cos_zenith_cut.
      // this to ensure misreconstructed events are not counted in 
      ///the signal efficiency

      for (int bx = 1; bx <= h_EFFAREA.GetNbinsX(); bx++)
      {
        for (int by = 1; by <= h_EFFAREA.GetNbinsY(); by++)
        {
          double costheta = h_EFFAREA.GetYaxis()->GetBinCenter(by);
          int bin = h_EFFAREA.GetBin( bx, by );
          if ( costheta > cos_zenith_cut )  h_EFFAREA.SetBinContent( bin , 0 );
        }
      }


      // now make the Aeff object

      Aeff = EffectiveArea(h_EFFAREA, det);

      // 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_EREC, 1.0, true);

      // compute the energy response as function on sindecl
      print("computing  h_EREC_sindecl", flavor_name(fl), rec_type);
      h_EREC_sindecl = convert_to_sindecl(det, h_EREC);
      print(h_EREC_sindecl.Integral());

      check(h_COSZENREC);
      check(h_EREC);
      check(h_ANGERR);
      check(h_BKGREC["ATMNU_conv"]);
      check(h_BKGREC["ATMNU_prompt"]);
    

      if (psf_mergebins != 0)
      {
        std::cout << "initializing PSF" << std::flush;
        psf.name = "PSF for " + name;
        psf.init(h_ANGERR, psf_mergebins, src_width, gaus);
      }

      h_BKGREC["ATM_all"] = h_BKGREC["ATMNU_conv"];
      h_BKGREC["ATM_all"].Add(&h_BKGREC["ATMNU_prompt"]);
      //  h_BKGREC["ATM_all"].Add(&h_BKGREC["ATMMUON"]);
    }
    else // muon flavor
    {

      cout << "getting "
           << "h_ATMMU_reco_" << rec_type << endl;
      TObject *tmp = f.Get(("h_ATMMU_reco_" + rec_type).c_str());
      bool use_3D_hist = tmp->InheritsFrom("TH3");

      for (auto fl_int_tmp : fl_int)
      {

        read_histograms(f, flavor_name(fl), fl_int_tmp, rec_type, use_3D_hist);
      }

      add_histograms(flavor_name(fl), fl_int, rec_type);

      f.Close();

      check(h_BKGREC["ATMMUON"]);
      h_BKGREC["ATM_all"] = h_BKGREC["ATMMUON"];
    }
    // compute_background vs sindecl
    h_BKGREC_sindecl = convert_to_sindecl(det, h_BKGREC["ATM_all"]);

    cout << "sum background: " << h_BKGREC["ATM_all"].Integral() << endl;
  }

  void init_det()
  {
    det.longitude = 0.278819;
    det.latitude = 0.633407;
    det.meridian_convergence_angle = 0.0100733;
  }

  /*!
  Function reads the histograms from the ingredients file
  3D histograms are projected to 2D histograms
  */
  void read_histograms(TFile &f, string flavor, flavor_interaction fl_int, string rec_type, bool use_3D_hist)
  {
    string idstr = flavor_interaction_name(fl_int) + "_" + rec_type;
    print("Reading IRFs for flavor=", flavor, " channel/rectype = ", rec_type, flavor_interaction_name(fl_int));

    if (fl_int != muonmuon)
    {
      string idstr = flavor_interaction_name(fl_int) + "_" + rec_type;

      read(f, h_COSZENREC_ALL[flavor_interaction_name(fl_int)], "h_COSZENREC_" + idstr);
      read(f, h_ANGERR_ALL[flavor_interaction_name(fl_int)], "h_ANGERR_" + idstr);
      read(f, h_EREC_ALL[flavor_interaction_name(fl_int)], "h_EREC_" + idstr);

      // extract 3D hists and store scaled 2D histograms
      if (use_3D_hist)
      {
        read(f, h_EFFAREA_3D_ALL[flavor_interaction_name(fl_int)], "h_EFFAREA_" + idstr);
        read(f, h_BKGREC_3D_ALL["ATMNU_conv_" + flavor_interaction_name(fl_int)], "h_ATMNU_reco_conv_" + idstr);
        read(f, h_BKGREC_3D_ALL["ATMNU_prompt_" + flavor_interaction_name(fl_int)], "h_ATMNU_reco_prompt_" + idstr);
        // read( f, h_BKGREC_3D_ALL["ATMMUON_" + flavor_interaction_name(fl_int)]      , "h_ATMMU_reco_"        + rec_type );

        // convert to 2d histograms
        h_EFFAREA_ALL[flavor_interaction_name(fl_int)] = *((TH2D *)h_EFFAREA_3D_ALL[flavor_interaction_name(fl_int)].Project3D("yx"));
        h_BKGREC_ALL["ATMNU_conv_" + flavor_interaction_name(fl_int)] = *((TH2D *)h_BKGREC_3D_ALL["ATMNU_conv_" + flavor_interaction_name(fl_int)].Project3D("yx"));
        h_BKGREC_ALL["ATMNU_prompt_" + flavor_interaction_name(fl_int)] = *((TH2D *)h_BKGREC_3D_ALL["ATMNU_prompt_" + flavor_interaction_name(fl_int)].Project3D("yx"));
        //  h_BKGREC_ALL["ATMMUON_" + flavor_interaction_name(fl_int)] =  *((TH2D*) h_BKGREC_3D_ALL["ATMMUON_" + flavor_interaction_name(fl_int)].Project3D("yx") );

        h_EFFAREA_ALL[flavor_interaction_name(fl_int)].Scale(1.0 / h_EFFAREA_3D_ALL[flavor_interaction_name(fl_int)].GetNbinsZ());
      }
      // assume 2D histograms
      else
      {
        read(f, h_EFFAREA_ALL[flavor_interaction_name(fl_int)], "h_EFFAREA_" + idstr);
        read(f, h_BKGREC_ALL["ATMNU_conv_" + flavor_interaction_name(fl_int)], "h_ATMNU_reco_conv_" + idstr);
        read(f, h_BKGREC_ALL["ATMNU_prompt_" + flavor_interaction_name(fl_int)], "h_ATMNU_reco_prompt_" + idstr);
        // read( f, h_BKGREC_ALL["ATMMUON_" + flavor_interaction_name(fl_int)]      , "h_ATMMU_reco_"        + rec_type );
      }
    }
    else
    {
      if (use_3D_hist)
      {
        read(f, h_BKGREC_3D_ALL["ATMMUON_" + flavor_interaction_name(fl_int)], "h_ATMMU_reco_" + rec_type);
        h_BKGREC_ALL["ATMMUON_" + flavor_interaction_name(fl_int)] = *((TH2D *)h_BKGREC_3D_ALL["ATMMUON_" + flavor_interaction_name(fl_int)].Project3D("yx"));
      }
      else
      {

        read(f, h_BKGREC_ALL["ATMMUON_" + flavor_interaction_name(fl_int)], "h_ATMMU_reco_" + rec_type);
      }
    }
  }

  /*!
  Function adds interaction contributions to one histogram
  */
  void add_histograms(string flavor, vector<flavor_interaction> interactions, string rec_type)
  {

    cout << "adding histograms" << endl;
    cout << "interaction " << flavor_interaction_name(interactions[0]) << endl;

    string idstr = flavor + "_" + rec_type;

    if (startswith(flavor, "a") || startswith(flavor, "n"))
    {

      h_COSZENREC = *((TH3D *)h_COSZENREC_ALL[flavor_interaction_name(interactions[0])].Clone(("h_COSZENREC_ALL_" + idstr).c_str()));
      h_EREC = *((TH3D *)h_EREC_ALL[flavor_interaction_name(interactions[0])].Clone(("h_EREC_ALL_" + idstr).c_str()));
      h_ANGERR = *((TH3D *)h_ANGERR_ALL[flavor_interaction_name(interactions[0])].Clone(("h_ANGERR_ALL_" + idstr).c_str()));
      h_EFFAREA = *((TH2D *)h_EFFAREA_ALL[flavor_interaction_name(interactions[0])].Clone(("h_EFFAREA_ALL_" + idstr).c_str()));

      h_BKGREC["ATMNU_conv"] = *((TH2D *)h_BKGREC_ALL["ATMNU_conv_" + flavor_interaction_name(interactions[0])].Clone(("h_ATMNU_reco_conv_ALL_" + idstr).c_str()));
      h_BKGREC["ATMNU_prompt"] = *((TH2D *)h_BKGREC_ALL["ATMNU_prompt_" + flavor_interaction_name(interactions[0])].Clone(("h_ATMNU_reco_prompt_ALL_" + idstr).c_str()));

      if (interactions.size() > 1)
      {
        for (unsigned i = 1; i < interactions.size(); i++)
        {
          cout << "interaction " << flavor_interaction_name(interactions[i]) << endl;

          h_COSZENREC.Add(&h_COSZENREC_ALL[flavor_interaction_name(interactions[i])]);
          h_EREC.Add(&h_EREC_ALL[flavor_interaction_name(interactions[i])]);
          h_ANGERR.Add(&h_ANGERR_ALL[flavor_interaction_name(interactions[i])]);

          h_EFFAREA.Add(&h_EFFAREA_ALL[flavor_interaction_name(interactions[i])]);
          h_BKGREC["ATMNU_conv"].Add(&h_BKGREC_ALL["ATMNU_conv_" + flavor_interaction_name(interactions[i])]);
          h_BKGREC["ATMNU_prompt"].Add(&h_BKGREC_ALL["ATMNU_prompt_" + flavor_interaction_name(interactions[i])]);
        }
      }
    }
    else if (startswith(flavor, "m"))
    {
      h_BKGREC["ATMMUON"] = *((TH2D *)h_BKGREC_ALL["ATMMUON_" + flavor_interaction_name(interactions[0])].Clone(("h_ATMMUON_" + idstr).c_str())); // make sure to not double add the muons!
    }
    else
    {

      fatal(" Error: flavour ", flavor, " not implemented");
    }
  }

  /*!
  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)
  {
    dbg("");
    double a = interpolate(&h_EREC, E_true, cos_theta_true, E_rec, true);
    double b = interpolate(&h_COSZENREC, E_true, cos_theta_true, cos_theta_rec, true);
    return a * b / (2 * pi);
  }

  /*! P(E_rec, cos_theta_rec | E_true, cos_theta_true ) */
  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_EREC.GetBin(bin_E_true, bin_cos_theta_true, bin_E_rec);
    double a = h_EREC.GetBinContent(bb);
    if (!isfinite(a))
      fatal("h_errr failed: ", &h_EREC, 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);
  }

  /*! P ( E_rec, rec_angle_from_true | E_true, cos_theta ) */
  // virtual double dP_dlogErec_dOmega_reco_error( double E_rec, double reco_error_deg, double E_true, double cos_theta_true)
  // {
  //   return dP_dlogErec_dOmega_reco_error_log( E_rec, log10(reco_error_deg), E_true, cos_theta_true );
  // }
  /*! P ( E_rec, log10(rec_angle_from_true) | E_true, cos_theta ) */
  virtual double dP_dlogErec_dOmega_reco_error_log(double E_rec, double log_reco_error_deg, double E_true, double cos_theta_true)
  {
    // MEMCHECK
    // print( E_rec, log_reco_error_deg, E_true, cos_theta_true );
    double dP_dOmega = psf.dP_dOmega(E_true, log_reco_error_deg);
    double r = interpolate(&h_EREC, E_true, cos_theta_true, E_rec, true) * dP_dOmega;
    if (::isnan(r))
    {
      fatal("nan detected in dP_dlogErec_dOmega_reco_error_log", E_true, cos_theta_true, E_rec, dP_dOmega);
    }
    return r;
  }

  /*! for a given true energy and declination, return the fraction
      of events with Erec > cutoff_erec */

  double erec_fraction_above_sindecl(double loge, double sindecl, double cutoff_logerec)
  {
    int bx = h_EREC_sindecl.GetXaxis()->FindBin(loge);
    int by = h_EREC_sindecl.GetYaxis()->FindBin(sindecl);

    double n(0), m(0);

    for (int bz = 1; bz <= h_EREC_sindecl.GetNbinsZ(); bz++)
    {
      // double logerec = h_EREC_sindecl.GetZaxis()->GetBinCenter( bz );

      double xlow  = h_EREC_sindecl.GetZaxis()->GetBinLowEdge(bz);
      double xup   = h_EREC_sindecl.GetZaxis()->GetBinUpEdge(bz);
      double width = h_EREC_sindecl.GetZaxis()->GetBinWidth(bz);

      int bin = h_EREC_sindecl.GetBin(bx, by, bz);
      double c = h_EREC_sindecl.GetBinContent(bin);

      if (xlow > cutoff_logerec)
        n += c; // the whole bin is above the cut
      else if (xup > cutoff_logerec)
        n += c * (xup - cutoff_logerec) / width;
      m += c;

      // print ("erec/n/m", bx,by,bz, bin, logerec, n, m, sindecl, h_EREC_sindecl.Integral()  );
    }
    if (m == 0)
      return 0;
    return n / m;
  }

  /*! for a given true energy and declination, return the fraction
      of events with Erec within an energy bin */
  double erec_fraction_range_above_sindecl(double loge, double sindecl, double cutoff_low_logerec, double cutoff_high_logerec)
  {
    int bx = h_EREC_sindecl.GetXaxis()->FindBin(loge);
    int by = h_EREC_sindecl.GetYaxis()->FindBin(sindecl);

    double n(0), m(0);

    for (int bz = 1; bz <= h_EREC_sindecl.GetNbinsZ(); bz++)
    {
      double xlow = h_EREC_sindecl.GetZaxis()->GetBinLowEdge(bz);
      double xup  = h_EREC_sindecl.GetZaxis()->GetBinUpEdge(bz);
      int bin     = h_EREC_sindecl.GetBin(bx, by, bz);
      double c    = h_EREC_sindecl.GetBinContent(bin);

      // which fraction of our bin overlaps with the requested enery interval?
      Interval ival( xlow, xup ); // from statutil.hh
      double weight = ival.fraction_in_range( cutoff_low_logerec,  cutoff_high_logerec );

      if (weight <0 || weight > 1 ) fatal ("programming error weight", weight);

      n += weight * c;
      m += c;
    }
    
    if (m == 0) return 0;
    return n / m;
  }



  /*! The number of background events in a cone of size cutoff_angle and with
   *  Erec > cutoff */

  double background_rate(double livetime,
                         double declination,
                         double cutoff_angle,
                         double cutoff_logerec)
  {
    int by = h_BKGREC_sindecl.GetYaxis()->FindBin(sin(declination));

    double r(0), f(0);

    int firstbin = h_BKGREC_sindecl.GetXaxis()->FindBin(cutoff_logerec);

    for (int bx = firstbin; bx <= h_BKGREC_sindecl.GetNbinsX(); bx++)
    {

      if (bx == firstbin)
        f = 1 - (cutoff_logerec - h_BKGREC_sindecl.GetXaxis()->GetBinLowEdge(bx)) / h_BKGREC_sindecl.GetXaxis()->GetBinWidth(bx);

      else
        f = 1;

      // print("f=",f);
      r += f * h_BKGREC_sindecl.GetBinContent(bx, by);
    }

    double g = 4 * pi * pow(sin(cutoff_angle / 2), 2) / (2 * pi * h_BKGREC_sindecl.GetYaxis()->GetBinWidth(by));
    return r * g * livetime / (365.0 * 24 * 60 * 60);
  }

  double background_rate_range(double livetime,
                               double declination,
                               double angle_min,
                               double angle_max,
                               double logerec_min,
                               double logerec_max)
  {
    int by = h_BKGREC_sindecl.GetYaxis()->FindBin(sin(declination));

    double r(0), f(0);

    int firstbin = h_BKGREC_sindecl.GetXaxis()->FindBin(logerec_min);
    int lastbin = h_BKGREC_sindecl.GetXaxis()->FindBin(logerec_max);

    if (firstbin == lastbin)
    {
      // print("oops");
      f = (logerec_max - logerec_min) / h_BKGREC_sindecl.GetXaxis()->GetBinWidth(firstbin);
      r += f * h_BKGREC_sindecl.GetBinContent(firstbin, by);
    }
    else
    {

      for (int bx = firstbin; bx <= lastbin; bx++)
      {

        if (bx == firstbin)
          f = 1 - (logerec_min - h_BKGREC_sindecl.GetXaxis()->GetBinLowEdge(bx)) / h_BKGREC_sindecl.GetXaxis()->GetBinWidth(bx);

        else if (bx == lastbin)
          f = (logerec_max - h_BKGREC_sindecl.GetXaxis()->GetBinLowEdge(bx)) / h_BKGREC_sindecl.GetXaxis()->GetBinWidth(bx);

        else
          f = 1;

        // print("bx", bx ,"f=",f, "bincont", h_BKGREC_sindecl.GetBinContent( bx, by ), "r", r);
        r += f * h_BKGREC_sindecl.GetBinContent(bx, by);
      }
    }

    double g = (4 * pi * pow(sin(angle_max / 2), 2) - 4 * pi * pow(sin(angle_min / 2), 2)) / (2 * pi * h_BKGREC_sindecl.GetYaxis()->GetBinWidth(by));
    return r * g * livetime / (365.0 * 24 * 60 * 60);
  }

  /*! return a human-readable table with rate information for flux */
  Table compute_rates_table(Flux &flux,
                            double livetime, // seconds
                            bool point,
                            double declination,
                            double cutoff_angle,
                            double cutoff_logerec); // in cc

  /*! Return total signal and background, optionally filling table */
  vector<double> compute_rates(Flux &flux,
                               double livetime, // seconds
                               bool point,
                               double declination,
                               double cutoff_angle,
                               double cutoff_logerec,
                               Table *T = nullptr); // in cc

  /*! Return total signal and background in energy/angle bin */
  vector<double> compute_rates_angle_energy_range(Flux &flux,
                                                  double livetime, // seconds
                                                  bool point,
                                                  double declination,
                                                  double angle_min,
                                                  double angle_max,
                                                  double logerec_min,
                                                  double logerec_max); // in cc

  virtual ~DefaultDetResponse()
  {
    // print ("dtor of DefaultDetResponse");
  }
  ClassDef(DefaultDetResponse, 1);
};

/*! A collection of DetResponse's for each channel and neutrino flavor */
struct FullDetResponse
{

  string name;
  // vector<flavor>  flavors={ muon, nue, anue, numu, anumu, nutau, anutau};
  vector<flavor> flavors; // let's build it from fl_int
  // we decide which ingredients are actually used by the list provided by user in fl_int
  vector<channel> channels;
  vector<flavor_interaction> fl_int;
  string filename;

  map<int, map<int, DefaultDetResponse>> response_map; // by pdg-type of incomming neutrino and by observed channeld (track, shower, ...)
  
  FullDetResponse() {}
  
  /*! 
   Initialize from File. 

  \param cos_zenith_cut : for bins above this cut (downgong = 1, upgoing = -1), the Aeff will be set to zero
  */


  FullDetResponse(string name,
                  string filename_,
                  vector<flavor_interaction> fl_int = {numuCC, numuNC, nueCC, nueNC, nutauNC, nutauCCshowerdecay, nutauCCmuondecay, anumuCC, anumuNC, anueCC, anueGLRES, anueNC, anutauNC, anutauCCshowerdecay, anutauCCmuondecay, muonmuon},
                  vector<channel> channels = {track, shower},
                  double src_width = 0.,
                  bool gaus = true,
                  double cos_zenith_cut = 999 )  : channels(channels)
  {
    flavors = get_flavor_from_flavor_interaction(fl_int);

    filename = filename_;
    int psf_mergebins = 2;

    for (auto chan : channels)
      for (auto flav : flavors)
      {
        if (get_flavor_interaction(flav, fl_int).size() > 0)
        {

          cout << "DEBUG in Fulldetresp " << flavor_name(flav) << " " << channel_name(chan) << endl;
          response_map[flav][chan] = DefaultDetResponse(filename,
                                                        flav,
                                                        get_flavor_interaction(flav, fl_int),
                                                        channel_name(chan),
                                                        src_width,
                                                        psf_mergebins,
                                                        gaus,
                                                        cos_zenith_cut );

          response_map[flav][chan].flav = flav;
          response_map[flav][chan].chan = chan;
          response_map[flav][chan].Aeff.flav = flav;
          response_map[flav][chan].Aeff.chan = chan;
        }
      }
  }

  Det &det() { return (response_map.begin()->second.begin()->second.det); }

  DefaultDetResponse &get(flavor nutype, channel chan)
  {
    if (response_map.find(nutype) == response_map.end())
    {
      fatal("FullDetResponse object has no DetResponse for flavor", nutype);
    }
    if (response_map[nutype].find(chan) == response_map[nutype].end())
    {
      fatal("FullDetResponse for flavor", nutype, "has no channel", chan);
    }
    return response_map[nutype][chan];
  }

  double dP_dlogErec_dOmega_zenith(flavor nutype, channel chan,
                                   double E_rec, double cos_theta_rec, double E_true, double cos_theta_true)
  {
    return get(nutype, chan).dP_dlogErec_dOmega_zenith(E_rec, cos_theta_rec, E_true, cos_theta_true);
  }

  double dP_dlogErec_dOmega_reco_error(flavor nutype, channel chan,
                                       double E_rec, double reco_error, double E_true, double cos_theta_true)
  {
    return get(nutype, chan).dP_dlogErec_dOmega_reco_error(E_rec, reco_error, E_true, cos_theta_true);
  }

  double dP_dlogErec_dOmega_zenith_bins(flavor nutype, channel chan,
                                        int bin_E_rec, int bin_cos_theta_rec, int bin_E_true, int bin_cos_theta_true)
  {
    return get(nutype, chan).dP_dlogErec_dOmega_zenith_bins(bin_E_rec, bin_cos_theta_rec, bin_E_true, bin_cos_theta_true);
  }

  double dP_dlogErec_dOmega_reco_error_log(flavor nutype, channel chan,
                                           double E_rec, double log_reco_error_deg, double E_true, double cos_theta_true)
  {
    return get(nutype, chan).dP_dlogErec_dOmega_reco_error_log(E_rec, log_reco_error_deg, E_true, cos_theta_true);
  }

  EffectiveArea &Aeff(flavor nutype, channel chan)
  {
    return get(nutype, chan).Aeff;
  }

  void show_rates(Flux &flux,
                  double livetime,
                  bool point,
                  double declination,
                  double cutoff_angle,
                  double cutoff_logerec);

  // void compute_rates( Flux& flux,         bool point, double declination = -99 , double cutoff_angle = pi , double cutoff_erec = 0  ); // in cc

  void compute_rates_summary(Flux &flux, bool point, double declination = -99, double max_angle = pi); // in cc
  virtual ~FullDetResponse(){};
  ClassDef(FullDetResponse, 1);
};
#endif