software/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 "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"

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

#include "JROOT/JManager.hh"

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

#include "JTools/JHashMap.hh"

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

#include "JAcoustics/JSoundVelocity.hh"
#include "JAcoustics/JEmitter.hh"
#include "JAcoustics/JAcousticsToolkit.hh"
#include "JAcoustics/JHit.hh"
#include "JAcoustics/JFitParameters.hh"
#include "JAcoustics/JKatoomba.hh"
#include "JAcoustics/JEvent.hh"
#include "JAcoustics/JEvt.hh"
#include "JAcoustics/JEvtToolkit.hh"
#include "JAcoustics/JSupport.hh"
#include "JAcoustics/JTransmission_t.hh"

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


/**
 * \file
 *
 * Application to make a global fit of the detector geometry to acoustic data.\n
 * \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;
  typedef JContainer< set<JTransmission_t> > disable_container;

  JMultipleFileScanner<>  inputFile;
  JFileRecorder<JTYPELIST<JEvt, JEvent, JMeta, TH1D, TH2D>::typelist>  outputFile;
  string                  detectorFile;
  JLimit_t&               numberOfEvents = inputFile.getLimit();
  JSoundVelocity          V              = getSoundVelocity;    // default sound velocity
  tripods_container       tripods;                              // tripods
  hydrophones_container   hydrophones;                          // hydrophones
  JFitParameters          parameters;                           // fit parameters
  int                     fit;                                  // type of fit
  bool                    unify;                                // unify weighing of pings
  disable_container       disable;                              // disable tansmissions
  int                     debug;

  try {

    JParser<> zap("Application to fit position calibration model to acoustic data.");

    zap['f'] = make_field(inputFile,       "output of JAcousticEventBuilder[.sh]");
    zap['o'] = make_field(outputFile);
    zap['n'] = make_field(numberOfEvents)                                    = JLimit::max();
    zap['a'] = make_field(detectorFile);
    zap['@'] = make_field(parameters)                                        = 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['M'] = make_field(getMechanics,    "mechanics data")                 = JPARSER::initialised();
    zap['F'] = make_field(fit,             "fit type")                       = linear_t, simplex_t, gandalf_t;
    zap['u'] = make_field(unify,           "unify weighing of pings");
    zap['!'] = make_field(disable,         "disable transmission")           = JPARSER::initialised();
    zap['d'] = make_field(debug)                                             = 1;

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


  cout.tie(&cerr);

  JDetector detector;

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

  const floor_range range = getRangeOfFloors(detector);

  JHashMap<int, JLocation>   receivers;
  JHashMap<int, JEmitter>    emitters;

  for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
    receivers[i->getID()] = i->getLocation();
  }

  for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) {
    emitters[i->getID()]  = JEmitter(i->getID(),
                                     i->getUTMPosition() - detector.getUTMPosition());
  }

  V.set(detector.getUTMZ());                                          // sound velocity at detector depth

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

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

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

  if (parameters.fixStrings) {
    JKatoomba_t::setOption(false);
  }
  
  JSimplex<JModel>::MAXIMUM_ITERATIONS = 10000;

  typedef JHit<JPDFGauss>  hit_type;

  TH1D h0("chi2/NDF", NULL, 50000,  0.0, 1000.0);
  TH1D h1("h1",       NULL,    51, -0.5,   50.5);

  JManager<int, TH2D> HA(new TH2D("ha[%]", NULL, 
                                  getNumberOfStrings(detector) + 0, -0.5, getNumberOfStrings(detector) - 0.5,
                                  range.getLength(), range.getLowerLimit() - 0.5, range.getLowerLimit() + 0.5));
          
  JManager<int, TH2D> HB(new TH2D("hb[%]", NULL, 
                                  getNumberOfStrings(detector) + 0, -0.5, getNumberOfStrings(detector) - 0.5,
                                  range.getLength(), range.getLowerLimit() - 0.5, range.getLowerLimit() + 0.5));

  for (Int_t i = 1; i <= HA->GetXaxis()->GetNbins(); ++i) {
    HA->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(geometry.at(i-1).first));
    HB->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(geometry.at(i-1).first));
  } 

  // input data

  string oid;                                    // detector identifier

  vector<JEvent> zbuf;

  JLimit_t     limit            = inputFile.getLimit();
  counter_type number_of_events = 0;

  for (JMultipleFileScanner<>::const_iterator file = inputFile.begin(); file != inputFile.end(); ++file) {

    JMultipleFileScanner<JEvent> in(*file);

    limit.setLowerLimit(limit.getLowerLimit() - in.skip(limit.getLowerLimit()));

    const size_t N = zbuf.size();
    
    for ( ; in.hasNext() && number_of_events != limit; ++number_of_events) {

      STATUS("input  " << setw(8) << number_of_events << '\r'); DEBUG(endl);

      const JEvent* evt = in.next();

      if (oid == "") {
        oid = evt->getOID();
      }

      if (oid != evt->getOID()) {                  // consistency check
        FATAL("Invalid detector identifier " << evt->getOID() << " != " << oid << endl);
      }
 
      zbuf.push_back(*evt);
    }

    vector<JEvent>::iterator q = zbuf.begin();

    advance(q, N);

    sort(q, zbuf.end());                           // sort according first time-of-emission

    if (q != zbuf.begin() && q != zbuf.end()) {

      vector<JEvent>::iterator p = q;

      --p;
      
      if (p->empty() || q->empty() || *q < *p) {      
        inplace_merge(zbuf.begin(), q, zbuf.end());
      }
    }
  }
  STATUS(endl);


  outputFile.open();

  outputFile.put(JMeta(argc, argv));

  int counter[] = { 0, 0 };

  for (vector<JEvent>::const_iterator p = zbuf.begin(), q; p != zbuf.end(); p = q) {
  
    STATUS("output " << setw(8) << distance(zbuf.cbegin(), p) << '\r'); DEBUG(endl);

    for (q = p; ++q != zbuf.end() && q->begin()->getToE() <= p->rbegin()->getToE() + parameters.Tmax_s; ) {}

    if (getNumberOfEmitters(p,q) >= parameters.Nmin) {

      h1.Fill((double) getNumberOfEmitters(p,q));

      map<int, int> numberOfPings;

      for (vector<JEvent>::const_iterator i = p; i != q; ++i) {
        numberOfPings[i->getID()] += 1;
      }

      for (map<int, int>::const_iterator i = numberOfPings.begin(); i != numberOfPings.end(); ++i) {
        DEBUG("Number of pings " << setw(2) << i->first << ' ' << setw(3) << i->second << endl);
      }

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

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

      vector<hit_type> data;

      for (vector<JEvent>::const_iterator i = p; i != q; ++i) {

        const JEmitter& emitter = emitters[i->getID()];

        const double    signal  = (unify ? (double) minimum_number_of_pings / (double) numberOfPings[i->getID()] : 1.0);

        // select transmission with highest quality

        map<int, vector<JTransmission> > buffer;

        for (JEvent::const_iterator hit = i->begin(); hit != i->end(); ++hit) {
          if (disable.count(JTransmission_t(i->getID(), hit->getID())) == 0 &&
              disable.count(JTransmission_t(-1,         hit->getID())) == 0) {
            buffer[hit->getID()].push_back(*hit);
          }
        }

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

          if (receivers.has(ps->first) && geometry.hasLocation(receivers[ps->first])) {

            sort(ps->second.begin(), ps->second.end(), make_comparator(&JTransmission::getQ, JComparison::gt()));

            if (ps->second.begin()->getQ() >= parameters.Qmin) {

              data.push_back(hit_type(emitter,
                                      i->getCounter(),
                                      receivers[ps->first],
                                      JPDFGauss(ps->second.begin()->getToA(), parameters.sigma_s, signal, parameters.background)));
            }
          }
        }
      }
      
      for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
        DEBUG("hit: " << *hit << endl);
      }



      JModel result = estimator(data.begin(), data.end());              // pre-fit

      DEBUG("prefit:"                           << endl
            << result                           << endl);


      for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
        HA[hit->getID()]->Fill(geometry.getIndex(hit->getString()), hit->getFloor(), 1.0);
      }


      // remove outliers

      if (parameters.stdev > 0.0) {

        double chi2_old =  evaluator(result, data.begin(), data.end());

        for ( ; ; ) {

          vector<hit_type>::iterator out = data.end();

          double xmax = 0.0;

          for (vector<hit_type>::iterator hit = data.begin(); hit != data.end(); ++hit) {

            const double x = fabs(hit->getValue() - estimator.getToA(result, *hit)) / hit->sigma;

            if (x > xmax) {
              xmax = x;
              out  = hit;
            }
          }

          if (xmax > parameters.stdev) {

            vector<hit_type>::iterator __end = data.end();

            iter_swap(out, --__end);

            result = estimator(data.begin(), __end);

            double chi2_new =  evaluator(result, data.begin(), __end);

            if (chi2_old - chi2_new  >  parameters.stdev * parameters.stdev) {

              DEBUG("Remove outlier " << __end->getLocation() << ' ' << xmax << endl);

              data.pop_back();

              chi2_old = chi2_new;

              continue;
            }
          }
           
          break;
        }
      }

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

      DEBUG("prefit:"                           << endl
            << result                           << endl);

      for (vector<hit_type>::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
        HB[hit->getID()]->Fill(geometry.getIndex(hit->getString()), hit->getFloor(), 1.0);
      }


      double chi2 = numeric_limits<double>::max();

      switch (fit) {

      case linear_t:

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

      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() - result.getN();

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

      h0.Fill(chi2 / (W - result.getN()));
      
      // store results

      if (chi2 / ndf <= parameters.chi2perNDF) {
      
        const JEvt evt = getEvt(JHead(oid, data.begin()->getValue(), data.rbegin()->getValue(), ndf, (W - result.getN()), chi2), result);

        outputFile.put(evt);

        for (vector<JEvent>::const_iterator i = p; i != q; ++i) {

          JEvent out(*i);                                                 // event with fitted times of emission

          const double toe = result.emitter[JEKey(i->getID(), i->getCounter())].t1;

          for (JEvent::iterator hit = out.begin(); hit != out.end(); ++hit) {
            hit->setToE(toe);
          }

          outputFile.put(out);
        }

        counter[0] += 1;
        
      } else {

        WARNING(endl << "Event not written: " << chi2 << '/' << ndf << endl);

        counter[1] += 1;
      }
    }
  }
  STATUS(endl);

  STATUS("Number of events written / rejected: " << counter[0] << " / " << counter[1] << endl);

  JMultipleFileScanner<JMeta> io(inputFile);

  io >> outputFile;

  outputFile.put(h0);
  outputFile.put(h1);

  for (JManager<int, TH2D>::const_iterator i = HA.begin(); i != HA.end(); ++i) { outputFile.put(*(i->second)); }
  for (JManager<int, TH2D>::const_iterator i = HB.begin(); i != HB.end(); ++i) { outputFile.put(*(i->second)); }

  outputFile.close();
}