JAcousticsTriggerProcessor.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
#include <map>
#include <algorithm>
 
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TGraph.h"
 
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JHydrophone.hh"
 
#include "JLang/JClonable.hh"
#include "JLang/JPredicate.hh"
#include "JLang/JComparator.hh"
 
#include "JTools/JRange.hh"
#include "JTools/JHashMap.hh"
 
#include "JROOT/JManager.hh"
#include "JROOT/JGraph.hh"
 
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMeta.hh"
 
#include "JAcoustics/JToA.hh"
#include "JAcoustics/JTransmission.hh"
#include "JAcoustics/JReceiver.hh"
#include "JAcoustics/JSoundVelocity.hh"
#include "JAcoustics/JConstantSoundVelocity.hh"
#include "JAcoustics/JAcousticsToolkit.hh"
#include "JAcoustics/JAcousticsSupportkit.hh"
#include "JAcoustics/JTriggerParameters.hh"
#include "JAcoustics/JEvent.hh"
#include "JAcoustics/JSupport.hh"
 
#include "JGeometry3D/JVector3D.hh"
#include "JGeometry3D/JVertex3D.hh"
#include "JGeometry3D/JPosition3D.hh"
#include "JGeometry3D/JCylinder3D.hh"
 
#include "Jeep/JContainer.hh"
#include "Jeep/JProperties.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::JVector3D;
  using JGEOMETRY3D::JVertex3D;
  using JGEOMETRY3D::JPosition3D;
  using JTRIGGER::JMatch;
  using JLANG::JClonable;
  using JTOOLS::JRange;
  using JACOUSTICS::JEvent;
 
 
  /**
   * 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 JAbstractSoundVelocity& 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 JAbstractSoundVelocity& V;
    const double Tmax_s;
  };
 
 
  /**
   * Event.
   */
  struct event_type :
    public JEvent,
    public JVertex3D
  {
    /**
     * Default constructor.
     */
    event_type() :
      JEvent(),
      JVertex3D()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  event          event
     * \param  vertex         vertex
     */
    event_type(const JEvent& event, const JVertex3D& vertex = JVertex3D()) :
      JEvent(event),
      JVertex3D(vertex)
    {}
 
 
    /**
     * Get event.
     *
     * \return                event
     */
    const JEvent& getEvent() const
    {
      return static_cast<const JEvent&>(*this);
    }
 
 
    /**
     * Get vertex.
     *
     * \return                vertex
     */
    const JVertex3D& getVertex() const
    {
      return static_cast<const JVertex3D&>(*this);
    }
  };
  
 
  /**
   * Match of two events considering overlap in time and position.
   */
  struct JEventOverlap {
    /**
     * Constructor.
     *
     * \param Tmax_s          maximal time difference between two events [s]
     * \param Dmax_m          maximal distance between two events [m]
     */
    JEventOverlap(const double Tmax_s, const double Dmax_m = std::numeric_limits<double>::max()) :
      Tmax_s(Tmax_s),
      Dmax_m(Dmax_m)
    {}
 
 
    /**
     * Match criterion.
     *
     * \param  first          first  event
     * \param  second         second event
     * \return                true if two events overlap; else false
     */
    bool operator()(const event_type& first, const event_type& second) const
    {
      using namespace JPP;
 
      return fabs(first.getT() - second.getT()) <= Tmax_s && getDistance(first.getPosition(), second.getPosition()) <= Dmax_m;
    }
 
    const double Tmax_s;
    const double Dmax_m;
  };
 
 
  /**
   * Vertex.
   */
  struct vertex_type :
    public JVertex3D
  {
    /**
     * Default constructor.
     */
    vertex_type() :
      JVertex3D(),
      N(0),
      Q(0.0)
    {}
 
 
    /**
     * Constructor.
     *
     * \param  p0               position
     * \param  t0               time
     * \param  n0               number of hits
     * \param  q0               total quality
     */
    vertex_type(const JVector3D& p0, const double t0, const int n0, const double q0) :
      JVertex3D(p0, t0),
      N(n0),
      Q(q0)
    {}
 
 
    /**
     * Get vertex.
     *
     * \return                vertex
     */
    const JVertex3D& getVertex() const
    {
      return static_cast<const JVertex3D&>(*this);
    }
 
    int    N;                   //!< number of hits
    double Q;                   //!< total quality
  };
 
 
  /**
   * Vertex locator.
   */
  struct JVelo :
    public std::vector<JPosition3D>
  {
    /**
     * Constructor.
     *
     * \param  V                sound velocity
     * \param  center           center for generation of vertices
     * \param  RMax_m           radius for generation of vertices [m]
     * \param  Xv_m             step size for generation of vertices [m]
     * \param  factor           multiplication factor for corresponding time window
     */
    JVelo(const JAbstractSoundVelocity& V,
          const JVector3D& center,
          const double     RMax_m,
          const double     Xv_m,
          const double     factor = 1.0) :
      V(V),
      Tmax_s(factor * Xv_m / V())
    {
      this->push_back(center);
 
      if (RMax_m > 0.0 && Xv_m > 0.0) {
        for (double x = 0.5*Xv_m; x <= RMax_m; x += Xv_m) {
          for (double y = 0.5*Xv_m; y <= RMax_m; y += Xv_m) {
            if (x*x + y*y <= RMax_m*RMax_m) {
              this->push_back(JVector3D(center.getX() - x, center.getY() + y, center.getZ()));
              this->push_back(JVector3D(center.getX() + x, center.getY() + y, center.getZ()));
              this->push_back(JVector3D(center.getX() - x, center.getY() - y, center.getZ()));
              this->push_back(JVector3D(center.getX() + x, center.getY() - y, center.getZ()));
            }
          }
        }
      }
    }
 
 
    /**
     * Check vertex.
     *
     * \param  vx               vertex
     * \param  __begin          begin of data
     * \param  __end            end   of data
     * \return                  number of hits
     */
    template<class T>
    inline int operator()(const JVertex3D& vx, T __begin, T __end) const
    {
      const JPosition3D& p0 = vx.getPosition();
      const double       t0 = vx.getT();
 
      int    n0 = 0;
      double q0 = 0.0;
 
      for (T p = __begin; p != __end; ++p) {
 
        const double t1 = p->getToA() - V.getTime(p0.getDistance(p->getPosition()), p0.getZ(), p->getZ());
 
        if (fabs(t1 - t0) <= Tmax_s) {
          n0 += 1;
          q0 += p->getQ();
        }
      }
 
      return n0;
    }
 
 
    /**
     * Update vertex at given position.
     *
     * \param  position         position
     * \param  root             root hit
     * \param  __begin          begin of data
     * \param  __end            end   of data
     * \return                  vertex
     */
    template<class T>
    inline vertex_type operator()(const JPosition3D& position, const hit_type& root, T __begin, T __end) const
    {
      const JPosition3D& p0 = position;
      const double       t0 = root.getToA() - V.getTime(p0.getDistance(root.getPosition()), p0.getZ(), root.getZ());
 
      int    n0 = 1;
      double q0 = root.getQ();
 
      for (T p = __begin; p != __end; ++p) {
 
        const double t1 = p->getToA() - V.getTime(p0.getDistance(p->getPosition()), p0.getZ(), p->getZ());
 
        if (fabs(t1 - t0) <= Tmax_s) {
          n0 += 1;
          q0 += p->getQ();
        }
      }
 
      return vertex_type(p0, t0, n0, q0);
    }
 
 
    /**
     * Locate vertex.
     *
     * \param  root             root hit
     * \param  __begin          begin of data
     * \param  __end            end   of data
     * \param  numberOfHits     number of hits
     * \return                  vertex
     */
    template<class T>
    inline vertex_type operator()(const hit_type& root, T __begin, T __end, const int numberOfHits = 0) const
    {
      vertex_type vertex;
 
      for (const_iterator i = this->cbegin(); i != this->cend(); ++i) {
 
        const JPosition3D& p0 = *i;
        const double       t0 = root.getToA() - V.getTime(p0.getDistance(root.getPosition()), p0.getZ(), root.getZ());
 
        int    n0 = 1;
        double q0 = root.getQ();
 
        for (T p = __begin; p != __end && n0 + distance(p, __end) >= numberOfHits && n0 + distance(p, __end) >= vertex.N; ++p) {
 
          const double t1 = p->getToA() - V.getTime(p0.getDistance(p->getPosition()), p0.getZ(), p->getZ());
 
          if (fabs(t1 - t0) <= Tmax_s) {
            n0 += 1;
            q0 += p->getQ();
          }
        }
 
        if (q0 > vertex.Q) {
          vertex = vertex_type(p0, t0, n0, q0);
        }
      }
 
      return vertex;
    }
 
    const JAbstractSoundVelocity& V;
    const double Tmax_s;
  };
 
 
  /**
   * Get average quality.
   *
   * \param  evt             event
   * \return                 quality
   */
  inline double getQuality(const JEvent& evt)
  {
    double Q = 0.0;
 
    if (!evt.empty()) {
 
      for (const auto& i : evt) {
        Q += i.getQ();
      }
 
      Q /= evt.size();
    }
 
    return Q;
  }
}
 
 
/**
 * \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, TH1D, TGraph>::typelist  typelist;
 
  typedef JContainer< vector<JHydrophone> >  hydrophones_container;
 
  JMultipleFileScanner<JToA> inputFile;
  JLimit_t&                numberOfEvents = inputFile.getLimit();
  JTriggerParameters       parameters;
  int                      factoryLimit   = 10000;
  double                   TMaxExtra_s    = 500.0e-6;
  double                   RMax_m         = 0.0;
  double                   DMax_m         = 0.0;
  double                   Xv_m           = 0.0;
  double                   factor         = 2.0;
  double                   Z_m            = 0.0;   // z-position of emitters
  double                   fraction       = 0.75;  // fraction of working modules
  JFileRecorder<typelist>  outputFile;
  JSoundVelocity           V = getSoundVelocity;   // default sound velocity
  string                   detectorFile;
  hydrophones_container    hydrophones;            // hydrophones
  double                   precision;
  set<int>                 ID;
  int                      debug;
 
  try {
 
    JProperties properties;
 
    properties.insert(gmake_property(parameters.Q));
    properties.insert(gmake_property(parameters.TMax_s));
    properties.insert(gmake_property(parameters.numberOfHits));
    properties.insert(gmake_property(fraction));
    properties.insert(gmake_property(factoryLimit));
    properties.insert(gmake_property(TMaxExtra_s));
    properties.insert(gmake_property(RMax_m));
    properties.insert(gmake_property(DMax_m));
    properties.insert(gmake_property(Xv_m));
    properties.insert(gmake_property(factor));
    properties.insert(gmake_property(Z_m));
    
    JParser<> zap("Main program to trigger acoustic data.");
    
    zap['f'] = make_field(inputFile,    "output of JToA");
    zap['n'] = make_field(numberOfEvents)                                    = JLimit::max();    
    zap['@'] = make_field(properties,   "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['W'] = make_field(ID)                                                = 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 JCylinder3D cylinder(detector.begin(), detector.end());
 
  const JConstantSoundVelocity V0 = V(cylinder.getZmax());
 
  const JVelo velo(V0, JVector3D(cylinder.getX(), cylinder.getY(), Z_m), RMax_m, Xv_m, factor);
 
  const JMatch3D      match(V0, TMaxExtra_s);
  const JEventOverlap overlap2D(parameters.TMax_s, DMax_m);
  const JEventOverlap overlap1D(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);
  }
 
  // monitoring
  
  JManager<int, TH1D> M1(new TH1D("M1[%]", NULL, 2001, -0.5, 2000.5));
  JManager<int, TH1D> M2(new TH1D("M2[%]", NULL, 2001, -0.5, 2000.5));
  JManager<int, TH1D> M3(new TH1D("M3[%]", NULL, 2001, -0.5, 2000.5));
  JManager<int, TH1D> M4(new TH1D("M4[%]", NULL, 2001, -0.5, 2000.5));
  JManager<int, TH1D> M5(new TH1D("M5[%]", NULL, 2001, -0.5, 2000.5));
 
  map<int, JGraph_t>  G2;
 
  outputFile.open();
 
  if (!outputFile.is_open()) {
    FATAL("Error opening file " << outputFile << endl);
  }
 
  outputFile.put(JMeta(argc, argv));
  outputFile.put(parameters);
 
 
  // input data
 
  typedef vector<hit_type> buffer_type;          // data type
 
  map<int, map<int, buffer_type> > f1;           // emitter -> receiver -> data
 
  set<int> ids;                                  // working receivers
 
  static double TOA_s = numeric_limits<double>::max();
 
  while (inputFile.hasNext()) {
 
    if (inputFile.getCounter()%100000 == 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);
    }
 
    if (receivers.has(parameters->DOMID)) {
 
      ids.insert(parameters->DOMID);
 
      const JReceiver& receiver = receivers[parameters->DOMID];
      const double     toa_s    = parameters->TOA_S();
 
      if (toa_s < TOA_s) {
        TOA_s = toa_s;
      }
 
      const JTransmission transmission(parameters->RUN,
                                       parameters->DOMID,
                                       parameters->QUALITYFACTOR,
                                       parameters->QUALITYNORMALISATION,
                                       receiver.getT(toa_s),
                                       receiver.getT(toa_s));
 
      if (ID.empty() || ID.count(getWaveformID(parameters->WAVEFORMID)) != 0) {
        f1[getWaveformID(parameters->WAVEFORMID)][receiver.getID()].push_back(hit_type(receiver.getPosition(), transmission));
      }
    }
  }
  STATUS(endl);
 
  if (parameters.numberOfHits == 0) {
 
    parameters.numberOfHits = (int) (ids.size() * fraction);
 
    STATUS("Number of hits " << parameters.numberOfHits << endl);
  }
 
  vector<event_type> queue;
  
  for (map<int, map<int, buffer_type> >::iterator i = f1.begin(); i != f1.end(); ++i) {
 
    buffer_type data;
 
    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)) {
          data.push_back(*p);
        }
      }
    }
 
    if (!data.empty()) {
 
      sort(data.begin(), data.end(), make_comparator(&hit_type::getToA, JComparison::lt()));          // sort according time-of-arrival
 
      double      t0 = TOA_s;
    
      buffer_type buffer(factoryLimit);               // local buffer of hits
      event_type  out[2];                             // FIFO of events
 
      TH1D* h1 = M1[i->first];
      TH1D* h2 = M2[i->first];
      TH1D* h3 = M3[i->first];
      TH1D* h4 = M4[i->first];
      TH1D* h5 = M5[i->first];
 
      for (buffer_type::const_iterator p = data.begin(); p != data.end(); ++p) {
 
        if (distance(data.cbegin(), p)%100 == 0) {
          STATUS("processed: " << setw(3) << i->first << ' ' << FIXED(9,3) << p->getToA() - data.begin()->getToA() << " [s]" << '\r' << flush); DEBUG(endl);
        }
 
        buffer_type::const_iterator q = p;
 
        while (++q != data.end() && q->getToA() - p->getToA() <= parameters.TMax_s) {}
 
        h1->Fill(distance(p,q));
 
        if (distance(p,q) >= parameters.numberOfHits) {
 
          if (distance(p,q) < factoryLimit) {
 
            if ((int) out[1].size() >= parameters.numberOfHits && velo(out[1].getVertex(), p, q) >= parameters.numberOfHits) {
 
              // ?
 
            } else {
            
              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;
                }
              }
 
              h2->Fill(distance(root,__q));
 
              if (distance(root,__q) >= parameters.numberOfHits) {
 
                __q = clusterize(__p, __q, match);
 
                h3->Fill(distance(root,__q));
 
                if (distance(root,__q) >= parameters.numberOfHits) {
 
                  const vertex_type vx = velo(*root, __p, __q, parameters.numberOfHits);
 
                  h4->Fill(vx.N);
 
                  if (vx.N >= parameters.numberOfHits) {
                    out[1] = event_type(JEvent(detector.getID(), out[0].getCounter() + 1, i->first, root, __q), vx.getVertex());
                  }
                }
              }
            }
 
          } else {
 
            out[1] = event_type(JEvent(detector.getID(), out[0].getCounter() + 1, i->first, p, q));
          }
 
          if (!out[1].empty()) {
 
            if (out[0].empty()) {
 
              out[0] = out[1];                           // shift
 
            } else if (overlap2D(out[0],out[1])) {
      
              out[0].merge(out[1]);                      // merge
 
            } else {
 
              const double t1 = out[0].begin()->getToA();
 
              STATUS(endl);
              STATUS("trigger:   "
                     << setw(3)    << i->first           << ' '
                     << FIXED(9,3) << t1 - t0 << " [s]"  << ' '
                     << setw(4)    << out[0].size()      << ' '
                     << FIXED(9,2) << getQuality(out[0]) << ' '
                     << out[0].getPosition()             << endl);
 
              //t0  = t1;
            
              h5->Fill(out[0].size());
 
              queue.push_back(out[0]);                   // store
              
              out[0] = out[1];                           // shift
            }
 
            out[1].clear();
          }
        }
      }
 
      if (!out[0].empty()) {
 
        queue.push_back(out[0]);                         // store
      }
      STATUS(endl);
    }
  }
 
  if (!queue.empty()) {
 
    struct time_t {
 
      void operator=(const double t0) { this->t0 = t0; }
 
      operator double() const { return this->t0; }
 
      double t0 = TOA_s;
    };
 
    map<int, time_t> ts;
 
    sort(queue.begin(), queue.end());
    
    for (vector<event_type>::iterator p = queue.begin(); p != queue.end(); ) {
 
      vector<event_type>::iterator q = p;
      vector<event_type>::iterator r = p;
 
      // select event with highest quality
 
      while (++q != queue.end() && overlap1D(*p,*q)) {
        if (getQuality(*q) > getQuality(*r)) {
          r = q;
        }
      }
 
      const double t0 = ts[r->getID()];
      const double t1 = r->begin()->getToA();
 
      STATUS("event:     "
             << setw(3)    << r->getID()         << ' '
             << FIXED(9,3) << t1 - t0 << " [s]"  << ' '
             << setw(4)    << r->size()          << ' '
             << FIXED(9,2) << getQuality(*r)     << ' '
             << r->getPosition()                 << endl);
 
      //ts[r->getID()] = t1;
 
      outputFile.put(*r);
 
      p = q;
 
      G2[r->getID()].put(r->getX(), r->getY());
    }
  }
  
  JMultipleFileScanner<JMeta> io(inputFile);
 
  io >> outputFile;
 
  for (const auto M : { &M1, &M2, &M3, &M4, &M5 }){ 
    for (const auto& i : *M) {
      outputFile.put(*i.second);
    }
  }
 
  for (map<int, JGraph_t>::const_iterator i = G2.begin(); i != G2.end(); ++i) {
    outputFile.put(JGraph(i->second, MAKE_CSTRING("G[" << FILL(2,'0') << i->first << "]")));
  }
 
  outputFile.close();
}