JSydney.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
#include <algorithm>
#include <limits>
#include <sys/stat.h>

#include <type_traits>
#include <functional>
#include <future>
#include <mutex>
#include <thread>
#include <vector>
#include <queue>

#include "TROOT.h"
#include "TFile.h"

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

#include "JSystem/JStat.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 "JDetector/JLocationRouter.hh"

#include "JFit/JGradient.hh"

#include "JTools/JHashMap.hh"
#include "JTools/JRange.hh"
#include "JTools/JQuantile.hh"

#include "JSupport/JMultipleFileScanner.hh"

#include "JAcoustics/JSoundVelocity.hh"
#include "JAcoustics/JEmitter.hh"
#include "JAcoustics/JAcousticsToolkit.hh"
#include "JAcoustics/JHit.hh"
#include "JAcoustics/JFitParameters.hh"
#include "JAcoustics/JKatoomba_t.hh"
#include "JAcoustics/JEvent.hh"
#include "JAcoustics/JSuperEvt.hh"
#include "JAcoustics/JSuperEvtToolkit.hh"
#include "JAcoustics/JSupport.hh"
#include "JAcoustics/JTransmission_t.hh"
#include "JAcoustics/JFremantle_t.hh"
#include "JAcoustics/JPlatypus_t.hh"

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


namespace JACOUSTICS {

  using JLANG::JValueOutOfRange;
  using JEEP::JContainer;
  using JFIT::JParameter_t;

  using namespace JDETECTOR;

  typedef JContainer< std::vector<JTripod> >       tripods_container;
  typedef JContainer< std::vector<JTransmitter> >  transmitters_container;
  typedef JContainer< std::vector<JHydrophone> >   hydrophones_container;
  typedef JContainer< std::set<JTransmission_t> >  disable_container;


  static const std::string initialise_t  =  "initialise";           //!< initialise
  static const std::string fix_t         =  "fix";                  //!< fix objects
  static const std::string string_t      =  "string";               //!< string
  static const std::string tripod_t      =  "tripod";               //!< tripod
  static const std::string stage_t       =  "stage";                //!< fit stage


  /**
   * Auxiliary data structure for handling file names.
   */
  struct JFilenames {
    std::string  detector;                                          //!< detector
    std::string  tripod;                                            //!< tripod
    std::string  hydrophone;                                        //!< hydrophone
    std::string  transmitter;                                       //!< transmitter
  };

    
  /**
   * Auxiliary data structure for setup of complete system.
   */
  struct JSetup {
    JDetector               detector;                               //!< detector
    tripods_container       tripods;                                //!< tripods
    struct :
      public hydrophones_container
    {
      /**
       * Check if there is a hydrophone on given string.
       *
       * \param  id                     string identifier
       * \return                        true if hydrophone present; else false
       */
      bool hasString(const int id) const
      {
        using namespace std;

        return (find_if(this->begin(), this->end(), make_predicate(&JHydrophone::getString, id)) != this->end());
      }
    }                       hydrophones;                            //!< hydrophones
    transmitters_container  transmitters;                           //!< transmitters
  };


  /**
   * Main class for pre-calibration using acoustics data.
   */
  struct JSydney {

    static constexpr double   RADIUS_M  =  1.0;                     //!< maximal horizontal distance between T-bar and emitter/hydrophone

