#ifndef ROOTUTIL_HH_INCLUDED
#define ROOTUTIL_HH_INCLUDED

#include "TTimeStamp.h"
#include "TCanvas.h"
#include "TStopwatch.h"
#include "TSystem.h"
#include "TMemberInspector.h"
#include "TDataMember.h"
#include "TClass.h"
#include "TTree.h"
#include "TFile.h"
#include "TH3.h"
#include "TH2D.h"
#include "TFunction.h"


#include "TStreamerElement.h"
#include "TClass.h"
#include "TVirtualStreamerInfo.h"
#include "TROOT.h"
#include "TClassRef.h"


#include "TDataType.h"
#include "stringutil.hh"
#include "foreach.hh"
#include "Table.hh"
#include "Timer.hh" // coudl remove

//#include "Evt.hh" // would like to get rid of (for dump())

#include <iostream>
#include <iomanip>
#include <cstdio>

using stringutil::simple_hash;
using stringutil::file_exists;
using stringutil::contains;
using stringutil::split;
using stringutil::str;



using namespace std;


// =============================================================
//  Two functions that are useful for generation numpy arrays
// =============================================================

inline void* _avoid(void* x) { return x; } // pun intended

inline void* _vector_array_address_and_size( const void* x , long& size )
{
  if ( x == 0 ) return 0;
  vector<int>* p = (vector<int>*) x;
  size = p->size();
  if (p->size() == 0 ) return 0;
  return &(p->operator[](0));
}

inline void* cling_jit( string function_definition )
{
  gInterpreter->Declare( function_definition.c_str() );
  void* p = gInterpreter->FindSym( gROOT->GetGlobalFunction("ffff")->GetMangledName() );
  return p;
}


/*! Sanity check a histogram - no bins should contain NaN */

inline void check( TH1& h , string ss ) {
	
	if ( h.GetNbinsX() < 1) 
	{
		print ( " ---- histogram is dummy ---- ", ss );
	}

	for (int b = 1 ; b < h.GetNbinsX(); b++ ) 
	{
		double c = h.GetBinContent( b );
		if ( ::isnan(c) ) {
			print ( " ---- histogram has NaN ---- ", ss );
			print ( h.GetName() );
			for (int bb = 1 ; bb < h.GetNbinsX(); bb++ ) {
				cout << h.GetBinContent(bb) << " ";
			}
			fatal("---------------------------------------");	
		}	
	}
}


// ==================================================
//  Convoluting and smearing
// ==================================================

inline void check (TH2D& src )
{
  bool bad = false;

  for (int bx = 0; bx <= src.GetNbinsX() + 1 ; bx ++ )
    for (int by = 0; by <= src.GetNbinsY() + 1 ; by ++ )
    {
      double v = src.GetBinContent(bx, by);
      if (!isfinite( v ) )
      {
        print (" bin ", by, bx , "of", src.GetName(), " not good ", v );
        bad = true;
      }
    }

  if (bad) fatal("histogram not okay");
}


inline void check (TH3D& src )
{
  bool bad = false;

  for (int bx = 0; bx <= src.GetNbinsX() + 1 ; bx ++ )
    for (int by = 0; by <= src.GetNbinsY() + 1 ; by ++ )
      for (int bz = 0; bz <= src.GetNbinsZ() + 1 ; bz ++ )
      {
        double v = src.GetBinContent(bx, by, bz);
        if (!isfinite( v ) )
        {
          print (" bin ", by, bx , bz, "of", src.GetName(), " not good ", v );
          bad = true;
        }
      }

  if (bad) fatal("histogram not okay");
}


/*! normalize each z-slices to one */

inline void normalize_z( TH3D& src , double integral = 1 , bool accept_zero = false )
{
  for (int bx = 1; bx <= src.GetNbinsX() ; bx ++ )
    for (int by = 1; by <= src.GetNbinsY() ; by ++ )
    {
      double c = 0;

      for (int bz = 1; bz <= src.GetNbinsZ() ; bz ++ )
      {
        int b = src.GetBin(bx, by, bz);
        c += src.GetBinContent( b ) * src.GetZaxis()->GetBinWidth(bz);
      }

      if ( c == 0 ) {
        if (accept_zero) continue;
        fatal("cannot normalize histogram with sum(bins) of zero ", src.GetName() );
      }


      for (int bz = 1; bz <= src.GetNbinsZ() ; bz ++ )
      {
        int b = src.GetBin(bx, by, bz);
        src.SetBinContent( b , src.GetBinContent(b) * integral / c );
      }
    }
}

/* normalize so that the sum_x sum_y value_i * withx_i * widthy_i = 1 */

