JTurbot.cc File Reference

Example program to search for correlations between triggered events and time slice data. More...

#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"

Go to the source code of this file.

Functions
int	main (int argc, char **argv)

Detailed Description

Example program to search for correlations between triggered events and time slice data.

This application can be used to identify optical modules that are out-of-sync with respect to the central White Rabbit clock.
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.
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 -N <number of time slices>.
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

Definition in file JTurbot.cc.

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 88 of file JTurbot.cc.

{
  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;
}