JPerth.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <type_traits>
#include <functional>
#include <future>
#include <mutex>
#include <thread>
#include <vector>
#include <set>
#include <memory>
#include <algorithm>
 
#include "km3net-dataformat/online/JDAQEvent.hh"
#include "km3net-dataformat/definitions/fitparameters.hh"
#include "JDAQ/JDAQEventIO.hh"
 
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JPMTParametersMap.hh"
#include "JDetector/JCalibration.hh"
 
#include "JDynamics/JDynamics.hh"
 
#include "JSupport/JMeta.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JParallelFileScanner.hh"
#include "JSupport/JSummaryFileRouter.hh"
 
#include "JReconstruction/JHitW0.hh"
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JReconstruction/JMuonGandalfParameters_t.hh"
 
#include "JTrigger/JHitL0.hh"
#include "JTrigger/JBuildL0.hh"
 
#include "JFit/JFitToolkit.hh"
#include "JFit/JLine1Z.hh"
#include "JFit/JLine3Z.hh"
#include "JFit/JModel.hh"
#include "JFit/JGandalf.hh"
#include "JFit/JLine3ZRegressor.hh"
#include "JFit/JGradient.hh"
 
#include "JLang/JSTDObjectReader.hh"
#include "JLang/JVectorize.hh"
 
#include "JGeometry3D/JRotation3D.hh"
 
#include "JMath/JMath.hh"
#include "JMath/JQuantile_t.hh"
#include "JMath/JMatrix2S.hh"
 
#include "JTools/JRange.hh"
 
#include "Jeep/JProperties.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
namespace JRECONSTRUCTION {
 
  using KM3NETDAQ::JDAQEvent;
  using JFIT::JParameter_t;
  using JFIT::JLine1Z;
  using JFIT::JLine3Z;
  using JDETECTOR::JDetector;
  using JTRIGGER::JHit;
  using JTRIGGER::getTimeSlewing;
 
  typedef std::vector<JHitW0>                  buffer_type;  //!< hits
  typedef std::map<int, buffer_type>           map_type;     //!< identifier -> hits
 
 
  /**
   * Get list of time-over-threshold identifiers.
   * The time-over-threshold identifier corresponds to the index of the tabulated time slewing correction.
   *
   * \return                       identifiers
   */
  inline std::set<int> getTimeOverThresholdIDs()
  {
    std::set<int> buffer;
 
    for (size_t i = 0; i != getTimeSlewing.size(); ++i) {
      buffer.insert(i);
    }
 
    return buffer;
  }
 
 
  /**
   * Auxiliary data structure to store data and fit in memory.
   */
  struct event_type {
    map_type  data;
    JLine1Z   value;
    bool      status  =  true;
  };
    
  typedef std::vector<event_type>                    data_type;
  typedef JLANG::JSTDObjectReader<const event_type>  input_type;
 
 
  /**
   * Auxiliary class for editing time offset.
   */
  struct JEditor :
    public JParameter_t
  {
    /**
     * Constructor.
     *
     * \param  data                  data
     * \param  id                    identifier
     */
    JEditor(data_type& data, const int id) :
      data(data),
      id(id),
      t0(0.0)
    {}
    
 
    /**
     * Apply step.
     *
     * \param  step                  step
     */
    virtual void apply(const double step) override
    {
      for (data_type::iterator evt = data.begin(); evt != data.end(); ++evt) {
 
        map_type::iterator p = evt->data.find(id);
 
        if ((evt->status = (p != evt->data.end()))) {
          for (buffer_type::iterator hit = p->second.begin(); hit != p->second.end(); ++hit) {
            static_cast<JHit&>(*hit) = JHit(hit->getT1() + step, hit->getToT());
          }
        }
      }
 
      this->t0 += step;
    }
 
    data_type& data;         //!< data
    int        id;           //!< identifier
    double     t0;           //!< time offset [ns]
  };
 
 
  typedef JFIT::JRegressorStorage<JFIT::JLine3Z, JFIT::JGandalf>   JRegressorStorage_t;
  typedef JFIT::JRegressor       <JFIT::JLine3Z, JFIT::JGandalf>   JRegressor_t;
 
 
  /**
   * Thread pool for fits to data.
   */
  struct JPerth {
    /**
     * Constructor.
     *
     * \param  storage            storage for PDFs
     * \param  data               data
     * \param  ns                 number of threads
     * \param  option             option
     */
    JPerth(const JRegressorStorage_t& storage,
           const data_type&           data,
           const size_t               ns,
           const int                  option) :
      input(data),
      stop(false)
    {
      using namespace std;
      using namespace JPP;
 
      Q.reset();
 
      for (size_t i = 0; i < ns; ++i) {
 
        thread worker([this, storage, option]() {
 
          JRegressor_t regressor(storage);
 
          regressor.parameters.resize(5);
 
          regressor.parameters[0] = JLine3Z::pX();
          regressor.parameters[1] = JLine3Z::pY();
          regressor.parameters[2] = JLine3Z::pT();
          regressor.parameters[3] = JLine3Z::pDX();
          regressor.parameters[4] = JLine3Z::pDY();
 
          buffer_type  data;
          JLine3Z      value;
          
          for ( ; ; ) {
 
            {
              unique_lock<mutex> lock(in);
 
              cv.wait(lock, [this]() { return stop || this->input.hasNext(); });
 
              if (stop && !this->input.hasNext()) {
                return;
              }
 
              const event_type* evt = this->input.next();
 
              // re-organise data
 
              data.clear();
 
              // 2 => evaluation of the derivative of the chi2 to each fit parameter.
 
              if (option != 2 || evt->status) {
                for (map_type::const_iterator p = evt->data.begin(); p != evt->data.end(); ++p) {
                  copy(p->second.begin(), p->second.end(), back_inserter(data));
                }
              }
 
              value = evt->value;
            }
 
            if (!data.empty()) {
 
              const double chi2 = regressor(value, data.begin(), data.end());
 
              {
                unique_lock<mutex> lock(out);
 
                Q.put(chi2);
              }
            }
          }
        });
 
        workers.emplace_back(std::move(worker));
      }
    }
 
 
    /**
     * Destructor.
     */   
    ~JPerth()
    {
      using namespace std;
 
      {
        unique_lock<mutex> lock(in);
 
        stop = true;
      }
      
      cv.notify_all();
 
      for (auto& worker : workers) {
        worker.join();
      }
    }
    
