JAcousticsTriggerProcessor.cc File Reference

Main program to trigger acoustic data. More...

#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/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"

Go to the source code of this file.

Classes
struct	JACOUSTICS::hit_type
	Transmission with position. More...
class	JACOUSTICS::JMatch3D
	3D match criterion for acoustic signals. More...
struct	JACOUSTICS::event_type
	Event with vertex. More...
class	JACOUSTICS::JEventOverlap
	Match of two events considering overlap in time and position. More...
struct	JACOUSTICS::vertex_type
	Vertex. More...
struct	JACOUSTICS::JVelo
	Vertex locator. More...

Namespaces
namespace	JACOUSTICS
	Auxiliary classes and methods for acoustic position calibration.

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

Detailed Description

Main program to trigger acoustic data.

An acoustic event is based on coincidences between times of arrival.
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

Definition in file JAcousticsTriggerProcessor.cc.

Function Documentation

main()

int main	(	int	argc,
		char **	argv )

Definition at line 423 of file JAcousticsTriggerProcessor.cc.

{
  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
  double                   grid           = 10.0;
  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));
    properties.insert(gmake_property(grid));
    
    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,           "waveform identifiers")              = 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 JVector3D   center(cylinder.getX(), cylinder.getY(), Z_m);
  
  //const JConstantSoundVelocity V0 = V(cylinder.getZmax());
  const JAbstractSoundVelocity& V0 = V;
 
  const JVelo Velo(V0, RMax_m, Xv_m*grid, factor);
  const JVelo velo(V0, RMax_m/grid, 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
 
  STATUS("reading data... " << flush); DEBUG(endl);
 
  while (inputFile.hasNext()) {
 
    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 < event_type::TOA_s) {
        event_type::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("OK" << 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) {
 
    STATUS("processing: " << setw(3) << i->first << endl);
 
    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
    
      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) {
 
        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) {
 
                  vertex_type vx = Velo(center, *root, __p, __q, parameters.numberOfHits);
 
                  if (vx.N >= parameters.numberOfHits) {
 
                    vx = velo(vx.getPosition(), *root, __p, __q, parameters.numberOfHits);
 
                    h4->Fill(vx.N);
 
                    if (vx.N >= parameters.numberOfHits) {
                      out[1] = event_type(vx.getVertex(), JEvent(detector.getID(), out[0].getCounter() + 1, i->first, root, __q)); 
                    }
                  }
                }
              }
            }
 
          } 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 {
 
              STATUS("trigger: " << out[0] << endl);
 
              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()) {
 
    sort(queue.begin(), queue.end());
    
    for (vector<event_type>::iterator p = queue.begin(); p != queue.end(); ) {
 
      // select event with highest quality within given time window
 
      vector<event_type>::iterator q = p;
      vector<event_type>::iterator i = p;
 
      STATUS(endl);
 
      for ( ; q != queue.end() && overlap1D(*p,*q); ++q) {
 
        STATUS("|        " << *q << endl);
 
        if (q->getQ() > i->getQ()) {
          i = q;
        }
      }
 
      STATUS("event:   " << *i << endl);
 
      outputFile.put(*i);
 
      G2[i->getID()].put(i->getX(), i->getY());
 
      p = q;
    }
  }
  
  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();
}