#include <string>
#include <iostream>
#include <iomanip>
#include <vector>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TRandom3.h"

#include "km3net-dataformat/online/JDAQ.hh"
#include "km3net-dataformat/online/JDAQClock.hh"
#include "km3net-dataformat/online/JDAQTimeslice.hh"
#include "JTimeslice/JRandomTimeslice.hh"
#include "JPhysics/JK40Rates.hh"
#include "JROOT/JManager.hh"

#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

namespace {

  /**
   * Check is data are time sorted.
   *
   * \param  frame        frame
   * \return              true if sorted; else false
   */
  inline bool is_sorted(const KM3NETDAQ::JDAQSuperFrame& frame)
  {
    using namespace KM3NETDAQ;

    if (frame.size() > 1) {
      for (JDAQSuperFrame::const_iterator q = frame.begin(), p = q++; q != frame.end(); ++p, ++q) {
	if (p->getT() > q->getT()) {
	  return false;
	}
      }
    }

    return true;
  }
}


/**
 * \file
 *
 * Auxiliary program to generate KM3NETDAQ::JDAQTimeslice with random data.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  string                 outputFile;
  Long64_t               numberOfSlices;
  JK40Rates              rates_Hz;
  pair<size_t, double>   recycling;
  UInt_t                 seed;
  int                    debug;

  try { 

    JParser<> zap("Auxiliary program to generate time slices with random data.");
    
    zap['o'] = make_field(outputFile,          "output file")                                          = "timeslice.root";
    zap['n'] = make_field(numberOfSlices);
    zap['B'] = make_field(rates_Hz,            "background rates [Hz]");
    zap['N'] = make_field(recycling,           "number of recycles / time interval for sampling data") = make_pair(0, 0.0);
    zap['S'] = make_field(seed,                "seed")        = 0;
    zap['d'] = make_field(debug,               "debug")       = 0;

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  gRandom->SetSeed(seed);

  const JDAQHit::JPMT_t  PMT_L0[] =  { 0, 1 };
  const JDAQHit::JPMT_t  PMT_L1   =    2;

  const double TMin_ns = 1.0;  // Minimal time between two consecutive L(0|1) hits [ns]

  JManager<size_t, TH1D> H0(new TH1D("h0[%]", NULL, 100, 0.0, getFrameTime()));

  for (Long64_t counter = 0; counter != numberOfSlices; ++counter) {

    STATUS("slice: " << setw(10) << counter << '\r'); DEBUG(endl);

    JRandomTimeslice timeslice;

    timeslice.resize(1);


    vector<JDAQHit> buffer;

    if (rates_Hz.getSinglesRate() > 0.0) {

      double period_ns = 1.0e9 / (NUMBER_OF_PMTS * rates_Hz.getSinglesRate());

      for (double t1 = 0.0 * getFrameTime() + gRandom->Exp(period_ns); t1 < 0.5 * getFrameTime(); t1 += TMin_ns + gRandom->Exp(period_ns)) {
	buffer.push_back(JDAQHit(PMT_L0[0], t1, 0));
      }

      period_ns *= 2.0;

      for (double t1 = 0.5 * getFrameTime() + gRandom->Exp(period_ns); t1 < 1.0 * getFrameTime(); t1 += TMin_ns + gRandom->Exp(period_ns)) {
	buffer.push_back(JDAQHit(PMT_L0[1], t1, 0));
      }
    }


    map<multiplicity_type, int> L1;

    for (multiplicity_type M = rates_Hz.getLowerL1Multiplicity(); M <= rates_Hz.getUpperL1Multiplicity(); ++M) {

      if (rates_Hz.getMultiplesRate(M) > 0.0) {

	const double period_ns = 1.0e9 / rates_Hz.getMultiplesRate(M);

	for (double t1 = gRandom->Exp(period_ns); t1 < getFrameTime(); t1 += TMin_ns + gRandom->Exp(period_ns)) {

	  L1[M] += 1;

	  for (multiplicity_type i = 0; i != M; ++i) {

	    const JDAQHit hit(PMT_L1, t1, 0);

	    vector<JDAQHit>::iterator p = lower_bound(buffer.begin(), buffer.end(), hit);

	    buffer.insert(p, hit);
	  }
	}
      }
    }

    timeslice[0].add(buffer.size(), buffer.data());


    TH1D* h0 = H0[0];

    for (JDAQSuperFrame::const_iterator hit = timeslice[0].begin(); hit != timeslice[0].end(); ++hit) {
      h0->Fill((Double_t) hit->getT());
    }


    for (size_t i = 1; i <= recycling.first; ++i) {

      
      timeslice.recycle(recycling.second);


      TH1D* h0 = H0[i];

      for (JDAQSuperFrame::const_iterator hit = timeslice[0].begin(); hit != timeslice[0].end(); ++hit) {
	h0->Fill((Double_t) hit->getT());
      }


      ASSERT(is_sorted(timeslice[0]), "Test if data are sorted after recycling.");

      map<multiplicity_type, int> N1;

      for (JDAQSuperFrame::const_iterator p = timeslice[0].begin(); p != timeslice[0].end(); ) {

	if (p->getPMT() == PMT_L1) {

	  JDAQSuperFrame::const_iterator q = p;

	  int n1 = 1;

	  for (++q; q != timeslice[0].end() && q->getT() == p->getT(); ++q) {
	    if (q->getPMT() == PMT_L1) {
	      ++n1;
	    }
	  }

	  N1[n1] += 1;
	  
	  p = q;

	} else {

	  ++p;
	}
      }

      
      for (multiplicity_type M = rates_Hz.getLowerL1Multiplicity(); M <= rates_Hz.getUpperL1Multiplicity(); ++M) {

	if (debug >= debug_t) {

	  cout << "Multiplicity " << setw(2) << M << ' ' << setw(4) << L1[M] << " ?= " << setw(4) << N1[M] << endl;

	  if (L1[M] != N1[M]) {

	    int count = 0;

	    for (JDAQSuperFrame::const_iterator i = timeslice[0].begin(); i != timeslice[0].end(); ++i) {

	      if (i->getPMT() == PMT_L1) {

		cout << setw(4) << count << ' ' << setw(2) << (int) i->getPMT() << ' ' << setw(10) << i->getT() << endl;

		++count;
	      }
	    }
	  }
	}

	ASSERT(L1[M] == N1[M], "Multiplicity " << setw(2) << M << ' ' << setw(4) << L1[M] << " ?= " << setw(4) << N1[M]);
      }
    }
  }
  STATUS(endl);
  

  TFile out(outputFile.c_str(), "recreate");

  out << H0;

  out.Write();
  out.Close();
}