#ifndef __JSUPERNOVA_JLIGHTCURVEBACKGROUNDGENERATOR__
#define __JSUPERNOVA_JLIGHTCURVEBACKGROUNDGENERATOR__

#include <vector>

#include "TFile.h"
#include "TKey.h"

#include "TH1D.h"
#include "TH2F.h"
#include "TF1.h"

#include "TRandom3.h"

#include "JRipple.hh"

/**
 * \file
 * 
 * This file provides a class to emulate a L0 (L1) background on the ms time scale for an arbitrarily sized detector.
 * The emulation is based on the sampling of real data (in principle from a small-sized detector) through the JRipple application.
 * Samples of given duration are taken in a contiguous sequence and summed bin-by-bin as they were coming from different DUs,
 * in this way time-correlations are used to emulate space-correlations which are not possible to probe due to the detector limited size.
 *
 * \author mlincett
 */

namespace JSUPERNOVA {

  using namespace std;

  typedef vector<vector <double> > bg_type; 

  /**
   * Load histogram values to vector, each bin is converted to an element.
   * No other information is preserved (binning, axis semantics, under/overflows).
   *
   * \param  in           input histogram
   * \return              vector with one value per bin              
   */
  vector<double> loadHistogram(TH1D* in) {
    int n = in->GetNbinsX();
    
    vector<double> out(n);
    
    for (int i = 0; i < n; i++) {
      out[i] = in->GetBinContent(i + 1);
    }

    return out;
  }

  /**
   * Class to emulate L0 background for an arbitrarily sized detector.
   */

  class JLightCurveBackgroundGenerator {
  private:
    vector<vector<double> > rt;
    vector<vector<double> > nc;

    vector<h2d_t*> rt2d;

    int    TWindow_ms;
    int    genRatio;  
    int    outputSize;

    double genScale;  
 
    TRandom* rnd;

  public:

    /**
     * Default constructor
     *
     * The input file contains a sequence of histograms with sampled L0 data;
     * Each histogram is loaded as a vector, which will be considered as an individual L0 dataset.
     *
     * \param  in           input file from JRipple
     * \param twoDim        make 2D histogram with distribution of hit time differences vs time
     */
    JLightCurveBackgroundGenerator(TFile* in, bool twoDim = false) {

      TWindow_ms = 100;
      genRatio   = 115;
      genScale   = 1;

      TString rt_tag = ".RT_DET_SUM";
      TString nc_tag = ".NC_DET_SUM";
      TString rt2d_tag = ".RT2D_DET";

      rnd = new TRandom();

      TIter iter(in->GetListOfKeys());

      for (TKey* key; (key = (TKey*) iter.Next()) != NULL; ) {

	const TString tag(key->GetName());

	// match 
	if (tag.EndsWith(rt_tag)) {

	  TString rt_hist_tag = tag;
	  TString run_tag = TString(tag, tag.Length() - rt_tag.Length());
	  
	  TString nc_hist_tag = run_tag + nc_tag;

	  TH1D* RT = (TH1D*) in->Get(rt_hist_tag);
	  TH1D* NC = (TH1D*) in->Get(nc_hist_tag);

	  vector<double> rt_buf = loadHistogram(RT);
	  vector<double> nc_buf = loadHistogram(NC);

	  int n = nc_buf.size();

	  rt.push_back( vector<double>() );
	  nc.push_back( vector<double>() );

	  // build parallel vectors 
	  // 1 NC bin <=> 100 RT bins
	  for (int i = 0; i < n; i++) {
	    // empty slices are discarded
	    if (nc_buf[i] > 0) {
	      nc.back().push_back( nc_buf[i] );
   	      for (int j = 100 * i; j < 100 * (i + 1); j++) {
		rt.back().push_back( rt_buf[j] );
	      }
	    }
	  }
	
	  if (twoDim) {
	    TString rt2d_hist_tag = run_tag + rt2d_tag;
	    h2d_t* RT2D = (h2d_t*) in->Get(rt2d_hist_tag);

	    if (RT2D != NULL) {
	      RT2D->SetDirectory(0);
	      rt2d.push_back(RT2D);
	    }
	  }

	}

      }

    }

