examples/JAcoustics/JKatoomba.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
#include <map>
#include <algorithm>
#include <chrono>
 
#include "TRandom3.h"
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
 
#include "JLang/JPredicate.hh"
#include "JLang/JSTDToolkit.hh"
 
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JTripod.hh"
#include "JDetector/JTransmitter.hh"
#include "JDetector/JHydrophone.hh"
#include "JDetector/JModule.hh"
 
#include "JROOT/JManager.hh"
#include "JROOT/JRootToolkit.hh"
 
#include "JAcoustics/JSoundVelocity.hh"
#include "JAcoustics/JEmitter.hh"
#include "JAcoustics/JAcousticsToolkit.hh"
#include "JAcoustics/JAcousticsSupportkit.hh"
#include "JAcoustics/JHit.hh"
#include "JAcoustics/JFitParameters.hh"
#include "JAcoustics/JGlobalfit.hh"
#include "JAcoustics/JTransmission_t.hh"
 
#include "Jeep/JContainer.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
namespace {
 
  using namespace JPP;
 
  static const int       WILDCARD         =  -1;          // wild card for emitter identifier
 
  static const int       NUMBER_OF_PINGS  =  11;          // default number of pings
  static const JWaveform EMITTER_WAVEFORM(2.5e6,          // default emitter power [quality]
                                          30.0);          // default emitter frequency [kHz]
}
 
 
/**
 * \file
 *
 * Example application to test fit of model to simulated acoustic data.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  typedef JContainer< vector<JTripod> >      tripods_container;
  typedef JContainer< vector<JTransmitter> > transmitters_container;
  typedef JContainer< vector<JHydrophone> >  hydrophones_container;
  typedef JContainer< map<int, JWaveform> >  waveform_container;
  typedef JContainer< set<JTransmission_t> > disable_container;
 
  string                  detectorFile;
  string                  outputFile;
  int                     numberOfEvents;
  map<int, int>           numberOfPings;
  JSoundVelocity          V              = getSoundVelocity;    // default sound velocity
  tripods_container       tripods;                              // tripods
  transmitters_container  transmitters;                         // transmitters
  hydrophones_container   hydrophones;                          // hydrophones
  JDetectorMechanics      mechanics;                            // mechanical model data
  JTimeRange              T_s(-1.0e-4, +1.0e-4);                // time range of background [s]
  JFitParameters          parameters;                           // fit parameters
  bool                    unify;                                // unify weighing of pings
  disable_container       disable;                              // disable tansmissions
  waveform_container      waveforms;                            // emitter power and frequency
  double                  background;
  ULong_t                 seed;
  int                     debug;
 
  numberOfPings[WILDCARD] = NUMBER_OF_PINGS;
  waveforms    [WILDCARD] = EMITTER_WAVEFORM;
 
  try {
    
    JParser<> zap("Example application to test fit of model to simulated acoustic data.");
 
    zap['a'] = make_field(detectorFile);
    zap['o'] = make_field(outputFile);
    zap['n'] = make_field(numberOfEvents);
    zap['@'] = make_field(parameters)                                        = JPARSER::initialised();
    zap['N'] = make_field(numberOfPings,   "number of pings per emitter")    = JPARSER::initialised();
    zap['V'] = make_field(V,               "sound velocity")                 = JPARSER::initialised();
    zap['T'] = make_field(tripods,         "tripod data");
    zap['Y'] = make_field(transmitters,    "transmitter data")               = JPARSER::initialised();
    zap['H'] = make_field(hydrophones,     "hydrophone data")                = JPARSER::initialised();
    zap['u'] = make_field(unify,           "unify weighing of pings");
    zap['M'] = make_field(mechanics,       "mechanics data")                 = JPARSER::initialised();
    zap['!'] = make_field(disable,         "disable transmission")           = JPARSER::initialised();
    zap['W'] = make_field(waveforms,       "emitter power and frequency")    = JPARSER::initialised();
    zap['B'] = make_field(background,      "background probability")         = 0.0;
    zap['S'] = make_field(seed)                                              = 0;
    zap['d'] = make_field(debug)                                             = 1;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  gRandom->SetSeed(seed);
 
  if (numberOfEvents <= 0) { FATAL("Invalid number of events " << numberOfEvents << endl); }
 
  JDetector detector;
 
  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }
 
  vector<JEmitter> emitters;
 
  for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
    emitters.push_back(JEmitter(i->getID(),
                                i->getUTMPosition() - detector.getUTMPosition()));
  }
 
  for (transmitters_container::const_iterator i = transmitters.begin(); i != transmitters.end(); ++i) {
    try {
      emitters.push_back(JEmitter(i->getID(),
                                  i->getPosition() + detector.getModule(i->getLocation()).getPosition()));
    }
    catch(const exception&) {
      continue;                                // if no string available, discard transmitter
    }
  }
 
  if (detector.empty()) { FATAL("No modules in detector." << endl); }
  if (emitters.empty()) { FATAL("No emitters in system."  << endl); }
 
 
  V.set(detector.getUTMZ());                                          // sound velocity at detector depth
 
  const double D_m = -detector.getUTMZ();                             // depth [m]
 
  JGeometry   geometry(detector, mechanics, hydrophones);
  JGlobalfit  katoomba(geometry, V, parameters);
 
  JKatoomba<JGandalf>::debug              = debug;
  JKatoomba<JGandalf>::MAXIMUM_ITERATIONS = 10000;
  
  DEBUG(geometry);
 
  typedef vector<JHit> data_type;
 
 
  TH1D h0("cpu",      NULL, 100, 0.0, 1.0e3);
  TH1D h1("chi2/NDF", NULL, 100, 0.0, 5.0);
 
  JManager<int, TH1D> H1(new TH1D("emitter[%]",  NULL, 100, -5.0e-5, +5.0e-5));
  JManager<int, TH2D> H2(new TH2D("string[%]",   NULL, 100, -2.0, +2.0,  100, -2.0, +2.0));
 
 
  for (int number_of_events = 0, count = 0; number_of_events != numberOfEvents; ++number_of_events) {
 
    STATUS("event: " << setw(10) << number_of_events << '\r'); DEBUG(endl);
 
    JModel model;                                                     // randomised model data
 
    for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
 
      if (model.string[module->getString()] == JMODEL::JString()) {
 
        model.string[module->getString()] = JMODEL::JString(gRandom->Uniform(-2.0e-2, +2.0e-2),
                                                            gRandom->Uniform(-2.0e-2, +2.0e-2));
      }
    }
 
    DEBUG("Model" << endl << model << endl);
 
    // set module positions according actual model data
 
    for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
      if (module->getFloor() != 0 && geometry.hasLocation(module->getLocation())) {
        module->set(geometry[module->getString()].getPosition(model.string[module->getString()], module->getFloor()));
      }
    }
 
 
    // generate hits
 
    int minimum_number_of_pings = numeric_limits<int>::max();
 
    for (vector<JEmitter>::const_iterator emitter = emitters.begin(); emitter != emitters.end(); ++emitter) {
      minimum_number_of_pings = min(minimum_number_of_pings, getValue(numberOfPings, emitter->getID(), numberOfPings[WILDCARD]));
    }
 
    data_type data;
 
    for (vector<JEmitter>::const_iterator emitter = emitters.begin(); emitter != emitters.end(); ++emitter) {
 
      const int number_of_pings = getValue(numberOfPings, emitter->getID(), numberOfPings[WILDCARD]);
 
      const double weight = (unify ? (double) minimum_number_of_pings / (double) number_of_pings : 1.0);
 
      for (int ping_counter = 0; ping_counter != number_of_pings; ++ping_counter) {
 
        ++count;
 
        const double toe_s = emitter->getID() * 1000.0  +  ping_counter * 10.0  +  gRandom->Uniform(-1.0, +1.0);
 
        model.emission[emitter->getID()][count] = toe_s;
 
        for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
 
          if (disable.count(JTransmission_t(emitter->getID(), module->getID())) == 0) {
 
            if (geometry.hasLocation(module->getLocation())) {
 
              const JWaveform& waveform = getValue(waveforms, emitter->getID(), waveforms[WILDCARD]);
 
              const double d_m   = getDistance(module->getPosition(), emitter->getPosition());
              const double toa_s = toe_s + V.getTime(d_m, emitter->getZ(), module->getZ());
              const double Q     = waveform.getQ(D_m, d_m);
 
              if (Q >= parameters.Qmin) {
 
                double t1_s = toa_s;
 
                if (gRandom->Rndm() >= background)
                  t1_s = gRandom->Gaus(toa_s, parameters.sigma_s);
                else
                  t1_s = gRandom->Uniform(toa_s + T_s.getLowerLimit(),
                                          toa_s + T_s.getUpperLimit());
 
                const JHit hit(*emitter,
                               count,
                               module->getLocation(),
                               t1_s,
                               parameters.sigma_s,
                               weight);
 
                DEBUG("hit: " << hit << ' ' << FIXED(7,1) << Q << endl);
 
                data.push_back(hit);
              }
            }
          }
        }
      }
    }
 
    const chrono::steady_clock::time_point t0 = chrono::steady_clock::now();
 
    const auto result = katoomba(data.begin(), data.end());
 
    const chrono::steady_clock::time_point t1 = chrono::steady_clock::now();
 
    if (debug >= debug_t) {
 
      cout << "result:"                                                     << ' '
           << FIXED(12,3) << result.chi2 << '/' << FIXED(6,1) << result.ndf << endl
           << result.value                                                  << endl;
      
      for (data_type::const_iterator hit = result.begin; hit != result.end; ++hit) {
        cout << "hit: " << *hit << ' ' << FIXED(9,3) << katoomba.evaluator(result.value, *hit) << endl;
      }
    }
 
    h0.Fill(chrono::duration_cast<chrono::milliseconds>(t1 - t0).count());
    h1.Fill(result.chi2 / result.ndf);
 
    for (JHashMap<int, JMODEL::JString>::const_iterator i = model.string.begin(); i != model.string.end(); ++i) {
 
      const double tx = (i->second.tx - result.value.string [i->first].tx) * 1.0e3;   // mrad
      const double ty = (i->second.ty - result.value.string [i->first].ty) * 1.0e3;  // mrad
 
      H2          ->Fill(tx, ty);
      H2[i->first]->Fill(tx, ty);                        
    }
 
    for (JHashMap<JEKey, JMODEL::JEmission>::const_iterator i = model.emission.begin(); i != model.emission.end(); ++i) {
 
      const double t1 = i->second.t1 - result.value.emission[i->first].t1;
 
      H1                  ->Fill(t1);
      H1[i->first.getID()]->Fill(t1);
    }
  }
  STATUS(endl);
 
 
  TFile out(outputFile.c_str(), "recreate");
 
  out << h0 
      << h1 
      << H1 << *H1
      << H2 << *H2;
 
  out.Write();
  out.Close();
}