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 "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 "JLang/JComparator.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"


/**
 * \file
 *
 * 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.\n
 * If the number of coincident times of emission exceeds a preset minimum, 
 * the event is triggered and subsequently stored in the output file.\n
 * 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.\n
 * In this, a super event refers to the coincidence of events from a variety of emitters. 
 * \author mdejong
 */
int main(int argc, char **argv)
{
  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();
}