    /**
     * Destructor
     */
    ~JLightCurveBackgroundGenerator() {
      delete rnd;
    }
    
    /**
     * Configure the duration of an output sample
     *
     * \param  T_ms sample duration in ms
     */
    void configureTimeWindow(const int T_ms) {
      TWindow_ms = T_ms;
    }

    /**
     * Set TRandom seed
     *
     * \param  uSeed seed
     */
    void setSeed(const UInt_t uSeed = 0) {
      rnd->SetSeed(uSeed);
    } 
    
    /**
     * Configure generation ratio
     *
     * \param  inputNumberOfLines  number of lines of input detector
     * \param  outputNumberOfLines  number of lines of output detector
     */
    void configureRatio(const int inputNumberOfLines, const int outputNumberOfLines) {
      outputSize   = inputNumberOfLines;
      // ratio is integer by design
      genRatio     = outputNumberOfLines / inputNumberOfLines;
      // scaling factor to compensate the ratio remainder
      genScale     = outputNumberOfLines / (1.0 * inputNumberOfLines * genRatio);
    }
    

    /**
     * Generate sample of L0 background
     * L0 data are randomly sampled from a single L0 dataset.
     *
     * \return 2D vector with hit count and number of active channels as a function of time
     */
    bg_type generate() {

      bg_type out(2, vector<double>(TWindow_ms));

      int k = rnd->Integer(rt.size());

      int limit = rt[k].size() - (genRatio * TWindow_ms);

      int start = rnd->Integer(limit);

      for (int i = 0; i < genRatio; i++) {

	for (int j = 0; j < TWindow_ms; j++) {

	  int bin = start + (i * TWindow_ms) +  j;

	  out[0][j] += genScale * rt[k].at(bin      );
	  out[1][j] += genScale * nc[k].at(bin / 100);

	}

      }

      return out;
    }

    /** 
     * Generate pure poissionian background
     */
    bg_type generate_poisson(const int domRate_Hz = 500) {

      bg_type out(2, vector<double>(TWindow_ms));

      double mu = outputSize * 18 * domRate_Hz * 1.0e-3; 

      for (unsigned i = 0; i < out[0].size(); i++) {

	double count  = rnd->Poisson(mu);
	
	out[0][i] = count;
	out[1][i] = outputSize * 18 * 31;
	
      } 

      return out;

    }

    
    /**
     * Generate sample of L0 background
     * The sampling of the L0 data is not sequential but random within the L0 dataset
     * \return 2D vector with hit count and number of active channels as a function of time
     */
    bg_type generate_shuffled(bool randomizeRun = false) {

      bg_type out(2, vector<double>(TWindow_ms));

      int nRuns = rt.size();

      int k = rnd->Integer(nRuns);

      // a series n = 'genRatio' elements is sampled randomly within the dataset
      // all the elements start with the same offset to the timeslice beginning

      int offset = rnd->Integer(100);

      for (int i = 0; i < genRatio; i++) {

	int start = 100 * rnd->Integer(nc[k].size() - ceil(TWindow_ms / 100.0)) + offset;

	for (int j = 0; j < TWindow_ms; j++) {

	  int bin = start + j;

	  out[0][j] += genScale * rt[k].at(bin      );
	  out[1][j] += genScale * nc[k].at(bin / 100);
	 	  
	}

	if (randomizeRun) {
	  k = rnd->Integer(nRuns);
	}

      }

      return out;
    }    
    
    
    /**
     * Generate 2D histogram
     */
    h2d_t* build_H2D (const int k, const int start) {

      h2d_t* proto = rt2d[k];

      int ny   = proto->GetNbinsY();
      int ymin = proto->GetYaxis()->GetXmin(); 
      int ymax = proto->GetYaxis()->GetXmax();

      h2d_t* out = new h2d_t("sample", NULL, TWindow_ms, 0, TWindow_ms, ny, ymin, ymax);	
      
      for (int i = 0; i < genRatio; i++) {

	for (int j = 0; j < TWindow_ms; j++) {

	  int xbin = 1 + start + (i * TWindow_ms) +  j;

	  for (int ybin = 1; ybin <= rt2d[k]->GetNbinsY(); ybin++) {
	      
	    double val = rt2d[k]->GetBinContent(xbin, ybin);	

	    double bcx = ((TAxis*)out->GetXaxis())->GetBinCenter(j + 1);
	    double bcy = ((TAxis*)out->GetYaxis())->GetBinCenter(ybin);
	      
	    out->Fill(bcx, bcy, val);

	  }

	}

      }

      return out;

    }


