JAcousticsTriggerProcessor.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/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/JTransmission.hh"
#include "JAcoustics/JReceiver.hh"
#include "JAcoustics/JSoundVelocity.hh"
#include "JAcoustics/JAcousticsToolkit.hh"
#include "JAcoustics/JAcousticsSupportkit.hh"
#include "JAcoustics/JTriggerParameters.hh"
#include "JAcoustics/JEvent.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"

#include "JTrigger/JMatch.hh"
#include "JTrigger/JAlgorithm.hh"


namespace JACOUSTICS {

  using JGEOMETRY3D::JPosition3D;
  using JTRIGGER::JMatch;
  using JLANG::JClonable;
  using JTOOLS::JRange;


  /**
   * Acoustic hit.
   */
  struct hit_type :
    public JPosition3D,
   public JTransmission
  {
    /**
     * Default constructor.
     */
    hit_type()
    {}

    
    /**
     * Constructor.
     *
     * \param  position         position
     * \param  transmission     transmission
     */
    hit_type(const JPosition3D& position, const JTransmission& transmission) :
      JPosition3D  (position),
      JTransmission(transmission)
    {}
  };


  /** 
   * 3D match criterion for acoustic signals.
   */
  class JMatch3D :
    public JClonable< JMatch<hit_type>, JMatch3D>
  {
  public:
    /**
     * Constructor.
     *
     * \param  V                sound velocity
     * \param  Tmax_s           maximal extra time [s]
     */
    JMatch3D(const JSoundVelocity& V, const double Tmax_s = 0.0) :
      V(V),
      Tmax_s(Tmax_s)
    {}
    

    /**
     * Match operator.
     *
     * \param  first    hit
     * \param  second   hit
     * \return          match result
     */
    virtual bool operator()(const hit_type& first, const hit_type& second) const override 
    {
      using namespace JPP;

      const double t1 = V.getTime(getDistance(first.getPosition(), second.getPosition()),
                                  first .getZ(),
                                  second.getZ());

      const double dt = fabs(first.getToA() - second.getToA());
        
      return dt <= t1  +  Tmax_s;
    }


    const JSoundVelocity V;
    const double         Tmax_s;
  };
  

  /**
   * Match of two events considering overlap in time.
   */
  struct JEventOverlap {

    typedef JRange<double> time_range;        //!< Type definition of time range.

    /**
     * Constructor.
     *
     * \param Tmax_s          maximal time difference between two consecutive events [s]
     */
    JEventOverlap(const double Tmax_s = 0.0) :
      range(-Tmax_s, +Tmax_s)
    {}


    /**
     * Match criterion.
     *
     * \param  first          first  event
     * \param  second         second event
     * \return                true if two events overlap in time; else false
     */
    bool operator()(const JEvent& first, const JEvent& second) const
    {
      using namespace JPP;

      if (first .empty()) return false;
      if (second.empty()) return false;

      time_range r1(make_array(first .begin(), first .end(), &JTransmission::getToA));
      time_range r2(make_array(second.begin(), second.end(), &JTransmission::getToA));

      r1.add(range);
      r2.add(range);
      
      return overlap(r1, r2);
    }

    const time_range range;
  };
}


