JPMTSignalProcessor.cc

Go to the documentation of this file.

 
#include <string>
#include <iostream>
#include <iomanip>
 
#include "JDetector/JPMTParametersMap.hh"
#include "JDetector/JPMTIdentifier.hh"
#include "JDetector/JPMTAnalogueSignalProcessor.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * \file
 *
 * Example program to test PMT signal processor probability.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  int               numberOfHits;
  JPMTParametersMap parameters;
  JPMTIdentifier    pmt;
  int               NPE;
  double            precision;
  int               debug;
 
  try {
 
    JParser<> zap("Example program to test PMT signal processor probability.");
 
    zap['n'] = make_field(numberOfHits)     = 10000;
    zap['P'] = make_field(parameters)       = JPARSER::initialised();
    zap['p'] = make_field(pmt)              = JPMTIdentifier(-1,0);
    zap['N'] = make_field(NPE)              = 1;
    zap['e'] = make_field(precision)        = 0.02;
    zap['d'] = make_field(debug)            = 3;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  const double QEmax = 1.0;
 
  JPMTParameters buffer = parameters.getPMTParameters(pmt);
 
  DEBUG("PMT parameters:" << endl);
  DEBUG(buffer.getProperties(JEquationParameters("=", "\n", "", "")) << endl);
 
  if (buffer.QE > QEmax) {
    
    WARNING("PMT QE set to " << QEmax << endl);
 
    buffer.QE = QEmax;
  }
 
  ASSERT(numberOfHits > 0);
 
  for (JPMTSignalProcessorInterface* cpu : {
      (JPMTSignalProcessorInterface*) new JPMTSignalProcessorInterface(),
      (JPMTSignalProcessorInterface*) new JPMTAnalogueSignalProcessor(buffer) }) {
 
    int number_of_hits = 0;
 
    for (int i = 0; i != numberOfHits; ++i) {
 
      int npe = 0;
    
      for (int __i = 0; __i != NPE; ++__i) {
        if (cpu->applyQE()) {
          ++npe;
        }
      }
 
      double x = cpu->getRandomCharge(npe);
 
      try {
 
        if (cpu->applyThreshold(x)) {
          ++number_of_hits;
        }
 
      } catch (const JValueOutOfRange& exception) {
 
        DEBUG(exception.what());
        continue;
      }
 
    }
 
    const double P = (double) number_of_hits / (double) numberOfHits; 
 
    DEBUG("Number of generated hits " << setw(8) << right << numberOfHits   << endl);
    DEBUG("Number of accepted  hits " << setw(8) << right << number_of_hits << endl);
    DEBUG("Probability (real) " << FIXED(8,5) << P                                  << endl);
    DEBUG("Probability (calc) " << FIXED(8,5) << cpu->getSurvivalProbability(NPE)   << endl);
 
    ASSERT(fabs(P - cpu->getSurvivalProbability(NPE)) < precision * cpu->getSurvivalProbability(NPE));
 
    delete cpu;
  }
 
  return 0;
}