JBillabong.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <algorithm>
#include <chrono>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleSupportkit.hh"

#include "JLang/JPredicate.hh"
#include "JLang/JComparator.hh"

#include "JTools/JHashMap.hh"
#include "JTools/JAbstractHistogram.hh"

#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMeta.hh"

#include "JGeometry3D/JCylinder3D.hh"
#include "JGeometry3D/JOmega3D.hh"
#include "JGeometry3D/JRotator3D.hh"

#include "JAcoustics/JToA.hh"
#include "JAcoustics/JTransmission.hh"
#include "JAcoustics/JReceiver.hh"
#include "JAcoustics/JSoundVelocity.hh"
#include "JAcoustics/JAcousticsToolkit.hh"
#include "JAcoustics/JEvent.hh"
#include "JAcoustics/JEventOverlap.hh"
#include "JAcoustics/JSupport.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"
#include "JTrigger/JBind2nd.hh"


namespace JACOUSTICS {

  using JGEOMETRY3D::JVector3D;
  using JGEOMETRY3D::JPosition3D;
  using JTRIGGER::JMatch;
  using JTRIGGER::JClonable;
  using JTRIGGER::JBind2nd;
  using JTOOLS::JAbstractHistogram;

  /**
   * Type definition for scan along axis
   */
  typedef JAbstractHistogram<double>         JHistogram_t;


  /**
   * 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 [m/s]
     * \param  Dmax_m           maximal distance [m]
     * \param  Tmax_s           maximal extra time [s]
     */
    JMatch3D(const double V, const double Dmax_m, const double Tmax_s) :
      V(V),
      Dmax_m(Dmax_m),
      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 D = getDistance(first.getPosition(), second.getPosition());

      if (D <= Dmax_m)
        return fabs(first.getToA() - second.getToA())  <=  D / V  +  Tmax_s;
      else
        return false;
    }

  private:
    const double V;
    const double Dmax_m;
    const double Tmax_s;
  };
  

  /** 
   * 1D match criterion for acoustic signals.
   */
  class JMatch1D :
    public JClonable< JMatch<hit_type>, JMatch1D>
  {
  public:
    /**
     * Constructor.
     *
     * \param  V                sound velocity [m/s]
     * \param  Rmax_m           maximal radial   distance [m]
     * \param  Zmax_m           maximal vertical distance [m]
     * \param  Tmax_s           maximal extra time [s]
     */
    JMatch1D(const double V, const double Rmax_m, const double Zmax_m, const double Tmax_s) :
      V(V),
      Rmax_m(Rmax_m),
      Zmax_m(Zmax_m),
      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 
    {
      /*
      const double x = first.getX() - second.getX();
      const double y = first.getY() - second.getY();
      const double z = first.getZ() - second.getZ();

      if (fabs(z) <= Zmax_m) {

        const double R = sqrt(x*x + y*y);

        if (R <= Rmax_m) {
          return fabs(first.getToA() - second.getToA())  <=  sqrt(R*R + z*z) / V  +  Tmax_s;
        }
      }

      return false;
      */
      return fabs(first.getZ() - second.getZ()) <= Zmax_m;
    }

  private:
    const double V;
    const double Rmax_m;
    const double Zmax_m;
    const double Tmax_s;
  };


  /**
   * Auxiliary data structure for final check on event.
   */
  struct JCheck {
    /**
     * Constructor.
     *
     * \param  X                x-scan
     * \param  Y                y-scan
     * \param  Z                z-scan
     * \param  V                sound velocity [m/s]
     * \param  Tmax_s           maximal time [s]
     * \param  Nmin             minimum number of hits
     */
    JCheck(const JHistogram_t& X,
           const JHistogram_t& Y,
           const JHistogram_t& Z,
           const double        V,
           const double        Tmax_s,
           const int           Nmin) :
      V(V),
      Dmax_m(Tmax_s * V),
      Nmin(Nmin)
    {
      for (double x = X.getLowerLimit(); x < X.getUpperLimit() + 0.5*X.getBinWidth(); x += X.getBinWidth()) {
        for (double y = Y.getLowerLimit(); y < Y.getUpperLimit() + 0.5*Y.getBinWidth(); y += Y.getBinWidth()) {
          for (double z = Z.getLowerLimit(); z < Z.getUpperLimit() + 0.5*Z.getBinWidth(); z += Z.getBinWidth()) {
            U.push_back(JVector3D(x,y,z));
          }
        }
      }
    }


    /**
     * Check validity.
     *
     * \return                  true if valid; else false
     */
    bool is_valid() const
    {
      return !U.empty();
    }


    /**
     * Check event.
     *
     * \param  root             root hit
     * \param  __begin          begin of data
     * \param  __end            end   of data
     * \return                  true if accepted; else false
     */
    template<class T>
    inline bool operator()(const hit_type& root, T __begin, T __end) const
    {
      using namespace std;
      using namespace JPP;

      for (vector<JVector3D>::const_iterator u = U.begin(); u != U.end(); ++u) {
          
        const JVector3D p0 = root.getPosition() + *u;                // vertex position
        const double    t0 = root.getToA() * V  -  u->getLength();   // vertex time
        
        int n = 1;

        for (T p = __begin; p != __end; ++p) {

          const double t1 = p->getToA() * V  -  p0.getDistance(*p);

          if (fabs(t1 - t0) <= Dmax_m) {
            n += 1;
          } else if (n + distance(p,__end) == Nmin) {
            break;
          }
        }

        if (n >= Nmin) {
          return true;
        }
      }

      return false;
    }

  private:
    std::vector<JVector3D> U;
    const double           V;
    const double           Dmax_m;
    const int              Nmin;
  };


  /**
   * Print event.
   *
   * \param  out              output stream
   * \param  event            event
   */
  inline void print(std::ostream& out, const JEvent& event)
  {
    using namespace std;
    
    out << "event: " << setw(6) << event.getCounter() << ' ' << FIXED(12,2) << 0.5 * (event.begin()->getToA() + event.rbegin()->getToA()) << " [s] " << setw(4) << event.size() << endl;
  }
}