    /**
     * List of identifiers.
     */
    struct ids_t :
      public std::set<int>
    {
      /**
       * Default constructor.
       */
      ids_t()
      {}


      /**
       * Copy constructor.
       *
       * \param  buffer                list of identifiers
       */
      ids_t(const std::vector<int>& buffer) :
        std::set<int>(buffer.begin(), buffer.end())
      {}


      /**
       * Difference constructor.
       * Make list of all elements in A that are not in B.
       *
       * \param  A                     list of identifiers
       * \param  B                     list of identifiers
       */
      ids_t(const ids_t& A,
            const ids_t& B)
      {
        std::set_difference(A.begin(), A.end(), B.begin(), B.end(), std::inserter(*this, this->begin()));
      }


      /**
       * Fix.
       * Keep list of elements that are not in B.
       *
       * \param  B                     list of identifiers
       */
      void fix(const ids_t& B)
      {
        ids_t A;

        this->swap(A);

        std::set_difference(A.begin(), A.end(), B.begin(), B.end(), std::inserter(*this, this->begin()));
      }


      /**
       * Read identifiers from input stream
       *
       * \param  in                    input stream
       * \param  object                identifiers
       * \return                       input stream
       */
      friend inline std::istream& operator>>(std::istream& in, ids_t& object)
      {
        for (int id; in >> id; ) {
          object.insert(id);
        }

        if (!in.bad()) {
          in.clear();
        }

        return in;
      }


      /**
       * Write identifiers to output stream
       *
       * \param  out                   output stream
       * \param  object                identifiers
       * \return                       output stream
       */
      friend inline std::ostream& operator<<(std::ostream& out, const ids_t& object)
      {
        for (const int id : object) {
          out << ' ' << id;
        }

        return out;
      }
    };


    /**
     * Auxiliary data structure for group of lists of identifiers of to-be-fitted objects.
     */
    struct fits_t {
      /**
       * Default constructor.
       */
      fits_t()
      {}


      /**
       * Initialise.
       *       
       * \param  setup                 setup
       */
      void initialise(const JSetup& setup)
      {
        using namespace JPP;

        strings       =  make_array(setup.detector    .begin(), setup.detector    .end(), &JModule     ::getString);
        tripods       =  make_array(setup.tripods     .begin(), setup.tripods     .end(), &JTripod     ::getID);
        hydrophones   =  make_array(setup.hydrophones .begin(), setup.hydrophones .end(), &JHydrophone ::getString);
        transmitters  =  make_array(setup.transmitters.begin(), setup.transmitters.end(), &JTransmitter::getString);
      }

      ids_t  strings;                                               //!< strings
      ids_t  tripods;                                               //!< tripods
      ids_t  hydrophones;                                           //!< hydrophones
      ids_t  transmitters;                                          //!< transmitters
    };


    /**
     * Auxiliary class to edit (z) position of module.
     */
    struct JModuleEditor :
      public JParameter_t
    {
      /**
       * Constructor.
       *
       * \param  module                module
       */
      JModuleEditor(JModule& module) :
        module(module)
      {}


      /**
       * Apply step.
       *
       * \param  step                  step
       */
      virtual void apply(const double step) override
      {
        using namespace JPP;

        module.add(JVector3Z_t * step);
      }

    private:
      JModule&  module;
    };


    /**
     * Auxiliary class to edit (x,y,z) position of string.
     */
    struct JStringEditor :
      public JParameter_t
    {
      /**
       * Constructor.
       *
       * The option <tt>true</tt> and <tt>false</tt> correspond to all modules and optical modules only, respectively.
       *
       * \param  setup                 setup
       * \param  id                    string number
       * \param  direction             direction
       * \param  option                option
       */
      JStringEditor(JSetup& setup, const int id, const JVector3D& direction, const bool option) :
        detector (setup.detector),
        direction(direction)
      {
        for (size_t i = 0; i != detector.size(); ++i) {
          if (detector[i].getString() == id && (detector[i].getFloor() != 0 || option)) {
            index.push_back(i);
          }
        }
      }


      /**
       * Apply step.
       *
       * \param  step                  step
       */
      virtual void apply(const double step) override
      {
        for (const auto i : index) {
          detector[i].add(direction * step);
        }
      }

    private:
      JDetector&           detector;
      JVector3D            direction;
      std::vector<size_t>  index;
    };


