JRateK40.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
 
#include "JPhysics/KM3NeT.hh"
#include "JPhysics/KM3NeT2D.hh"
#include "JPhysics/Antares.hh"
#include "JPhysics/JPDF.hh" //for module radius
#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);
 
  //For Geant4 inputs see https://git.km3net.de/vkulikovskiy/g4ly
  const double wmin = 280.0;   // minimal wavelength [nm]
  const double wmax = 700.0;   // maximal wavelength [nm]
  double ng = 47.11;           // Geant4 simulations of K40 decays in (wmin,wmax) window 
  const int npe  =   1;        // number of photo-electrons for each decay
  double cpow = 2.156;         // 1/lambda^cpow shape of the Cherenkov photons distribution is a fit from Geant4
 
  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]
 
  double norm[] = { 0.0, 0.0 };
 
  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) / pow(w,cpow);
    R[1] += W[1] * dw * getAbsorptionLength(w) / pow(w,cpow);
    norm[0] += W[0] * dw * getQE(w) / pow(w,cpow);
    norm[1] += W[1] * dw / pow(w,cpow);
  }
 
  double avabs[] = { 0.0, 0.0 };
 
  cout << "survival probability: " << cpu.getSurvivalProbability(npe) << endl;
 
  for (int i = 0; i != sizeof(R)/sizeof(R[0]); ++i) {
    avabs[i] = R[i]/norm[i];
 
    R[i] *= cpow-1.0;
    R[i] /= pow(wmin, 1.0-cpow) - pow(wmax,1.0-cpow);
 
    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;
    double corr = R[i]*(1.0-exp(-MODULE_RADIUS_M/avabs[i])); //rate due to coincidences within the DOM
    cout << "Rcor "<<  setw(6) << left << option[i] << " [kHz]: " << R[i]-corr << " <abs> [m]: "<< avabs[i] << endl;
  }
}