/**
 * \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, TH1D, JMeta>::typelist  typelist;

  JMultipleFileScanner<JToA> inputFile;
  JLimit_t&                  numberOfEvents = inputFile.getLimit();

  struct {
    int                      factoryLimit   = 100000;
    int                      numberOfHits   =   0;
    double                   Q              =   0.0;
    double                   TMax_s         =   0.0;
    double                   DMax_m         = 500.0;
    double                   ZMax_m         = 200.0;
    double                   TMaxExtra_s    = 100.0e-6; 
    double                   TMaxVertex_s   =   1.0e-2;
    double                   gridAngle_deg  =  10.0;
  } parameters;

  JHistogram_t               X;
  JHistogram_t               Y;
  JHistogram_t               Z;
  JFileRecorder<typelist>    outputFile;
  JSoundVelocity             V = getSoundVelocity;   // default sound velocity
  string                     detectorFile;
  int                        waveform;
  int                        debug;

  try {

    JProperties properties;

    properties.insert(gmake_property(parameters.factoryLimit));
    properties.insert(gmake_property(parameters.numberOfHits));
    properties.insert(gmake_property(parameters.Q));
    properties.insert(gmake_property(parameters.TMax_s));
    properties.insert(gmake_property(parameters.DMax_m));
    properties.insert(gmake_property(parameters.ZMax_m));
    properties.insert(gmake_property(parameters.TMaxExtra_s));
    properties.insert(gmake_property(parameters.TMaxVertex_s));
    properties.insert(gmake_property(parameters.gridAngle_deg));
    
    JParser<> zap("Main program to trigger acoustic data.");
    
    zap['f'] = make_field(inputFile);
    zap['n'] = make_field(numberOfEvents)                                    = JLimit::max();    
    zap['@'] = make_field(properties,   "trigger parameters")                = JPARSER::initialised();
    zap['X'] = make_field(X,            "x-scan");
    zap['Y'] = make_field(Y,            "y-scan");
    zap['Z'] = make_field(Z,            "z-scan");
    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['W'] = make_field(waveform,     "waveform identifier");
    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 double v0 = V(cylinder.getZmax());

  if (parameters.TMax_s < parameters.DMax_m / v0) {

    parameters.TMax_s = parameters.DMax_m / v0;

    NOTICE("Set maximal time [s] " << FIXED(7,3) << parameters.TMax_s << endl);
  }
  
  const JMatch3D      match3D(v0, parameters.DMax_m, parameters.TMaxExtra_s);
  const JMatch1D      match1D(v0, parameters.DMax_m, parameters.ZMax_m, parameters.TMaxExtra_s);

  const JEventOverlap overlap(parameters.TMax_s);

  const JCheck        check(X, Y, Z, v0, parameters.TMaxVertex_s, parameters.numberOfHits);

  const JOmega3D      omega(JAngle3D(0.0, 0.0), JOmega3D::range_type(0.0, 0.5*PI), parameters.gridAngle_deg * PI / 180.0);
  const JRotator3D    rotator(omega);
  

  JHashMap<int, JReceiver>   receivers;

  const JPosition3D   center(cylinder.getX(), cylinder.getY(), 0.5 * (cylinder.getZmin() + cylinder.getZmax()));

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

    if (module->getFloor() != 0) {

      receivers[module->getID()] = JReceiver(module->getID(), 
                                             module->getPosition() + getPiezoPosition() - center,
                                             module->getT0() * 1.0e-9);
    }
  }

  TH1D h0("h0", NULL, 21, -0.5, 20.5);

  vector<TH1D*> M1(10, NULL);

  for (size_t i = 2; i <= 5; ++i) {
    M1[i] = new TH1D(MAKE_CSTRING("M[" << i << "]"), NULL, 1001, -0.5, 1000.5);
  }
      
  outputFile.open();

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

  outputFile.put(JMeta(argc, argv));
  
  // input data

  typedef vector<hit_type> buffer_type;          // data type

  buffer_type data;

  while (inputFile.hasNext()) {

    if (inputFile.getCounter()%1000 == 0) {
      STATUS("counter: " << setw(8) << inputFile.getCounter() << '\r' << flush); DEBUG(endl);
    }

    const JToA* p = inputFile.next();

    if (detector.getID() != p->DETID) {          // consistency check
      FATAL("Invalid detector identifier " << p->DETID << " != " << detector.getID() << endl);
    }

    if (p->WAVEFORMID == waveform && receivers.has(p->DOMID)) {

      if (p->QUALITYFACTOR >= parameters.Q * (parameters.Q <= 1.0 ? p->QUALITYNORMALISATION : 1.0)) {

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

    const double toa = p->TOA_S();

        const JTransmission transmission(p->RUN,
                                         p->DOMID,
                                         p->QUALITYFACTOR,
                                         p->QUALITYNORMALISATION,
                                         receiver.getT(toa),
                                         receiver.getT(toa));
        
        data.push_back(hit_type(receiver.getPosition(), transmission));
      }
    }
  }
  STATUS(endl);


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


  buffer_type buffer(parameters.factoryLimit);    // local buffer of hits
  JEvent      out[2];                             // FIFO of events
 
  const chrono::steady_clock::time_point t0 = chrono::steady_clock::now();

  size_t trigger_counter = 0;

  for (buffer_type::const_iterator p = data.begin(), q = p; p != data.end(); ++p) {

    if (distance(data.cbegin(),p)%10000 == 0) {
      //STATUS("processed: " << FIXED(5,1) << (double) (100 * distance(data.cbegin(),p)) / (double) data.size() << "%" << ' ' << FIXED(12,2) << p->getToA() << " [s]" << '\r' << flush); DEBUG(endl);

      const chrono::steady_clock::time_point t1 = chrono::steady_clock::now();

      STATUS("processed: "
             << FIXED(5,1) << (double) (100 * distance(data.cbegin(),p)) / (double) data.size() << "%"    << ' '
             << FIXED(8,1) << p->getToA()                                                       << " [s]" << ' '
             << setw(6)    << chrono::duration_cast<chrono::seconds>(t1 - t0).count()           << " [s]" << ' '
             << setw(6)    << trigger_counter                                                             << endl);
    }

    h0.Fill(0.0);

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

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

      h0.Fill(1.0);

      if (distance(p,q) < parameters.factoryLimit) {

        h0.Fill(2.0);

        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 (match3D(*root,*i)) {
            *(__q++) = *i;
          }
        }

        M1[2]->Fill(distance(root,__q));

        if (distance(root,__q) >= parameters.numberOfHits) { 

          h0.Fill(3.0);

          __q = clusterize(__p, __q, match3D);

          M1[3]->Fill(distance(root,__q));

          if (distance(root,__q) >= parameters.numberOfHits) {
            
            h0.Fill(4.0);

            const double W = 1.0 / (double) rotator.size();

            for (const JRotation3D& R : rotator) {

              for (buffer_type::iterator i = root; i != __q; ++i) {
                i->rotate(R);
              }

              buffer_type::iterator __z = partition(__p, __q, JBind2nd(match1D,*root));
              
              M1[4]->Fill(distance(root,__z), W);

              if (distance(root,__z) >= parameters.numberOfHits) {

                h0.Fill(5.0, W);
                
                __z = clusterize(__p, __z, match1D);

                M1[5]->Fill(distance(root,__z), W);

                if (distance(root,__z) >= parameters.numberOfHits) {

                  h0.Fill(6.0, W);

                  if (!check.is_valid() || check(*root, __p, __z)) {

                    trigger_counter += 1;

                    h0.Fill(7.0, W);

                    out[1] = JEvent(detector.getID(), out[0].getCounter() + 1, waveform, p, q);

                    break;
                  }
                }
              }
            }
          }
        }
          
      } else {

        trigger_counter += 1;

        out[1] = JEvent(detector.getID(), out[0].getCounter() + 1, waveform, p, q);
      }

      if (!out[1].empty()) {

        if (out[0].empty()) {
          
          out[0] = out[1];                         // shift

        } else if (overlap(out[0],out[1])) {

          out[0].merge(out[1]);                    // merge

        } else {

          if (debug >= notice_t) { print(cout, out[0]); }

          outputFile.put(out[0]);                  // write

          out[0] = out[1];                         // shift
        }
      }

      out[1].clear();
    }
  }

  if (!out[0].empty()) {

    if (debug >= notice_t) { print(cout, out[0]); }

    outputFile.put(out[0]);                        // write
  }
  STATUS(endl);

  outputFile.put(h0);

  for (TH1D* h1 : M1) {
    if (h1 != NULL) {
      outputFile.put(*h1);
    }
  }
  
  JMultipleFileScanner<JMeta> io(inputFile);

  io >> outputFile;

  outputFile.close();
}