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 "JDAQ/JDAQEventIO.hh"
#include "JDAQ/JDAQTimesliceIO.hh"
#include "JDAQ/JDAQEvaluator.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JDAQHitRouter.hh"
#include "JDetector/JPMTParametersMap.hh"
#include "JTrigger/JBuildL0.hh"
#include "JTrigger/JBuildL1.hh"
#include "JLang/JSinglePointer.hh"
#include "JGizmo/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 45 of file JTurbot.cc.

                 {
    DEFAULT       =  0,      //!< module exists in detector
    NODATA        = -2,      //!< module not present in data or no entries
    OUT_OF_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 70 of file JTurbot.cc.

{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  JSingleFileScanner<> inputFile;
  JLimit_t&          numberOfEvents = inputFile.getLimit();
  string             outputFile;
  string             detectorFile;
  double             TMaxLocal_ns;
  string             pmtFile;
  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['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['T'] = make_field(TMaxLocal_ns,        "time window for local coincidences (L1),")                 = 20.0;
    zap['P'] = make_field(pmtFile,             "optional PMT file (to set QE of PMTs).")                   = "";
    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) {

    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_OF_SYNC : ERROR);
    }
  }


  if (pmtFile != "") {

    JPMTParametersMap parameters;
    
    try {
      parameters.load(pmtFile.c_str());
    }
    catch(const JException& error) {}

    for (map<int, JStatus_t>::const_iterator i = status.begin(); i != status.end(); ++i) {

      if (i->second != IN_SYNC) {

        NOTICE("Module " << setw(8) << i->first << " set QEs of all PMTs to zero." << endl);
        
        for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
          parameters[JPMTIdentifier(i->first, pmt)].QE = 0.0;
        }
      }
    }

    ofstream out(pmtFile.c_str());

    parameters.comment.add(JMeta(argc, argv));
    
    out << parameters << endl;

    out.close();
  }

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