JAcousticsEventBuilder.cc File Reference

Main program to trigger events in acoustic data. More...

#include <iostream>
#include <iomanip>
#include <vector>
#include <map>
#include <algorithm>
#include "TROOT.h"
#include "TFile.h"
#include "km3net-dataformat/definitions/module_status.hh"
#include "JLang/JComparator.hh"
#include "JDB/JToAshort.hh"
#include "JDB/JSupport.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JDetectorSupportkit.hh"
#include "JDetector/JTripod.hh"
#include "JDetector/JTransmitter.hh"
#include "JDetector/JHydrophone.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 "JSupport/JMeta.hh"
#include "JAcoustics/JEmitterID.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"

Go to the source code of this file.

Functions
int	main (int argc, char **argv)

Detailed Description

Main program to trigger events in acoustic data.

An acoustic event is based on a coincidence of estimated times of emission from the same emitter.
If the number of coincident times of emission exceeds a preset minimum, the event is triggered and subsequently stored in the output file.
Note that an event counter is added which can be used to disambiguate the different cycles from the same emitter within one so-called "super" event that is used for the global fit.
In this, a super event refers to the coincidence of events from a variety of emitters.

Author

mdejong

Definition in file JAcousticsEventBuilder.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 63 of file JAcousticsEventBuilder.cc.

