#ifndef __JROOT__JROOTTESTKIT__
#define __JROOT__JROOTTESTKIT__

#include <cstddef>
#include <vector>
#include <algorithm>

#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wall"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TRandom3.h"
#pragma GCC diagnostic pop


/**
 * \author mdejong
 */

namespace JROOT {}
namespace JPP { using namespace JROOT; }

namespace JROOT {

  /**
   * Fill 1D histogram according Poisson statistics with expectation values from given 1D function.
   *
   * \param  h1            histogram
   * \param  f1            function
   */
  template<class T>
  inline void FillRandom(TH1* h1, const T& f1)
  {
    for (Int_t ix = 1; ix <= h1->GetXaxis()->GetNbins(); ++ix){

      const double x = h1->GetXaxis()->GetBinCenter(ix);

      h1->SetBinContent(ix, gRandom->Poisson(f1(x)));

      if (h1->GetSumw2N() != 0) {
	h1->SetBinError(ix, sqrt(h1->GetBinContent(ix)));
      }
    }
  }


  /**
   * Fill 2D histogram according Poisson statistics with expectation values from given 2D function.
   *
   * \param  h2            histogram
   * \param  f2            function
   */
  template<class T>
  inline void FillRandom(TH2* h2, const T& f2)
  {
    for (Int_t ix = 1; ix <= h2->GetXaxis()->GetNbins(); ++ix){
      for (Int_t iy = 1; iy <= h2->GetYaxis()->GetNbins(); ++iy){

	const double x = h2->GetXaxis()->GetBinCenter(ix);
	const double y = h2->GetYaxis()->GetBinCenter(iy);
	
	h2->SetBinContent(ix, iy, gRandom->Poisson(f2(x,y)));

	if (h2->GetSumw2N() != 0) {
	  h2->SetBinError(ix, iy, sqrt(h2->GetBinContent(ix,iy)));
	}
      }
    }
  }


  /**
   * Fill 3D histogram according Poisson statistics with expectation values from given 3D function.
   *
   * \param  h3            histogram
   * \param  f3            function
   */
  template<class T>
  inline void FillRandom(TH3* h3, const T& f3)
  {
    for (Int_t ix = 1; ix <= h3->GetXaxis()->GetNbins(); ++ix){
      for (Int_t iy = 1; iy <= h3->GetYaxis()->GetNbins(); ++iy){
	for (Int_t iz = 1; iz <= h3->GetZaxis()->GetNbins(); ++iz){

	  const double x = h3->GetXaxis()->GetBinCenter(ix);
	  const double y = h3->GetYaxis()->GetBinCenter(iy);
	  const double z = h3->GetZaxis()->GetBinCenter(iz);
	
	  h3->SetBinContent(ix, iy, iz, gRandom->Poisson(f3(x,y,z)));

	  if (h3->GetSumw2N() != 0) {
	    h3->SetBinError(ix, iy, iy, sqrt(h3->GetBinContent(ix,iy,iz)));
	  }
	}
      }
    }
  }


  /**
   * Fill 1D histogram according PDF as given 1D function.
   *
   * \param  h1            histogram
   * \param  f1            function
   * \param  ns            number of entries
   */
  template<class T>
  inline void FillRandom(TH1* h1, const T& f1, const size_t ns)
  {
    using namespace std;

    struct tuple {

      Double_t x;
      double   W;

      inline bool operator<(const double W) { return this->W < W; }
    };

    vector<tuple> buffer;

    double W = 0.0;

    for (Int_t ix = 1; ix <= h1->GetXaxis()->GetNbins(); ++ix){

      const Double_t x = h1->GetXaxis()->GetBinCenter(ix);
      
      W += f1(x);

      buffer.push_back({x, W});
    }

    for (size_t i = 0; i != ns; ++i) {

      typename vector<tuple>::const_iterator p = lower_bound(buffer.begin(), buffer.end(), W * gRandom->Rndm());

      h1->Fill(p->x);
    }
  }


  /**
   * Fill 2D histogram according PDF as given 2D function.
   *
   * \param  h2            histogram
   * \param  f2            function
   * \param  ns            number of entries
   */
  template<class T>
  inline void FillRandom(TH2* h2, const T& f2, const size_t ns)
  {
    using namespace std;

    struct tuple {

      Double_t x;
      Double_t y;
      double   W;

      inline bool operator<(const double W) { return this->W < W; }
    };

    vector<tuple> buffer;

    double W = 0.0;

    for (Int_t ix = 1; ix <= h2->GetXaxis()->GetNbins(); ++ix){
      for (Int_t iy = 1; iy <= h2->GetYaxis()->GetNbins(); ++iy){

	const Double_t x = h2->GetXaxis()->GetBinCenter(ix);
	const Double_t y = h2->GetYaxis()->GetBinCenter(iy);
      
	W += f2(x,y);

	buffer.push_back({x, y, W});
      }
    }

    for (size_t i = 0; i != ns; ++i) {

      typename vector<tuple>::const_iterator p = lower_bound(buffer.begin(), buffer.end(), W * gRandom->Rndm());

      h2->Fill(p->x, p->y);
    }
  }


  /**
   * Fill 3D histogram according PDF as given 3D function.
   *
   * \param  h3            histogram
   * \param  f3            function
   * \param  ns            number of entries
   */
  template<class T>
  inline void FillRandom(TH3* h3, const T& f3, const size_t ns)
  {
    using namespace std;

    struct tuple {

      Double_t x;
      Double_t y;
      Double_t z;
      double   W;

      inline bool operator<(const double W) { return this->W < W; }
    };

    vector<tuple> buffer;

    double W = 0.0;

    for (Int_t ix = 1; ix <= h3->GetXaxis()->GetNbins(); ++ix){
      for (Int_t iy = 1; iy <= h3->GetYaxis()->GetNbins(); ++iy){
	for (Int_t iz = 1; iz <= h3->GetYaxis()->GetNbins(); ++iz){

	  const Double_t x = h3->GetXaxis()->GetBinCenter(ix);
	  const Double_t y = h3->GetYaxis()->GetBinCenter(iy);
	  const Double_t z = h3->GetZaxis()->GetBinCenter(iz);
      
	  W += f3(x,y,z);

	  buffer.push_back({x, y, z, W});
	}
      }
    }

    for (size_t i = 0; i != ns; ++i) {

      typename vector<tuple>::const_iterator p = lower_bound(buffer.begin(), buffer.end(), W * gRandom->Rndm());

      h3->Fill(p->x, p->y, p->z);
    }
  }
}

#endif