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 "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleSupportkit.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/JToA.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 "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

main()

int main	(	int	argc,
		char **	argv )

Definition at line 61 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<JToA> 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 JToA");
    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 = getPiezoPosition();
 
      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_misses;
  map_type number_of_triggers;
  map_type number_of_events;
 
  // input data
 
  typedef vector<JTransmission> buffer_type;     // data type
 
  map<int, map< int, buffer_type > >  f1;        // emitter -> receiver -> data
 
  while (inputFile.hasNext()) {
 
    if (inputFile.getCounter()%10000 == 0) {
      STATUS("counter: " << setw(8) << inputFile.getCounter() << '\r' << flush); DEBUG(endl);
    }
 
    JToA* parameters = inputFile.next();
 
    if (detector.getID() != parameters->DETID) { // consistency check
      FATAL("Invalid detector identifier " << parameters->DETID << " != " << detector.getID() << endl);
    }
 
    try {
 
      const int id = getEmitterID(parameters->WAVEFORMID);
 
      number_of_entries[parameters->WAVEFORMID].put(1);
 
      if (emitters.has(id)) {
 
        if (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));
        }
        
      } else {
        number_of_misses[id].put(1);
      }
    }
    catch(const exception&) {
      number_of_aliens[parameters->WAVEFORMID].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);
  }
 
  for (map_type::const_iterator i = number_of_misses.begin();  i != number_of_misses.end();  ++i) {
    STATUS("number of misses:   " << setw(3) << i->first << ' ' << setw(8) << i->second.getCount() << endl);
  }
 
  const JEventOverlap match(parameters.TMax_s);
 
  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) {
 
      // filter similar hits
 
      sort(receiver->second.begin(), receiver->second.end(), JTransmission::compare(precision));
 
      buffer_type::iterator __end = unique(receiver->second.begin(), receiver->second.end(), JTransmission::equals(precision));
 
      // selection based on quality
 
      for (buffer_type::const_iterator p = receiver->second.begin(); p != __end; ++p) {
        if (p->getQ() >= parameters.Q * (parameters.Q <= 1.0 ? p->getW() : 1.0)) {
          buffer.push_back(*p);
        }
      }
    }
 
    sort(buffer.begin(), buffer.end());          // sort according time-of-emisson
 
    JQuantile      Q;                            // monitor of trigger
    JEvent         out;                          // temporary storage of event(s)
 
    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) {
 
        JEvent event(detector.getID(), number_of_triggers.sum(), i->first, p, q);
        
        number_of_triggers[i->first].put(event.size());
 
        DEBUG("trigger: " << endl << event);
 
        if        (out.empty()) {
 
          out = event;
 
        } else if (match(out, event)) {
 
          out.merge(event);
 
        } else {
 
          outputFile.put(out);
 
          number_of_events[out.getID()].put(out.size());
 
          out = event;
        }
      }
    }
 
    if (!out.empty()) {
 
      outputFile.put(out);                       // pending event
 
      number_of_events[out.getID()].put(out.size());
    }
 
    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);
  }
  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() << endl);
  }
 
  JMultipleFileScanner<JMeta> io(inputFile);
 
  io >> outputFile;
 
  outputFile.close();
 
  return 0;
}