inline void normalize( TH2D& src , double integral = 1 , bool accept_zero = false)
{
  double c = src.Integral();

  if ( c == 0 ) {
    if (accept_zero) return;
    fatal("cannot normalize histogram with sum(bins) of zero ", src.GetName() );
  }

  for (int bx = 1; bx <= src.GetNbinsX() ; bx ++ )
    for (int by = 1; by <= src.GetNbinsY() ; by ++ )
    {
      int b = src.GetBin(bx, by);
      double w = src.GetXaxis()->GetBinWidth(bx) * src.GetYaxis()->GetBinWidth(by);
      src.SetBinContent( b , src.GetBinContent(b) * integral / c / w );
    }
}




/*! Take a src histogram of x,y and turn it into a histogram of x' y'
    using a function f that is is pdf( x'| y ) */

template<typename F>
inline TH1D smeared( TH1D& src, F f )
{
  TH1D dst = src;
  dst.Reset();


  for (int bx = 1; bx <= src.GetNbinsX() ; bx ++ )
  {
    double x = src.GetXaxis()->GetBinCenter( bx );
    double v = src.GetBinContent( bx );

    for (int bxx = 1; bxx <= dst.GetNbinsX() ; bxx ++ )
    {
      double xx = dst.GetXaxis()->GetBinCenter( bxx );
      double w = dst.GetXaxis()->GetBinWidth( bxx );
      dst.AddBinContent( bxx, v * f( xx, x ) * w );
    }
  }
  return dst;
}




/*! Take a src histogram of x,y and turn it into a histogram of x' y'
    using a function f that is is pdf( x', y' | x, y ) */

template<typename F>
inline TH2D smeared( TH2D& src, F f )
{
  TH2D dst = src;
  dst.SetName("dst");
  dst.Reset();

  check(dst);

  for (int bx = 1; bx <= src.GetNbinsX() ; bx ++ )
    for (int by = 1; by <= src.GetNbinsY() ; by ++ )
    {
      double x = src.GetXaxis()->GetBinCenter( bx );
      double y = src.GetYaxis()->GetBinCenter( by );
      double b = src.GetBin(bx, by);
      double v = src.GetBinContent( b );


      for (int bxx = 1; bxx <= dst.GetNbinsX() ; bxx ++ )
        for (int byy = 1; byy <= dst.GetNbinsY() ; byy ++ )
        {
          double xx = dst.GetXaxis()->GetBinCenter( bxx );
          double yy = dst.GetYaxis()->GetBinCenter( byy );
          double bb = dst.GetBin(bxx, byy);
          double ww = dst.GetXaxis()->GetBinWidth( bxx ) * dst.GetYaxis()->GetBinWidth( byy );
          dst.AddBinContent( bb, v * f( xx, yy, x, y ) * ww );
        }
    }


  check(dst);
  return dst;
}


/*! version of smeared, where the function f takes bin-numbers
    directly -- only works if f uses histograms with identical binnin to src */

template<typename F>
inline TH2D smeared_bin( TH2D& src, F f  )
{
  TH2D dst = src;
  dst.SetName("dst");
  dst.Reset();

  for (int bx = 1; bx <= src.GetNbinsX() ; bx ++ )
    for (int by = 1; by <= src.GetNbinsY() ; by ++ )
    {
      double b = src.GetBin(bx, by);
      double v = src.GetBinContent( b );
      double w = src.GetXaxis()->GetBinWidth( bx ) * src.GetYaxis()->GetBinWidth( by );

      if (!isfinite(v) || !isfinite(w)) {
        fatal("non-finite value in", src.GetName() , bx, by , b, v, w );
      }


      for (int bxx = 1; bxx <= dst.GetNbinsX() ; bxx ++ )
        for (int byy = 1; byy <= dst.GetNbinsY() ; byy ++ )
        {
          double bb = dst.GetBin(bxx, byy);

          double vv = f(bxx, byy, bx, by);

          if ( !isfinite( vv )) {
            fatal("non-finite value in pdf function f", bxx, byy, bx, by );
          }

          dst.AddBinContent( bb, v * vv * w );
        }
    }

  check(dst);
  return dst;
}


/*! split along the x-axis */
inline vector<TH1D>  split_x( TH2D& src )
{
  vector<TH1D> dst;
  for ( int bx = 0 ; bx <= src.GetXaxis()->GetNbins() + 1 ; bx ++ ) // include overflow
  {
    dst.push_back( * src.ProjectionY( ("slicex" + str(bx)).c_str() , bx, bx) );
  }
  return dst;
}