    /**
     * Auxiliary class to edit length of Dyneema ropes.
     */
    struct JDyneemaEditor :
      public JParameter_t
    {
      /**
       * Constructor.
       *
       * \param  setup                 setup
       * \param  id                    string number
       * \param  z0                    reference position
       */
      JDyneemaEditor(JSetup& setup, const int id, const double z0 = 0.0) :
        detector (setup.detector),
        z0       (z0)
      {
        for (size_t i = 0; i != detector.size(); ++i) {
          if (detector[i].getString() == id && detector[i].getFloor() != 0) {
            index.push_back(i);
          }
        }
      }


      /**
       * Apply step.
       *
       * \param  step              step
       */
      virtual void apply(const double step) override
      {
        for (const auto i : index) {

          JModule& module = detector[i];
          
          if      (step > 0.0)
            module.set(JVector3D(module.getX(), module.getY(), z0 + (module.getZ() - z0) * (1.0 + step)));
          else if (step < 0.0)
            module.set(JVector3D(module.getX(), module.getY(), z0 + (module.getZ() - z0) / (1.0 - step)));
        }
      }

    private:
      JDetector&           detector;
      double               z0;
      std::vector<size_t>  index;
    };


    /**
     * Auxiliary class to edit (x,y,z) position of tripod.
     */
    struct JTripodEditor :
      public JParameter_t
    {
      /**
       * Constructor.
       *
       * \param  setup                 setup
       * \param  id                    tripod identifier
       * \param  direction             direction
       */
      JTripodEditor(JSetup& setup, const int id, const JVector3D& direction) :
        tripods  (setup.tripods),
        direction(direction)
      {
        using namespace std;
        using namespace JPP;

        index = distance(tripods.begin(), find_if(tripods.begin(), tripods.end(), make_predicate(&JTripod::getID, id)));
      }


      /**
       * Apply step.
       *
       * \param  step                  step
       */
      virtual void apply(const double step) override
      {
        tripods[index].add(direction * step);
      }

    private:
      std::vector<JTripod>& tripods;
      JVector3D             direction;
      size_t                index;
    };


    /**
     * Auxiliary class to edit rotation of anchor.
     */
    struct JAnchorEditor :
      public JParameter_t
    {
      /**
       * Constructor.
       *
       * \param  setup                 setup
       * \param  id                    string identifier
       */
      JAnchorEditor(JSetup& setup, const int id) :
        hydrophones (setup.hydrophones),
        transmitters(setup.transmitters)
      {
        using namespace std;
        using namespace JPP;

        index[0] = distance(hydrophones .begin(), find_if(hydrophones .begin(), hydrophones .end(), make_predicate(&JHydrophone ::getString, id)));
        index[1] = distance(transmitters.begin(), find_if(transmitters.begin(), transmitters.end(), make_predicate(&JTransmitter::getString, id)));
      }


      /**
       * Apply step.
       *
       * \param  step                  step
       */
      virtual void apply(const double step) override
      {
        using namespace JPP;

        const JRotation3Z R(step);
          
        if (index[0] != hydrophones .size()) { hydrophones [index[0]].rotate(R); }
        if (index[1] != transmitters.size()) { transmitters[index[1]].rotate(R); }
      }

    private:
      std::vector<JHydrophone>&  hydrophones;
      std::vector<JTransmitter>& transmitters;
      size_t                     index[2];
    };


