JAcousticsEventMonitor.cc

Go to the documentation of this file.

#include <iostream>
#include <iomanip>
 
#include "JAcoustics/JEvent.hh"
#include "JAcoustics/JSupport.hh"
 
#include "JTools/JQuantile.hh"
 
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JTreeScanner.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
namespace {
 
  using JACOUSTICS::JTransmission;
 
  /**
   * Auxiliary data structure to compare transmissions.
   */
  struct JCompare {
    /**
     * Compare transmissions.
     *
     * \param  first             first  transmission
     * \param  second            second transmission
     * \return                   true if first transmission before second; else false
     */
    inline bool operator()(const JTransmission& first, const JTransmission& second) const
    {
      if (first.getID() == second.getID())
        if (fabs(first.getToA() - second.getToA()) <= waveform_s)
          return false;
        else
          return first.getToA() < second.getToA();
      else
        return first.getID()  < second.getID();
    }
 
    double waveform_s;
    
  } compare;
}
 
 
/**
 * \file
 *
 * Example program to monitor acoustic events.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  JMultipleFileScanner<>  inputFile;
  JLimit_t&               numberOfEvents = inputFile.getLimit();
  double                  Tmax_s;
  int                     debug;
 
  try {
 
    JParser<> zap("Example program to monitor acoustic events.");
 
    zap['f'] = make_field(inputFile,           "output of JAcousticEventBuilder[.sh]");
    zap['n'] = make_field(numberOfEvents)       = JLimit::max();
    zap['p'] = make_field(compare.waveform_s,  "wavefrom duration [s]")                =  5.0e-3;
    zap['D'] = make_field(Tmax_s,              "time difference between events [s]")   = 20.0e-3;
    zap['d'] = make_field(debug)                = 2;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  typedef JTreeScanner<JEvent, JEvent::JEvaluator>  JTreeScanner_t;
 
  JEvent* p0 = new JEvent();
 
  counter_type events = 0;
 
  JQuantile Q;
  
  for (JMultipleFileScanner<>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
 
    JTreeScanner_t in(*i);
 
    for (JTreeScanner_t::iterator p1 = in.begin(); p1 != in.end() && events != inputFile.getLimit(); ++p1, ++events) {
 
      STATUS("event " << setw(8) << events << '\r'); DEBUG(endl);
 
      if (!p1->empty() && !p0->empty()) {
 
        if (p1->begin()->getToE() < p0->rbegin()->getToE() + Tmax_s) {
 
          if (debug >= debug_t) {
 
            cout << endl;
 
            cout << p0->getDetectorID()                   << endl
                 << FIXED(20,6) << p0-> begin()->getToE() << ' '
                 << FIXED(20,6) << p0->rbegin()->getToE() << endl
                 << FIXED(20,6) << p1-> begin()->getToE() << ' '
                 << FIXED(20,6) << p1->rbegin()->getToE() << endl;
          }
 
          JEvent evt[] = { *p0, *p1 };
 
          for (int i = 0; i != 2; ++i) {
            sort(evt[i].begin(), evt[i].end(), compare);
          }
 
          size_t n = 0;
          
          for (JEvent::const_iterator
                 i0 = evt[0].begin(),
                 i1 = evt[1].begin(); i0 != evt[0].end() && i1 != evt[1].end(); ) {
 
            if      (compare(*i0, *i1))
              ++i0;
            else if (compare(*i1, *i0))
              ++i1;
              
            else {
 
              n += 1;
              
              if (debug >= debug_t) {
 
                cout << setw(2)     << p0->getID()         << ' '
                     << setw(10)    << i0->getID()         << ' '
                     << FIXED(20,6) << i0->getToA()        << ' '
                     << FIXED( 8,0) << i0->getQ()          << ' ';
 
                cout << setw(2)     << p1->getID()         << ' '
                     << setw(10)    << i1->getID()         << ' '
                     << FIXED(20,6) << i1->getToA()        << ' '
                     << FIXED( 8,0) << i1->getQ()          << ' ';
 
                cout << FIXED( 9,6) << (i1->getToA() - i0->getToA()) << endl;
              }
 
              if      (i0->getToA() < i1->getToA())
                ++i0;
              else if (i1->getToA() < i0->getToA())
                ++i1;
              else {
                ++i0;
                ++i1;
              }
            }
          }
 
          if (n != 0) {
            Q.put(n);
          }
        }
      }
      
      new (p0) JEvent(*p1);
    }
  }
  STATUS(endl);
 
  delete p0;
 
  cout << "Number of errors / events: " << Q.getCount() << " / " << events << endl;
 
  if (Q.getCount() > 0) {
    Q.print(cout);
  }
}