JCalibrateToT.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <limits>

#include "TROOT.h"
#include "TFile.h"
#include "TH2D.h"

#include "km3net-dataformat/online/JDAQ.hh"
#include "km3net-dataformat/online/JDAQHeader.hh"

#include "JDAQ/JDAQTimesliceIO.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JTrigger/JHitR0.hh"
#include "JTrigger/JMatchL0.hh"
#include "JTrigger/JSuperFrame1D.hh"
#include "JTrigger/JSuperFrame2D.hh"
#include "JTrigger/JPreprocessor.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"
#include "JROOT/JManager.hh"
#include "JMath/JMathToolkit.hh"

#include "JROOT/JROOTClassSelector.hh"
#include "JROOT/JRootToolkit.hh"
#include "JROOT/JRootFileWriter.hh"
#include "JLang/JObjectMultiplexer.hh"
#include "JTrigger/JTriggerParameters.hh"
#include "JSupport/JTriggerParametersSupportkit.hh"
#include "JCalibrate/JCalibrateToT.hh"

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

/**
 * \file 
 *
 * Monitoring of PMT time-over-threshold distributions.
 * \author mkarel
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  JMultipleFileScanner<JDAQTimesliceTypes_t> inputFile;
  JLimit_t&          numberOfEvents = inputFile.getLimit();
  string             outputFile;
  string             detectorFile;
  JRange<double>     T_ns;
  double             ctMax;
  JRange<int>        multiplicity;
  double             deadTime_us;
  JROOTClassSelector selector;
  JPreprocessor      option;
  int                debug;

  try { 

    JParser<> zap("Monitoring of PMT time-over-threshold distributions.");

    zap['f'] = make_field(inputFile,      "input file.");
    zap['o'] = make_field(outputFile,     "output file.")                                 = "calibrate_tot.root";
    zap['n'] = make_field(numberOfEvents)                                                 = JLimit::max();
    zap['a'] = make_field(detectorFile,   "detector file.");
    zap['T'] = make_field(T_ns,           "time window [ns].")                            = JRange<double>(10.0, 25.0);
    zap['c'] = make_field(ctMax,          "maximal cosine space angle between PMT axes.") =   0.0;
    zap['M'] = make_field(multiplicity,   "multiplicity range of hits on DOM.")           = JRange<int>(1, 2);
    zap['D'] = make_field(deadTime_us,    "L0 dead time (us)")                            = 10;
    zap['C'] = make_field(selector,       "timeslice selector, e.g. JDAQTimesliceL1.")    = getROOTClassSelection<JDAQTimesliceTypes_t>();
    zap['O'] = make_field(option,         "hit pre-processing option.")                   = JPreprocessor::getOptions();
    zap['d'] = make_field(debug,          "debug flag.")                                  = 1;

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }

  cout.tie(&cerr);
  
  //-----------------------------------------------------------
  // check the input parameters
  //-----------------------------------------------------------
  
  if (!T_ns.is_valid()) {
    FATAL("Invalid time window [ns] " << T_ns << endl);
  }

  if (selector == JDAQTimeslice  ::Class_Name() ||
      selector == JDAQTimesliceL1::Class_Name()) {

    JTriggerParameters parameters;

    try {
      parameters = getTriggerParameters(inputFile);
    }
    catch(const JException& error) {
      FATAL("No trigger parameters from input." << endl);
    }

    if ((selector == JDAQTimeslice  ::Class_Name() && parameters.writeL1.prescale > 0) ||
        (selector == JDAQTimesliceL1::Class_Name())) {

      if (parameters.TMaxLocal_ns < T_ns.getUpperLimit()) {
        FATAL("Option -T <T_ns> = " << T_ns.getUpperLimit() << " is larger than in the trigger " << parameters.TMaxLocal_ns << endl);
      }
    }
  }
  
  if (!multiplicity.is_valid())          { FATAL("Invalid multiplicity " << multiplicity << endl); }
  if ( multiplicity.getLowerLimit() < 1) { FATAL("Invalid multiplicity " << multiplicity << endl); }
  
  //-----------------------------------------------------------
  // load detector file
  //-----------------------------------------------------------
  
  JDetector detector;

  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }

  if (detector.empty()) {
    FATAL("Empty detector." << endl);
  }

  const JModuleRouter router(detector);

  //-----------------------------------------------------------
  // initialise histograms
  //-----------------------------------------------------------

  const double zmin =   -0.5;   // [ns]
  const double zmax =  256.5;   // [ns]
  const int    nz   = (int) ((zmax-zmin) / 1.0);

  JManager<int, TH2D> manager(new TH2D(MAKE_CSTRING("%" << _2SToT), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS-0.5, nz, zmin, zmax));

  typedef JHitR0        hit_type;
  typedef JSuperFrame2D<hit_type>   JSuperFrame2D_t;
  typedef JSuperFrame1D<hit_type>   JSuperFrame1D_t;

  const JMatchL0<hit_type> match(T_ns.getUpperLimit());     // time window self-coincidences [ns]

  const double deadTime_ns  =  deadTime_us * 1e3;

  JObjectMultiplexer<JDAQTimesliceTypes_t, JDAQTimeslice> in(inputFile, selector);


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

  putObject(out, JMeta(argc, argv));

  JDAQHeader header;
  
  for (counter_type counter = 0; in.hasNext() && counter != inputFile.getLimit(); ++counter) {

    STATUS("event: " << setw(10) << counter << '\r'); DEBUG(endl);

    const JDAQTimeslice* timeslice = in.next();

    if (header.getDetectorID() != timeslice->getDetectorID() ||
        header.getRunNumber () != timeslice->getRunNumber ()) {

      header = timeslice->getDAQHeader();

      putObject(out, header);
    }

    for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {

      if (router.hasModule(frame->getModuleID())) {

        TH2D*          h2     = manager[frame->getModuleID()];
        const JModule& module = router.getModule(frame->getModuleID());
        
        JSuperFrame2D_t& buffer = JSuperFrame2D_t::demultiplex(*frame, module);

        buffer.preprocess(option, match);

        JSuperFrame1D_t& data   = JSuperFrame1D_t::multiplex(buffer);

        vector<double> t0(NUMBER_OF_PMTS, numeric_limits<double>::lowest());

        for (JSuperFrame1D_t::const_iterator p = data.begin(); p != data.end(); ) {

          JSuperFrame1D_t::const_iterator q = p;

          // reject genuine coincidences

          double ct_max = -1.0;
          double dt_min =  numeric_limits<double>::max();
          
          while (++q != data.end() && q->getT() - p->getT() < T_ns.getUpperLimit()) {

            const double ct = getDot(module.getPMT(p->getPMT()), module.getPMT(q->getPMT()));
            const double dt = q->getT() - p->getT();

            if (ct > ct_max) {
              ct_max = ct;
            }

            if (dt < dt_min) {
              dt_min = dt;
            }
          }
          
          if (multiplicity(distance(p,q)) && ct_max < ctMax && dt_min > T_ns.getLowerLimit()) {

            for (JSuperFrame1D_t::const_iterator i = p; i != q; ++i) {
            
              if (i->getT() > t0[i->getPMT()] + deadTime_ns) {
                h2->Fill(i->getPMT(), i->getToT());
              }
              
              t0[i->getPMT()] = i->getT();
            }
          }
          
          p = q;
        }
      }
    }
  }
  
  STATUS(endl);
  
  out << manager;

  out.Close();
  out.Write();
}