JAcousticsEventReader.cc

Go to the documentation of this file.

#include <iostream>
#include <iomanip>
#include <limits>
#include <map>
 
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
 
#include "JDAQ/JDAQEventIO.hh"
#include "JDAQ/JDAQSummarysliceIO.hh"
#include "JDAQ/JDAQEvaluator.hh"
 
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JTreeScanner.hh"
#include "JSupport/JSupport.hh"
 
#include "JROOT/JRootToolkit.hh"
 
#include "JAcoustics/JEvent.hh"
#include "JAcoustics/JSupport.hh"
#include "JAcoustics/JAcousticsToolkit.hh"
#include "JAcoustics/JEvtToolkit.hh"
 
#include "Jeep/JProperties.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * \file
 *
 * Example program to process DAQ and acoustic events.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
 
  JMultipleFileScanner<JDAQEvent>  inputFile;
  JMultipleFileScanner<>  acousticsFile;
  JLimit_t&               numberOfEvents = inputFile.getLimit();
  string                  outputFile;
  double                  Tmax_s         = 300.0;               // time window [s]
  size_t                  Nmin           =   3;                 // minimum number of emitters
  int                     debug;
 
  try {
 
    JProperties properties;
 
    properties.insert(gmake_property(Tmax_s));
    properties.insert(gmake_property(Nmin));
 
    JParser<> zap("Example program to process DAQ and acoustic events.");
 
    zap['f'] = make_field(inputFile);
    zap['A'] = make_field(acousticsFile);
    zap['o'] = make_field(outputFile)       = "example.root";
    zap['n'] = make_field(numberOfEvents)   = JLimit::max();
    zap['@'] = make_field(properties)       = JPARSER::initialised();
    zap['d'] = make_field(debug)            = 2;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  TH1D h1("h1", NULL, 101, -0.5,   +100.5);
  TH1D hn("hn", NULL, 101, -0.5,   +100.5);
  TH1D ht("ht", NULL, 100, -1.0e3, +1.0e3);
 
 
  map<double, double> data;     // time-of-emission -> range of times-of-emission
 
  typedef JTreeScanner<JEvent, JEvent::JEvaluator>  JTreeScanner_t;
 
  JTreeScanner_t in(acousticsFile);
 
  for (JTreeScanner_t::iterator p = in.begin(), q; p != in.end(); p = q) {
 
    STATUS("event " << FIXED(20,6) << p->begin()->getToE() << '\r'); DEBUG(endl);
 
    for (q = p; ++q != in.end() && q->begin()->getToE() <= p->rbegin()->getToE() + Tmax_s; ) {}
 
    h1.Fill((double) getNumberOfEmitters(p,q));
    hn.Fill((double) distance(p,q));
 
    if (getNumberOfEmitters(p,q) >= Nmin) {
 
      double range = 0.0;
 
      for (JTreeScanner_t::iterator i = p; i != q; ++i) {
        if (i->rbegin()->getToE() - p->begin()->getToE() > range) {
          range = i->rbegin()->getToE() - p->begin()->getToE();
        }
      }
 
      data[p->begin()->getToE()] = range;
    }
  }
  STATUS(endl);
 
 
  while (inputFile.hasNext()) {
 
    const JDAQEvent* evt = inputFile.next();
    const double     t0  = getUNIXTime(*evt);
 
    STATUS("DAQ event " << FIXED(20,6) << t0 << '\r'); DEBUG(endl); 
 
    map<double, double>::const_iterator p = data.lower_bound(t0);
 
    double t1 = (p != data.end() ? p->first : 0.0);
 
    if (p != data.begin()) {
 
      --p;
 
      if (fabs(t1 - t0) > fabs(p->first             - t0)) { t1 = p->first; }
      if (fabs(t1 - t0) > fabs(p->first + p->second - t0)) { t1 = p->first + p->second; }
    }
 
    ht.Fill(t0 - t1);
  }
  STATUS(endl);
 
 
  TFile out(outputFile.c_str(), "recreate");
 
  out << h1 << hn << ht;
 
  out.Write();
  out.Close();
}