JRateK40.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <map>
 
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TFormula.h"
 
#include "JROOT/JRootToolkit.hh"
 
#include "JTools/JFunction1D_t.hh"
 
#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/JContainer.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * \file
 *
 * Example program to calculate singles rate.\n
 * The calculation is based on Antares internal note ANTARES-PHYS-2006-005 by Juergen Brunner.\n
 * According Antares internal note ANTARES-PHYS-2012-013, the absorption length is used (and not the attenuation length).
 * \author mdejong, vkulikovskiy
 */
int main(int argc, char* argv[])
{
  using namespace std;
  using namespace JPP;
 
  typedef JContainer< map<double, double> > container_type;
 
  string         outputFile;
  double         bequerel;
  JPMTParameters parameters;
  container_type fd;
  TFormula       f1;
  int            debug;
 
  try {
 
    JProperties properties = parameters.getProperties();
 
    JParser<> zap("Example program to calculate singles rate.");
 
    zap['o'] = make_field(outputFile)                        = "k40.root";
    zap['b'] = make_field(bequerel,    "radioactivity")      = 13.75e3;     // [m^-3 s^-1]
    zap['P'] = make_field(properties,  "PMT parameters")     = JPARSER::initialised();
    zap['f'] = make_field(fd,          "filter data")        = JPARSER::initialised();
    zap['F'] = make_field(f1,          "filter function")    = JPARSER::initialised();
    zap['d'] = make_field(debug)                    = 3;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  using namespace NAMESPACE;
 
  JPolint1Function1D_t g1;
 
  for (const auto& i : fd) {
    g1[i.first] = i.second;
  }
 
  g1.compile();
  g1.setExceptionHandler(new typename JPolint1Function1D_t::JDefaultResult(JMATH::zero));  
  
  const 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]
  const double ng   = 47.11;    // Geant4 simulations of K40 decays in (wmin,wmax) window 
  const int    npe  =  1;       // number of photo-electrons for each decay
  const double cpow = 2.156;    // 1/lambda^cpow shape of the Cherenkov photons distribution is a fit from Geant4
 
  TFile out(outputFile.c_str(), "recreate");
 
  const double dx = 1.5;   // [nm]
  const int    nx = (int) ((wmax - wmin) / dx);
  
  TH1D h0("h0", NULL, nx, wmin, wmax);
  TH1D h1("h1", NULL, nx, wmin, wmax);
 
 
  double R[] = { 0.0, 0.0 };
  double W[] = { 0.0, 0.0 };
 
  const string option[] = { "1Dx1D", "2D" };
  
  for (double x = -1.0, dx = 0.02; x <= +1.0; x += dx) {
    W[0] += getPhotocathodeArea() * getAngularAcceptance(x) * dx;
  }
 
  double Y = 1.0;
 
  Y *= cpow - 1.0;
  Y /= pow(wmin, 1.0 - cpow)  -  pow(wmax, 1.0 - cpow);
 
  Y *= bequerel * ng;
 
  Y *= cpu.getSurvivalProbability(npe);
  Y *= 0.5e-3;
 
 
  double A[] = { 0.0, 0.0 };
 
  for (int ix = 1; ix <= h0.GetXaxis()->GetNbins(); ++ix) {
 
    const double w  = h0.GetXaxis()->GetBinCenter(ix);
    const double dw = h0.GetXaxis()->GetBinWidth (ix);
 
    W[1] = 0.0;
 
    for (double x = -1.0, dx = 0.02; x <= +1.0; x += dx) {
      W[1] += KM3NET2D::getPhotocathodeArea2D(x, w) * dx;
    }
 
    double U = Y / pow(w,cpow);
 
    if (!g1.empty()) {
      U *= g1(w);
    }
    if (f1.IsValid()) {
      U *= f1.Eval(w);
    }
 
    R[0] += W[0] * U * dw * getQE(w) * getAbsorptionLength(w);
    R[1] += W[1] * U * dw * getAbsorptionLength(w);
    A[0] += W[0] * U * dw * getQE(w);
    A[1] += W[1] * U * dw;
 
    h0.SetBinContent(ix, W[0] * U * getQE(w) * getAbsorptionLength(w));
    h1.SetBinContent(ix, W[1] * U * getAbsorptionLength(w));
  }
 
  double labs[] = { 0.0, 0.0 };
 
  cout << "PMT survival probability: " << cpu.getSurvivalProbability(npe) << endl;
 
  for (int i = 0; i != sizeof(R)/sizeof(R[0]); ++i) {
 
    labs[i] = R[i]/A[i];
 
    const double corr = R[i]*(1.0-exp(-MODULE_RADIUS_M/labs[i])); //rate due to coincidences within the DOM
 
    cout << setw(6) << left << option[i] + ":" << right
         << " rate "      << FIXED(7,3) << R[i]        << " [kHz]"
         << " corrected " << FIXED(7,3) << R[i] - corr << " [kHz]"
         << " <l_abs> "   << FIXED(7,3) << labs[i]     << " [m]" << endl;
  }
 
  out.Write();
  out.Close();
}