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 "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 "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/JParser.hh"
#include "Jeep/JMessage.hh"


namespace {

  typedef std::multimap<int, int>             JTDC_t;

  typedef std::pair<JTDC_t::const_iterator,
                    JTDC_t::const_iterator>   range_type;

  static const int WILD_CARD = -1;            //!< Wild catd for module identifier asnd PMT channel


  /**
   * Get range of TDC identifiers for given module identifier.
   *
   * \param  TDC       TDC identifiers
   * \param  id        module identifier
   * \return           range of TDC identifiers
   */
  inline range_type get_range(const JTDC_t& TDC, const int id)
  {
    range_type range = TDC.equal_range(id);

    if (range.first == range.second) {
      range = TDC.equal_range(WILD_CARD);
    }

    return range;
  }


  /**
   * 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           tim
   */
  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.
 * \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;
  JRange_t           T_ns;
  int                debug;

  try { 

    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['C'] = make_field(selector,          "timeslice selector, e.g. JDAQTimesliceL1.")        = getROOTClassSelection<JDAQTimesliceTypes_t>();
    zap['T'] = make_field(T_ns)              = JRange_t(-10.0, +10.0);
    zap['d'] = make_field(debug,             "debug.")                                           = 1;
    
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }


  for (JTDC_t::iterator tdc = TDC.begin(); tdc != TDC.end(); ) {

    if (tdc->second == WILD_CARD) {

      const int id = tdc->first;

      TDC.erase(tdc);

      for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
        tdc = TDC.insert(make_pair(id,pmt));
      }
      
    } else {

      ++tdc;
    }
  }

  for (JTDC_t::const_iterator tdc = TDC.begin(); tdc != TDC.end(); ++tdc) {
    if (tdc->second < 0 || tdc->second >= NUMBER_OF_PMTS) {
      FATAL("Illegal TDC (= readout channel) identifier: " << tdc->second << " {0, .., " << NUMBER_OF_PMTS - 1 << "}" << endl);
    }
  }

  if (laserFrequency_Hz <= 0.0) {
    FATAL("Invalid laser frequency " << laserFrequency_Hz << endl);
  }


  const double laserPeriod_ns = 1.0e9 / laserFrequency_Hz;

  cout.tie(&cerr);

  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 int    nx   =  2   * (int) laserPeriod_ns;
  const double xmin = -0.5 * laserPeriod_ns;
  const double xmax =  0.5 * laserPeriod_ns;

  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 range_type range = get_range(TDC, 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 range_type range = get_range(TDC, 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])) + [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);
    f1.SetParameter(1, x0);
    f1.SetParameter(2, sigma);
    f1.SetParameter(3, ymin);
    
    
    // fit
    
    TFitResultPtr result = h1->Fit(&f1, option.c_str(), "same", x0 - 3 * sigma, x0 + 3 * sigma);

    
    if (debug >= notice_t || !result->IsValid()) {
      cout << "Histogram " << h1->GetName() << " fit " << (result->IsValid() ? "" : "failed") << endl;
      cout << "\tx0 = " << FIXED(12,3) << x0                 << endl;
      cout << "\tt0 = " << FIXED(12,3) << f1.GetParameter(1) << endl;
    }

    t0[pmt] = f1.GetParameter(1);    
  }


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

    NOTICE("Store calibration data on file " << detectorFile << endl);

    detector.comment.add(JMeta(argc, argv));

    for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {

      const range_type range = get_range(TDC, module->getID());

      if (range.first != range.second) {

        const double t1 = t0[JPMTIdentifier(range.first->first,
                                            range.first->second)];   // t0 of the first reference PMT in this optical module

        for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {

          map<JPMTIdentifier, double>::const_iterator p = t0.find(JPMTIdentifier(module->getID(),pmt));
          
          if (p != t0.end())
            module->getPMT(pmt).subT0(p->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
        }
      }
    }

    store(detectorFile, detector);
  }

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