/*! split along the y-axis */
inline vector<TH1D>  split_y( TH2D& src )
{
  vector<TH1D> dst;
  for ( int by = 0 ; by <= src.GetYaxis()->GetNbins() + 1 ; by ++ ) // include overflow
  {
    dst.push_back( * src.ProjectionX( ("slicey" + str(by)).c_str() , by, by) );
  }
  return dst;
}


inline vector<TH2D> _split_hist(  TH3D& src , TAxis* A, string opt )
{
  vector<TH2D> dst;

  int r0 = A->GetFirst();
  int r1 = A->GetLast();

  for ( int b = 0 ; b <= A->GetNbins() + 1 ; b ++ ) // include overflow
  {
    A->SetRange(b, b);
    dst.push_back( * (TH2D*) src.Project3D(opt.c_str() ) );
    dst.rbegin()->SetName( (opt + "slice" +str(b)).c_str() );
  }

  A->SetRange(r0, r1);
  return dst;
}



inline vector<TH2D>  split_x( TH3D& src )
{
  string opt = "ZY";
  TAxis* A   = src.GetXaxis();
  return _split_hist( src, A, opt );
}

inline vector<TH2D>  split_y( TH3D& src )
{
  string opt = "ZX";
  TAxis* A   = src.GetYaxis();
  return _split_hist( src, A, opt );
}

inline vector<TH2D>  split_z( TH3D& src )
{
  string opt = "YX";
  TAxis* A   = src.GetZaxis();
  return _split_hist( src, A, opt );
}

/*! Return a th2d that consists of one z-layer at the specified value.
    The x, and y axes of the 2-distogram correspond to the x and y of the 3d- one. */

inline TH2D slice_at_z( TH3D& src, double z )
{
  int b = src.GetZaxis()->FindBin(z);

  if ( b < 1 || b > src.GetZaxis()->GetNbins() )
  {
    print ("z out of range", z);
  }

  src.GetZaxis() -> SetRange(b, b);
  TH2D* rp =  ( TH2D* ) src.Project3D("YX");
  rp->SetDirectory(0);
  TH2D r = *rp;
  delete rp;
  return r;
}


inline TH1D slice_at_xy( TH3D& src, double x, double y )
{
  int bx = src.GetXaxis()->FindBin(x);
  int by = src.GetYaxis()->FindBin(y);
  src.GetXaxis() -> SetRange(bx, bx);
  src.GetYaxis() -> SetRange(by, by);

  TH1D *rp = ( TH1D* ) src.Project3D("Z");
  rp->SetDirectory(0);

  TH1D r = *rp;
  delete rp;
  return r;
}



/*! stack a vector of 2-d histograms in z-direction to form a 3d histogram */
inline TH3D join( vector<TH2D>& v , string ztitle, double zmin, double zmax , bool overflow = true )
{
  TH2D& h = v[0];

  TH3D r("joined", "joined",
         h.GetXaxis()->GetNbins(), h.GetXaxis()->GetXmin(), h.GetXaxis()->GetXmax() ,
         h.GetYaxis()->GetNbins(), h.GetYaxis()->GetXmin(), h.GetYaxis()->GetXmax() ,
         v.size() - overflow * 2 , zmin,                    zmax );

  r.GetXaxis()->SetTitle( h.GetXaxis()->GetTitle() );
  r.GetYaxis()->SetTitle( h.GetYaxis()->GetTitle() );
  r.GetXaxis()->SetName( h.GetXaxis()->GetName() );
  r.GetYaxis()->SetName( h.GetYaxis()->GetName() );
  r.GetZaxis()->SetTitle(ztitle.c_str());
  r.GetZaxis()->SetTitle(ztitle.c_str());


  for (unsigned i = 0 ; i < v.size() ; i++ )
  {
    int bz = i + (!overflow); // if v conttains overflow start at 0, otherwise at 1
    for (int bx = 1; bx <= v[i].GetNbinsX() ; bx ++ )
      for (int by = 1; by <= v[i].GetNbinsX() ; by ++ )
      {
        int b = r.GetBin(bx, by, bz);
        r.SetBinContent( b, v[i].GetBinContent( v[i].GetBin( bx, by ) ) );
      }
  }
  return r;
}





// ==================================================
//  Get an object from a root file, reasonably
// ==================================================

template<typename T>
inline void getobj( TFile& f, string name, T& dst )
{
  getobj( f, dst, name );
}
  
template<typename T>
inline void getobj( TFile& f , T& dst, string name )
{
  T* p = (T*) f.Get( name.c_str() );
  if (!p)
  {
    f.ls();
    fatal("failed to get object", name );
  }

  dst = *p;
  dst.SetDirectory(0);
}


