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