    /**
     * Extended data structure for parameters of stage.
     */
    struct JParameters_t :
      public JFitParameters
    {
      /**
       * Default constuctor.
       */
      JParameters_t() :
        Nmax   (std::numeric_limits<size_t>::max()),
        Nextra (0),
        epsilon(1.0e-4),
        debug  (3)
      {}


      /**
       * Read parameters from input stream
       *
       * \param  in                    input stream
       * \param  object                parameters
       * \return                       input stream
       */
      friend inline std::istream& operator>>(std::istream& in, JParameters_t& object)
      {
        object = JParameters_t();

        in >> object.option
           >> object.mestimator
           >> object.sigma_s
           >> object.stdev
           >> object.Nextra;

        if (in) {

          for (double value; in >> value; ) {
            object.steps.push_back(value);
          }

          if (!in.bad()) {
            in.clear();
          }
        }
          
        return in;
      }


      /**
       * Write parameters to output stream
       *
       * \param  out                   output stream
       * \param  object                parameters
       * \return                       output stream
       */
      friend inline std::ostream& operator<<(std::ostream& out, const JParameters_t& object)
      {
        using namespace std;

        out << setw(2)         << object.option      << ' '
            << setw(2)         << object.mestimator  << ' '
            << SCIENTIFIC(9,3) << object.sigma_s     << ' '
            << SCIENTIFIC(9,3) << object.stdev       << ' '
            << setw(3)         << object.Nextra;

        for (const double value : object.steps) {
          out << ' ' << FIXED(9,5) << value;
        }

        return out;
      }

      size_t  Nmax; 
      size_t  Nextra;
      double  epsilon;
      int     debug;
      std::vector<double> steps;
    };


    /**
     * Constructor.
     *
     * \param  filenames              file names
     * \param  V                      sound velocity
     * \param  jobs                   jobs
     * \param  debug                  debug
     */
    JSydney(const JFilenames&     filenames,
            const JSoundVelocity& V,
            const size_t          jobs,
            const int             debug) :
      filenames(filenames),
      V(V),
      jobs(jobs),
      debug(debug)
    {
      load(filenames.detector, setup.detector);

      setup.tripods.load(filenames.tripod.c_str());

      if (filenames.hydrophone  != "") { setup.hydrophones .load(filenames.hydrophone .c_str()); }
      if (filenames.transmitter != "") { setup.transmitters.load(filenames.transmitter.c_str()); }

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

      // detach PMTs

      detector = setup.detector;

      for (JDetector::iterator module = setup.detector.begin(); module != setup.detector.end(); ++module) {
        module->clear();
      }

      router.reset(new JLocationRouter(setup.detector));

      fits.initialise(setup);

      this->V.set(setup.detector.getUTMZ());                              // sound velocity at detector depth

      JKatoomba<JGandalf>::debug              = debug;
      JKatoomba<JGandalf>::MAXIMUM_ITERATIONS = 10000;

      JKatoomba<JEstimator>::MATRIX_INVERSION = LDU_t;

      ROOT::EnableThreadSafety();
    }


    /**
     * Fit procedure to determine the positions of tripods using strings that are fixed.
     *
     * \param  parameters            parameters
     */
    void stage_0(const JParameters_t& parameters)
    {
      using namespace std;
      using namespace JPP;

      this->parameters = parameters;

      JDetector buffer[2];

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

        if (fits.strings.count(module->getString()) != 0)
          buffer[0].push_back(*module);                            // new strings
        else
          buffer[1].push_back(*module);                            // old strings
      }

      setup.detector.swap(buffer[1]);

      JGradient fit(parameters.Nmax, parameters.Nextra, parameters.epsilon, parameters.debug);

      for (const int i : fits.tripods) {
        fit.push_back(JModifier_t(MAKE_STRING(LEFT(12) << "tripod.x: " << RIGHT(4) << i), new JTripodEditor(setup, i, JVector3X_t), parameters.steps[0]));
        fit.push_back(JModifier_t(MAKE_STRING(LEFT(12) << "tripod.y: " << RIGHT(4) << i), new JTripodEditor(setup, i, JVector3Y_t), parameters.steps[0]));
        fit.push_back(JModifier_t(MAKE_STRING(LEFT(12) << "tripod.z: " << RIGHT(4) << i), new JTripodEditor(setup, i, JVector3Z_t), parameters.steps[0]));
      }

      const double chi2 = fit(*this);

      copy(buffer[0].begin(), buffer[0].end(), back_inserter(setup.detector));