// ==================================================
//  compile and load some c++ code on the fly
// ==================================================

bool inline cxx_file( string filename,
                      string gcc_opts       = "-Wno-ignored-qualifiers -std=c++11 ",
                      string aclic_options = "" )
{
  string S =
    "cd $BuildDir; "
    "g++ -c $Opt -pipe -m64 -W -fPIC -DR__HAVE_CONFIG -pthread $IncludePath $SourceFiles; "
    "g++ $ObjectFiles -shared -Wl,-soname,$LibName.so -m64 -O3  -Wl,--no-undefined -Wl,--as-needed $LinkedLibs -o $SharedLib ";

  cout << S << endl;

  gSystem->SetBuildDir( gSystem->GetBuildDir(), true );
  gSystem->SetMakeSharedLib( S.c_str() );
  gSystem->SetFlagsOpt( gcc_opts.c_str() );
  gSystem->SetIncludePath("-I$ROOTSYS/include -I$AADIR -I$AADIR/evt -I$AADIR/util");
  return gSystem->CompileMacro ( filename.c_str() , aclic_options.c_str() );
}

bool inline cxx( string code,
                 string filename       = "",
                 string gcc_opts       = "-Wno-ignored-qualifiers",
                 string aclic_options = "O",
                 bool use_cling = true)
{
  if ( use_cling )
  {
    return gInterpreter->Declare( code.c_str() );
  }


  // for newer compilers, it may help to uncomment the following
  // gcc_opts += " -std=c++11";


  if (!contains( code, "Evt.hh" )) code = "#include \"Evt.hh\"\n" + code;

  bool writefile = true;

  if ( filename == "" )
  {
    filename = form("tmpsrc%u.cc", simple_hash( code ) );
    if ( file_exists( filename ) ) writefile = false;
  }

  if ( writefile )
  {
    ofstream outfile( filename.c_str() );
    outfile << code << endl;
  }

  return cxx_file( filename );
}

// ==================================================
//  Iterate over the bins of histogram
// ==================================================

struct Looper
{
  TH1D* h;
  bool include_overflow;

  Looper(TH1D& h_, bool include_overflow = false ) :
    h(&h_),
    include_overflow(include_overflow) {}

  struct bin
  {
    TH1* h;
    int  b;

    double center, low, up;

    bin( TH1D& h, int bin_ ): h(&h), b(bin_),
      center( h.GetBinCenter ( b ) ),
      low   ( h.GetBinLowEdge( b ) ),
      up    ( h.GetBinLowEdge( b ) + h.GetBinWidth( b ) ) {}

    bin() {}

    void set( TH1D& h_, int bin_ )
    {
      h = &h_;
      b = bin_;
      center =  h->GetBinCenter( b );
      low =   h->GetBinLowEdge( b );
      up =  h->GetBinLowEdge( b ) + h->GetBinWidth( b );
    }

    operator double()
    {
      return h->GetBinContent(b);
    }

    void operator=(double v)
    {
      h->SetBinContent( b, v );
      // return *this;
    }

    void prnt()
    {
      print("bin object", h, b, center );
    }

  };


  struct iterator
  {
    TH1D* h;
    int b;

    bin _mybin;

    iterator(TH1D& h, int binnr ) : h(&h), b(binnr) {}

    iterator operator++(int) //it++
    {
      iterator me = *this;
      b++;
      return me;
    }

    iterator operator++() // ++it;
    {
      b++;
      return *this;
    }

    bin& operator*()
    {
      _mybin.set(*h, b );
      return _mybin;
    }

    bool operator==(const iterator& i )
    {
      return h == i.h && b == i.b;
    }

    bool operator!=(const iterator& i )
    {
      // dbg ( !(operator==(i)) , b, i.b );
      return !(operator==(i));
    }

    ClassDefNV( iterator, 2 );
  };


  iterator begin()
  {
    return iterator( *h,
                     include_overflow ? 0 : 1 );
  }

  iterator end()
  {
    return iterator( *h,
                     h->GetNbinsX() + (include_overflow ? 2 : 1) );
  }

  ClassDefNV( Looper , 1 );
};


inline void print_ascii( TH1D& h ,
                         int nlines = 20 )
{
  for (int iline = nlines; iline >= 0; iline--)
  {
    double ymin = iline * h.GetMaximum() / nlines;

    string c = form("%8.2f |", ymin);

    foreach ( bin, Looper( h ) )
    {
      c += iline ? ( (bin < ymin ) ? " " : "#" ) : "-";
    }

    print(c);
  }
}



// ==================================================
//  access to TTree draw results from c++ and pyroot
// ==================================================

