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/JDAQEvent.hh"
#include "JDAQ/JDAQTimeslice.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/JHitL0.hh"
#include "JTrigger/JBuildL0.hh"

#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"

#include "JLang/JObjectMultiplexer.hh"
#include "JROOT/JROOTClassSelector.hh"

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


/**
 * \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 multimap<int, int>  JTDC_t;
  typedef std::pair<JTDC_t::const_iterator, JTDC_t::const_iterator>  range_type;

  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;
  int                debug;

  try { 

    JParser<> zap("Application for dark room time calibration.");
    
    zap['f'] = make_field(inputFile,         "input file (output from JCalibrateK40).");
    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.")                   = "L";
    zap['C'] = make_field(selector,          "timeslice selector, e.g. JDAQTimesliceL1.")        = getROOTClassSelection<JDAQTimesliceTypes_t>();
    zap['d'] = make_field(debug,             "debug.")                                           = 1;
    
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }


  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;


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

    const 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);
      
      ostringstream os;
      
      os  << "s"  << module->getString()
          << "_f" << module->getFloor() 
          << "_m" << module->getID() 
          << "_c" << setfill('0') << setw(2) << i->second;
      
      zmap.insert(make_pair(JPMTIdentifier(module->getID(),i->second), new TH1D(os.str().c_str(), NULL, nx, xmin, xmax)));
    }
  }


  typedef vector<JHitL0>     JDataL0_t;

  const JBuildL0<JHitL0>     buildL0;


  JObjectMultiplexer<JDAQTimesliceTypes_t, JDAQTimeslice> in(inputFile, selector);
  
  for (counter_type counter = 0; 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 int id = frame->getModuleID();

      DEBUG("Module id =" << id << endl);

      double t0 = getTimeOfFrame(frame->getFrameIndex());

      while (t0 > xmax) {
        t0 -= laserPeriod_ns;  
      }

      JDataL0_t dataL0;

      buildL0(*frame, moduleRouter.getModule(id), back_inserter(dataL0));

      for (JDataL0_t::const_iterator hit = dataL0.begin(); hit != dataL0.end(); ++hit) {
        
        map_type::const_iterator p = zmap.find(JPMTIdentifier(id,hit->getPMTAddress()));
        
        if (p != zmap.end()) {
          
          double t1 = t0 + hit->getT();
          
          while (t1 > xmax) {
            t1 -= laserPeriod_ns;  
          }
          while (t1 < xmin) {
            t1 += laserPeriod_ns;  
          }
          
          p->second->Fill(t1);
        }
      }
    }
  }

  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 (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 = TDC.equal_range(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();  
}