      STATUS("detector: " << FIXED(9,4) << chi2 << endl);
    }


    /**
     * Fit procedure to determine the positions of the strings and tripods.
     *
     * \param  parameters            parameters
     */
    void stage_a(const JParameters_t& parameters)
    {
      using namespace std;
      using namespace JPP;

      this->parameters = parameters;

      JGradient fit(parameters.Nmax, parameters.Nextra, parameters.epsilon, parameters.debug);

      for (const int i : fits.strings) {
        fit.push_back(JModifier_t(MAKE_STRING(LEFT(12) << "string.x: " << RIGHT(4) << i), new JStringEditor(setup, i, JVector3X_t, true),  parameters.steps[0]));
        fit.push_back(JModifier_t(MAKE_STRING(LEFT(12) << "string.y: " << RIGHT(4) << i), new JStringEditor(setup, i, JVector3Y_t, true),  parameters.steps[0]));
        fit.push_back(JModifier_t(MAKE_STRING(LEFT(12) << "string.z: " << RIGHT(4) << i), new JStringEditor(setup, i, JVector3Z_t, false), parameters.steps[0]));
      }

      for (const int i : fits.tripods) {
        fit.push_back(JModifier_t(MAKE_STRING(LEFT(12) << "tripod.x: " << RIGHT(4) << i), new JTripodEditor(setup, i, JVector3X_t), parameters.steps[1]));
        fit.push_back(JModifier_t(MAKE_STRING(LEFT(12) << "tripod.y: " << RIGHT(4) << i), new JTripodEditor(setup, i, JVector3Y_t), parameters.steps[1]));
        fit.push_back(JModifier_t(MAKE_STRING(LEFT(12) << "tripod.z: " << RIGHT(4) << i), new JTripodEditor(setup, i, JVector3Z_t), parameters.steps[1]));
      }

      for (const int i : ids_t(fits.hydrophones, fits.transmitters)) {
        fit.push_back(JModifier_t(MAKE_STRING(LEFT(12) << "anchor.R: " << RIGHT(4) << i), new JAnchorEditor(setup, i), parameters.steps[0] / RADIUS_M));
      }

      for (const int i : fits.transmitters) {

        JModule& module = setup.detector.getModule(router->getAddress(JLocation(i,0)));

        fit.push_back(JModifier_t(MAKE_STRING(LEFT(12) << "anchor.R: " << RIGHT(4) << i), new JAnchorEditor(setup, i), parameters.steps[0] / RADIUS_M));
        fit.push_back(JModifier_t(MAKE_STRING(LEFT(12) << "anchor.z: " << RIGHT(4) << i), new JModuleEditor(module),   parameters.steps[0]));
      }

      const double chi2 = fit(*this);

      STATUS("detector: " << FIXED(9,4) << chi2 << endl);
    }


    /**
     * Fit procedure to determine the stretching and z-positions of individual strings.
     *
     * \param  parameters            parameters
     */
    void stage_b(const JParameters_t& parameters)
    {
      using namespace std;
      using namespace JPP;

      this->parameters = parameters;

      map<int, JDetector> buffer;
      
      double z0 = 0.0;
      
      for (JDetector::iterator module = setup.detector.begin(); module != setup.detector.end(); ++module) {

        buffer[module->getString()].push_back(*module);
        
        if (module->getZ() > z0) {
          z0 = module->getZ();
        }
      }

      setup.detector.clear();      

      for (const int i : fits.strings) {

        setup.detector.swap(buffer[i]);

        router.reset(new JLocationRouter(setup.detector));

        JGradient fit(parameters.Nmax, parameters.Nextra, parameters.epsilon, parameters.debug);

        fit.push_back(JModifier_t(MAKE_STRING(LEFT(12) << "string.M: " << RIGHT(4) << i), new JDyneemaEditor(setup, i, z0),                 parameters.steps[0]));
        fit.push_back(JModifier_t(MAKE_STRING(LEFT(12) << "string.z: " << RIGHT(4) << i), new JStringEditor (setup, i, JVector3Z_t, false), parameters.steps[1]));

        const double chi2 = fit(*this);

        STATUS("string: " << setw(4) << i << ' ' << FIXED(9,4) << chi2 << endl);

        setup.detector.swap(buffer[i]);
      }

      setup.detector.clear();      
      
      for (const auto& detector : buffer) {
        copy(detector.second.begin(), detector.second.end(), back_inserter(setup.detector));
      }

      router.reset(new JLocationRouter(setup.detector));
    }