struct TreeVals
{
  TTree* tree;
  TreeVals(TTree* tree) : tree(tree) {}

  struct iterator
  {
    TTree* tree;
    int index;

    bool operator==( const iterator& a ) { return index == a.index; }
    bool operator!=( const iterator& a ) { return index != a.index; }

    iterator( TTree* tree , int index ) : tree(tree), index(index) {}

    iterator& operator++() {index++; return *this;}
    iterator operator++(int) // it++
    {
      iterator me = *this;
      operator++();
      return me;
    }

    vector<double> operator*()
    {
      vector<double> r;
      for (int j = 0; j < tree->GetPlayer()->GetDimension(); j++)
      {
        r.push_back( tree-> GetVal( j )[index] );
      }
      return r;
    }
  };

  iterator begin() { return iterator( tree, 0); }
  iterator end()   { return iterator( tree, tree->GetSelectedRows());}
};



inline void clamp( int& x, int xmin, int xmax )
{
  if (x < xmin ) x = xmin;
  if (x > xmax ) x = xmax;
}

inline void clamp( double& x, double xmin, double xmax )
{
  if (x < xmin ) x = xmin;
  if (x > xmax ) x = xmax;
}

inline void clamp_bincenter(double& x , TAxis* A )
{
  const double eps = 1e-8;
  const double xmin = A->GetBinCenter( 1 )              + eps;
  const double xmax = A->GetBinCenter( A->GetNbins() )  - eps;
  clamp( x, xmin, xmax );
}


inline double interpolate( TH1D* h , double x )
{
  return h->Interpolate( x );
}

inline void interpolate_derivative( TH1D* h , double x, double& f, double& df )
{

  if ( x != x )
  {
    fatal("x = nan");
  }

  int xbin = h->FindBin(x);

  if ( x < h->GetBinCenter(1) )
  {
    f = h->GetArray()[1];
    df = 0;
    return;
  }
  else if ( x >= h->GetBinCenter( h->GetNbinsX() ))
  {
    f = h->GetArray()[ h->GetNbinsX() ];
    df = 0;
    return;
  }

  double y0, y1, x0, x1;

  if ( x <= h->GetBinCenter( xbin ) )
  {
    y0 = h->GetArray()[ xbin - 1];
    y1 = h->GetArray()[ xbin    ];
    x0 = h->GetBinCenter( xbin - 1);
    x1 = h->GetBinCenter( xbin  );
  }
  else
  {
    y0 = h->GetArray()[ xbin    ];
    y1 = h->GetArray()[ xbin + 1 ];
    x0 = h->GetBinCenter( xbin  );
    x1 = h->GetBinCenter( xbin + 1);
  }


  df = (y1 - y0) / (x1 - x0);
  f  = y0 + (x - x0) * df;


  //  print ("inter", x, x0, x1, y0,y1, df, f );
}

inline string str_axis( TAxis& a )
{
  ostringstream out;
  out << "TAxis " << a.GetName() << "bins " << a.GetNbins() << " from " << a.GetXmin() << " to " << a.GetXmax();
  return out.str();
}