{
  using namespace std;
  using namespace JPP;
  
  typedef JTYPELIST<JEvent, JTriggerParameters, JMeta>::typelist  typelist;

  typedef JContainer< vector<JTripod> >      tripods_container;
  typedef JContainer< vector<JTransmitter> > transmitters_container;
  typedef JContainer< vector<JHydrophone> >  hydrophones_container;

  JMultipleFileScanner<JToAshort> inputFile;
  JLimit_t&                numberOfEvents = inputFile.getLimit();
  JTriggerParameters       parameters;
  JFileRecorder<typelist>  outputFile;
  JSoundVelocity           V = getSoundVelocity;   // default sound velocity
  string                   detectorFile;
  tripods_container        tripods;                // tripods
  transmitters_container   transmitters;           // transmitters
  hydrophones_container    hydrophones;            // hydrophones
  double                   precision;
  int                      debug;

  try {

    JParser<> zap("Main program to trigger events in acoustic data.");
    
    zap['f'] = make_field(inputFile,    "output of JConvertDB -q toashort");
    zap['n'] = make_field(numberOfEvents)                                    = JLimit::max();    
    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['Y'] = make_field(transmitters,    "transmitter data")               = JPARSER::initialised();
    zap['H'] = make_field(hydrophones,  "hydrophone data")                   = JPARSER::initialised();
    zap['W'] = make_field(getEmitterID, "waveform identification 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);
  }

  V.set(detector.getUTMZ());

  JHashMap<int, JReceiver>   receivers;
  JHashMap<int, JEmitter>    emitters;

  for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {

    if ((i->getFloor() == 0 && !i->has(HYDROPHONE_DISABLE))  ||
        (i->getFloor() != 0 && !i->has(PIEZO_DISABLE))) {

      JPosition3D pos(0.0, 0.0, 0.0);

      if (i->getFloor() == 0) {                    // get relative position of hydrophone

        try {
          pos = getPosition(hydrophones.begin(),
                            hydrophones.end(),
                            make_predicate(&JHydrophone::getString, i->getString()));
        }
        catch(const exception&) {
          continue;                                // if no position available, discard hydrophone
        }
      }

      receivers[i->getID()] = JReceiver(i->getID(), 
                                        i->getPosition() + pos, 
                                        i->getT0() * 1.0e-9);
    }
  }

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

  for (transmitters_container::const_iterator i = transmitters.begin(); i != transmitters.end(); ++i) {
    try {
      emitters[i->getID()]  = JEmitter(i->getID(),
                                       i->getPosition() + detector.getModule(i->getLocation()).getPosition());
    }
    catch(const exception&) {
      continue;                                // if no string available, discard transmitter
    }
  }


  outputFile.open();

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

  outputFile.put(JMeta(argc, argv));
  outputFile.put(parameters);

  // statistics

  struct map_type :
   public map<int, JQuantile>
  {
    long long int sum() const 
    {
      long long int count = 0;

      for (const_iterator i = this->begin(); i != this->end(); ++i) {
        count += i->second.getCount();
      }

      return count;
    }
  };

  map_type number_of_entries;
  map_type number_of_aliens;
  map_type number_of_triggers;
  map_type number_of_events;
  map_type number_of_outliers;

  // input data

  string oid;                                    // detector identifier

  typedef vector<JTransmission> buffer_type;     // data type

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

  while (inputFile.hasNext()) {

    if (inputFile.getCounter()%1000 == 0) {
      STATUS("counter: " << setw(8) << inputFile.getCounter() << '\r' << flush); DEBUG(endl);
    }

    JToAshort* parameters = inputFile.next();

    if (oid == "") {
      oid = parameters->DETID;
    }

    if (oid != parameters->DETID) {              // consistency check
      FATAL("Invalid detector identifier " << parameters->DETID << " != " << oid << endl);
    }

    try {

      const int id = getEmitterID(parameters->EMITTERID);

      number_of_entries[parameters->EMITTERID].put(1);

      if (emitters.has(id) && receivers.has(parameters->DOMID)) {

        const JTransceiver transceiver(emitters[id], receivers[parameters->DOMID]);

        f1[transceiver.emitter.getID()][transceiver.receiver.getID()].push_back(transceiver.getTransmission(*parameters, V));
      }
    }
    catch(const exception&) {
      number_of_aliens[parameters->EMITTERID].put(1);
    }    
  }
  STATUS(endl);

  for (map_type::const_iterator i = number_of_entries.begin();  i != number_of_entries.end();  ++i) {
    STATUS("number of entries:  " << setw(3) << i->first << ' ' << setw(8) << i->second.getCount() << endl);
  }

  for (map_type::const_iterator i = number_of_aliens.begin();  i != number_of_aliens.end();  ++i) {
    STATUS("number of aliens:   " << setw(3) << i->first << ' ' << setw(8) << i->second.getCount() << 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(), JTransmission::compare(precision));

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


  // output events

  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;
      double          Qmin = parameters.Q;

      if (parameters.Q > 0.0 &&                  // use quantile to evaluate corresponding quality
          parameters.Q < 1.0) {                  // 

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

        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(buffer.begin(), buffer.end());          // sort according time-of-emisson

    JQuantile Q;

    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) {}

      Q.put(distance(p,q));

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

        STATUS("trigger: " << setw(8) << number_of_triggers.sum() << '\r' << flush); DEBUG(endl);

        JEvent event(oid, number_of_triggers.sum(), i->first, p, q);
        
        number_of_triggers[i->first].put(event.size());

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

        data.push_back(event);
      }
    }

    STATUS("number of toes:     " << setw(3) << i->first << ' ' << setw(8) << buffer.size());
    if (Q.getCount() > 0) {
      STATUS(' ' << FIXED(6,1) << Q.getMean() << 
             ' ' << FIXED(6,1) << Q.getXmin() << 
             ' ' << FIXED(6,1) << Q.getXmax());
    }
    STATUS(endl);

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

    for (JTriggerOutput::const_iterator event = data.begin(); event != data.end(); ++event) {

      number_of_events[event->getID()].put(event->size());

      if (parameters.quantile < 1.0) {

        const int M = (int) (event->size() * 0.5 * (1.0 - parameters.quantile));

        JEvent::const_iterator         p = event-> begin();  advance(p, M);
        JEvent::const_reverse_iterator q = event->rbegin();  advance(q, M);

        const double Tmin = p->getToE();
        const double Tmax = q->getToE();

        const double W = (1.0 - parameters.quantile) / parameters.quantile;
        
        while (p != event-> begin() && (p - 1)->getToE() >= Tmin - 0.5 * (Tmax - Tmin) * W) { --p; }
        while (q != event->rbegin() && (q - 1)->getToE() <= Tmax + 0.5 * (Tmax - Tmin) * W) { --q; }
        
        number_of_outliers[event->getID()].put(event->size() - distance(p, q.base()));

        outputFile.put(JEvent(event->getOID(),
                              event->getCounter(),
                              event->getID(),
                              p, q.base()));

      } else {
      
        number_of_outliers[event->getID()].put(0);

        outputFile.put(*event);
      }
    }
  }
  STATUS(endl);

  for (map_type::const_iterator i = number_of_triggers.begin(); i != number_of_triggers.end(); ++i) {
    STATUS("number of triggers: " << setw(3) << i->first << ' ' << setw(8) << i->second.getCount() << endl);
  }

  for (map_type::const_iterator i = number_of_events.begin();   i != number_of_events.end();   ++i) {
    STATUS("number of events:   " << setw(3) << i->first << ' ' << setw(8) << i->second.getCount() << ' ' << FIXED(7,3) << i->second.getMean() << ' ' << FIXED(7,3) << number_of_outliers[i->first].getMean() << endl);
  }

  JMultipleFileScanner<JMeta> io(inputFile);

  io >> outputFile;

  outputFile.close();
}