/**
 * \file
 *
 * Main program to trigger acoustic data.
 *
 * An acoustic event is based on coincidences between times of arrival.\n
 * If the number of coincident times of arrival exceeds a preset minimum, 
 * the event is triggered and subsequently stored in the output file.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  
  typedef JTYPELIST<JEvent, JTriggerParameters, JMeta>::typelist  typelist;

  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;
  hydrophones_container    hydrophones;            // hydrophones
  double                   precision;
  int                      ID;
  set<int>                 waveforms;
  int                      debug;

  try {

    JParser<> zap("Main program to trigger 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['H'] = make_field(hydrophones,  "hydrophone data")                   = JPARSER::initialised();
    zap['p'] = make_field(precision,    "precision time-of-arrival")         = 1.0e-6;
    zap['D'] = make_field(ID)                                                = 10000;
    zap['w'] = make_field(waveforms)                                         = JPARSER::initialised();
    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());


  const int           FACTORY_LIMIT  =  10000;
  const double        TMaxExtra_s    =  1.0e-3;

  const JMatch3D      match(V, TMaxExtra_s);
  const JEventOverlap overlap;                   // (parameters.TMax_s);

    
  JHashMap<int, JReceiver>   receivers;

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

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


  outputFile.open();

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

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


  // input data

  string oid;                                    // detector identifier

  typedef vector<hit_type> buffer_type;          // data type

  map<int, buffer_type>  f1;                     // 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);
    }

    if (waveforms.empty() || waveforms.count(parameters->EMITTERID)) {
      
      if (receivers.has(parameters->DOMID)) {

        const JReceiver&    receiver = receivers[parameters->DOMID];

        const JTransmission transmission(parameters->RUNNUMBER,
                                         parameters->DOMID,
                                         parameters->QUALITYFACTOR,
                                         parameters->QUALITYNORMALISATION,
                                         parameters->UNIXTIMEBASE  +  receiver.getT(parameters->TOA_S),
                                         parameters->UNIXTIMEBASE  +  receiver.getT(parameters->TOA_S));

        f1[receiver.getID()].push_back(hit_type(receiver.getPosition(), transmission));
      }
    }
  }
  STATUS(endl);


  // filter similar hits

  for (map< int, buffer_type>::iterator receiver = f1.begin(); receiver != f1.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());
  }


  buffer_type data;

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

    // selection based on quality

    buffer_type&    f2   = receiver->second;
    double          Qmin = parameters.Q;

    if (parameters.Q > 0.0 && 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(data));
  }


  sort(data.begin(), data.end(), make_comparator(&hit_type::getToA, JComparison::lt()));          // sort according time-of-arrival


  buffer_type buffer(FACTORY_LIMIT);               // local buffer of hits
  JEvent      out[2];                              // FIFO of events
  int         number_of_events = 0;
 
  for (buffer_type::const_iterator p = data.begin(); p != data.end(); ++p) {

    if (distance(data.cbegin(),p)%1000 == 0) {
      STATUS("processed: " << FIXED(5,1) << (double) (100 * distance(data.cbegin(),p)) / (double) data.size() << "%" << '\r' << flush); DEBUG(endl);
    }
    
    buffer_type::const_iterator q = p;

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

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

      if (distance(p,q) < FACTORY_LIMIT) {

        buffer_type::iterator root = buffer.begin();
        buffer_type::iterator __p  = buffer.begin();
        buffer_type::iterator __q  = buffer.begin();

        *root = *p;
        
        ++__p;
        ++__q;

        for (buffer_type::const_iterator i = p; ++i != q; ) {
          if (match(*p,*i)) {
            *__q = *i;
            ++__q;
          }
        }

        if (distance(root,__q) >= parameters.numberOfHits) { 
        
          __q = clusterize(__p, __q, match, parameters.numberOfHits - 1);

          if (distance(root,__q) >= parameters.numberOfHits) {
            out[1] = JEvent(oid, ++number_of_events, ID, root, __q);
          }
        }

      } else {

        out[1] = JEvent(oid, ++number_of_events, ID, p, q);
      }

      if (!overlap(out[0],out[1])) {

        if (!out[0].empty()) {
          outputFile.put(out[0]);                  // write
        }

        out[0] = out[1];                           // shift
              
      } else {
              
        out[0].merge(out[1]);                      // merge
      }

      out[1].clear();
    }
  }

  if (!out[0].empty()) {
    outputFile.put(out[0]);                        // write
  }
  STATUS(endl);

  STATUS("triggers: " << setw(7) << number_of_events << endl);
  
  JMultipleFileScanner<JMeta> io(inputFile);

  io >> outputFile;

  outputFile.close();
}