JAcoustics/JToAshortMonitor.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 "TH2D.h"

#include "JTools/JHashCollection.hh"
#include "JTools/JHashMap.hh"
#include "JTools/JRange.hh"
#include "JTools/JAbstractHistogram.hh"

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"

#include "JSupport/JMultipleFileScanner.hh"

#include "JLang/JComparator.hh"
#include "JLang/JComparison.hh"

#include "JROOT/JManager.hh"

#include "JAcoustics/JEvent.hh"
#include "JAcoustics/JSupport.hh"

#include "Jeep/JContainer.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


/**
 * \file
 *
 * Example program to monitor acoustic events.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  typedef JRange<double>              JRange_t;
  typedef JAbstractHistogram<double>  JHistogram_t;

  JMultipleFileScanner<JEvent>  inputFile;
  JLimit_t&               numberOfEvents = inputFile.getLimit();
  string                  outputFile;
  string                  detectorFile;
  JRange_t                T_s;
  JHistogram_t            x;
  int                     id;                                   // emitter identifier
  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['o'] = make_field(outputFile)       = "toashort.root";
    zap['a'] = make_field(detectorFile);
    zap['T'] = make_field(T_s)              = JRange_t(0.0, 10.0);
    zap['x'] = make_field(x,               "histogram x-binning")            = JHistogram_t(100000, -2.0e-2, +2.0e-2);
    zap['E'] = make_field(id,              "emitter identifier (-1 = all)")  = -1;
    zap['d'] = make_field(debug)            = 2;

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


  JDetector detector;

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

  JHashMap<int, JLocation>   receivers;

  for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
    receivers[i->getID()] = i->getLocation();
  }


  JManager<JLocation, TH1D> H1(new TH1D("[%].t1",  NULL,  x.getNumberOfBins(), x.getLowerLimit(), x.getUpperLimit()));


  map<int, vector<JEvent> > zbuf;

  while (inputFile.hasNext()) {

    STATUS("input  " << setw(6) << inputFile.getCounter() << '\r'); DEBUG(endl);

    const JEvent* evt = inputFile.next();

    zbuf[evt->getID()].push_back(*evt);
  }
  STATUS(endl);


  for (map<int, vector<JEvent> >::iterator emitter = zbuf.begin(); emitter != zbuf.end(); ++emitter) {

    STATUS("Emitter " << setw(2) << emitter->first << ' ' << setw(6) << emitter->second.size() << endl);

    if (emitter->first == id || id == -1) {
    
      if (emitter->second.size() > 1) {

        sort(emitter->second.begin(), emitter->second.end());                         // sort by time

        for (vector<JEvent>::iterator q = emitter->second.begin(), p = q++; q != emitter->second.end(); ++p, ++q) {

          if (!p->empty() && !q->empty()) {

            DEBUG("Time difference " << FIXED(7,3) << q->begin()->getToE() - p->begin()->getToE() << endl);

            if (T_s(q->begin()->getToE() - p->begin()->getToE())) {

              sort(p->begin(), p->end(), make_comparator(&JTransmission::getID));     // sort by identifier
              sort(q->begin(), q->end(), make_comparator(&JTransmission::getID));     // sort by identifier

              map<JLocation, double> buffer;

              JEvent::const_iterator __p = p->begin();
              JEvent::const_iterator __q = q->begin();

              while (__p != p->end() && __q != q->end()) {

                while (__p != p->end() && __q != q->end() && __p->getID() < __q->getID()) { ++__p; }
                while (__p != p->end() && __q != q->end() && __q->getID() < __p->getID()) { ++__q; }

                if (__p != p->end() && __q != q->end() && __q->getID() == __p->getID()) {

                  if (receivers.has(__p->getID())) {

                    const JLocation& location = receivers[__p->getID()];
                    const double     t1       = __q->getToA() - __p->getToA();

                    buffer[location] = t1;
                  }

                  ++__p;
                  ++__q;
                }
              }

              if (buffer.size() > 1) {

                for (map<JLocation, double>::const_iterator q = buffer.begin(), p = q++; q != buffer.end(); ++p, ++q) {

                  if (p->first.getString() == q->first.getString() && p->first.getFloor() + 1 == q->first.getFloor()) {

                    const double t1 = q->second - p->second;

                    H1[p->first]->Fill(t1);
                    H1          ->Fill(t1);
                  }
                }
              }
            }
          }
        }
      }
    }
  }


  TFile out(outputFile.c_str(), "recreate");

  out << H1 << *H1;

  out.Write();
  out.Close();
}