examples/JAcoustics/JKatoomba.cc

Go to the documentation of this file.

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

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

#include "JLang/JPredicate.hh"

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JTripod.hh"
#include "JDetector/JModule.hh"
#include "JDetector/JHydrophone.hh"

#include "JROOT/JManager.hh"
#include "JROOT/JRootToolkit.hh"

#include "JAcoustics/JSoundVelocity.hh"
#include "JAcoustics/JEmitter.hh"
#include "JAcoustics/JAcousticsToolkit.hh"
#include "JAcoustics/JHit.hh"
#include "JAcoustics/JKatoomba.hh"

#include "Jeep/JProperties.hh"
#include "Jeep/JContainer.hh"
#include "Jeep/JTimer.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


/**
 * \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<JHydrophone> >  hydrophones_container;

  string                  detectorFile;
  string                  outputFile;
  int                     numberOfEvents;
  map<int, int>           numberOfPings;
  JSoundVelocity          V              = getSoundVelocity;    // default sound velocity
  tripods_container       tripods;                              // tripods
  hydrophones_container   hydrophones;                          // hydrophones
  int                     mestimator     = EM_LINEAR;           // M-estimator
  double                  sigma_s        =  50.0e-6;            // time-of-arrival resolution [s]
  double                  background     =   1.0e-4;            // probability background
  JTimeRange              T_s(-1.0e-4, +1.0e-4);                // time range of background [s]
  int                     fit;                                  // type of fit
  bool                    unique;                               // one ping per cycle
  UInt_t                  seed;
  int                     debug;

  numberOfPings[-1] = 11; // default

  try {
    
    JProperties properties;

    properties.insert(gmake_property(mestimator));
    properties.insert(gmake_property(sigma_s));
    properties.insert(gmake_property(background));
    properties.insert(gmake_property(T_s));

    JParser<> zap("Example application to test fit of model to simukated acoustic data.");

    zap['a'] = make_field(detectorFile);
    zap['o'] = make_field(outputFile);
    zap['n'] = make_field(numberOfEvents);
    zap['@'] = make_field(properties)                                        = 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['H'] = make_field(hydrophones,     "hydrophone data")                = JPARSER::initialised();
    zap['F'] = make_field(fit)                                               = linear_t, simplex_t, gandalf_t;
    zap['u'] = make_field(unique,          "one ping per cycle");
    zap['S'] = make_field(seed)                                              = 0;
    zap['d'] = make_field(debug)                                             = 1;

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  cout.tie(&cerr);

  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()));
  }

  if (detector.empty()) { FATAL("No modules in detector." << endl); }
  if (emitters.empty()) { FATAL("No emitters in system." << endl); }


  const JGeometry geometry(detector, hydrophones);

  DEBUG(geometry);


  JKatoomba<JEstimator>         estimator(geometry, V);
  JKatoomba<JAbstractMinimiser> evaluator(geometry, V);
  JKatoomba<JSimplex>           simplex  (geometry, V);
  JKatoomba<JGandalf>           gandalf  (geometry, V);

  simplex.estimator.reset(getMEstimator(mestimator));
  gandalf.estimator.reset(getMEstimator(mestimator));

  simplex.debug = debug;
  gandalf.debug = debug;

  JSimplex<JModel>::MAXIMUM_ITERATIONS = 10000;


  typedef JHit<JPDFGauss>  hit_type;


  TH1D h0("cpu",      NULL, 100, 1.0, 7.0);
  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, -1.0, +1.0,  100, -1.0, +1.0));


  JTimer timer;

  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);

    // set module positions according actual model data

    for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
      if (geometry.hasLocation(module->getLocation())) {
        module->set(estimator.detector[module->getString()].getPosition(model.string[module->getString()], module->getFloor()));
      }
    }


    // generate hits

    vector<hit_type> data;

    for (vector<JEmitter>::const_iterator emitter = emitters.begin(); emitter != emitters.end(); ++emitter) {

      int number_of_pings = numberOfPings[-1];

      if (numberOfPings.find(emitter->getID()) != numberOfPings.end()) {
        number_of_pings = numberOfPings[emitter->getID()];
      }

      int weight = numeric_limits<int>::max();

      for (map<int, int>::const_iterator i = numberOfPings.begin(); i != numberOfPings.end(); ++i) {
        if (i->second < weight) {
          weight = i->second;
        }
      }

      const double signal = (unique ? (double) weight / (double) number_of_pings : 1.0);

      for (int ping_counter = 0; ping_counter != number_of_pings; ++ping_counter) {

        ++count;

        const double toe_s = ping_counter * 10.0  +  gRandom->Uniform(-1.0, +1.0);

        model.emitter[emitter->getID()][count] = toe_s;

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

          if (geometry.hasLocation(module->getLocation())) {

            const double toa_s = toe_s + V.getTime(getDistance(module->getPosition(), emitter->getPosition()), emitter->getZ(), module->getZ());

            double t1_s = toa_s;

            if (gRandom->Rndm() >= background)
              t1_s = gRandom->Gaus(toa_s, sigma_s);
            else
              t1_s = gRandom->Uniform(toa_s + T_s.getLowerLimit(),
                                    toa_s + T_s.getUpperLimit());

            data.push_back(hit_type(*emitter,
                                    count,
                                    module->getLocation(),
                                    JPDFGauss(t1_s, sigma_s, signal, background)));
          }
        }
      }
    }

    DEBUG("Model" << endl << model << endl);

    for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
      DEBUG("hit: " << *hit << endl);
    }


    timer.reset();
    timer.start();

    JModel result = estimator(data.begin(), data.end());              // pre-fit
    double chi2   = numeric_limits<double>::max();

    switch (fit) {

    case linear_t:

      chi2   = evaluator(result, data.begin(), data.end());
      break;

    case simplex_t:

      simplex.value = result;                                         // start value

      chi2   = simplex(data.begin(), data.end()) / simplex.estimator->getRho(1.0);
      result = simplex.value;
      break;

    case gandalf_t:

      gandalf.value = result;                                         // start value

      chi2   = gandalf(data.begin(), data.end()) / gandalf.estimator->getRho(1.0);
      result = gandalf.value;
      break;

    default:
      break;
    }

    timer.stop();

    double W = 0.0;

    for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
      W += hit->getWeight();
    }

    const int ndf = data.size() - model.getN();

    DEBUG("Final values"                      << endl
          << FIXED(9,3) << chi2 << '/' << ndf << endl
          << result                           << endl);


    h0.Fill(log10((double) timer.usec_wall));
    h1.Fill(chi2 / (double) (W - model.getN()));

    for (JHashMap<int, JMODEL::JString>::const_iterator i = model.string.begin(); i != model.string.end(); ++i) {
      H2[i->first]->Fill((i->second.tx - result.string [i->first].tx) * 1.0e3,   // mrad
                         (i->second.ty - result.string [i->first].ty) * 1.0e3);  // mrad
    }

    for (JHashMap<JEKey, JMODEL::JEmitter>::const_iterator i = model.emitter.begin(); i != model.emitter.end(); ++i) {
      H1[i->first.getID()]->Fill(i->second.t1 - result.emitter[i->first].t1);
    }
  }
  STATUS(endl);


  TFile out(outputFile.c_str(), "recreate");

  out << h0 
      << h1 
      << H1
      << H2;

  out.Write();
  out.Close();
}