JK40Rates.cc

Go to the documentation of this file.

 
#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <cmath>
 
#include "TRandom3.h"
 
#include "JPhysics/JK40Rates.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * \file
 *
 * Example program to check background calculation.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  JK40Rates      rates_Hz;
  double         QE;
  double         livetime_s;
  double         precision;
  ULong_t        seed;
  int            debug;
 
  try {
 
    JParser<> zap("Example program to check background calculation.");
 
    zap['B'] = make_field(rates_Hz);
    zap['Q'] = make_field(QE)                  = 1.0;
    zap['T'] = make_field(livetime_s)          = 1.0e5;
    zap['e'] = make_field(precision)           = 0.02;
    zap['S'] = make_field(seed)                = 0;
    zap['d'] = make_field(debug)               = 1;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
  gRandom->SetSeed(seed);
 
  NOTICE("K40 rates (original)  " << rates_Hz << endl);
 
  JK40Rates corrected_rates_Hz = rates_Hz;
  
  corrected_rates_Hz.correct(QE);
  
  NOTICE("K40 rates (corrected) " << corrected_rates_Hz << endl);
 
  if (livetime_s > 0.0) {
  
    vector<double> buffer;
 
    buffer.resize(corrected_rates_Hz.getUpperL1Multiplicity() + 1);
    
    const double W = 1.0 / livetime_s;
      
    for (size_t M = corrected_rates_Hz.getLowerL1Multiplicity(); M <= corrected_rates_Hz.getUpperL1Multiplicity(); ++M) {
      
      const int numberOfEvents = (int) (livetime_s * corrected_rates_Hz.getMultiplesRate(M)); 
      
      for (int number_of_events = 0; number_of_events != numberOfEvents; ++number_of_events) {
        
        int n = 0;
        
        for (size_t i = 0; i != M; ++i) {
          if (gRandom->Rndm() <= QE) {
            ++n;
          }
        }
        
        buffer[n] += W;
      }
    }
 
    NOTICE("simulated" << endl);
    
    for (size_t M = corrected_rates_Hz.getLowerL1Multiplicity(); M <= corrected_rates_Hz.getUpperL1Multiplicity(); ++M) {
      NOTICE(setw(2) << M << ' ' << FIXED(8,4) << buffer[M] << endl);
    }
 
    // test
    
    for (size_t M = rates_Hz.getLowerL1Multiplicity(); M <= rates_Hz.getUpperL1Multiplicity(); ++M) {
      if (fabs(buffer[M] - rates_Hz.getMultiplesRate(M)) > precision * rates_Hz.getMultiplesRate(M)) {
        FATAL("R[" << setw(2) << M << "] [Hz] " << FIXED(8,4) << buffer[M] << " != " << rates_Hz.getMultiplesRate(M) << endl);
      }
    }
  }
}