    /**
     * Fit procedure to determine the z-positions of the modules.
     *
     * \param  parameters            parameters
     */
    void stage_c(const JParameters_t& parameters)
    {
      using namespace std;
      using namespace JPP;

      this->parameters = parameters;

      JGradient fit(parameters.Nmax, parameters.Nextra, parameters.epsilon, parameters.debug);

      for (JDetector::iterator module = setup.detector.begin(); module != setup.detector.end(); ++module) {
        if (fits.strings.count(module->getString()) != 0 && module->getFloor() != 0) {
          fit.push_back(JModifier_t(MAKE_STRING(LEFT(12) << "module.z: " << right << module->getLocation()), new JModuleEditor(*module), parameters.steps[0]));
        }
      }

      const double chi2 = fit(*this);

      STATUS("detector: " << FIXED(9,4) << chi2 << endl);
    }


    /**
     * Get chi2.
     *
     * \param  option                option
     * \return                       chi2/NDF
     */
    double operator()(const int option) const
    {
      using namespace std;
      using namespace JPP;

      const JGeometry geometry(setup.detector, setup.hydrophones);

      JHashMap<int, JEmitter>    emitters;

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

      for (transmitters_container::const_iterator i = setup.transmitters.begin(); i != setup.transmitters.end(); ++i) {
        try {
          emitters[i->getID()]  = JEmitter(i->getID(), i->getPosition()    + router->getModule(i->getLocation()).getPosition());
        }
        catch(const exception&) {}                                        // if no module available, discard transmitter
      }

      if        (option == 0 ||                                           // step wise improvement of the chi2
                 option == 1) {                                           // evaluation of the chi2 before the determination of the gradient of the chi2

        JSTDObjectWriter<JSuperEvt>  out(this->odat);                     // storage of data for subsequent use

        this->odat.clear();

        JFremantle::output = (option == 1 ? &out : NULL);

        {
          JFremantle fremantle(geometry, V, parameters, jobs);

          for (idat_type::const_iterator superevt = idat.begin(); superevt != idat.end(); ++superevt) {

            const JWeight getWeight(superevt->begin(), superevt->end());

            set<int>      counter;
            vector<JHit>  data;

            for (ievt_type::const_iterator evt = superevt->begin(); evt != superevt->end(); ++evt) {

              const JEmitter& emitter = emitters [evt->getID()];
              const double    weight  = getWeight(evt->getID());

              for (JEvent::const_iterator i = evt->begin(); i != evt->end(); ++i) {

                if (geometry.hasLocation(receivers[i->getID()])) {

                  data.push_back(JHit(emitter,
                                      distance(superevt->begin(), evt),
                                      receivers[i->getID()],
                                      i->getToA(),
                                      parameters.sigma_s,
                                      weight));

                  counter.insert(evt->getID());
                }
              }
            }

            if (counter.size() >= parameters.Nmin) {
              fremantle.enqueue(data);
            }
          }
        }

        return JFremantle::Q.getMean(numeric_limits<float>::max());

      } else if (option == 2) {                                           // evaluation of the derivative of the chi2 to each fit parameter

        {
          JPlatypus platypus(geometry, emitters, V, parameters, jobs);

          for (odat_type::const_iterator superevt = odat.begin(); superevt != odat.end(); ++superevt) {
            platypus.enqueue(*superevt);
          }
        }

        return JPlatypus::Q.getMean(numeric_limits<float>::max());

      } else {

        return numeric_limits<float>::max();
      }
    }


