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

#include "TRandom3.h"

#include "JPhysics/JK40Rates.hh"

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


/**
 * \file
 *
 * Example program to check background calculation.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  JK40Rates      rates_Hz;
  double         QE;
  double         livetime_s;
  double         precision;
  UInt_t         seed;
  int            debug;

  try {

    JParser<> zap("Example program to check background calculation.");

    zap['B'] = make_field(rates_Hz);
    zap['Q'] = make_field(QE)                  = 1.0;
    zap['T'] = make_field(livetime_s)          = 1.0e5;
    zap['e'] = make_field(precision)           = 0.02;
    zap['S'] = make_field(seed)                = 0;
    zap['d'] = make_field(debug)               = 1;

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

  gRandom->SetSeed(seed);

  NOTICE("K40 rates (original)  " << rates_Hz << endl);

  JK40Rates corrected_rates_Hz = rates_Hz;
  
  corrected_rates_Hz.correct(QE);
  
  NOTICE("K40 rates (corrected) " << corrected_rates_Hz << endl);

  if (livetime_s > 0.0) {
  
    vector<double> buffer;

    buffer.resize(corrected_rates_Hz.getUpperL1Multiplicity() + 1);
    
    const double W = 1.0 / livetime_s;
      
    for (size_t M = corrected_rates_Hz.getLowerL1Multiplicity(); M <= corrected_rates_Hz.getUpperL1Multiplicity(); ++M) {
      
      const int numberOfEvents = (int) (livetime_s * corrected_rates_Hz.getMultiplesRate(M)); 
      
      for (int number_of_events = 0; number_of_events != numberOfEvents; ++number_of_events) {
	
	int n = 0;
	
	for (size_t i = 0; i != M; ++i) {
	  if (gRandom->Rndm() <= QE) {
	    ++n;
	  }
	}
	
	buffer[n] += W;
      }
    }

    NOTICE("simulated" << endl);
    
    for (size_t M = corrected_rates_Hz.getLowerL1Multiplicity(); M <= corrected_rates_Hz.getUpperL1Multiplicity(); ++M) {
      NOTICE(setw(2) << M << ' ' << FIXED(8,4) << buffer[M] << endl);
    }

    // test
    
    for (size_t M = rates_Hz.getLowerL1Multiplicity(); M <= rates_Hz.getUpperL1Multiplicity(); ++M) {
      if (fabs(buffer[M] - rates_Hz.getMultiplesRate(M)) > precision * rates_Hz.getMultiplesRate(M)) {
	FATAL("R[" << setw(2) << M << "] [Hz] " << FIXED(8,4) << buffer[M] << " != " << rates_Hz.getMultiplesRate(M) << endl);
      }
    }
  }
}