software/JAcoustics/JKatoomba.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
#include <map>
#include <deque>
#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 "JROOT/JROOTClassSelector.hh"
 
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JTripod.hh"
#include "JDetector/JTransmitter.hh"
#include "JDetector/JModule.hh"
#include "JDetector/JHydrophone.hh"
 
#include "JTools/JHashMap.hh"
#include "JTools/JRange.hh"
 
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JSingleFileScanner.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/JGlobalfit.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<JTransmitter> > transmitters_container;
  typedef JContainer< vector<JHydrophone> >  hydrophones_container;
  typedef JContainer< set<JTransmission_t> > disable_container;
 
  JMultipleFileScanner<JEvent>  inputFile;
  JFileRecorder<JTYPELIST<JEvt, JEvent, JFitParameters, JMeta, TH1D, TH2D>::typelist>  outputFile;
  string                  detectorFile;
  JLimit_t&               numberOfEvents = inputFile.getLimit();
  JSoundVelocity          V              = getSoundVelocity;    // default sound velocity
  tripods_container       tripods;                              // tripods
  transmitters_container  transmitters;                         // transmitters
  hydrophones_container   hydrophones;                          // hydrophones
  JFitParameters          parameters;                           // fit parameters
  disable_container       disable;                              // disable tansmissions
  JROOTClassSelection     selection      = getROOTClassSelection<JTYPELIST<JEvent>::typelist>();
  JRange<double>          utc;
  bool                    squash;
  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['Y'] = make_field(transmitters,    "transmitter data")               = JPARSER::initialised();
    zap['H'] = make_field(hydrophones,     "hydrophone data")                = JPARSER::initialised();
    zap['M'] = make_field(getMechanics,    "mechanics data")                 = JPARSER::initialised();
    zap['!'] = make_field(disable,         "disable transmission")           = JPARSER::initialised();
    zap['C'] = make_field(selection,
                          "Precede name of data structure by a '+' or '-' "
                          "to add or remove data types in the output, respectively."
                          "\nROOT wildcards are accepted.") = JPARSER::initialised();
    zap['u'] = make_field(utc,             "select UTC time range")          = JRange<double>();
    zap['q'] = make_field(squash,          "squash meta data");
    zap['d'] = make_field(debug)                                             = 1;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  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());
  }
 
  for (transmitters_container::const_iterator i = transmitters.begin(); i != transmitters.end(); ++i) {
    try {
      emitters[i->getID()]  = JEmitter(i->getID(), i->getPosition() + detector.getModule(i->getLocation()).getPosition());
    }
    catch(const exception&) {}                       // if no module available, discard transmitter
  }
 
  V.set(detector.getUTMZ());                         // sound velocity at detector depth
 
  JGeometry   geometry(detector, hydrophones);
  JGlobalfit  katoomba(geometry, V, parameters);
 
  JKatoomba<JGandalf>::debug              = debug;
  JKatoomba<JGandalf>::MAXIMUM_ITERATIONS = 10000;
  
  DEBUG(geometry);
  
  typedef vector<JHit> data_type;
 
  TH1D h0("chi2/NDF", NULL, 50000,  0.0, 1000.0);
  TH1D h1("h1",       NULL,    51, -0.5,   50.5);
  TH1D hn("hn",       NULL,   100,  0.0,    6.0);
 
  JManager<int, TH2D> HA(new TH2D("ha[%]", NULL, 
                                  getNumberOfStrings(detector) + 0, -0.5, getNumberOfStrings(detector) - 0.5,
                                  range.getLength() + 1, range.getLowerLimit() - 0.5, range.getUpperLimit() + 0.5));
          
  JManager<int, TH2D> HB(new TH2D("hb[%]", NULL, 
                                  getNumberOfStrings(detector) + 0, -0.5, getNumberOfStrings(detector) - 0.5,
                                  range.getLength() + 1, range.getLowerLimit() - 0.5, range.getUpperLimit() + 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));
  } 
 
  if (inputFile.size() > 1u) {                       // sort input files
 
    map<double, string> zmap;
 
    for (const string& file_name : inputFile) {
 
      STATUS(file_name << '\r'); DEBUG(endl);
 
      for (JSingleFileScanner<JEvent> in(file_name, 1); in.hasNext(); ) {
 
        const JEvent* evt = in.next();
 
        if (detector.getID() != evt->getDetectorID()) {
          FATAL("Invalid detector identifier " << evt->getDetectorID() << " != " << evt->getDetectorID() << endl);
        }
 
        if (!evt->empty()) {
          zmap[evt->begin()->getToE()] = file_name;
        }
      }
    }
    STATUS(endl);
 
    inputFile.clear();
 
    for (map<double, string>::const_iterator i = zmap.begin(); i != zmap.end(); ++i) {
      inputFile.push_back(i->second);
    }
  }
 
  outputFile.open();
 
  outputFile.put(JMeta(argc, argv));
  outputFile.put(parameters);
 
  try {                                              // limit RAM usage
    dynamic_cast<JTreeWriterObjectOutput<JEvent>&>(*outputFile).getTreeWriter().SetAutoFlush(5000);
  }
  catch(const exception&) {}
    
  int counter[] = { 0, 0 };
 
  typedef deque<JEvent> buffer_type;
 
  for (buffer_type zbuf; inputFile.hasNext(); ) {
 
    STATUS(inputFile.getFilename() << '\r'); DEBUG(endl);
 
    // read one file at a time
    
    for (const string file_name = inputFile.getFilename(); inputFile.hasNext() && file_name == inputFile.getFilename(); ) {
 
      const JEvent* evt = inputFile.next();
 
      if (!evt->empty() && emitters.has(evt->getID())) {
        if (utc(evt->begin()->getToA()) && utc(evt->rbegin()->getToA())) {
          zbuf.push_back(*evt);
        }
      }
    }
 
    sort(zbuf.begin(), zbuf.end());                  // sort according first time-of-emission
 
    for (buffer_type::iterator p = zbuf.begin(), q; p != zbuf.end(); p = q) {
  
      for (q = p; ++q != zbuf.end() && q->begin()->getToE() <= p->rbegin()->getToE() + parameters.Tmax_s; ) {}
 
      if (q == zbuf.end()) {
        
        if (inputFile.hasNext()) {
          
          zbuf.erase(zbuf.begin(), p);               // remove processed data and continue reading
 
          break;
        }
      }
 
      JEvent::overlap(p, q, parameters.deadTime_s);  // empty overlapping events
 
      if (getNumberOfEmitters(p,q) >= parameters.Nmin) {
 
        h1.Fill((double) getNumberOfEmitters(p,q));
 
        const JWeight getWeight(p, q);
 
        data_type data;
 
        for (buffer_type::iterator evt = p; evt != q; ++evt) {
 
          sort(evt->begin(), evt->end(), JKatoomba<>::compare);
 
          JEvent::iterator __end  = unique(evt->begin(), evt->end(), make_comparator(&JTransmission::getID, JComparison::eq()));
 
          const JEmitter& emitter = emitters [evt->getID()];
          const double    weight  = getWeight(evt->getID());
 
          for (JEvent::const_iterator i = evt->begin(); i != __end; ++i) {
 
            if (disable.count(JTransmission_t(evt->getID(), i->getID())) == 0 &&
                disable.count(JTransmission_t(-1,           i->getID())) == 0) {
 
              if (receivers.has(i->getID()) && geometry.hasLocation(receivers[i->getID()]) && i->getQ() >= parameters.Qmin * (parameters.Qmin <= 1.0 ? i->getW() : 1.0)) {
 
                data.push_back(JHit(emitter,
                                    distance(p,evt),
                                    receivers[i->getID()],
                                    i->getToA(),
                                    parameters.sigma_s,
                                    weight));
              }
            }
          }
        }
 
        if (getMinimumNumberOfEmitters(data.begin(), data.end()) >= parameters.Nmin) {
        
          for (data_type::const_iterator hit = data.begin(); hit != data.end(); ++hit) {
            HA[hit->getID()]->Fill(geometry.getIndex(hit->getString()), hit->getFloor(), 1.0);
          }
 
          // fit
        
          const auto result = katoomba(data.begin(), data.end());
 
          for (data_type::const_iterator hit = result.begin; hit != result.end; ++hit) {
            HB[hit->getID()]->Fill(geometry.getIndex(hit->getString()), hit->getFloor(), 1.0);
          }
 
          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;
            }
          }
 
          hn.Fill(log10(katoomba.gandalf.numberOfIterations));
          h0.Fill(result.chi2 / result.ndf);
 
          // store results
 
          if ((parameters.chi2perNDF > 0.0 && result.chi2 / result.ndf <= +parameters.chi2perNDF) ||
              (parameters.chi2perNDF < 0.0 && result.chi2 / result.ndf >= -parameters.chi2perNDF)) {
 
            const JEvt evt = getEvt(JHead(detector.getID(),
                                          result.getTimeRange(),
                                          data  .size(),
                                          result.size(),
                                          result.value.getN(),
                                          result.ndf,
                                          result.chi2),
                                    result.value);
 
            outputFile.put(evt);
 
            if (selection.is_valid<JEvent>()) {
 
              for (buffer_type::iterator i = p; i != q; ++i) {
 
                JEvent out(*i);                      // event with fitted times of emission
 
                const double toe = result.value.emission[JEKey(i->getID(), distance(p,i))].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 " << FIXED(12,3) << result.chi2 << '/' << FIXED(6,1) << result.ndf << endl);
 
            counter[1] += 1;
          }
 
        } else {
 
          WARNING(endl << "Event not accepted - minimum number of emitters " << getMinimumNumberOfEmitters(data.begin(), data.end()) << endl);
 
          counter[1] += 1;
        }
      }
    }
  }
  STATUS(endl);
 
  STATUS("Number of events written / rejected: " << counter[0] << " / " << counter[1] << endl);
 
  if (!squash) {
 
    JMultipleFileScanner<JMeta> io(inputFile);
 
    io >> outputFile;
  }
 
  outputFile.put(h0);
  outputFile.put(h1);
  outputFile.put(hn);
 
  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();
}