#ifndef HISTOCONV_INCD
#define HISTOCONV_INCD

#include "Det.hh"
#include "foreach.hh"
#include "rootutil.hh"
#include "TH2D.h"

struct convert_decl
{
  TH2D pivot_matrix;



};


/*! for every sin-decl, you have a slice of cos-theta with the bin-content
 *  indication how much time the source spends there */

inline TH2D conversion_table_costheta_sindecl( Det& det,
                                               int nbinsxy,
                                               int Nsamples )
{
  static TH2D R;

  if ( R.GetNbinsY() == nbinsxy ) return R; // todo: check the rest of the parameters

  print ("generating new conversion table for nbins=", nbinsxy );

  R = TH2D("h_conversion_table","conversion from costheta_to_sindecl", 
          nbinsxy, -1,1, nbinsxy, -1,1 );

  TTimeStamp T(2020, 1, 1, 2, 2, 2, 2 );

  for (int bx = 1; bx <= R.GetNbinsX() ; bx++ )
  {
    double sindecl = R.GetXaxis()->GetBinCenter(bx);
    double decl = asin( sindecl );  

    for ( int i = 0; i < Nsamples ; i++ )
    {
      double ra = i * 2 * pi / Nsamples;
      EquatorialCoords coords( ra, decl );

      Vec v = -coords.track_direction( det, T );
      double costheta = v.z;
      R.Fill( sindecl,costheta, 1.0/Nsamples);
    }
  } 
  return R;
}


/*! take a 2d histogram with cos-zenith on the y-axis and turn it into a histogram 
 *  with sin-decl on the y-axis
 * */

inline TH2D* h_declzen_conversion( TH2D* p = nullptr )
{
  static TH2D* x;
  if ( p ) x = p;
  return x;
}



inline TH2D convert_to_sindecl( Det& det, TH2D& H )
{
  static TH2D C; 

  if ( H.GetNbinsY() != C.GetNbinsY() )
  {
    print ( "gen new conv table", H.GetNbinsY() , C.GetNbinsY());
    C = conversion_table_costheta_sindecl( det, H.GetNbinsY(), 1000 );
    h_declzen_conversion(&C);
  }

  TH2D R = H; // copy
  R.Reset();
  R.SetName( ( string( H.GetName()) + "__sindecl" ).c_str() );
  R.GetYaxis()->SetTitle("sin(decl)");

  for (int hy = 1; hy <= H.GetNbinsY(); hy++ )
    {
      for (int ry = 1; ry <= R.GetNbinsY(); ry++ )
      {
        double w = C.GetBinContent( ry, hy );

        for (int hx=1; hx <= H.GetNbinsX(); hx++ )
        {
          double u = R.GetBinContent( hx, ry );
          R.SetBinContent( hx, ry,  u + w * H.GetBinContent( hx, hy ) );
        }
      }
    }

  return R;
}

/*! Given a 3d histogram who's y-axis contains cos-zenith, Return one
 *  with y-axis representing sin(decl) */


inline TH3D convert_to_sindecl( Det& det, TH3D& H )
{
  const TH2D& C = conversion_table_costheta_sindecl( det, H.GetNbinsY(), 1000 );

  TH3D R = H; // copy
  R.Reset();
  R.SetName( ( string( H.GetName()) + "__sindecl" ).c_str() );
  R.GetYaxis()->SetTitle("sin(decl)");

  if ( H.Integral() == 0 ) fatal("empty histogram in convert_to_sindecl");


  for (int hy = 1; hy <= H.GetNbinsY(); hy++ )
    {
      for (int ry = 1; ry <= R.GetNbinsY(); ry++ )
      {
        double w = C.GetBinContent( ry, hy );

        if ( w == 0 ) {
          continue;
        }

        for (int hz=1; hz <= H.GetNbinsZ(); hz++ )
        {  
          for (int hx=1; hx <= H.GetNbinsX(); hx++ )
          {
            double u = R.GetBinContent( hx, ry , hz );
            R.SetBinContent( hx, ry,  hz, u + w * H.GetBinContent( hx, hy, hz ) );
          }
        }
      }
    }

  return R;
}




#endif