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/JBase.hh"
#include "JDetector/JTripod.hh"
#include "JDetector/JModule.hh"

#include "JTools/JQuantile.hh"
#include "JGizmo/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/JContainer.hh"
#include "Jeep/JTimer.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


namespace {
  
  using namespace JPP;

  /**
   * Auxiliary data structure for module location and position.
   */
  struct module_type :
    public JLocation, 
    public JPosition3D
  {
    /**
     * Constructor.
     *
     * \paramater  location       module location
     * \paramater  position       module position
     */
    module_type(const JLocation&   location,
                const JPosition3D& position) :
      JLocation  (location),
      JPosition3D(position)
    {}


    /**
     * Less-than operator.
     *
     * \parameter  first          first  module
     * \parameter  second         second module
     * \return                    true if floor of first module less than that of second; else false
     */
    friend inline bool operator<(const module_type& first, const module_type& second)
    {
      return first.getFloor() < second.getFloor();
    }
  };
}


/**
 * \file
 *
 * Example application to test fit of model to acoustic data.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  string                 detectorFile;
  string                 outputFile;
  int                    numberOfEvents;
  int                    numberOfCycles;
  JSoundVelocity         V = getSoundVelocity;   // default sound velocity
  JContainer<JTripod>    tripods;                // tripods
  JContainer<JBase>      bases;                  // base modules
  int                    mestimator;
  UInt_t                 seed;
  int                    debug;

  try {

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

    zap['a'] = make_field(detectorFile);
    zap['o'] = make_field(outputFile);
    zap['n'] = make_field(numberOfEvents);
    zap['N'] = make_field(numberOfCycles,  "number of cycles per emitter")   = 10;
    zap['V'] = make_field(V,               "sound velocity")                 = JPARSER::initialised();
    zap['T'] = make_field(tripods,         "tripod data");
    zap['B'] = make_field(bases,           "base module data")               = JPARSER::initialised();
    zap['M'] = make_field(mestimator)                                        = EM_LINEAR, EM_LORENTZIAN, EM_NORMAL;
    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);

  const double     sigma_s = 1.0e-5;                                  // resolution  of time-of_arrival [s]
  const double     P_bg    = 1.0e-4;                                  // probability of random background
  const JTimeRange T_s(-1.0e-4, +1.0e-4);                             // time range  of random background [s]

  if (numberOfEvents <= 0) { FATAL("Invalid number of events " << numberOfEvents << endl); }
  if (numberOfCycles <= 0) { FATAL("Invalid number of cycles " << numberOfCycles << endl); }

  JDetector detector;

  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }

  vector<JEmitter> emitters;

  for (JContainer<JTripod>::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); }


  JGeometry::JDetector abc;

  {
    map<int, vector<module_type> > buffer;                            // string -> module

    for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
      buffer[module->getString()].push_back(module_type(module->getLocation(), module->getPosition()));
    }

    for (map<int, vector<module_type> >::iterator i = buffer.begin(); i != buffer.end(); ++i) {

      JPosition3D p0;                                                 // reference position of string

      JContainer<JBase>::const_iterator p = find_if(bases.begin(), bases.end(), make_predicate(&JBase::getString, i->first));
      
      if (p != bases.end()) {

        p0 = p->getPosition();

      } else if (i->second.size() >= 2) {

        sort(i->second.begin(), i->second.end());

        const double       z0 = TBARZ_M;                              // position of T-bar
        const module_type& p1 = i->second[0];

        JVector3D slope;

        for (vector<module_type>::const_iterator p2 = i->second.begin(); ++p2 != i->second.end(); ) {
          slope += (*p2 - p1) / p2->getDistance(p1);
        }

        slope /= (double) (i->second.size() - 1);

        // linear extrapolation
        
        p0  =  p1  +  (z0 - p1.getZ()) * slope;

      } else {

        FATAL("Invalid data for string " << i->first << endl);

        continue;
      }

      NOTICE("Position string " 
             << setw(4)    << i->first  << ' ' 
             << FIXED(7,2) << p0.getX() << ' ' 
             << FIXED(7,2) << p0.getY() << ' ' 
             << FIXED(7,2) << p0.getZ() << ' ' 
             << (p != bases.end() ? "user": "calculated") << endl);

      abc[i->first] = JGeometry::JString(p0, i->second.begin(), i->second.end());

      abc[i->first].hydrophone = JPosition3D();          // ?
    }

    DEBUG(abc);
  }


  JKatoomba<JEstimator>         estimator(abc, V);
  JKatoomba<JAbstractMinimiser> evaluator(abc, V);

  evaluator.estimator.reset(getMEstimator(mestimator));


  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, -2.0e-5, +2.0e-5));
  JManager<int, TH2D> H2(new TH2D("string[%]",   NULL, 100, -2.0e-4, +2.0e-4,  100, -2.0e-4, +2.0e-4));


  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(-1.0e-2, +1.0e-2),
                                                            gRandom->Uniform(-1.0e-2, +1.0e-2));
      }
    }


    // set module positions according actual model data

    for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
      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) {

      for (int number_of_cycles = 0; number_of_cycles != numberOfCycles; ++number_of_cycles) {

        ++count;

        const double toe_s = number_of_cycles * 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) {

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

          double t1_s = toa_s;

          if (gRandom->Rndm() >= P_bg)
            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, P_bg)));
        }
      }
    }

    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   = evaluator(result, data.begin(), data.end());

    timer.stop();


    DEBUG("Final values"                << endl
          << FIXED(9,3) << chi2         << '/'
          << data.size() - model.getN() << endl
          << result                     << endl);


    h0.Fill(log10((double) timer.usec_wall));
    h1.Fill(2.0 * chi2 / (double) (data.size() - 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,
                         i->second.ty - result.string [i->first].ty);
    }

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