/* copy of TH3::Interpolate, but with decent error handling */
inline Double_t TH3_Interpolate(TH3D& H ,
                                Double_t x, Double_t y, Double_t z,
                                bool clamp = false )
{
  TAxis& fXaxis = * H.GetXaxis();
  TAxis& fYaxis = * H.GetYaxis();
  TAxis& fZaxis = * H.GetZaxis();

  // work aound bug wehn we are at the bincenter of the last bin
  const double eps = 1e-10;
  double x1 = fXaxis.GetBinCenter( fXaxis.GetNbins() );
  double y1 = fYaxis.GetBinCenter( fYaxis.GetNbins() );
  double z1 = fZaxis.GetBinCenter( fZaxis.GetNbins() );

  if ( x == x1 ) x -= eps;
  if ( y == y1 ) y -= eps;
  if ( z == z1 ) z -= eps;

  if ( clamp ) {
    clamp_bincenter( x, &fXaxis );
    clamp_bincenter( y, &fYaxis );
    clamp_bincenter( z, &fZaxis );
  }


  Int_t ubx = fXaxis.FindFixBin(x);
  if ( x < fXaxis.GetBinCenter(ubx) ) ubx -= 1;
  Int_t obx = ubx + 1;

  Int_t uby = fYaxis.FindFixBin(y);
  if ( y < fYaxis.GetBinCenter(uby) ) uby -= 1;
  Int_t oby = uby + 1;

  Int_t ubz = fZaxis.FindFixBin(z);
  if ( z < fZaxis.GetBinCenter(ubz) ) ubz -= 1;
  Int_t obz = ubz + 1;


//    if ( IsBinUnderflow(GetBin(ubx, uby, ubz)) ||
//         IsBinOverflow (GetBin(obx, oby, obz)) ) {
  if ( ubx <= 0 || uby <= 0 || ubz <= 0 ||
       obx > fXaxis.GetNbins() || oby > fYaxis.GetNbins() || obz > fZaxis.GetNbins() ) {
    print("Interpolate", "Cannot interpolate outside domain of hist", H.GetName());

    print ( str_axis(fXaxis) , x , obx );
    print ( str_axis(fYaxis) , y , oby );
    print ( str_axis(fZaxis) , z , obz );
    exit(1);
  }

  Double_t xw = fXaxis.GetBinCenter(obx) - fXaxis.GetBinCenter(ubx);
  Double_t yw = fYaxis.GetBinCenter(oby) - fYaxis.GetBinCenter(uby);
  Double_t zw = fZaxis.GetBinCenter(obz) - fZaxis.GetBinCenter(ubz);

  Double_t xd = (x - fXaxis.GetBinCenter(ubx)) / xw;
  Double_t yd = (y - fYaxis.GetBinCenter(uby)) / yw;
  Double_t zd = (z - fZaxis.GetBinCenter(ubz)) / zw;


  Double_t v[] = { H.GetBinContent( ubx, uby, ubz ), H.GetBinContent( ubx, uby, obz ),
                   H.GetBinContent( ubx, oby, ubz ), H.GetBinContent( ubx, oby, obz ),
                   H.GetBinContent( obx, uby, ubz ), H.GetBinContent( obx, uby, obz ),
                   H.GetBinContent( obx, oby, ubz ), H.GetBinContent( obx, oby, obz )
                 };


  Double_t i1 = v[0] * (1 - zd) + v[1] * zd;
  Double_t i2 = v[2] * (1 - zd) + v[3] * zd;
  Double_t j1 = v[4] * (1 - zd) + v[5] * zd;
  Double_t j2 = v[6] * (1 - zd) + v[7] * zd;


  Double_t w1 = i1 * (1 - yd) + i2 * yd;
  Double_t w2 = j1 * (1 - yd) + j2 * yd;


  Double_t result = w1 * (1 - xd) + w2 * xd;

  return result;
}



inline double interpolate( TH3D* h , double x, double y, double z , bool clamp = false )
{
  return TH3_Interpolate( *h, x, y, z , clamp );

}

inline bool err( int e = -999 )
{
  static bool r = false ;
  if (e == -999) return r;
  else return r = e;
}

#include "TError.h"
inline void errorhandler (int level, Bool_t abort, const char *location,
                          const char *msg)
{
  cout << " error " << level << endl;
  err(1);
  return;
}


inline void get_interpolation_consts( TAxis* ax, double x,
                                      int& b1, int& b2, double& w1, double& w2)
{
  b1 = ax->FindBin( x );
  const double c1 = ax->GetBinCenter( b1 );
  b2 = x > c1 ? b1 + 1 : b1 - 1;
  const double c2 =  ax->GetBinCenter( b2 );
  const double delta = c2 - c1;
  w2 = (x - c1) / delta;
  w1 = 1.0 - w2;
}

/*! returned bin numbers run from zero to nbins (not including) */

inline void get_interpolation_consts( int nbins, double xmin, double xmax,
                                      double x,
                                      int& b1, int& b2, double& w1, double& w2)
{
  const double bw = ( xmax - xmin ) / nbins;
  b1 = int ( ( x - xmin ) / bw ) ;
  const double c1 = xmin + ( b1 + 0.5 ) * bw;
  b2 = x > c1 ? b1 + 1 : b1 - 1;
  double delta = x > c1 ? bw : -bw;
  w2 = (x - c1) / delta;
  w1 = 1.0 - w2;
}



