JRateK40.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
 
#include "JPhysics/KM3NeT.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 KM3NET;
  //using namespace ANTARES;
 
  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] += 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;
  }
}