JTurbot.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <limits>
#include <map>
#include <vector>
 
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TF1.h"
#include "TMath.h"
 
#include "JROOT/JMinimizer.hh"
 
#include "km3net-dataformat/online/JDAQ.hh"
#include "km3net-dataformat/definitions/pmt_status.hh"
#include "km3net-dataformat/definitions/module_status.hh"
 
#include "JDAQ/JDAQEventIO.hh"
#include "JDAQ/JDAQTimesliceIO.hh"
#include "JDAQ/JDAQEvaluator.hh"
 
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JDAQHitRouter.hh"
#include "JTrigger/JBuildL0.hh"
#include "JTrigger/JBuildL1.hh"
#include "JLang/JSinglePointer.hh"
#include "JROOT/JManager.hh"
#include "JTools/JWeight.hh"
#include "JROOT/JROOTClassSelector.hh"
#include "JSupport/JSingleFileScanner.hh"
#include "JSupport/JTreeScannerInterface.hh"
#include "JSupport/JTreeScanner.hh"
#include "JSupport/JAutoTreeScanner.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
namespace {
 
  using JTOOLS::JWeight;
 
  // status of module
 
  enum JStatus_t {
    DEFAULT       =  0,      //!< module exists in detector
    NODATA        = -2,      //!< module not present in data or no entries
    OUT_SYNC      = -3,      //!< module is out-of-sync
    ERROR         = -4,      //!< module is out-of-sync, possibly multiple time offsets
    IN_SYNC       = +1       //!< module is in-sync
  };
 
  /**
   * Get estimated number of background events.
   *
   * \param  sn          weight of signal
   * \param  bg          weight of background
   * \return             background
   */
  double getBackground(const JWeight& sn, const JWeight& bg)
  {
    return std::max(bg.getTotal(), 1.0) * (double) sn.getCount() / (double) bg.getCount();
  }
}
 
