JMultiplicityK40.cc File Reference

Example program to calculate multiples rate. More...

#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/JRange.hh"
#include "JTools/JWeight.hh"
#include "JTools/JQuantile.hh"
#include "JPhysics/KM3NeT.hh"
#include "JPhysics/KM3NeT2D.hh"
#include "JPhysics/JConstants.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 "JROOT/JManager.hh"
#include "JROOT/JRootToolkit.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions
int	main (int argc, char *argv[])

Detailed Description

Example program to calculate multiples rate.

Author

mdejong

Definition in file JMultiplicityK40.cc.

Function Documentation

int main	(	int	argc,
		char *	argv[]
	)

Definition at line 682 of file JMultiplicityK40.cc.

{
  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 NAMESPACE;


  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 = KM3NET2D::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(8,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();
}