inline double cached_interpolate( TH3D* h, double x, double y, double z ) // not tested (much)
{
  static double oldx(-1e100), oldy(-1e100), oldz(-1e100);

  static int b1x, b2x; static double  w1x, w2x;
  static int b1y, b2y; static double  w1y, w2y;
  static int b1z, b2z; static double  w1z, w2z;

  if ( x != oldx ) get_interpolation_consts( h->GetXaxis(), x, b1x, b2x, w1x, w2x );
  if ( y != oldy ) get_interpolation_consts( h->GetYaxis(), y, b1y, b2y, w1y, w2y );
  if ( z != oldz ) get_interpolation_consts( h->GetZaxis(), z, b1z, b2z, w1z, w2z );

  const double c11 = h -> GetBinContent( b1x, b1y, b1z ) * w1x  +  h -> GetBinContent( b2x, b1y, b1z ) * w2x;
  const double c12 = h -> GetBinContent( b1x, b2y, b1z ) * w1x  +  h -> GetBinContent( b2x, b2y, b1z ) * w2x;
  const double c21 = h -> GetBinContent( b1x, b1y, b2z ) * w1x  +  h -> GetBinContent( b2x, b1y, b2z ) * w2x;
  const double c22 = h -> GetBinContent( b1x, b2y, b2z ) * w1x  +  h -> GetBinContent( b2x, b2y, b2z ) * w2x;

  const double c1  = c11 * w1y + c12 * w2y;
  const double c2  = c21 * w1y + c22 * w2y;

  return c1 * w1z + c2 * w2z;
}



// get interpolation for y-value at the center of a y-bin
// interpolated along x and z

inline double interpol_xz( TH3D* h, double x, double z , int ybin )
{
  static double oldx(-1e100), oldz(-1e100);

  static int b1x, b2x; static double  w1x, w2x;
  static int b1z, b2z; static double  w1z, w2z;

  if ( x != oldx ) get_interpolation_consts( h->GetXaxis(), x, b1x, b2x, w1x, w2x );
  if ( z != oldz ) get_interpolation_consts( h->GetZaxis(), z, b1z, b2z, w1z, w2z );

  const double c1 = h -> GetBinContent( b1x, ybin, b1z ) * w1x  +  h -> GetBinContent( b2x, ybin, b1z ) * w2x;
  const double c2 = h -> GetBinContent( b1x, ybin, b2z ) * w1x  +  h -> GetBinContent( b2x, ybin, b2z ) * w2x;
  return c1 * w1z + c2 * w2z;
}

// get interpolation for y-value at the center of a y-bin
// interpolated along x and z

inline void interpol_xz( TH3D* h, double x, vector<double> vz , TH1D* dst , bool doclamp = true )
{
  static double oldx(-1e100);

  if (doclamp)
  {
    clamp_bincenter( x, h->GetXaxis() );
  }

  static int b1x, b2x; static double  w1x, w2x;
  if ( x != oldx ) get_interpolation_consts( h->GetXaxis(), x, b1x, b2x, w1x, w2x );

  foreach ( z, vz )
  {
    int b1z, b2z;
    double  w1z, w2z;
    if (doclamp) clamp_bincenter( z, h->GetZaxis() );
    get_interpolation_consts( h->GetZaxis(), z, b1z, b2z, w1z, w2z );

    for (int ybin = 1; ybin <= dst->GetNbinsX(); ybin++)
    {
      const double c1 = h -> GetBinContent( b1x, ybin, b1z ) * w1x  +  h -> GetBinContent( b2x, ybin, b1z ) * w2x;
      const double c2 = h -> GetBinContent( b1x, ybin, b2z ) * w1x  +  h -> GetBinContent( b2x, ybin, b2z ) * w2x;
      dst -> SetBinContent( ybin, c1 * w1z + c2 * w2z + dst -> GetBinContent( ybin ) );
    }
  }

}

inline double intr( TH1D* h, double x )
{
  int b1 = h->FindBin( x );
  double c1 = h->GetBinCenter( b1 );

  int b2 = x > c1 ? b1 + 1 : b1 - 1;
  double c2 =  h->GetBinCenter( b2 );

  double delta = c2 - c1;

  double f2 = (x - c1) / delta;
  double f1 = 1 - f2;


  // the idea is to remember f1, f1, b1, and b2 for fast evaluation of
  // many histograms at the same x.

  // value
  print ( x, " b=", b1, b2, " f= " ,  f1, f2, " y=", h->GetBinContent( b1 ),  h->GetBinContent( b2 ) );

  double f = f1 * h->GetBinContent( b1 ) + f2 * h->GetBinContent( b2 );
  return f;
}

inline TH1D h()
{
  TH1D h("h", "h", 10, 0, 20);
  for (int i = 0; i < 10; i++) h.SetBinContent( i, (i * i * i) % 3 );
  return h;
}

inline TH1D test_intr()
{
  TH1D h0 = h();
  TH1D h1("h1", "h1", 100, 0, 20);

  for (int i = 0; i < 100; i++) h1.SetBinContent( i, intr( &h0, double(i) / 5) );
  return h1;
}


