JTurbot2D.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 "TH2D.h"
#include "TF1.h"

#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 "JROOT/JManager.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 "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


/**
 * \file
 *
 * Example program to search for out of sync shifts around integral timeslices evolving over time.
 * \author shallmann
 */
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;
  double             TMax_ns;
  int                numberOfTimeslices;
  double             binWidth_ns;
  double             binWidth_min_timeEvolution;
  JROOTClassSelector selector;
  int                debug;

  try { 

    JParser<> zap("Example program to search for out of sync shifts around integral timeslices evolving over time.");
    
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)          = "";
    zap['a'] = make_field(detectorFile);
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['T'] = make_field(TMax_ns)             = 20.0;
    zap['N'] = make_field(numberOfTimeslices)  = 40;
    zap['W'] = make_field(binWidth_ns)         = 10e3;
    zap['M'] = make_field(binWidth_min_timeEvolution) = 2;
    zap['C'] = make_field(selector)            = "", getROOTClassSelection<JDAQTimesliceTypes_t>();
    zap['d'] = make_field(debug)               =  2;

    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }



  cout.tie(&cerr);

  
  map<int, int> MASK;                // mask PMTs with permament high-rate veto.

  const int WR     =  0x80000000;    // White Rabbit
  
  MASK[808969848]  =  0x00000020;    // floor 03, pmt 05
  MASK[809544061]  =  0x00000080;    // floor 05, pmt 07
  MASK[808432835]  =  0x00004000;    // floor 18, pmt 14

  
  JDetector detector;

  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }
  
  const JDAQHitRouter router(detector);


  typedef  double    hit_type;
  
  JBuildL0<hit_type> buildL0;
  JBuildL1<hit_type> buildL1(JBuildL1Parameters(TMax_ns, true));
  vector  <hit_type> data;  
  
  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);
  }

  typedef JManager<int, TH2D>  JManager_t;

  // shift in integer time slices
  const Double_t ymin = -(numberOfTimeslices + 0.5);
  const Double_t ymax = +(numberOfTimeslices + 0.5);
  const Int_t    ny   = (Int_t) (ymax - ymin);
  // time in M minute binning
  const Double_t xmin = inputFile.getLimit().min()*getFrameTime()*1e-9;// time of first event - offset
  const Double_t xmax = min(inputFile.getLimit().max(), ps->getEntries())*getFrameTime()*1e-9; // time of last event
  const Int_t    nx   = (Int_t) ((xmax - xmin) / (60*binWidth_min_timeEvolution)) + 1;
  
  JManager_t manager(new TH2D("M2D_%", ";time [s];shift [timeslices]", nx, xmin, xmax, ny, ymin, ymax));
  
  typedef map<int, vector<double> > map_type;    // frame index -> event times

  map_type buffer;
    
  for (JTreeScanner<JDAQEvent> in(inputFile.getFilename(), inputFile.getLimit()); 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) {

      // Test that none of the PMTs has high-rate veto.
          
      if ((frame->getStatus() & ~MASK[frame->getModuleID()] & ~WR) == 0) {
          
        data.clear();

        buildL1(*frame, router.getModule(frame->getModuleID()), back_inserter(data));
        
        TH2D* h2 = 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;

          // fill only shift is near integer time slice
              const double timeFromInt = (t1 - t0)/getFrameTime() - round((t1 - t0)/getFrameTime());
          if ( abs(timeFromInt)*getFrameTime() < binWidth_ns/2. ){
                h2->Fill(timeslice->getFrameIndex() * getFrameTime() * 1e-9, round(t1 - t0)/getFrameTime());
          }
            }
          }
        }
      }
    }
  }
  STATUS(endl);

  if (outputFile != "") {
    manager.Write(outputFile.c_str());
  }
}