JPulsar.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 "TF1.h"
#include "TFitResult.h"
 
#include "JROOT/JMinimizer.hh"
 
#include "JDAQ/JDAQTimesliceIO.hh"
 
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JPMTIdentifier.hh"
 
#include "JTrigger/JHit.hh"
#include "JTrigger/JFrame.hh"
#include "JTrigger/JTimeslice.hh"
#include "JTrigger/JHitR0.hh"
#include "JTrigger/JBuildL0.hh"
 
#include "JCalibrate/JTDC_t.hh"
 
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"
 
#include "JTools/JRange.hh"
#include "JLang/JObjectMultiplexer.hh"
#include "JROOT/JROOTClassSelector.hh"
 
#include "Jeep/JProperties.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
namespace {
 
  /**
   * Get time within given period.
   *
   * The returned time is constraint to the interval <tt>[-T/2,+T/2]</tt>.
   *
   * \param  t1        time
   * \param  T         period
   * \return           time
   */
  inline double get_time(double t1, const double T)
  {
    for (int buffer[] = { 1000, 100, 10, 1, 0 }, *i = buffer; *i != 0; ++i) {
 
      const double xmin = -0.5 * (*i) * T;
      const double xmax = +0.5 * (*i) * T;
 
      while (t1 > xmax) { t1 -= (*i) * T; }
      while (t1 < xmin) { t1 += (*i) * T; }
    }
 
    return t1;
  }
}
 
 
/**
 * \file
 *
 * Application for dark room time calibration using external laser.
 *
 * The time calibration is made for the specified reference PMTs (option <tt>-! "<module identifier> <PMT channel>"</tt>.\n
 * In this, the wild card value <tt>-1</tt> can be used.
 * A ROOT histogram corresponding to each reference PMT is created, filled, fitted and stored to the specified output file.\n
 * The detector can be updated using option <tt>-A</tt>\n
 * The laser should be synchronised with the same clock system as the detector and 
 * produce single pulses at a fixed frequency (option <tt>-l <laserFrequency_Hz></tt>.
 *
 * In addition, the transit-time distribution of the reference PMTs can be monitored (see JLegolas.cc).\n
 * For this, the expectation values should be less than one and the time-over-threshold distribution nominal.
 *
 * \author mbouwhuis
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
 
  typedef   JRange<double>  JRange_t;
 
  JMultipleFileScanner<JDAQTimesliceTypes_t> inputFile;
  JLimit_t&          numberOfEvents = inputFile.getLimit();
  string             outputFile;
  string             detectorFile;
  JTDC_t             TDC;
  double             laserFrequency_Hz;
  bool               overwriteDetector;
  JROOTClassSelector selector;
  string             option;
  double             Wmin              =   1e3;    // Minimal signal intensity
  JRange_t           T_ns;
  JRange_t           X;
  int                debug;
 
  try { 
 
    JProperties properties;
 
    properties.insert(gmake_property(Wmin));
 
    JParser<> zap("Application for dark room time calibration.");
    
    zap['f'] = make_field(inputFile,         "input file (time slice data from laser calibration).");
    zap['o'] = make_field(outputFile,        "output file.")                                     = "pulsar.root";
    zap['a'] = make_field(detectorFile,      "detector file.");
    zap['!'] = make_field(TDC,
                          "Set reference PMTs, e.g."
                          "\n-! \"808969848 0 808982077 23\" sets PMT 0 of module 808969848 and PMT 23 of module 808982077 as references.");
    zap['n'] = make_field(numberOfEvents)    = JLimit::max(); 
    zap['l'] = make_field(laserFrequency_Hz, "laser frequency [Hz]")                             = 10000;
    zap['A'] = make_field(overwriteDetector, "overwrite detector file provided through '-a' with correct time offsets.");
    zap['O'] = make_field(option,            "ROOT fit option, see TH1::Fit.")                   = "LS";
    zap['@'] = make_field(properties)                                                            = JPARSER::initialised();
    zap['C'] = make_field(selector,          "timeslice selector, e.g. JDAQTimesliceL1.")        = getROOTClassSelection<JDAQTimesliceTypes_t>();
    zap['T'] = make_field(T_ns,              "time window for time-over-threshold monitor")      = JRange_t(-10.0, +10.0);
    zap['x'] = make_field(X,                 "time window for histogram")                        = JRange_t();
    zap['d'] = make_field(debug,             "debug.")                                           = 1;
    
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
 
  TDC.is_valid(true);
 
  if (laserFrequency_Hz <= 0.0) {
    FATAL("Invalid laser frequency " << laserFrequency_Hz << endl);
  }
 
  const double laserPeriod_ns = 1.0e9 / laserFrequency_Hz;
 
  if (option.find('R') == string::npos) { option += 'R'; }
  if (option.find('S') == string::npos) { option += 'S'; }
  if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; }
 
 
  JDetector detector;
 
  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }
 
  const JModuleRouter moduleRouter(detector);
 
 
  // ROOT I/O
 
  TFile out(outputFile.c_str(), "recreate");
 
  putObject(out, JMeta(argc, argv));
 
  typedef map<JPMTIdentifier, TH1D*> map_type;
 
  map_type     zmap;
 
  const double xmin = (X == JRange_t() ? -0.5 * laserPeriod_ns : X.getLowerLimit());
  const double xmax = (X == JRange_t() ? +0.5 * laserPeriod_ns : X.getUpperLimit());
  const int    nx   =  2 * (int) (xmax - xmin);
 
  TH1D h0("h0", NULL, nx, xmin, xmax);
  TH1D h1("h1", NULL, 256, -0.5, +255.5);
 
  map<JPMTIdentifier, int> counts;
 
  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
 
    const JTDC_t::range_type range = TDC.equal_range(module->getID());
 
    for (JTDC_t::const_iterator i = range.first; i != range.second; ++i) {
 
      NOTICE("Reference PMT at"
             << " string  " << setw(3) << module->getString() 
             << " floor   " << setw(2) << module->getFloor()
             << " module  " << setw(8) << module->getID() 
             << " channel " << setw(2) << i->second << endl);
 
      const JPMTIdentifier id(module->getID(),i->second);
      
      ostringstream os;
      
      os << getLabel(module->getLocation()) << ' ' << getLabel(id);
      
      zmap.insert(make_pair(id, new TH1D(os.str().c_str(), NULL, nx, xmin, xmax)));
    }
  }
 
 
  typedef vector<JHitR0>     JDataL0_t;
 
  const JBuildL0<JHitR0>     buildL0;
 
 
  JObjectMultiplexer<JDAQTimesliceTypes_t, JDAQTimeslice> in(inputFile, selector);
  
  counter_type counter = 0;
 
  for ( ; in.hasNext() && counter != inputFile.getLimit(); ++counter) {
 
    STATUS("event: " << setw(10) << counter << '\r'); DEBUG(endl);
 
    const JDAQTimeslice* timeslice = in.next();
 
    for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
 
      const JTDC_t::range_type range = TDC.equal_range(frame->getModuleID());
 
      if (range.first != range.second) {
 
        const double t0 = get_time(getTimeOfFrame(frame->getFrameIndex()), laserPeriod_ns);
 
        JDataL0_t dataL0;
 
        buildL0(*frame, moduleRouter.getModule(frame->getModuleID()), back_inserter(dataL0));
 
        for (JDataL0_t::const_iterator hit = dataL0.begin(); hit != dataL0.end(); ++hit) {
 
          const JPMTIdentifier id(frame->getModuleID(), hit->getPMT());
 
          map_type::const_iterator p = zmap.find(id);
 
          if (p != zmap.end()) {
 
            const double t1 = get_time(t0 + hit->getT(), laserPeriod_ns);
 
            p->second->Fill(t1);
 
            h0.Fill(t1);
 
            if (T_ns(t1)) {
 
              h1.Fill(hit->getToT());
 
              counts[id] += 1;
            }
          }
        }
      }
    }
  }
  STATUS(endl);
 
 
  map<JPMTIdentifier, double> t0;     // fitted time offsets
 
 
  TF1 f1("f1", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2]))/(TMath::Sqrt(2.0*TMath::Pi())*[2]) + [3]");
 
  for (map_type::iterator it = zmap.begin(); it != zmap.end(); ++it) {
 
    JPMTIdentifier pmt = it->first;
    TH1D*          h1  = it->second;
 
    if (h1->GetEntries() == 0) {
      WARNING("Histogram " << h1->GetName() << " empty" << endl);
      continue;
    }
 
    STATUS("--- PMT = " << pmt << "; histogram " << h1->GetName() <<  endl);
 
    // start values
    
    Double_t ymax  = 0.0;
    Double_t x0    = 0.0;    // [ns]
    Double_t sigma = 2.0;    // [ns]
    Double_t ymin  = 0.0;
    
    for (int i = 1; i != h1->GetNbinsX(); ++i) {
      
      const Double_t x = h1->GetBinCenter (i);
      const Double_t y = h1->GetBinContent(i);
      
      if (y > ymax) {
        ymax = y;
        x0   = x;
      }
    }
 
    f1.SetParameter(0, ymax/sigma);
    f1.SetParameter(1, x0);
    f1.SetParameter(2, sigma);
    f1.SetParameter(3, ymin);
 
    for (Int_t i = 0; i != f1.GetNpar(); ++i) {
      f1.SetParError(i, 0.0);
    }
    
    
    // fit
    
    TFitResultPtr result = h1->Fit(&f1, option.c_str(), "same", x0 - 3 * sigma, x0 + 3 * sigma);
 
    if (result.Get() == NULL) {
      FATAL("Invalid TFitResultPtr " << h1->GetName() << endl);
    }
    
    if (debug >= notice_t || !result->IsValid() || f1.GetParameter(0) < Wmin) {
      cout << "Histogram " << h1->GetName() << " fit " << (result->IsValid() ? "ok" : "failed") << endl;
      cout << "\tw  = " << FIXED(12,3) << f1.GetParameter(0) << endl;
      cout << "\tx0 = " << FIXED(12,3) << x0                 << endl;
      cout << "\tt0 = " << FIXED(12,3) << f1.GetParameter(1) << endl;
    }
 
 
    // check for multiple peaks
 
    int number_of_peaks = 0;
 
    Double_t dx = 2.0 * f1.GetParameter(2);
    Double_t W  = 0.0;
    Double_t Y  = f1.GetParameter(3);
 
    if (dx < (xmax - xmin) / nx) {
      dx = (xmax - xmin) / nx;                    // minimal width
    }
 
    for (Int_t il = 1, ir = il; ir <= nx; ) {
 
      for ( ; ir <= nx && h1->GetBinCenter(ir) <= h1->GetBinCenter(il) + dx; ++ir) {
        W += h1->GetBinContent(ir) - Y;
      }
 
      if (W >= Wmin) {
 
        number_of_peaks += 1;
 
        W   = 0.0;                                // reset
        il  = ir;
        ir  = il;
 
      } else {
 
        W  -= h1->GetBinContent(il) - Y;          // shift
        il += 1;
      }
    }
 
    if (number_of_peaks != 1) {
      WARNING("Number of peaks " << h1->GetName() << ' ' << number_of_peaks << " != 1" << endl);
    }
 
    if (result->IsValid() && f1.GetParameter(0) >= Wmin) {
      t0[pmt] = f1.GetParameter(1);
    }
  }
 
  out.Write();
  out.Close();
 
 
  if (counter != 0) {
 
    const double W = laserFrequency_Hz * counter * getFrameTime() * 1.0e-9;
 
    NOTICE("Expection values [npe]" << endl);
 
    for (map<JPMTIdentifier, int>::const_iterator i = counts.begin(); i != counts.end(); ++i) {
      NOTICE(i->first << ' ' << FIXED(7,3) << i->second / W << endl);
    }
  }
 
 
  if (overwriteDetector) { 
 
    int errors = 0;
 
    for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
 
      const JTDC_t::range_type range = TDC.equal_range(module->getID());
 
      if (range.first != range.second) {
 
        const JPMTIdentifier id(module->getID(), range.first->second); // select t0 of the first reference PMT in this optical module
 
        map<JPMTIdentifier, double>::const_iterator p0 = t0.find(id);
 
        if (p0 != t0.end()) {
 
          const double t1 = p0->second;
 
          for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
 
            map<JPMTIdentifier, double>::const_iterator p1 = t0.find(JPMTIdentifier(module->getID(),pmt));
          
            if (p1 != t0.end())
              module->getPMT(pmt).subT0(p1->second);                   // offset with t0 of this reference PMT
            else
              module->getPMT(pmt).subT0(t1);                           // offset with t0 of first reference PMT in same optical module
          }
 
        } else {
 
          if (!module->getPMT(id.getPMTAddress()).has(PMT_DISABLE)) {
            ++errors;
          }
 
          ERROR("Module/PMT " 
                << setw(10) << module->getID() << "/" << FILL(2,'0') << id.getPMTAddress() << FILL() << ' '
                << "missing or insufficient signal." << endl); 
        }
      }
    }
 
    if (errors == 0) {
 
      NOTICE("Store calibration data on file " << detectorFile << endl);
 
      detector.comment.add(JMeta(argc, argv));
 
      store(detectorFile, detector);
 
    } else {
 
      FATAL("Number of errors " << errors << " != 0" << endl);
    }
  }
}