#include <string>
#include <iostream>

#include "JPhysics/KM3NeT.hh"
#include "JPhysics/KM3NeT2D.hh"
#include "JPhysics/Antares.hh"
#include "JDetector/JPMTAnalogueSignalProcessor.hh"

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


/**
 * \file
 *
 * Example program to calculate singles rate.
 * The calculation is based on the Antares internal note ANTARES-PHYS-2006-005 by Jurgen Brunner.
 * According Antares internal note ANTARES-PHYS-2012-013, the absorption length is used (and not the attenuation length).
 * \author mdejong
 */
int main(int argc, char* argv[])
{
  using namespace std;
  using namespace JPP;

  double         bequerel;
  JPMTParameters parameters;
  int            debug;

  try {

    JProperties properties = parameters.getProperties();

    JParser<> zap("Example program to calculate singles rate.");

    zap['b'] = make_field(bequerel)     = 13.75e3;     // [m^-3 s^-1]
    zap['P'] = make_field(properties)   = JPARSER::initialised();
    zap['d'] = make_field(debug)        = 3;

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


  using namespace NAMESPACE;

  JPMTAnalogueSignalProcessor cpu(parameters);

  const double wmin = 280.0;   // minimal wavelength [nm]
  const double wmax = 700.0;   // maximal wavelength [nm]
  double       ng   =  41.0;   // average number of photons per decay in given wavelength range 
  const int    npe  =   1;     // number of photo-electrons for each decay

  ng *=  (wmax-wmin) / (wmin*wmax) * (300.0*600.0)/(600.0-300.0);

  double R[] = { 0.0, 0.0 };
  double W[] = { 0.0, 0.0 };

  const char* option[] = { "1Dx1D", "2D" };
  
  for (double x = -1.0, dx = 0.02; x <= +1.0; x += dx) {
    W[0] += getPhotocathodeArea() * getAngularAcceptance(x) * dx;
  }


  const double dw = 1.5;   // [nm]

  for (double w = wmin; w <= wmax; w += dw) {

    W[1] = 0.0;

    for (double x = -1.0, dx = 0.02; x <= +1.0; x += dx) {
      W[1] += KM3NET2D::getPhotocathodeArea2D(x, w) * dx;
    }

    R[0] += W[0] * dw * getQE(w) * getAbsorptionLength(w) / (w*w);
    R[1] += W[1] * dw * getAbsorptionLength(w) / (w*w);
  }
 
  for (int i = 0; i != sizeof(R)/sizeof(R[0]); ++i) {

    R[i] *= wmax*wmin / (wmax - wmin);
    R[i] *= bequerel * ng;

    R[i] *= cpu.getSurvivalProbability(npe);
    R[i] *= 0.5e-3;

    cout << "Rate " << setw(6) << left << option[i] << " [kHz]: " << R[i] << endl;
  }
}