JMultiplicityK40.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <map>
#include <cmath>
 
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TRandom3.h"
 
#include "JTools/JConstants.hh"
#include "JTools/JRange.hh"
#include "JTools/JWeight.hh"
#include "JTools/JQuantile.hh"
#include "JPhysics/KM3NeT.hh"
#include "JGeometry3D/JPosition3D.hh"
#include "JGeometry3D/JDirection3D.hh"
#include "JMath/JMathToolkit.hh"
#include "JDetector/JModule.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JDetectorAddressMapToolkit.hh"
#include "JGizmo/JManager.hh"
#include "JROOT/JRootToolkit.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
namespace {
 
  using JTOOLS::PI;
 
  typedef JTOOLS::JRange<double> JRange_t;
 
  /**
   * Generation result.
   */
  struct JResult {
    double D;           //!< generated distance [m]
    double V;           //!< generation volume [m^3]
  };
 
 
  /**
   * Generation interface.
   */
  struct JGenerator {
    virtual const JResult& next() = 0;
  };
 
 
  /**
   * Template class for generation of distance between surface of optical module and decay vertex.\n
   * The template value corresponds to the power for the generation of this distance.
   */
  template<int N> struct JEnigma;
 
 
  /**
   * Generator according D^2.
   */
  template<>
  struct JEnigma<+2> :
    public  JGenerator,
    private JResult
  {
    /*
     * Constructor.
     *
     * \param  D_m        distance range [m]
     */
    JEnigma(const JRange_t& D_m) :
      D3_min(D_m.getLowerLimit() * D_m.getLowerLimit() * D_m.getLowerLimit()),
      D3_max(D_m.getUpperLimit() * D_m.getUpperLimit() * D_m.getUpperLimit())
    {
      this->V = (4.0/3.0)*PI*(D3_max - D3_min);
    }
 
    virtual const JResult& next() 
    {
      const double D3 = gRandom->Uniform(D3_min, D3_max);
 
      this->D = cbrt(D3);
      
      return *this;
    }
 
    const double D3_min;
    const double D3_max;
  };
 
 
  /**
   * Generator according D^0.
   */
  template<>
  struct JEnigma<0> :
    public  JGenerator,
    private JResult
  {
    /*
     * Constructor.
     *
     * \param  D_m        distance range [m]
     */
    JEnigma(const JRange_t& D_m) :
      D_min(D_m.getLowerLimit()),
      D_max(D_m.getUpperLimit())
    {}
 
    virtual const JResult& next() 
    {
      this->D = gRandom->Uniform(D_min, D_max);
      this->V = 4.0*PI*D*D * (D_max - D_min);
      
      return *this;
    }
 
    const double D_min;
    const double D_max;
  };
 
 
  /**
   * Generator according D^-2.
   */
  template<>
  struct JEnigma<-2> :
    public  JGenerator,
    private JResult
  {
    /*
     * Constructor.
     *
     * \param  D_m        distance range [m]
     */
    JEnigma(const JRange_t& D_m) :
      Dinv_min(1.0/D_m.getUpperLimit()),
      Dinv_max(1.0/D_m.getLowerLimit())
    {}
 
    virtual const JResult& next() 
    {
      const double Dinv = gRandom->Uniform(Dinv_min, Dinv_max);
 
      this->D = 1.0/Dinv;
      this->V = 4.0*PI*D*D * D*D * (Dinv_max - Dinv_min);
      
      return *this;
    }
 
    const double Dinv_min;
    const double Dinv_max;
  };
 
 
  /**
   * Auxiliary base class for random number generation.
   */
  struct JCDF_t :
    public std::map<double, double>
  {
    /**
     * Get function value.
     *
     * \param  x       random number
     * \return         function value
     */
    double operator()(const double x) const
    {
      const_iterator i = this->lower_bound(x);
 
      if (i != this->begin()) {
        --i;
      }
 
      return i->second;
    }
  };
 
 
  /**
   * Auxiliary class to generate number of Cherenkov photons from K40 beta decay.\n
   * Reference: Antares internal note, ANTARES-PHYS-2006-005, J. Brunner.
   */
  class JK40BetaDecay :
    public JCDF_t
  {
  public:
    /**
     * Default constructor.
     */
    JK40BetaDecay()
    {
      (*this)[0.000000] =   0.000;
      (*this)[0.140820] =   1.000;
      (*this)[0.165515] =   2.000;
      (*this)[0.185075] =   3.000;
      (*this)[0.202024] =   4.000;
      (*this)[0.216980] =   5.000;
      (*this)[0.231162] =   6.000;
      (*this)[0.244796] =   7.000;
      (*this)[0.257734] =   8.000;
      (*this)[0.270294] =   9.000;
      (*this)[0.282401] =  10.000;
      (*this)[0.294095] =  11.000;
      (*this)[0.305618] =  12.000;
      (*this)[0.316936] =  13.000;
      (*this)[0.328064] =  14.000;
      (*this)[0.338988] =  15.000;
      (*this)[0.349740] =  16.000;
      (*this)[0.360311] =  17.000;
      (*this)[0.370834] =  18.000;
      (*this)[0.381182] =  19.000;
      (*this)[0.391350] =  20.000;
      (*this)[0.401436] =  21.000;
      (*this)[0.411309] =  22.000;
      (*this)[0.421253] =  23.000;
      (*this)[0.430885] =  24.000;
      (*this)[0.440576] =  25.000;
      (*this)[0.450050] =  26.000;
      (*this)[0.459358] =  27.000;
      (*this)[0.468670] =  28.000;
      (*this)[0.477844] =  29.000;
      (*this)[0.487045] =  30.000;
      (*this)[0.496031] =  31.000;
      (*this)[0.505032] =  32.000;
      (*this)[0.513970] =  33.000;
      (*this)[0.522793] =  34.000;
      (*this)[0.531627] =  35.000;
      (*this)[0.540303] =  36.000;
      (*this)[0.548891] =  37.000;
      (*this)[0.557355] =  38.000;
      (*this)[0.565805] =  39.000;
      (*this)[0.574089] =  40.000;
      (*this)[0.582449] =  41.000;
      (*this)[0.590614] =  42.000;
      (*this)[0.598731] =  43.000;
      (*this)[0.606772] =  44.000;
      (*this)[0.614729] =  45.000;
      (*this)[0.622647] =  46.000;
      (*this)[0.630523] =  47.000;
      (*this)[0.638299] =  48.000;
      (*this)[0.646020] =  49.000;
      (*this)[0.653625] =  50.000;
      (*this)[0.661176] =  51.000;
      (*this)[0.668671] =  52.000;
      (*this)[0.675974] =  53.000;
      (*this)[0.683343] =  54.000;
      (*this)[0.690566] =  55.000;
      (*this)[0.697753] =  56.000;
      (*this)[0.704796] =  57.000;
      (*this)[0.711778] =  58.000;
      (*this)[0.718749] =  59.000;
      (*this)[0.725594] =  60.000;
      (*this)[0.732284] =  61.000;
      (*this)[0.738976] =  62.000;
      (*this)[0.745591] =  63.000;
      (*this)[0.752171] =  64.000;
      (*this)[0.758597] =  65.000;
      (*this)[0.764951] =  66.000;
      (*this)[0.771278] =  67.000;
      (*this)[0.777509] =  68.000;
      (*this)[0.783618] =  69.000;
      (*this)[0.789677] =  70.000;
      (*this)[0.795685] =  71.000;
      (*this)[0.801588] =  72.000;
      (*this)[0.807409] =  73.000;
      (*this)[0.813077] =  74.000;
      (*this)[0.818736] =  75.000;
      (*this)[0.824251] =  76.000;
      (*this)[0.829650] =  77.000;
      (*this)[0.835090] =  78.000;
      (*this)[0.840348] =  79.000;
      (*this)[0.845611] =  80.000;
      (*this)[0.850801] =  81.000;
      (*this)[0.855803] =  82.000;
      (*this)[0.860776] =  83.000;
      (*this)[0.865648] =  84.000;
      (*this)[0.870441] =  85.000;
      (*this)[0.875074] =  86.000;
      (*this)[0.879596] =  87.000;
      (*this)[0.884102] =  88.000;
      (*this)[0.888459] =  89.000;
      (*this)[0.892717] =  90.000;
      (*this)[0.896964] =  91.000;
      (*this)[0.901105] =  92.000;
      (*this)[0.905168] =  93.000;
      (*this)[0.909097] =  94.000;
      (*this)[0.912907] =  95.000;
      (*this)[0.916656] =  96.000;
      (*this)[0.920329] =  97.000;
      (*this)[0.923924] =  98.000;
      (*this)[0.927393] =  99.000;
      (*this)[0.930810] = 100.000;
      (*this)[0.934111] = 101.000;
      (*this)[0.937303] = 102.000;
      (*this)[0.940378] = 103.000;
      (*this)[0.943379] = 104.000;
      (*this)[0.946364] = 105.000;
      (*this)[0.949168] = 106.000;
      (*this)[0.951929] = 107.000;
      (*this)[0.954558] = 108.000;
      (*this)[0.957099] = 109.000;
      (*this)[0.959577] = 110.000;
      (*this)[0.961934] = 111.000;
      (*this)[0.964206] = 112.000;
      (*this)[0.966378] = 113.000;
      (*this)[0.968517] = 114.000;
      (*this)[0.970553] = 115.000;
      (*this)[0.972528] = 116.000;
      (*this)[0.974349] = 117.000;
      (*this)[0.976112] = 118.000;
      (*this)[0.977841] = 119.000;
      (*this)[0.979426] = 120.000;
      (*this)[0.980973] = 121.000;
      (*this)[0.982446] = 122.000;
      (*this)[0.983772] = 123.000;
      (*this)[0.985085] = 124.000;
      (*this)[0.986341] = 125.000;
      (*this)[0.987526] = 126.000;
      (*this)[0.988621] = 127.000;
      (*this)[0.989622] = 128.000;
      (*this)[0.990565] = 129.000;
      (*this)[0.991450] = 130.000;
      (*this)[0.992296] = 131.000;
      (*this)[0.993054] = 132.000;
      (*this)[0.993792] = 133.000;
      (*this)[0.994464] = 134.000;
      (*this)[0.995083] = 135.000;
      (*this)[0.995633] = 136.000;
      (*this)[0.996158] = 137.000;
      (*this)[0.996639] = 138.000;
      (*this)[0.997088] = 139.000;
      (*this)[0.997470] = 140.000;
      (*this)[0.997795] = 141.000;
      (*this)[0.998098] = 142.000;
      (*this)[0.998372] = 143.000;
      (*this)[0.998612] = 144.000;
      (*this)[0.998825] = 145.000;
      (*this)[0.999005] = 146.000;
      (*this)[0.999162] = 147.000;
      (*this)[0.999314] = 148.000;
      (*this)[0.999441] = 149.000;
      (*this)[0.999543] = 150.000;
      (*this)[0.999620] = 151.000;
      (*this)[0.999695] = 152.000;
      (*this)[0.999755] = 153.000;
      (*this)[0.999812] = 154.000;
      (*this)[0.999847] = 155.000;
      (*this)[0.999885] = 156.000;
      (*this)[0.999915] = 157.000;
      (*this)[0.999939] = 158.000;
      (*this)[0.999956] = 159.000;
      (*this)[0.999966] = 160.000;
      (*this)[0.999975] = 161.000;
      (*this)[0.999984] = 162.000;
      (*this)[0.999986] = 163.000;
      (*this)[0.999992] = 164.000;
      (*this)[0.999996] = 165.000;
      (*this)[0.999997] = 166.000;
      (*this)[0.999998] = 167.000;
      (*this)[0.999999] = 168.000;
      (*this)[0.999999] = 169.000;
      (*this)[0.999999] = 170.000;
      (*this)[1.000000] = 171.000;
      /*
      (*this)[1.000000] = 172.000;
      (*this)[1.000000] = 173.000;
      (*this)[1.000000] = 174.000;
      (*this)[1.000000] = 175.000;
      (*this)[1.000000] = 176.000;
      (*this)[1.000000] = 177.000;
      (*this)[1.000000] = 178.000;
      (*this)[1.000000] = 179.000;
      (*this)[1.000000] = 180.000;
      (*this)[1.000000] = 181.000;
      (*this)[1.000000] = 182.000;
      (*this)[1.000000] = 183.000;
      (*this)[1.000000] = 184.000;
      (*this)[1.000000] = 185.000;
      (*this)[1.000000] = 186.000;
      (*this)[1.000000] = 187.000;
      (*this)[1.000000] = 188.000;
      (*this)[1.000000] = 189.000;
      (*this)[1.000000] = 190.000;
      (*this)[1.000000] = 191.000;
      (*this)[1.000000] = 192.000;
      (*this)[1.000000] = 193.000;
      (*this)[1.000000] = 194.000;
      (*this)[1.000000] = 195.000;
      (*this)[1.000000] = 196.000;
      (*this)[1.000000] = 197.000;
      (*this)[1.000000] = 198.000;
      (*this)[1.000000] = 199.000;
      (*this)[1.000000] = 200.000;
      */
      //compile();
    }
 
    /**
     * Get branching ratio.
     *
     * \return         branching ratio
     */
    static double getBranchingRatio()
    {
      return 0.8928;
    }
  };
 
 
  /**
   * Auxiliary class to generate number of Cherenkov photons from K40 electron capture.\n
   * Reference: Antares internal note, ANTARES-PHYS-2006-005, J. Brunner.
   */
  class JK40ElectronCapture :
    public JCDF_t
  {
  public:
    /**
     * Default constructor.
     */
    JK40ElectronCapture()
    {
      (*this)[0.000000] =   0.000;
      (*this)[0.000964] =   1.000;
      (*this)[0.002391] =   2.000;
      (*this)[0.004139] =   3.000;
      (*this)[0.005983] =   4.000;
      (*this)[0.008092] =   5.000;
      (*this)[0.010555] =   6.000;
      (*this)[0.013245] =   7.000;
      (*this)[0.016417] =   8.000;
      (*this)[0.019601] =   9.000;
      (*this)[0.023333] =  10.000;
      (*this)[0.027037] =  11.000;
      (*this)[0.031226] =  12.000;
      (*this)[0.035502] =  13.000;
      (*this)[0.040167] =  14.000;
      (*this)[0.045124] =  15.000;
      (*this)[0.050383] =  16.000;
      (*this)[0.055996] =  17.000;
      (*this)[0.061656] =  18.000;
      (*this)[0.067462] =  19.000;
      (*this)[0.073877] =  20.000;
      (*this)[0.080352] =  21.000;
      (*this)[0.086917] =  22.000;
      (*this)[0.093883] =  23.000;
      (*this)[0.100961] =  24.000;
      (*this)[0.108294] =  25.000;
      (*this)[0.116022] =  26.000;
      (*this)[0.123821] =  27.000;
      (*this)[0.131823] =  28.000;
      (*this)[0.139890] =  29.000;
      (*this)[0.147991] =  30.000;
      (*this)[0.156341] =  31.000;
      (*this)[0.164641] =  32.000;
      (*this)[0.173283] =  33.000;
      (*this)[0.181640] =  34.000;
      (*this)[0.189983] =  35.000;
      (*this)[0.198629] =  36.000;
      (*this)[0.207442] =  37.000;
      (*this)[0.216109] =  38.000;
      (*this)[0.225215] =  39.000;
      (*this)[0.234364] =  40.000;
      (*this)[0.243404] =  41.000;
      (*this)[0.252236] =  42.000;
      (*this)[0.261426] =  43.000;
      (*this)[0.270904] =  44.000;
      (*this)[0.280399] =  45.000;
      (*this)[0.290229] =  46.000;
      (*this)[0.300124] =  47.000;
      (*this)[0.309784] =  48.000;
      (*this)[0.319648] =  49.000;
      (*this)[0.329575] =  50.000;
      (*this)[0.339769] =  51.000;
      (*this)[0.350323] =  52.000;
      (*this)[0.360685] =  53.000;
      (*this)[0.371175] =  54.000;
      (*this)[0.381689] =  55.000;
      (*this)[0.392278] =  56.000;
      (*this)[0.402836] =  57.000;
      (*this)[0.413363] =  58.000;
      (*this)[0.424129] =  59.000;
      (*this)[0.435091] =  60.000;
      (*this)[0.445833] =  61.000;
      (*this)[0.456412] =  62.000;
      (*this)[0.466519] =  63.000;
      (*this)[0.477130] =  64.000;
      (*this)[0.487492] =  65.000;
      (*this)[0.497307] =  66.000;
      (*this)[0.507566] =  67.000;
      (*this)[0.517633] =  68.000;
      (*this)[0.527133] =  69.000;
      (*this)[0.536413] =  70.000;
      (*this)[0.545223] =  71.000;
      (*this)[0.554018] =  72.000;
      (*this)[0.562511] =  73.000;
      (*this)[0.570839] =  74.000;
      (*this)[0.579232] =  75.000;
      (*this)[0.587318] =  76.000;
      (*this)[0.595435] =  77.000;
      (*this)[0.603502] =  78.000;
      (*this)[0.611422] =  79.000;
      (*this)[0.619045] =  80.000;
      (*this)[0.626384] =  81.000;
      (*this)[0.633823] =  82.000;
      (*this)[0.641268] =  83.000;
      (*this)[0.648831] =  84.000;
      (*this)[0.656397] =  85.000;
      (*this)[0.663693] =  86.000;
      (*this)[0.671029] =  87.000;
      (*this)[0.678402] =  88.000;
      (*this)[0.685922] =  89.000;
      (*this)[0.693255] =  90.000;
      (*this)[0.700336] =  91.000;
      (*this)[0.707653] =  92.000;
      (*this)[0.714999] =  93.000;
      (*this)[0.721974] =  94.000;
      (*this)[0.728990] =  95.000;
      (*this)[0.736015] =  96.000;
      (*this)[0.742894] =  97.000;
      (*this)[0.750246] =  98.000;
      (*this)[0.757448] =  99.000;
      (*this)[0.764563] = 100.000;
      (*this)[0.771738] = 101.000;
      (*this)[0.778704] = 102.000;
      (*this)[0.785757] = 103.000;
      (*this)[0.793025] = 104.000;
      (*this)[0.800100] = 105.000;
      (*this)[0.807125] = 106.000;
      (*this)[0.814274] = 107.000;
      (*this)[0.821156] = 108.000;
      (*this)[0.828160] = 109.000;
      (*this)[0.834846] = 110.000;
      (*this)[0.841731] = 111.000;
      (*this)[0.848563] = 112.000;
      (*this)[0.855346] = 113.000;
      (*this)[0.862256] = 114.000;
      (*this)[0.868982] = 115.000;
      (*this)[0.875899] = 116.000;
      (*this)[0.882461] = 117.000;
      (*this)[0.888889] = 118.000;
      (*this)[0.895478] = 119.000;
      (*this)[0.901776] = 120.000;
      (*this)[0.908026] = 121.000;
      (*this)[0.914094] = 122.000;
      (*this)[0.920233] = 123.000;
      (*this)[0.926076] = 124.000;
      (*this)[0.931717] = 125.000;
      (*this)[0.937147] = 126.000;
      (*this)[0.942499] = 127.000;
      (*this)[0.947630] = 128.000;
      (*this)[0.952460] = 129.000;
      (*this)[0.956957] = 130.000;
      (*this)[0.961478] = 131.000;
      (*this)[0.965667] = 132.000;
      (*this)[0.969667] = 133.000;
      (*this)[0.973330] = 134.000;
      (*this)[0.976881] = 135.000;
      (*this)[0.980044] = 136.000;
      (*this)[0.982943] = 137.000;
      (*this)[0.985614] = 138.000;
      (*this)[0.987847] = 139.000;
      (*this)[0.990126] = 140.000;
      (*this)[0.991874] = 141.000;
      (*this)[0.993441] = 142.000;
      (*this)[0.994695] = 143.000;
      (*this)[0.995898] = 144.000;
      (*this)[0.996831] = 145.000;
      (*this)[0.997633] = 146.000;
      (*this)[0.998305] = 147.000;
      (*this)[0.998762] = 148.000;
      (*this)[0.999114] = 149.000;
      (*this)[0.999362] = 150.000;
      (*this)[0.999534] = 151.000;
      (*this)[0.999705] = 152.000;
      (*this)[0.999801] = 153.000;
      (*this)[0.999876] = 154.000;
      (*this)[0.999919] = 155.000;
      (*this)[0.999953] = 156.000;
      (*this)[0.999966] = 157.000;
      (*this)[0.999972] = 158.000;
      (*this)[0.999978] = 159.000;
      (*this)[0.999978] = 160.000;
      (*this)[0.999984] = 161.000;
      (*this)[0.999988] = 162.000;
      (*this)[0.999988] = 163.000;
      (*this)[0.999988] = 164.000;
      (*this)[0.999988] = 165.000;
      (*this)[0.999988] = 166.000;
      (*this)[0.999988] = 167.000;
      (*this)[0.999988] = 168.000;
      (*this)[0.999994] = 169.000;
      (*this)[0.999994] = 170.000;
      (*this)[0.999994] = 171.000;
      (*this)[0.999994] = 172.000;
      (*this)[0.999994] = 173.000;
      (*this)[0.999994] = 174.000;
      (*this)[0.999994] = 175.000;
      (*this)[0.999994] = 176.000;
      (*this)[0.999997] = 177.000;
      (*this)[0.999997] = 178.000;
      (*this)[0.999997] = 179.000;
      (*this)[0.999997] = 180.000;
      (*this)[1.000000] = 181.000;
      /*
      (*this)[1.000000] = 182.000;
      (*this)[1.000000] = 183.000;
      (*this)[1.000000] = 184.000;
      (*this)[1.000000] = 185.000;
      (*this)[1.000000] = 186.000;
      (*this)[1.000000] = 187.000;
      (*this)[1.000000] = 188.000;
      (*this)[1.000000] = 189.000;
      (*this)[1.000000] = 190.000;
      (*this)[1.000000] = 191.000;
      (*this)[1.000000] = 192.000;
      (*this)[1.000000] = 193.000;
      (*this)[1.000000] = 194.000;
      (*this)[1.000000] = 195.000;
      (*this)[1.000000] = 196.000;
      (*this)[1.000000] = 197.000;
      (*this)[1.000000] = 198.000;
      (*this)[1.000000] = 199.000;
      (*this)[1.000000] = 200.000;
      */
      //compile();
    }
 
    /**
     * Get branching ratio.
     *
     * \return         branching ratio
     */
    static double getBranchingRatio()
    { 
      return 0.1072;
    }
  };
 