    /**
     * Run.
     *
     * \param  script                steering script
     */
    void run(const std::string& script)
    {
      using namespace std;
      using namespace JPP;

      ifstream in(script.c_str());

      while (in) {

        string key, value;

        if (in >> key) {

          if        (key == initialise_t) {            // set object identifiers

            fits.initialise(setup);
            
          } else if (key == fix_t) {                   // fix object identifiers

            getline(in, value);

            istringstream is(value);

            string type;                               // type of object
            ids_t  id;                                 // identifiers
        
            if (is >> type >> id) {
              if        (type == string_t) {
                fits.strings     .fix(id);
                fits.hydrophones .fix(id);
                fits.transmitters.fix(id);
              } else if (type == tripod_t) {
                fits.tripods     .fix(id);
              } else {
                THROW(JValueOutOfRange, "Invalid type <" << type << ">");
              }
            }

          } else if (key == stage_t) {                 // stage

            getline(in, value);

            istringstream is(value);

            string        stage;
            JParameters_t input;

            if (is >> stage >> input) {

              STATUS("stage " << setw(3) << stage << " {" << input << "}" << endl);

              JTimer timer;

              timer.start();

              ofstream out(MAKE_CSTRING("stage-" << stage << ".log"));

              {
                JRedirectStream redirect(cout, out);
              
                switch (stage[stage.size() - 1]) {

                case '0':
                  stage_0(input);
                  break;

                case 'a':
                case 'A':
                  stage_a(input);
                  break;

                case 'b':
                case 'B':
                  stage_b(input);
                  break;

                case 'c':
                case 'C':
                  stage_c(input);
                  break;

                default:
                  THROW(JValueOutOfRange, "Invalid stage <" << stage << ">");
                  break;
                }
              }
            
              out.close();

              store(stage);
              store();

              timer.stop();

              STATUS("Elapsed time " << FIXED(12,3) << timer.usec_wall * 1.0e-6 << " s." << endl);
            }

          } else {
            THROW(JValueOutOfRange, "Invalid key <" << key << ">");
          }
        }
      }

      in.close();
    }


    /**
     * Store data in given directory.
     *
     * \param  dir                    directory
     */
    void store(const std::string& dir = ".")
    {
      using namespace JPP;

      if (getFileStatus(dir.c_str()) || (mkdir(dir.c_str(), S_IRWXU | S_IRWXG) != -1)) {

        // attach PMTs

        for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
          module->set(router->getModule(module->getLocation()).getPosition());
        }

        JDETECTOR::store(getFilename(dir, filenames.detector), detector);

        setup.tripods.store(getFilename(dir, filenames.tripod).c_str());

        if (filenames.hydrophone  != "") { setup.hydrophones .store(getFilename(dir, filenames.hydrophone) .c_str()); }
        if (filenames.transmitter != "") { setup.transmitters.store(getFilename(dir, filenames.transmitter).c_str()); }

      } else {

        THROW(JValueOutOfRange, "Invalid directory <" << dir << ">");
      }
    }


    JFilenames       filenames;      
    JSoundVelocity   V;
    size_t           jobs;
    int              debug;
    JSetup           setup;
    fits_t           fits;
    
    JTOOLS::JHashMap<int, JLocation>    receivers;
    std::unique_ptr<JLocationRouter>    router;

    typedef std::vector<JEvent>         ievt_type;
    typedef std::vector<ievt_type>      idat_type;

    idat_type                           idat;

  private:
    typedef std::vector<JSuperEvt>      odat_type;
    mutable odat_type                   odat;

    JDetector        detector;                         //!< PMTs
    JFitParameters   parameters;
  };
}


