toa.cc

Go to the documentation of this file.

#include <iostream>
#include <iomanip>
#include <map>
 
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
 
#include "JGizmo/JManager.hh"
#include "JUTM/JUTMPosition.hh"
 
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
 
#include "JDB/JDB.hh"
#include "JDB/JSelector.hh"
#include "JDB/JSelectorSupportkit.hh"
#include "JDB/JToAshort.hh"
#include "JDB/JDBToolkit.hh"
 
#include "JLang/JLangToolkit.hh"
#include "JLang/JPredicate.hh"
#include "JTools/JRange.hh"
#include "JTools/JQuantile.hh"
#include "JROOT/JRootPrinter.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
namespace {
 
  /**
   * Auxiliary data structure for acoustic data.
   */
  struct JElement_t {
    /**
     * Constructor.
     *
     * \param  t1         time stamp [s]
     * \param  q1         quality
     */
    JElement_t(double t1,
               double q1) :
      t1(t1),
      q1(q1)
    {}
 
    /**
     * Less-than operator.
     *
     * \param  first      first  measurement
     * \param  second     second measurement
     * \return            true if first measurement earlier than second; esle false
     */
    friend inline bool operator<(const JElement_t& first, const JElement_t& second)
    {
      return first.t1 < second.t1;
    }
 
    double t1;   //!< time stamp [s]
    double q1;   //!< quality
  };
}
 
 
/**
 * \file
 *
 * Example program to analyse acoustic data from data base.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  typedef JRange<int> JRange_t; 
 
  string    usr;
  string    pwd;
  string    cookie;
  string    detectorFile;
  string    outputFile;
  string    detid;
  JRange_t  runs;
  double    TMax_s;
  double    Q;
  int       debug;
 
  try {
 
    JParser<> zap("Example program to analyse acoustic data from data base.");
    
    zap['u'] = make_field(usr)              = "";
    zap['!'] = make_field(pwd)              = "";
    zap['C'] = make_field(cookie)           = "";
    zap['a'] = make_field(detectorFile);
    zap['o'] = make_field(outputFile)       = "toa.root";
    zap['D'] = make_field(detid);
    zap['R'] = make_field(runs);
    zap['T'] = make_field(TMax_s);
    zap['Q'] = make_field(Q)                = 0.0;
    zap['d'] = make_field(debug)            = 1;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  const double SOUND_VELOCITY  =  1481.0;   // [m/s]
 
 
  JDetector detector;
 
  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }
 
  const JModuleRouter router(detector);
 
  const int MODULE[] = {
    detector.getModule(JModuleLocation(1,  1)).getID(),
    detector.getModule(JModuleLocation(1, 17)).getID()
  };
 
  /*
    Waveform #12 (sweep 26 kHz -28 kHz. 5 ms long)--> Tripod 1
    Waveform #14 (sweep 30 kHz -32 kHz. 5 ms long)--> Tripod 2
    Waveform #16 (sweep 34 kHz -36 kHz. 5 ms long)--> Tripod 3 
  */
 
  map<int, JPosition3D> tripod;
 
  tripod[12] = (JUTMPosition(587888, 4017079, detector.getUTMZ()) - detector.getUTMPosition()).getPosition();
  tripod[14] = (JUTMPosition(587436, 4016992, detector.getUTMZ()) - detector.getUTMPosition()).getPosition();
  tripod[16] = (JUTMPosition(587843, 4016727, detector.getUTMZ()) - detector.getUTMPosition()).getPosition();
 
  DEBUG("Tripods" << endl);
 
  for (map<int, JPosition3D>::const_iterator i = tripod.begin(); i != tripod.end(); ++i) {
    
    const JPosition3D pos = i->second  -  JPosition3D(detector.getModule(JModuleLocation(1, 1)).getX(),
                                                      detector.getModule(JModuleLocation(1, 1)).getY(),
                                                      0.0);
 
    DEBUG(setw(2)    << i->first   << ' ' 
          << FIXED(7,2) << pos.getX() << ' ' 
          << FIXED(7,2) << pos.getY() << ' ' 
          << FIXED(7,2) << pos.getZ() << ' ' 
          << FIXED(7,2) << sqrt(pos.getX()*pos.getX() + pos.getY()*pos.getY())
          << endl);
  }
 
  JManager<int, TH1D> H1(new TH1D("M[%]", NULL, 10000, -0.25, +0.25));
  
  typedef vector<JElement_t> buffer_type;
 
  map<int,                    // emitter
    map< int,                 // module
    buffer_type> > f1;        // data
 
  try {
 
    JDB::reset(usr, pwd, cookie);
 
    if (is_integer(detid))
      detid = getDetector(to_value<int>(detid));
    else
      getDetector(detid);
 
    JSelector selector = getSelector<JToAshort>(detid,
                                                runs.getLowerLimit(),
                                                runs.getUpperLimit());
 
    ResultSet& rs  = getResultSet(getTable<JToAshort>(), selector);
 
    int number_of_rows = 0;
 
    for (JToAshort parameters; rs >> parameters; ++number_of_rows) {
 
      STATUS(setw(8) << number_of_rows << '\r' << flush); DEBUG(endl);
 
      if (debug >= debug_t) {
        JRootPrinter::print(cout, parameters);
      }
 
      const double d1 = (router.getModule(parameters.DOMID).getPosition() - tripod[parameters.EMITTERID]).getLength();
      const double t1 =  parameters.UNIXTIMEBASE  +  parameters.TOA_S  -  d1 / SOUND_VELOCITY;
      const double q1 =  parameters.QUALITYFACTOR;
 
      f1[parameters.EMITTERID][parameters.DOMID].push_back(JElement_t(t1,q1));
    }
    STATUS(endl);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
 
  for (map<int, map< int, buffer_type> >::iterator emitter = f1.begin(); emitter != f1.end(); ++emitter) {
 
    TH1D* h1 = H1[emitter->first];
 
    // data selection
 
    for (int i = 0; i != sizeof(MODULE)/sizeof(MODULE[0]); ++i) {
 
      buffer_type& buffer = emitter->second[MODULE[i]];
 
      JQuantile Q1("quality", buffer.begin(), buffer.end(), &JElement_t::q1);
 
      const double Qmin = Q1.getQuantile(Q);
 
      buffer_type::iterator __end = partition(buffer.begin(), buffer.end(), make_predicate(&JElement_t::q1, Qmin, JComparison::gt()));
 
      buffer.erase(__end, buffer.end());
 
      sort(buffer.begin(), buffer.end());
    }
 
    // correlator
 
    buffer_type::const_iterator p0   = emitter->second[MODULE[0]].begin();
    buffer_type::const_iterator p1   = emitter->second[MODULE[1]].begin();
 
    buffer_type::const_iterator end0 = emitter->second[MODULE[0]].end();
    buffer_type::const_iterator end1 = emitter->second[MODULE[1]].end();
 
    if (p0 != end0 && p1 != end1) {
 
      for ( ; ; ++p0, ++p1) {
 
        while (p0 != end0  &&  p0->t1 < p1->t1 - TMax_s) { ++p0; }
        while (p1 != end1  &&  p1->t1 < p0->t1 - TMax_s) { ++p1; }
 
        if (p0 != end0 && p1 != end1)
          h1->Fill(p1->t1 - p0->t1);
        else
          break;
      }
    }
  }
 
 
  TFile out(outputFile.c_str(), "recreate");
 
  out << H1;
 
  out.Write();
  out.Close();
}