    static JMATH::JQuantile_t Q;
    
  private:
    std::vector<std::thread>  workers;
    input_type                input;
    std::mutex                in;
    std::mutex                out;
    std::condition_variable   cv;
    bool                      stop;    
  };
 
  /**
   */
  JMATH::JQuantile_t JPerth::Q;
 
 
  /**
   * Auxiliary data structure for chi2 function object.
   */
  struct JPerth_t {
    /**
     * Get chi2.
     *
     * \param  option                option
     * \return                       chi2/NDF
     */
    double operator()(const int option) const
    {
      JPerth(storage, data, ns, option);
 
      return JPerth::Q.getMean(0.0);
    }
 
    const JRegressorStorage_t& storage;
    const data_type&           data;
    const size_t               ns;
  };
}
 
 
/**
 * \file
 *
 * Program to determine string or optical module time calibrations.
 *
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
 
  typedef JParallelFileScanner< JTypeList<JDAQEvent, JFIT::JEvt> > JParallelFileScanner_t;
  typedef JParallelFileScanner_t::multi_pointer_type               multi_pointer_type;
  typedef JMultipleFileScanner<calibration_types>                  JCalibration_t;
 
  typedef JRange<double>   range_type;
 
 
  JMultipleFileScanner<>   inputFile;
  JLimit_t&                numberOfEvents = inputFile.getLimit();
  string                   detectorFile;
  JCalibration_t           calibrationFile;
  double                   Tmax_s;
  string                   pdfFile;
  JMuonGandalfParameters_t parameters;
  bool                     overwriteDetector;
  set<int>                 strings;                              // string identifiers
  set<int>                 modules;                              // module identifiers
  set<int>                 indices;                              // time-over-threshold identifiers
  range_type               X;                                    // time-over-threshold values
  int                      number_of_iterations    =  1000;
  int                      number_of_extra_steps   =     0;
  double                   epsilon                 =   1.0e-4;
  double                   T_ns                    =   2.5;      // step size
  JPMTParametersMap        pmtParameters;
  size_t                   threads;                              // number of threads
  int                      debug;
 
  parameters.ZMin_m = -1.0e4;                                    // set range hit selection
  parameters.ZMax_m = +1.0e4;
 
  const int DEFAULT_ID = -1;
 
  try { 
 
    JProperties properties;
 
    properties.insert(gmake_property(number_of_iterations));
    properties.insert(gmake_property(number_of_extra_steps));
    properties.insert(gmake_property(epsilon));
    properties.insert(gmake_property(T_ns));
 
    JParser<> zap("Program to determine string or optical module time calibrations.");
    
    zap['f'] = make_field(inputFile);
    zap['a'] = make_field(detectorFile);
    zap['+'] = make_field(calibrationFile)     = JPARSER::initialised();
    zap['T'] = make_field(Tmax_s)              = 100.0;
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['F'] = make_field(pdfFile);
    zap['P'] = make_field(pmtParameters,       "PMT simulation data (or corresponding file name)")        = JPARSER::initialised();
    zap['@'] = make_field(parameters)          = JPARSER::initialised();
    zap['A'] = make_field(overwriteDetector, "overwrite detector file provided through '-a' with fitted time offsets.");
    zap['S'] = make_field(strings,           "string identifiers          (" << DEFAULT_ID << " => all)") = JPARSER::initialised();
    zap['M'] = make_field(modules,           "module identifiers          (" << DEFAULT_ID << " => all)") = JPARSER::initialised();
    zap['I'] = make_field(indices,           "time-over-threshold indices (" << DEFAULT_ID << " => all)") = JPARSER::initialised();
    zap['x'] = make_field(X,                 "fit range for time slewing correction")                     = range_type(4.0, 175.0);
    zap['%'] = make_field(properties,        "fit parameters")                                            = JPARSER::initialised();
    zap['N'] = make_field(threads,           "number of threads")                                         = 1;
    zap['d'] = make_field(debug)               = 1;
    
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
 
  if ((strings.empty() ? 0 : 1) +
      (modules.empty() ? 0 : 1) +
      (indices.empty() ? 0 : 1) != 1) {
    FATAL("Set either strings (option -S) or modules (option -M) or indices (option -I)." << endl);
  }
 
 
  typedef JDAQHit::JTOT_t                      JTOT_t;
 
  const JRange<double> range(0.0, numeric_limits<JTOT_t>::max());
 
  JDetector detector;
 
  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }
 
  unique_ptr<JDynamics> dynamics;
 
  if (!calibrationFile.empty()) {
 
    try {
 
      dynamics.reset(new JDynamics(detector, Tmax_s));
 
      dynamics->load(calibrationFile);
    }
    catch(const exception& error) {
      FATAL(error.what());
    }
  }
 
  const JModuleRouter router(dynamics ? dynamics->getDetector() : detector);
 
  NOTICE("Reading PDFs... " << flush);
 
  const JRegressorStorage_t storage(pdfFile, JTimeRange(parameters.TMin_ns, parameters.TMax_ns), parameters.TTS_ns);
 
  NOTICE("OK" << endl);
 
  JRegressor_t::debug              = debug;
  JRegressor_t::Vmax_npe           = parameters.VMax_npe;
  JRegressor_t::MAXIMUM_ITERATIONS = parameters.NMax;
 
  // preprocess input data
  
  NOTICE("Reading data" << endl);
 
  const JBuildL0<JHitL0>  buildL0;
 
  data_type    data;
 
  counter_type skip    = inputFile.getLowerLimit();
  counter_type counter = inputFile.getLowerLimit();
  
  for (JMultipleFileScanner<>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
 
    JSummaryFileRouter  summary(*i);
 
    for (JParallelFileScanner_t in(*i); (skip -= in.skip(skip)) == 0 && in.hasNext() && counter != numberOfEvents; ++counter) {
 
      STATUS("event: " << setw(10) << counter << '\r'); DEBUG(endl);
 
      multi_pointer_type ps  = in.next();
    
      const JDAQEvent*   tev = ps;
      JFIT::JEvt*        evt = ps;
 
      summary.update(*tev);
 
      if (dynamics) {
        dynamics->update(*tev);
      }
 
      JFIT::JEvt::iterator __end = partition(evt->begin(), evt->end(), JHistory::is_application(JMUONGANDALF));
 
      if (evt->begin() != __end) {
 
        sort(evt->begin(), __end, qualitySorter);
 
        vector<JHitL0>  dataL0;
              
        buildL0(*tev, router, true, back_inserter(dataL0));
 
        const JFIT::JFit&  fit = (*evt)[0];
      
        const JRotation3D  R (getDirection(fit));
        const JLine1Z      tz(getPosition (fit).rotate(R), fit.getT());
        JRange<double>     Z_m;
 
        if (fit.getW(JSTART_LENGTH_METRES, 0.0) > 0.0) {
          Z_m = JZRange(fit.getW(JSTART_ZMIN_M) + parameters.ZMin_m,
                        fit.getW(JSTART_ZMAX_M) + parameters.ZMax_m);
        }
 
        const JFIT::JModel<JLine1Z> match(tz, parameters.roadWidth_m, JTimeRange(parameters.TMin_ns, parameters.TMax_ns), Z_m);
 
        // hit selection based on start value
 
        buffer_type buffer;
 
        for (vector<JHitL0>::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
 
          const JPMTIdentifier id(i->getModuleID(), i->getPMTAddress());
 
          const JPMTParameters& wip = pmtParameters.getPMTParameters(id);
 
          const int    type = wip.getType();
          const double QE   = wip.QE;
          const double R_Hz = summary.getRate(i->getPMTIdentifier(), parameters.R_Hz);
 
          JHitW0 hit(*i, type, QE, R_Hz);
 
          hit.rotate(R);
 
          if (match(hit)) {
            buffer.push_back(hit);
          }
        }
 
        // select first hit
 
        sort(buffer.begin(), buffer.end(), JHitL0::compare);
 
        buffer_type::const_iterator __end = unique(buffer.begin(), buffer.end(), equal_to<JDAQPMTIdentifier>());
 
        // re-organise data
        
        map_type map;
 
        for (buffer_type::const_iterator hit = buffer.begin(); hit != __end; ++hit) {
          
          const JModule& module = router.getModule(hit->getModuleID());
          const JTOT_t   index  = (JTOT_t) range.constrain(hit->getToT());
 
          if (!strings.empty()) { map[module.getString()].push_back(*hit); }
          if (!modules.empty()) { map[module.getID()]    .push_back(*hit); }
          if (!indices.empty()) { map[index]             .push_back(*hit); }
        }
        
        data.push_back({map, tz, true});
      }
    }
  }
  STATUS(endl);
  NOTICE("OK" << endl);
 
 
  // fit
 
  JGradient fit(number_of_iterations, number_of_extra_steps, epsilon, 3);
 
  if (strings.count(DEFAULT_ID) != 0) { strings = getStringIDs(detector); }
  if (modules.count(DEFAULT_ID) != 0) { modules = getModuleIDs(detector); }
  if (indices.count(DEFAULT_ID) != 0) { indices = getTimeOverThresholdIDs(); }
 
  for (int id : strings) { fit.push_back(JModifier_t(MAKE_STRING("string: " << FILL(4,'0') << id << FILL()), new JEditor(data, id), T_ns)); }
  for (int id : modules) { fit.push_back(JModifier_t(MAKE_STRING("module: " << setw(8)     << id),           new JEditor(data, id), T_ns)); }
  for (int id : indices) { fit.push_back(JModifier_t(MAKE_STRING("index:  " << setw(3)     << id),           new JEditor(data, id), T_ns)); }
 
  const JPerth_t perth = {storage, data, threads};
 
  const double   chi2  = fit(perth);
 
  STATUS("result: " << FIXED(12,6) << chi2 << ' ' << setw(6) << fit.numberOfIterations << endl);
 
  for (size_t i = 0; i != fit.size(); ++i) {
    {
      const JEditor* p = dynamic_cast<const JEditor*>(fit[i].get());
 
      if (p != NULL) { STATUS(fit[i].name << ' ' << FIXED(9,3) << p->t0 << " [ns]" << endl); }
    }
  }
 
 
  if (!indices.empty()) {
 
    // solve y = ax + b
 
    JMatrix2S V;
    double    Y[2] = { 0.0, 0.0 };
 
    for (size_t i = 0; i != fit.size(); ++i) {
 
      const JEditor* p = dynamic_cast<const JEditor*>(fit[i].get());
 
      if (p != NULL) {
 
        const double x = p->id;
        const double y = p->t0;
 
        if (X(x)) {
 
          V.a00 += x*x;  V.a01 += x;
          V.a10 += x;    V.a11 += 1;
 
          Y[0] += x*y;
          Y[1] += y;
        }
      }
    }
 
    try {
 
      V.invert();
 
      const double a = V.a00 * Y[0] + V.a01 * Y[1];
      //const double b = V.a10 * Y[0] + V.a11 * Y[1];
 
      cout << "// Produced by JPerth.cc; to be included in JTimeSlewing.hh" << endl;
 
      for (size_t i = 0; i != getTimeSlewing.size(); ++i) {
 
        const double x = i;
        const double y = a * (X.constrain(x) - TIME_OVER_THRESHOLD_NS);
 
        cout << "this->push_back(" << FIXED(6,2) << getTimeSlewing[i] - y << ");"  << endl;
      }
    }
    catch(const exception& error) {
      FATAL(error.what());
    }
  }
 
  if (overwriteDetector) {
 
    detector.comment.add(JMeta(argc, argv));
 
    map<int, double> calibration;
 
    for (size_t i = 0; i != fit.size(); ++i) {
 
      const JEditor* p = dynamic_cast<const JEditor*>(fit[i].get());
 
      if (p != NULL) {
        calibration[p->id] = p->t0;
      }
    }
 
    const double t0 = getAverage(get_values(calibration), 0.0);
 
    for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
 
      if (!module->empty()) {
 
        map<int, double>::const_iterator p = (!strings.empty() ? calibration.find(module->getString()) :
                                              !modules.empty() ? calibration.find(module->getID())     :
                                              calibration.end());
        if (p != calibration.end()) {
          for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
            module->getPMT(pmt).addT0(p->second - t0);
          }
        }
      }
    }
 
    NOTICE("Store calibration data on file " << detectorFile << endl);
 
    store(detectorFile, detector);
  }
}