/**
 * \file
 *
 * Application to perform acoustic pre-calibration.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  JMultipleFileScanner<JEvent>  inputFile;
  JLimit_t&               numberOfEvents = inputFile.getLimit();
  JFilenames              filenames;                            // file names
  JFitParameters          parameters;                           // fit parameters
  string                  script;                               // script file
  JSoundVelocity          V              = getSoundVelocity;    // default sound velocity
  disable_container       disable;                              // disable tansmissions
  size_t                  jobs;                                 // number of parallel jobs
  int                     debug;

  try {

    JParser<> zap("Application to perform acoustic pre-calibration.");

    zap['f'] = make_field(inputFile,              "output of JAcousticEventBuilder[.sh]");
    zap['n'] = make_field(numberOfEvents)                                           = JLimit::max();
    zap['a'] = make_field(filenames.detector);
    zap['@'] = make_field(parameters)                                               = JPARSER::initialised();
    zap['s'] = make_field(script,                 "steering script");
    zap['V'] = make_field(V,                      "sound velocity")                 = JPARSER::initialised();
    zap['T'] = make_field(filenames.tripod,       "tripod file");
    zap['Y'] = make_field(filenames.transmitter,  "transmitter file")               = JPARSER::initialised();
    zap['H'] = make_field(filenames.hydrophone,   "hydrophone file")                = JPARSER::initialised();
    zap['M'] = make_field(getMechanics,           "mechanics data")                 = JPARSER::initialised();
    zap['!'] = make_field(disable,                "disable transmission")           = JPARSER::initialised();
    zap['N'] = make_field(jobs,                   "number of parallel jobs")        = 1;
    zap['d'] = make_field(debug)                                                    = 1;

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

  if (jobs == 0) {
    FATAL("Invalid number of jobs " << jobs << endl);
  }

  JSydney   sydney(filenames, V, jobs, debug);

  set<int>  receivers;
  set<int>  emitters;

  for (JDetector::const_iterator              i = sydney.setup.detector    .begin(); i != sydney.setup.detector    .end(); ++i) {
    if (i->getFloor() != 0 || sydney.setup.hydrophones.hasString(i->getString())) {
      receivers.insert(i->getID());
    }
  }
  
  for (tripods_container::const_iterator      i = sydney.setup.tripods     .begin(); i != sydney.setup.tripods     .end(); ++i) {
    {
      emitters.insert(i->getID());
    }
  }

  for (transmitters_container::const_iterator i = sydney.setup.transmitters.begin(); i != sydney.setup.transmitters.end(); ++i) {
    if (sydney.setup.detector.hasModule(i->getLocation())) {
      emitters.insert(i->getID());
    }
  }

  typedef vector<JEvent> buffer_type;

  buffer_type zbuf;
  
  while (inputFile.hasNext()) {

    const JEvent* evt = inputFile.next();

    if (emitters.count(evt->getID())) {
      zbuf.push_back(*evt);
    }
  }

  const JRange<double> unit(0.0, 1.0);

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

    JEvent::overlap(p, q, parameters.deadTime_s);  // empty overlapping events

    JSydney::ievt_type buffer;

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

      for (JEvent::iterator i = evt->begin(); i != __end; ) {

        if (disable.count(JTransmission_t(evt->getID(), i->getID())) == 0 &&
            disable.count(JTransmission_t(-1,           i->getID())) == 0) {

          if (receivers.count(i->getID()) && i->getQ() >= parameters.Qmin * (unit(parameters.Qmin) ? i->getW() : 1.0)) {
            ++i; continue;
          }
        }

        iter_swap(i, --__end);
      }

      buffer.push_back(JEvent(evt->getOID(), buffer.size(), evt->getID(), evt->begin(), __end));
    }

    if (getNumberOfEmitters(buffer.begin(), buffer.end()) >= parameters.Nmin) {
      sydney.idat.push_back(buffer);
    }
  }

  sydney.run(script);
}