/**
 * \file
 *
 * Example program to search for correlations between triggered events and time slice data.\n
 * 
 * This application can be used to identify optical modules that are out-of-sync with respect to the central White Rabbit clock.\n
 * To this end, time correlations are searched for between triggered events and local coincidences, referred to as L1 hits.
 * The local coincidences are independently obtained from time slice data and can thereby cross the time slice border of a triggered event.\n
 * The MODULE_OUT_OF_SYNC of status the modules that are found to be out-of-sync are set.
 *
 * The number of time slices to look for L1 hits away from the triggered event can be set via command line option <tt>-N <number of time slices></tt>.\n
 * The larger this value, the bigger time offsets can be covered. 
 * To identify the modules that are out-of-sync as fast as possible (without determining the time offset), 
 * this value can be set to 0 with the caveat that a out-of-sync occurring during the data taking run may then be missed.
 *
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
 
  JSingleFileScanner<> inputFile;
  JLimit_t&          numberOfEvents = inputFile.getLimit();
  string             outputFile;
  string             detectorFile;
  bool               overwriteDetector;
  double             TMaxLocal_ns;
  int                numberOfTimeslices;
  double             binWidth_ns;
  double             deadTime_us;
  double             Pmin;
  JROOTClassSelector selector;
  int                qaqc;
  int                debug;
 
  try { 
 
    JParser<> zap("Example program to search for correlations between triggered events and timeslice data.");
    
    zap['f'] = make_field(inputFile,           "input file.");
    zap['o'] = make_field(outputFile,          "optional output file containing ROOT histograms.")         = "";
    zap['a'] = make_field(detectorFile,         "detector file.");
    zap['A'] = make_field(overwriteDetector,    "overwrite detector file provided through '-a' with module (PMT) status.");
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['T'] = make_field(TMaxLocal_ns,        "time window for local coincidences (L1),")                 = 20.0;
    zap['N'] = make_field(numberOfTimeslices,  "time slice difference between triggered event and L1.")    = 100;
    zap['W'] = make_field(binWidth_ns,         "bin width for output histograms." )                        = 10e3;
    zap['D'] = make_field(deadTime_us,         "L1 dead time (us)")                                        = 200.0;
    zap['p'] = make_field(Pmin,                "minimal probability for background to be signal.")         = 1.0e-7;
    zap['C'] = make_field(selector,            "data type selection (default is all).")                    = "", getROOTClassSelection<JDAQTimesliceTypes_t>();
    zap['Q'] = make_field(qaqc)                =  0;
    zap['d'] = make_field(debug)               =  2;
 
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
 
  JDetector detector;
 
  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }
 
  if (detector.empty()) {
    FATAL("Empty detector " << detectorFile << endl);
  }
  
  const JDAQHitRouter router(detector);
 
  const double TMax_ns      =  max(binWidth_ns, getMaximalTime(detector));  // signal time window
  const double BOOST        =  20.0;                                        // scaling of time window for background
  const double deadTime_ns  =  deadTime_us * 1e3;
 
  NOTICE("Time window " << FIXED(7,1) << TMax_ns << " [ns]" << endl);
 
  typedef  double    hit_type;
  
  JBuildL0<hit_type> buildL0;
  JBuildL1<hit_type> buildL1(JBuildL1Parameters(TMaxLocal_ns, true));
  vector  <hit_type> data;  
  
  typedef JManager<int, TH1D>  JManager_t;
 
  const Double_t xmin = -(numberOfTimeslices + 1) * getFrameTime() - 0.5*binWidth_ns;
  const Double_t xmax = +(numberOfTimeslices + 1) * getFrameTime() + 0.5*binWidth_ns;
  const Int_t    nx   = (Int_t) ((xmax - xmin) / binWidth_ns);
  
  JManager_t manager(new TH1D("M_%", NULL, nx, xmin, xmax));
 
  JTreeScanner<>::debug = debug;
  
  JSinglePointer< JTreeScannerInterface<JDAQTimeslice> > ps;
 
  JAutoTreeScanner<JDAQTimeslice> zmap = JType<JDAQTimesliceTypes_t>();
  
  if (selector == "") {
 
    if ((ps = new JTreeScanner< JAssertConversion<JDAQTimesliceSN, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
        (ps = new JTreeScanner< JAssertConversion<JDAQTimesliceL2, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
        (ps = new JTreeScanner< JAssertConversion<JDAQTimesliceL1, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
        (ps = new JTreeScanner< JAssertConversion<JDAQTimesliceL0, JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
        (ps = new JTreeScanner< JAssertConversion<JDAQTimeslice,   JDAQTimeslice> >(inputFile))->getEntries() != 0) {
    } else {
      FATAL("No timeslice data." << endl);
    }
 
    NOTICE("Selected class " << ps->getClassName() << endl);
 
  } else {
 
    ps = zmap[selector];
    
    ps->configure(inputFile);
  }
 
  ps->setLimit(inputFile.getLimit());
 
  typedef map<int, vector<double> > map_type;    // frame index -> event times
 
  map_type buffer;
    
  for (JTreeScanner<JDAQEvent> in(inputFile.getFilename()); in.hasNext(); ) {
 
    if (in.getCounter()%1000 == 0) {
      STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);
    }
 
    JDAQEvent* event = in.next();
 
    // Average time of triggered hits
    
    double t0 = 0.0;
    
    for (JDAQEvent::const_iterator<JDAQTriggeredHit> hit = event->begin<JDAQTriggeredHit>(); hit != event->end<JDAQTriggeredHit>(); ++hit) {
      t0 += getTime(*hit, router.getPMT(*hit));
    }
    
    t0 /= event->size<JDAQTriggeredHit>();
    t0 += event->getFrameIndex() * getFrameTime();
    
    buffer[event->getFrameIndex()].push_back(t0);
  }
  STATUS(endl);
 
 
  while (ps->hasNext()) {
    
    if (ps->getCounter()%1000 == 0) {
      STATUS("event: " << setw(10) << ps->getCounter() << '\r'); DEBUG(endl);
    }
 
    JDAQTimeslice* timeslice = ps->next();
 
    map_type::const_iterator p = buffer.lower_bound(timeslice->getFrameIndex() - numberOfTimeslices);
    map_type::const_iterator q = buffer.upper_bound(timeslice->getFrameIndex() + numberOfTimeslices);
 
    int number_of_events = 0;
 
    for (map_type::const_iterator i = p; i != q; ++i) {
      number_of_events += i->second.size();
    }
 
    if (number_of_events == 0) {
      continue;
    }
    
    for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
 
      data.clear();
 
      buildL1(*frame, router.getModule(frame->getModuleID()), back_inserter(data));
 
      TH1D*  h1 = manager[frame->getModuleID()];
 
      double t1 = numeric_limits<double>::lowest();
 
      for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
          
        const double t2 = *hit + frame->getFrameIndex() * getFrameTime();
          
        if (t2 > t1 + deadTime_ns) {
 
          for (map_type::const_iterator i = p; i != q; ++i) {
            for (map_type::mapped_type::const_iterator t0 = i->second.begin(); t0 != i->second.end(); ++t0) {
              h1->Fill(t2 - *t0);
            }
          }
 
          t1 = t2;
        }
      }
    }
  }
  STATUS(endl);
 
 
  // fit function
 
  TF1 f1("f1", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2])) + [3]");
    
  string option = "L";
 
  if (debug < debug_t) {
    option += "Q";
  }
 
 
  map<int, JStatus_t> status;     
 
  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
 
    if (module->getFloor() != 0) {
 
      status[module->getID()] = DEFAULT;
    
      JManager_t::iterator p = manager.find(module->getID());
 
      if (p == manager.end() || p->second->GetEntries() == 0) { 
      
        status[module->getID()] = NODATA;
 
        continue;
      }
 
      TH1D* h1 = p->second;
 
      // start values
          
      Double_t ymax  =   0.0;
      Double_t x0    =   0.0;    // [ns]
      Double_t sigma = 250.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;
        }
 
        ymin += y;
      }
 
      ymin /= h1->GetNbinsX();
        
      f1.SetParameter(0, ymax);
      f1.SetParameter(1, x0);
      if (binWidth_ns < sigma) 
        f1.SetParameter(2, sigma);
      else
        f1.FixParameter(2, binWidth_ns/sqrt(12.0));
      f1.SetParameter(3, ymin);
          
      for (Int_t i = 0; i != f1.GetNpar(); ++i) {
        f1.SetParError(i, 0.0);
      }
 
      // fit
    
      h1->Fit(&f1, option.c_str(), "same", x0 - 5 * sigma, x0 + 5 * sigma);
 
      if (fabs(f1.GetParameter(1)) <= 0.5*getFrameTime() &&    // position
          f1.GetParameter(0)       >= f1.GetParameter(3)) {    // S/N
 
        status[module->getID()] = IN_SYNC;
      }
 
      double fm = fmod(f1.GetParameter(1), getFrameTime());
 
      if (fm < -0.5 * getFrameTime()) { fm += getFrameTime(); }
      if (fm > +0.5 * getFrameTime()) { fm -= getFrameTime(); }
 
      STATUS("Module/fit "
             << setw(10)    << module->getID()    << ' '
             << FIXED(15,3) << f1.GetParameter(1) << ' '
             << FIXED(12,3) << f1.GetParameter(0) << ' '
             << FIXED(12,3) << f1.GetParameter(3) << ' '
             << FIXED(12,3) << fm                 << ' '
             << status[module->getID()] << endl);
 
      // look for peaks at time offsets equal to a multiple of frame time
          
      map<int, JWeight> bg;    // background
      map<int, JWeight> sn;    // signal(s)
 
      for (Int_t i = 1, ns = -(numberOfTimeslices + 1); i <= h1->GetNbinsX(); ++i) {
 
        const Double_t x = h1->GetBinCenter (i);
        const Double_t y = h1->GetBinContent(i);
 
        while (x > (ns + 1) * getFrameTime() - BOOST * TMax_ns) {
          ++ns;
        }
            
        if      (fabs(x - ns * getFrameTime()) < TMax_ns)
          sn[ns].put(y);
        else if (fabs(x - ns * getFrameTime()) < BOOST * TMax_ns)
          bg[ns].put(y);
      }
 
      for (map<int, JWeight>::iterator i = sn.begin(); i != sn.end(); ) {
        
        const Double_t y = getBackground(i->second, bg[i->first]);
        const Double_t P = TMath::PoissonI(i->second.getTotal(), y);
 
        if (debug >= debug_t || status[module->getID()] != IN_SYNC) {
          cout << "Module/peak " << setw(10) << module->getID() << ' '
               << setw(4) << i->first << ' '
               << "["
               << FIXED(15,1) << i->first * getFrameTime() - TMax_ns
               << ","
               << FIXED(15,1) << i->first * getFrameTime() + TMax_ns
               << "]"
               << " S/N = " 
               << FIXED(7,1) << i->second.getTotal() << " / "
               << FIXED(7,1) << y                    << ' ' 
               << SCIENTIFIC(12,3) << P              << ' '
               << (i->second.getTotal() > y && P < Pmin && i->first != 0 ? "***" : "") << endl;
        }
        
        if (i->second.getTotal() > y && P < Pmin)
          ++i;               // keep
        else
          sn.erase(i++);     // remove
      }
 
      if (!(sn.size()         == 1 && 
            sn.begin()->first == 0)) {
 
        status[module->getID()] = (sn.size() == 1 ? OUT_SYNC : ERROR);
 
        ERROR("Module/error "
              << setw(10)    << module->getID()    << ' '
              << "number of peaks " << sn.size()   << ' '
              << "peak " << (sn.size() == 1 ? MAKE_CSTRING(sn.begin()->first) : "?") << ' ' 
              << status[module->getID()] << endl);
      }
    }
  }
 
  if (overwriteDetector) {
 
    detector.comment.add(JMeta(argc, argv));
 
    for (map<int, JStatus_t>::const_iterator i = status.begin(); i != status.end(); ++i) {
 
      if (i->second != IN_SYNC && 
          i->second != NODATA) {
 
        NOTICE("Module " << setw(8) << i->first << " set out-of-sync." << endl);
 
        JModule& module = detector.getModule(router.getAddress(i->first));
 
        module.getStatus().set(MODULE_OUT_OF_SYNC);
        
        for (JModule::iterator pmt = module.begin(); pmt != module.end(); ++pmt) {
          pmt->getStatus().set(OUT_OF_SYNC);
        }
      }
    }
 
    store(detectorFile.c_str(), detector);
  }
 
  if (outputFile != "") {
    manager.Write(outputFile.c_str());
  }
 
  int nin  = 0;
  int nout = 0;
 
  for (map<int, JStatus_t>::const_iterator i = status.begin(); i != status.end(); ++i) {
    if (i->second == IN_SYNC) {
      ++nin;
    }
    if (i->second != IN_SYNC && 
        i->second != NODATA) {
      ++nout;
    }
  }
 
  NOTICE("Number of modules out/in sync " << nout << '/' << nin << endl);
 
  QAQC(nin << ' ' << nout << endl);
 
  return 0;
}