JCharge.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/JPMTParametersMap.hh"
#include "JDetector/JPMTAnalogueSignalProcessor.hh"

#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


/**
 * \file
 *
 * Example program to histogram charge distributions.
 * \author mdejong & bjjung 
 */
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 charge distributions.");

    zap['o'] = make_field(outputFile)       = "charge.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.01;
    zap['d'] = make_field(debug)            = 0;

    zap(argc, argv);

  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }

  const JPMTParameters parameters = parametersMap.getPMTParameters(pmt);

  const JPMTAnalogueSignalProcessor cpu      (parameters);

  ASSERT(numberOfHits > 0);


  TFile out(outputFile.c_str(), "recreate");

  const double dQ   = (cpu.threshold - cpu.thresholdBand) / 5.5;
  const double Qmin =   0.0 - 0.5 * dQ;
  const double Qmax =  10.0 + 0.5 * dQ;
  const int    nq   = (int) ((Qmax - Qmin) / dQ);

  TH1D h0("h0", NULL, nq, Qmin, Qmax);
  TH1D h1("h1", NULL, nq, Qmin, Qmax);

  h1.Sumw2();

  for (int i=1; i <= h0.GetNbinsX(); ++i) {
    
    const double npe = h0.GetBinCenter(i);
    const double p   = cpu.getChargeProbability(npe,NPE);

    h0.SetBinContent(i,p);
  }


  if (numberOfHits > 0) {

    int number_of_hits = 0;
    
    for (int i=0; i != numberOfHits; ++i) {

      const double npe = cpu.getRandomCharge(NPE);

      try {

        if (cpu.applyThreshold(npe) > cpu.BELOW_THRESHOLDBAND) {

          ++number_of_hits;
          
          h1.Fill(npe);
        }

      } catch (const JValueOutOfRange& exception) {

        DEBUG(exception.what());
        continue;
      }

    }

    h1.Scale(1.0 / (double) number_of_hits / dQ);
  }
  
  out.Write();
  out.Close();


  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);

    DEBUG("[" << setw(3) << i << "]"  << ' ' 
          << FIXED(5,3) << x  << ' '
          << FIXED(6,4) << y0 << ' '
          << FIXED(6,4) << y1 << endl);

    ASSERT(fabs(y0 - y1) < precision);
  }
}