JCharge.cc

Go to the documentation of this file.

 
#include <string>
#include <iostream>
#include <iomanip>
#include <map>
 
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
 
#include "JDetector/JPMTParametersMap.hh"
#include "JDetector/JPMTAnalogueSignalProcessor.hh"
 
#include "JROOT/JRootToolkit.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);
 
  map<int, TH1D*> H0;
  map<int, TH1D*> H1;
 
  {
    const double dx   = JPMTSignalProcessorInterface::getQmin();
    const double xmin = 0.0  -  0.5*dx;
    const int    nx   = (int) ((10.0 - xmin) / dx);
    
    H0[1] = new TH1D("1.0", NULL, nx, xmin, xmin + nx * dx);
    H1[1] = new TH1D("1.1", NULL, nx, xmin, xmin + nx * dx);
  }
  {
    const double dx   = 0.02;
    const double xmin = parameters.threshold - parameters.thresholdBand;
    const int    nx   = (int) ((10.0 - xmin) / dx);
    
    H0[2] = new TH1D("2.0", NULL, nx, xmin, xmin + nx * dx);
    H1[2] = new TH1D("2.1", NULL, nx, xmin, xmin + nx * dx);
  }
 
  H1[1]->Sumw2();
  H1[2]->Sumw2();
 
  int test = 1;
 
  for (JPMTSignalProcessorInterface* cpu : {
      (JPMTSignalProcessorInterface*) new JPMTSignalProcessorInterface(),
      (JPMTSignalProcessorInterface*) new JPMTAnalogueSignalProcessor(parameters)
        }) {
 
    for (int i=1; i <= H0[test]->GetNbinsX(); ++i) {
    
      const double npe = H0[test]->GetBinCenter(i);
      const double p   = cpu->getChargeProbability(npe,NPE);
 
      H0[test]->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)) {
 
            ++number_of_hits;
          
            H1[test]->Fill(npe);
          }
 
        } catch (const JValueOutOfRange& exception) {
 
          DEBUG(exception.what());
          continue;
        }
      }
 
      const double dx = (H1[test]->GetXaxis()->GetXmax() - H1[test]->GetXaxis()->GetXmin()) / H1[test]->GetXaxis()->GetNbins();
      
      H1[test]->Scale(1.0 / (double) number_of_hits / dx);
 
      delete cpu;
 
      ++test;
    }
  }
  
  TFile out(outputFile.c_str(), "recreate");
 
  for (auto& i : H0) { out << *i.second; }
  for (auto& i : H1) { out << *i.second; }
  
  out.Write();
  out.Close();
 
  // test
 
  ASSERT(numberOfHits > 0);
 
  for (int i = 1; i != test; ++i) {
 
    DEBUG("Test " << i << endl);
    
    for (int ix = 1; ix <= H0[i]->GetNbinsX(); ++ix) {
 
      const Double_t x  = H0[i]->GetBinCenter (ix);
      const Double_t y0 = H0[i]->GetBinContent(ix);
      const Double_t y1 = H1[i]->GetBinContent(ix);
 
      DEBUG("[" << setw(3) << ix << "]"  << ' ' 
            << FIXED(5,3) << x  << ' '
            << FIXED(6,4) << y0 << ' '
            << FIXED(6,4) << y1 << endl);
 
      ASSERT(fabs(y0 - y1) < precision);
    }
  }
 
  return 0;
}