    /**
     * Build NC sequence
     */
    vector<double> build_NC (const int k, const int start) {
      
      vector<double> out(TWindow_ms);

      for (int i = 0; i < genRatio; i++) {

	for (int j = 0; j < TWindow_ms; j++) {

	  int bin = start + (i * TWindow_ms) +  j;

	  out[j] += genScale * nc[k].at(bin / 100);

	}

      }

      return out;
    }

    /** 
     * Rebin vector
     */
    static vector<double> rebin(const vector<double> &in, const int rb) {

      const int n = in.size();

      vector<double> out(n/rb);
     
      for (int i = 0; i < n; i++) {
	out[i/rb] += in[i];
      }

      return out;

    }

    /**
     * Generate 2D sample
     */
    h2d_t* generate_H2D(int rb = 1) {

      h2d_t* out = NULL;
      
      if (rt2d.size() > 0) {

	int k = rnd->Integer(rt2d.size());

	int limit = rt2d[k]->GetNbinsX() - (genRatio * TWindow_ms);

	int start = rnd->Integer(limit);

	out = build_H2D(k, start);

      }

      if (rb != 1) { out->RebinX(rb); }

      return out;
    
    }

    
    /**
     * Fit 2D sample
     */
    static bg_type fit_H2D(h2d_t* in) {

      int n = in->GetNbinsX();

      // 0 -> signal
      // 1 -> fitted background
      // 2 -> background fit error

      bg_type out(3, vector<double>(n));

      for (int bin = 1; bin <= n; bin++) {

	double sg = 0, bg = 0, er = 0;

	TH1D* h = in->ProjectionY("timeBin", bin, bin);

	if (h->GetEntries() > 0) {
	  
	  TF1 f("f", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2]))/(TMath::Sqrt(2*TMath::Pi())*[2]) + [3]");

	  f.SetParameter(0, h->GetMaximum());
	  f.SetParameter(1, h->GetMean());
	  f.SetParameter(2, h->GetRMS() * 0.25);
	  f.SetParameter(3, h->GetMinimum());
	  
	  h->Fit(&f, "Q", "same");

	  bg = f.GetParameter(3) * h->GetNbinsX();
	  sg = h->GetSumOfWeights() - bg;
	  er = f.GetParError(3) * h->GetNbinsX();

	}

	int i = bin - 1;
      
	out[0][i] = sg;
	out[1][i] = bg;
	out[2][i] = er;

	delete h;

      }
      
      return out;

    }


    /**
     * Generate fitted L1 sample 
     */

    bg_type generate_fitted(int rb = 1) {

      bg_type out;

      int k = rnd->Integer(rt.size());

      int limit = rt[k].size() - (genRatio * TWindow_ms);

      int start = rnd->Integer(limit);

      h2d_t* sample = build_H2D(k, start);
      
      if (rb != 1) { sample->RebinX(rb); }

      bg_type fit = fit_H2D(sample);
      
      out.push_back(fit[0]);

      delete sample;
      
      if (rb == 1) {
	out.push_back(build_NC(k, start));
      }	else {
	out.push_back(rebin(build_NC(k, start), rb));
      }

      out.push_back(fit[1]);
      out.push_back(fit[2]);

      return out;      
    
    }	

  };
}

#endif