  /**
   * Function objects for K40 beta.
   */
  static const JK40BetaDecay        k40_beta_decay;
  static const JK40ElectronCapture  k40_electron_capture;
 
  /**
   * Parametrised angular acceptence of KM3NeT PMT.\n
   * Reference: Antares internal note, ANTARES-PHYS-2006-005, J. Brunner.
   */
  double a;
 
  /**
   * Parametrised angular acceptence of KM3NeT PMT.
   *
   * \param  ct     cosine of angle of incidence
   * \return        probability
   */
  inline double get_angular_acceptance(const double ct)
  {
    const double w = 0.25*a*(1.0 + ct)*(1.0 + ct) - ct;
 
    if (w > 0.0)
      return 1.4 * w;
    else
      return 0.0;
  }
}
 
 
/**
 * \file
 *
 * Example program to calculate multiples rate.
 * \author mdejong
 */
int main(int argc, char* argv[])
{
  using namespace std;
  using namespace JPP;
 
  typedef long long int  counter_type;
 
  string        outputFile;
  double        bequerel;
  JRange_t      D_m;
  counter_type  numberOfEvents;
  double        QE;
  UInt_t        seed;
  int           generator;
  double        focus;
  int           Aefftype;
  bool          exclusive;
  int           debug;
 
  try {
 
    JParser<> zap("Example program to calculate multiples rate.");
 
    zap['o'] = make_field(outputFile)      =  "k40.root";
    zap['b'] = make_field(bequerel)        =  13.75e3;     // [m^-3 s^-1]
    zap['n'] = make_field(numberOfEvents)  =  1e7;
    zap['D'] = make_field(D_m)             =  JRange_t(0.216, 10);
    zap['Q'] = make_field(QE)              =  1.0;
    zap['S'] = make_field(seed)            =  0;
    zap['G'] = make_field(generator)       =  0, +2, -2;
    zap['F'] = make_field(focus)           =  1.0;
    zap['A'] = make_field(Aefftype)        =  1,  2,  3;
    zap['U'] = make_field(exclusive);
    zap['a'] = make_field(a)               =  0.0;
    zap['d'] = make_field(debug)           =  2;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
  gRandom->SetSeed(seed);
 
  using namespace KM3NET;
  using namespace JPP;
 
  const int     id     = 1;
  const JModule module = getModule(getDetectorAddressMap(14).get(id), id);
 
  DEBUG(module << endl);
 
  bequerel /= focus;           // correct decay rate for focussing of photons
 
  const double  R_m = 17.0 * 2.54 * 0.5e-2;                // radius 17''   optical module [m]
  const double  A   = PI * R_m * R_m;                      // cross section optical module [m^2]
 
  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 double WAVELENGTH_EXPANSION =  (wmax-wmin) / (wmin*wmax) * (300.0*600.0)/(600.0-300.0);
 
  JGenerator* enigma = NULL;         // distance generation
 
  switch (generator) {
 
  case +2: enigma = new JEnigma<+2>(D_m); break;
  case  0: enigma = new JEnigma< 0>(D_m); break;
  case -2: enigma = new JEnigma<-2>(D_m); break;
 
  default: FATAL("No generator type " << generator << endl);
  }
 
  const double vmin = 1.0 / wmax;    // wavelength generation
  const double vmax = 1.0 / wmin;    // as lambda^-2
 
  double QEmax = 0.0;
 
  for (double w = wmin; w <= wmax; w += 1.0) {
    if (getQE(w) > QEmax) {
      QEmax = getQE(w);
    }
  }
 
  NOTICE("Maximal QE                  " << FIXED(5,3) << QEmax                << endl);
  NOTICE("Wavelength expansion        " << FIXED(5,3) << WAVELENGTH_EXPANSION << endl);
  NOTICE("Number of photons per decay " << FIXED(5,2) << ng                   << endl);
 
  typedef JManager<int, TH1D>  JManager_t;
 
  JManager_t H1(new TH1D("M[%]", NULL, 100, D_m.getLowerLimit(), D_m.getUpperLimit()));
 
  H1->Sumw2();
 
  TH1D pmt("pmt", NULL, 1000, -1.0, +1.0);
 
  for (Int_t i = 1; i != pmt.GetNbinsX(); ++i) {
 
    const double dot = pmt.GetBinCenter(i);
    
    double y = 0.0;
 
    switch (Aefftype) {
      
    case 1:
      y = getAngularAcceptance(dot);
      break;
      
    case 3:
      y = get_angular_acceptance(dot);
      break;
    }
 
    pmt.SetBinContent(i, y);
  }
 
 
  map<int, JWeight>   h1;
  map<int, JQuantile> P2;
 
 
  vector<double>  pi(module.size());
  vector<int>     buffer;
  
  for (counter_type event_count = 0; event_count != numberOfEvents; ++event_count) {
 
    if (event_count%10000 == 0) {
      STATUS("event: " << setw(10) << event_count << '\r'); DEBUG(endl);
    }
 
    const JResult& result = enigma->next();
 
    const double D = result.D;
    const double V = result.V;
 
    // check first photons on cross section of optical module at its face 
    // it is assumed that the light from the K40 decay is purely isotropic
 
    double W = A / (4*PI*(D-R_m)*(D-R_m));
 
    if (W > 0.5) {
      W = 0.5;
    }
 
    //
    double x = gRandom->Rndm();      // decay mode
    double y = 0.0;                  // number of photons
 
    if      ((x -= k40_beta_decay      .getBranchingRatio()) <= 0.0)
      y = k40_beta_decay      (gRandom->Rndm());
    else if ((x -= k40_electron_capture.getBranchingRatio()) <= 0.0)
      y = k40_electron_capture(gRandom->Rndm());
    
    const int N = gRandom->Poisson(y  * WAVELENGTH_EXPANSION * QE * W * QEmax * focus);
    //
    /*
    const int N = gRandom->Poisson(ng * WAVELENGTH_EXPANSION * QE * W * QEmax * focus);
    */
    if (N >= 2) {
 
      // generate vertex
 
      const double ct  = gRandom->Uniform(-1.0, +1.0);
      const double phi = gRandom->Uniform(-PI,  +PI);
      
      const double st  = sqrt((1.0 - ct) * (1.0 + ct));    
 
      const JPosition3D vertex(D * st * cos(phi),
                               D * st * sin(phi),
                               D * ct);
      
      buffer.clear();
 
      for (int i = 0; i != N; ++i) {
 
        // generate wavelength of photon
 
        const double v     = gRandom->Uniform(vmin, vmax);
        const double w     = 1.0 / v;
        
        const double l_abs = getAbsorptionLength(w);
 
        double P = 0.0;
 
        for (size_t pmt = 0; pmt != module.size(); ++pmt) {
 
          JPosition3D pos(module[pmt].getPosition());
 
          pos -= vertex;
        
          const double d   = pos.getLength();
 
          pos /= d;                // direction of photon
 
          if (d < D - R_m) {
            ERROR("Distance " << d << " < " << D << endl);
          }
 
          const double dot = getDot(pos, module[pmt].getDirection());
 
          const double U   = 0.5 * (1.0 - d / sqrt(d*d + getPhotocathodeArea() / PI));
 
          double p = 0.0;
 
          switch (Aefftype) {
 
          case 1:
            p = getAngularAcceptance(dot)   * getQE(w);
            break;
 
          case 2:
            p = getPhotocathodeArea2D(dot, w) / getPhotocathodeArea();
            break;
 
          case 3:
            p = get_angular_acceptance(dot) * getQE(w);
            break;
          }
 
          P += pi[pmt] = U * p * exp(-d/l_abs);
        }
 
        if (P > W) {
          ERROR("Probability " << P << " > " << W << endl);
        }
         
        if (W * QEmax * gRandom->Rndm() < P) {
 
          int    pmt = 0;
          double y   = gRandom->Uniform(P);
          
          for (vector<double>::const_iterator i = pi.begin(); i != pi.end() && (y -= *i) > 0.0; ++i, ++pmt) {}
 
          buffer.push_back(pmt);
        }
      }
 
      if (!buffer.empty()) {
 
        int M = buffer.size();
 
        if (exclusive) {
 
          sort(buffer.begin(), buffer.end());
 
          M = distance(buffer.begin(), unique(buffer.begin(), buffer.end()));
        }
 
        h1[M].put(V);
        H1[M]->Fill(D, V);
 
        for (int i = 2; i <= M; ++i) {
          P2[i].put((double) (buffer.size() - M) / (double) M, V);
        }
 
      }
    }
  }
  STATUS(endl);
 
  for (JManager_t::iterator i = H1.begin(); i != H1.end(); ++i) {
    i->second->Scale(bequerel / (double) numberOfEvents);
  }
 
  for (size_t M = 2; M != 7; ++M) {
    cout << "Rate[" << M << "] = "
         << FIXED(7,3) << bequerel * h1[M].getTotal() / (double) numberOfEvents
         << " +/- "
         << FIXED(7,3) << bequerel * h1[M].getError() / (double) numberOfEvents
         << " Hz" << endl;
  }
 
  for (size_t M = 2; M != 7; ++M) {
    if (P2[M].getCount() != 0) {
      cout << "P2[" << M << "] = " << P2[M].getMean() << endl;
    }
  }
 
  TFile out(outputFile.c_str(), "recreate");
 
  out << H1 << pmt;
 
  out.Write();
  out.Close();
}