#include <string>
#include <iostream>
#include <iomanip>
#include <limits>

#include "JDAQ/JDAQTimesliceIO.hh"

#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JSupport.hh"

#include "JTools/JGrid.hh"
#include "JTools/JHashMap.hh"

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


/**
 * \file
 *
 * Example program to determine multiplicity probabilities from data.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  JMultipleFileScanner<JDAQTimesliceL0> inputFile;
  counter_type       numberOfEvents;
  double             T_ns;
  int                debug;

  try { 

    JParser<> zap("Example program to determine multiplicity probabilities from data.");
    
    zap['f'] = make_field(inputFile);
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['T'] = make_field(T_ns);
    zap['d'] = make_field(debug)               = 2;
    
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }


  const int nx = (int) (getFrameTime() / T_ns);

  const JGrid<double> grid(nx, 0.0, nx * T_ns);
  
  JHashMap<int, int> buffer;

  map<int, double> zmap;

  const int N = 10;
  double    W = 0.0;

  for (int i = 0; i != N; ++i) {
    zmap[i] = 0.0;
  }


  while (inputFile.hasNext()) {

    STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);

    const JDAQTimeslice* timeslice = inputFile.next();

    for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {

      buffer.clear();

      for (JDAQSuperFrame::const_iterator hit = frame->begin(); hit != frame->end(); ++hit) {
	
	const int index = grid.getIndex(hit->getT());
	
	buffer[index] += 1;
      }

      W       += grid.getSize();                        // total number of measurements
      zmap[0] += grid.getSize() - buffer.size();        // number of measurements with multiplicity zero

      for (JHashMap<int, int>::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
	zmap[i->second] += 1.0;                         // number of measurements with given multiplicity
      }
    }
  }
  STATUS(endl);

  for (int i = 0; i != N; ++i) {
    cout << setw(2) << i << ' ' << SCIENTIFIC(12,3) << zmap[i] / W << " +/- " << SCIENTIFIC(12,3) << sqrt(zmap[i]) / W << endl;
  }
}