inline TH1D histogram_slice_y( TH3D* src, double x, double z ,
                               bool clamp = true )
{
  // SetErrorHandler( errorhandler );

  if (clamp)
  {
    clamp_bincenter( x, src->GetXaxis() );
    clamp_bincenter( z, src->GetZaxis() );
  }

  TAxis& ax = * ( src -> GetYaxis() );
  string name = string(src->GetName()) + "slice_y";

  TH1D r( name.c_str(), name.c_str(), ax.GetNbins(), ax.GetXbins()->GetArray() );

  for (int b = 1 ; b <= r.GetNbinsX(); b++ )
  {
    double y = r.GetBinCenter( b );
    double v;

    clamp_bincenter( y, r.GetXaxis() );
    {
      //Timer _("slice interpolate");
      v = interpolate( src, x, y, z );
    }

    r.SetBinContent( b, v );
  }
  return r;
}
//------------------------------------------------------------
// hammer projection, also known as hammer-aitoff projection
// https://en.wikipedia.org/wiki/Hammer_projection
//------------------------------------------------------------

#include "TGraph.h"
#include "TGraph2D.h"
#include "TMarker.h"
#include "TMultiGraph.h"


inline vector<double> hammer_projection( double lon, double lat )
{
  // you can't touch this!

  while ( lon < -M_PI ) lon += 2 * M_PI;
  while ( lon > M_PI  ) lon -= 2 * M_PI;

  if ( lat < -M_PI / 2 || lat > M_PI / 2 )
  {
    fatal("hammer_projection: lattitute out of range (-pi,pi) ", lat );
  }

  double d = sqrt( 1 + cos(lat) * cos(lon / 2) );

  vector<double> r(2);
  r[0] = 2 * sqrt(2) * cos( lat ) * sin( lon / 2 ) / d ;
  r[1] = sqrt(2) * sin( lat ) / d;
  return r;
}

inline void hammer_projection( TGraph* g )
{
  for (int i = 0; i < g->GetN(); i++)
  {
    vector<double> c = hammer_projection( g->GetX()[i], g->GetY()[i] );
    g->SetPoint(i, c[0], c[1] );
  }
}

inline void hammer_projection( TGraph2D* g )
{
  for (int i = 0; i < g->GetN(); i++)
  {
    vector<double> c = hammer_projection( g->GetX()[i], g->GetY()[i] );
    g->SetPoint(i, c[0], c[1] , g->GetZ()[i] );
  }
}

inline void hammer_projection( TMarker* m )
{
  vector<double> c = hammer_projection( m->GetX(), m->GetY() );
  m->SetX( c[0] );
  m->SetY( c[1] );
}

inline void hammer_projection( TMultiGraph* g)
{
  TList* l = g->GetListOfGraphs();
  for (int i = 0; i < l->GetSize(); i++) hammer_projection( (TGraph*) l->At(i) );
}



// A sequence of n evenly spaced points, the first
// being xmin, and the last xmax

inline vector<double> floatrange( double xmin, double xmax, int n )
{
  vector<double> r(n);
  r[0]   = xmin;
  r[n - 1] = xmax;
  for (int i = 1; i < n - 1; i++) r[i] = xmin + i * (xmax - xmin) / (n - 1);
  return r;
}

inline TMultiGraph* grid( int nlines_lon = 7, int nlines_lat = 5 )
{
  TMultiGraph* g = new TMultiGraph;

  const double pi = acos(-1);
  const double xsize = 2 * pi;
  const double ysize = pi;

  const int nsteps = 100;

  vector<double> lon_lines = floatrange( -xsize / 2 , xsize / 2 , nlines_lon );
  vector<double> lon_steps = floatrange( -xsize / 2 , xsize / 2 , nsteps     );

  const double dy = (ysize) / nlines_lat;
  vector<double> lat_lines = floatrange( -ysize / 2 + dy , ysize / 2 - dy  , nlines_lat );
  vector<double> lat_steps = floatrange( -ysize / 2 , ysize / 2 , nsteps     );


  foreach ( lon, lon_lines )
  {
    vector<double> x, y;
    foreach ( lat, lat_steps )
    {
      x.push_back( lon );
      y.push_back( lat );
    }
    g->Add( new TGraph( std::min( x.size(), y.size()) , &(x[0]), &(y[0]) ) );
  }

  foreach ( lat, lat_lines )
  {
    vector<double> x, y;
    foreach ( lon, lon_steps )
    {
      x.push_back( lon );
      y.push_back( lat );
    }
    g->Add( new TGraph( std::min( x.size(), y.size()) , &(x[0]), &(y[0]) ) );
  }

  return g;
}

inline TMultiGraph* sky_grid( int nlines_lon = 7, int nlines_lat = 5 )
{
  TMultiGraph* g = grid( nlines_lon, nlines_lat );
  hammer_projection( g );
  return g;
}




#endif