JAcousticsEventBuilder.cc

Go to the documentation of this file.

#include <iostream>
#include <iomanip>
#include <vector>
#include <map>
#include <algorithm>

#include "TROOT.h"
#include "TFile.h"

#include "JLang/JComparator.hh"

#include "JDB/JToAshort.hh"
#include "JDB/JSupport.hh"

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

#include "JLang/JPredicate.hh"
#include "JTools/JRange.hh"
#include "JTools/JQuantile.hh"
#include "JTools/JHashMap.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JFileRecorder.hh"

#include "JAcoustics/JTransceiver.hh"
#include "JAcoustics/JTransmission.hh"
#include "JAcoustics/JAcousticsToolkit.hh"
#include "JAcoustics/JAcousticsSupportkit.hh"
#include "JAcoustics/JTriggerParameters.hh"
#include "JAcoustics/JEvent.hh"
#include "JAcoustics/JEventOverlap.hh"
#include "JAcoustics/JTriggerOutput.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 trigger acoustic data.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  
  JMultipleFileScanner<JToAshort> inputFile;     // ROOT file (after JConvertDB -q toashort)
  JTriggerParameters     parameters;
  JFileRecorder<JEvent>  outputFile;
  JSoundVelocity         V = getSoundVelocity;   // default sound velocity
  string                 detectorFile;
  JContainer<JTripod>    tripods;                // tripods
  JContainer<JBase>      bases;                  // base modules
  double                 precision;
  int                    debug;

  try {

    JParser<> zap("Example program to trigger acoustic data.");
    
    zap['f'] = make_field(inputFile,    "output of JConvertDB -q toashort")  = JPARSER::initialised();
    zap['@'] = make_field(parameters,   "trigger parameters");
    zap['o'] = make_field(outputFile,   "output file")                       = "event.root";
    zap['V'] = make_field(V,            "sound velocity")                    = JPARSER::initialised();
    zap['a'] = make_field(detectorFile, "detector file");
    zap['T'] = make_field(tripods,      "tripod data");
    zap['B'] = make_field(bases,        "base module data")                  = JPARSER::initialised();
    zap['p'] = make_field(precision,    "precision time-of-arrival")         = 1.0e-6;
    zap['d'] = make_field(debug)                                             = 1;

    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, JReceiver>   receivers;
  JHashMap<int, JEmitter>    emitters;

  for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
    receivers[i->getID()] = JReceiver(i->getID(), i->getPosition(), getT0(*i) * 1.0e-9, PIEZO_DELAYTIME_US * 1.0e-6);
  }

  for (JContainer<JBase>::const_iterator i = bases.begin(); i != bases.end(); ++i) {
    receivers[i->getID()] = JReceiver(i->getID(), i->getPosition(), getT0(*i) * 1.0e-9, HYDROPHONE_DELAYTIME_US * 1.0e-6);
  }

  for (JContainer<JTripod>::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
    emitters[i->getID()]  = JEmitter(i->getID(), i->getUTMPosition() - detector.getUTMPosition());
  }


  outputFile.open();

  if (!outputFile.is_open()) {
    FATAL("Error opening file " << outputFile << endl);
  }


  // input data

  typedef vector<JTransmission> buffer_type;     // acoustic data type

  map<int, map< int, buffer_type > >  f1;        // emitter -> receiver -> data

  for (int counter = 0; inputFile.hasNext(); ++counter) {

    STATUS("counter: " << setw(8) << counter << '\r' << flush); DEBUG(endl);

    const JToAshort*   parameters = inputFile.next();
    const JTransceiver transceiver(emitters[parameters->EMITTERID], receivers[parameters->DOMID]);

    f1[transceiver.emitter.getID()][transceiver.receiver.getID()].push_back(transceiver.getTransmission(*parameters, V[detector.getUTMZ()]));
  }
  STATUS(endl);


  // filter similar hits

  for (map<int, map< int, buffer_type> >::iterator i = f1.begin(); i != f1.end(); ++i) {

    for (map< int, buffer_type>::iterator receiver = i->second.begin(); receiver != i->second.end(); ++receiver) {

      buffer_type& buffer = receiver->second;

      sort(buffer.begin(), buffer.end(), make_comparator(&JTransmission::getToA, JComparison::lt()));

      buffer.erase(unique(buffer.begin(), buffer.end(), JTransmission::compare(precision)), buffer.end());
    }
  }


  // output events

  int counter = 0;

  for (map<int, map< int, buffer_type> >::iterator i = f1.begin(); i != f1.end(); ++i) {

    buffer_type buffer;

    for (map< int, buffer_type>::iterator receiver = i->second.begin(); receiver != i->second.end(); ++receiver) {

      // selection based on quality

      buffer_type& f2 = receiver->second;

      const JQuantile Q1("quality", f2.begin(), f2.end(), &JTransmission::getQ, true);

      const double Qmin = Q1.getQuantile(parameters.Q);

      buffer_type::iterator __end = partition(f2.begin(), f2.end(), make_predicate(&JTransmission::getQ, Qmin, JComparison::gt()));

      copy(f2.begin(), __end, back_inserter(buffer));
    }


    // sort according time-of-emisson

    sort(buffer.begin(), buffer.end());


    JTriggerOutput data;

    for (buffer_type::const_iterator p = buffer.begin(); p != buffer.end(); ++p) {

      buffer_type::const_iterator q = p;

      // basic correlator

      while (++q != buffer.end() && q->getToE() - p->getToE() <= parameters.TMax_s) {}

      if (distance(p,q) >= parameters.numberOfHits) {

        ++counter;

        STATUS("event:   " << setw(8) << counter << '\r' << flush); DEBUG(endl);

        JEvent event(counter, i->first, p, q);

        DEBUG("event: " << endl << event);

        data.push_back(event);
      }
    }

    data.merge(JEventOverlap(parameters.TMax_s));

    for (JTriggerOutput::const_iterator event = data.begin(); event != data.end(); ++event) {
      outputFile.put(*event);
    }
  }
  STATUS(endl);

  outputFile.close();
}