JToT.cc

Go to the documentation of this file.

 
#include <string>
#include <iostream>
#include <iomanip>
 
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
 
#include "JDetector/JPMTParameters.hh"
#include "JDetector/JPMTAnalogueSignalProcessor.hh"
#include "JDetector/JPMTDefaultSimulator.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
#include "JSupport/JMeta.hh"
 
 
/**
 * \file
 *
 * Example program to histogram time over threshold distributions.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  string            outputFile;
  JPMTParametersMap parametersMap;
  JPMTIdentifier    pmt;
  int               NPE;
  int               numberOfHits;
  double            precision;
  int               debug;
 
  try {
 
    JParser<> zap("Example program to histogram time over threshold distributions.");
 
    zap['o'] = make_field(outputFile)       = "tot.root";
    zap['P'] = make_field(parametersMap)    = JPARSER::initialised();
    zap['p'] = make_field(pmt)              = JPMTIdentifier(1,0);
    zap['N'] = make_field(NPE)              = 1;
    zap['n'] = make_field(numberOfHits)     = 1000000;
    zap['e'] = make_field(precision)        = 0.005;
    zap['d'] = make_field(debug)            = 0;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
  const JPMTParameters parameters = parametersMap.getPMTParameters(pmt);
 
  const JCalibration                calibration;
  const JStatus                     status;
  const JPMTDefaultSimulator        simulator(parameters, pmt);
  const JPMTAnalogueSignalProcessor cpu      (parameters);
 
  JPMTData<JPMTSignal> input;
  JPMTData<JPMTPulse>  output;
 
 
  TFile out(outputFile.c_str(), "recreate");
 
  const double dx   =   1.0;
  const double xmin =   0.0 - 0.5 * dx;
  const double xmax = round(parameters.saturation) - 0.5*dx;
  const int    nx   = (int) ((xmax - xmin) / dx);
 
  TH1D h0("h0", NULL, nx, xmin, xmax);
  TH1D h1("h1", NULL, nx, xmin, xmax);
  TH1D h2("h2", NULL, nx, xmin, xmax);
  TH1D h3("h3", NULL, 500, 0.0, 100.0);
 
  h1.Sumw2();
  h2.Sumw2();
 
  DEBUG("    ToT      Probability " << endl);
 
  for (int i = 1; i <= h0.GetNbinsX(); ++i) {
 
    const double x        = h0.GetBinCenter(i);
    const double y        = cpu.getTimeOverThresholdProbability(x,NPE);
 
    DEBUG("  "
          << FIXED(7,3) << x        << "   "
          << FIXED(7,5) << y << endl);
 
    h0.SetBinContent(i, y);
  }
 
 
  if (numberOfHits > 0) {
 
    int number_of_hits = 0;
 
    for (int i = 0; i != numberOfHits; ++i) {
 
      const double npe = cpu.getRandomCharge(NPE);
 
      double tot_ns;
 
      try {
        tot_ns = cpu.getTimeOverThreshold(npe);
      } catch (const JValueOutOfRange& exception) {
        continue;
      }
      
      if (tot_ns > 0.0) {
        
        ++number_of_hits;
 
        h1.Fill(tot_ns);
 
      }
    }
 
    h1.Scale(1.0 / (double) number_of_hits / dx);
  }
 
 
  if (numberOfHits > 0) {
 
    int number_of_hits = 0;
 
    const double t_ns = 25.0;
    
    for (int i = 0; i != numberOfHits; ++i) {
      
      input .clear();
      output.clear();
 
      input.push_back(JPMTSignal(t_ns, NPE));
 
      simulator.processHits(pmt, calibration, status, input, output);
 
      for (JPMTData<JPMTPulse>::const_iterator hit = output.begin(); hit != output.end(); ++hit) {
 
        ++number_of_hits;
 
        h2.Fill(hit->tot_ns);
        h3.Fill(hit->t_ns);
      }
    }
 
    h2.Scale(1.0 / (double) number_of_hits / dx);
  }
 
  out.Write();
  out.Close();
 
  ASSERT(numberOfHits > 0);
 
  for (int i = 1; i <= h0.GetNbinsX(); ++i) {
 
    const Double_t x  = h0.GetBinCenter (i);
    const Double_t y0 = h0.GetBinContent(i);
    const Double_t y1 = h1.GetBinContent(i);
    const Double_t y2 = h2.GetBinContent(i);
 
    DEBUG("[" << setw(3) << i << "]"  << ' ' 
          << FIXED(5,1) << x  << ' '
          << FIXED(6,4) << y0 << ' '
          << FIXED(6,4) << y1 << ' '
          << FIXED(6,4) << y2 << endl);
 
    ASSERT(fabs(y0 - y1) < precision);
    ASSERT(fabs(y0 - y2) < precision);
  }
}