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

Enumerations
enum	JStatus_t

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 QE of the PMTs of the modules that are out-of-sync are set to 0.
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 maybe missed.

Author

mdejong

Definition in file JTurbot.cc.

Enumeration Type Documentation

enum JStatus_t

Definition at line 46 of file JTurbot.cc.

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

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 71 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             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['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);
  }


  cout.tie(&cerr);

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

  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<JDAQTimeslice,   JDAQTimeslice> >(inputFile))->getEntries() != 0 ||
        (ps = new JTreeScanner< JAssertConversion<JDAQTimesliceL0, 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(); ) {

    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()) {
    
    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()];
        
      for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
          
        const double t1 = *hit + frame->getFrameIndex() * getFrameTime();
          
        for (map_type::const_iterator i = p; i != q; ++i) {
          for (map_type::mapped_type::const_iterator j = i->second.begin(); j != i->second.end(); ++j) {
            
            const double t0 = *j;
            
            h1->Fill(t1 - t0);
          }
        }
      }
    }
  }
  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);
          
      // 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;
      }


      // 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() - TMax_ns) { ++ns; }
            
        if      (fabs(x - ns * getFrameTime()) < 1.0  * TMax_ns)
          sn[ns].put(y);
        else if (fabs(x - ns * getFrameTime()) < 10.0 * TMax_ns)
          bg[ns].put(y);
      }

      DEBUG("Module " << setw(8) << module->getID() << endl);

      for (map<int, JWeight>::iterator i = sn.begin(); i != sn.end(); ) {
        
        const Double_t y = bg[i->first].getTotal() * i->second.getCount() / bg[i->first].getCount();
        const Double_t P = TMath::PoissonI(i->second.getTotal(), y);

        DEBUG("Peak at " << setw(4) << i->first << " [frame time] " 
              << "S/N = " 
              << FIXED(7,1) << i->second.getTotal() << " / "
              << FIXED(7,1) << y                    << ' ' 
              << SCIENTIFIC(12,3) << P              << endl);
            
        if (i->second.getTotal() > y && P < Pmin)
          ++i;               // keep
        else
          sn.erase(i++);     // remove
      }

      if (!(sn.size()         == 1 && 
            sn.begin()->first == 0)) {

        ERROR("Module " << setw(8) << module->getID() << " number of peaks " << sn.size() << endl);

        for (map<int, JWeight>::const_iterator i = sn.begin(); i != sn.end(); ++i) {
          ERROR("Peak at " << setw(4) << i->first << " [frame time] "
                << "S/N = " 
                << FIXED(7,1) << i->second.getTotal()                                                     << " / "
                << FIXED(7,1) << bg[i->first].getTotal() * i->second.getCount() / bg[i->first].getCount() << endl);
        }

        status[module->getID()] = (sn.size() == 1 ? OUT_SYNC : ERROR);
      }
    }
  }

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