JPerth.cc File Reference

Program to determine inter-string time calibration. More...

#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 "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/JRegressorHelper.hh"
#include "JFit/JLine3ZRegressor.hh"
#include "JFit/JGradient.hh"
#include "JLang/JSTDObjectReader.hh"
#include "JGeometry3D/JRotation3D.hh"
#include "JTools/JRange.hh"
#include "JTools/JQuantile.hh"
#include "Jeep/JProperties.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Classes
struct	JRECONSTRUCTION::event_type
	Auxiliary data structure to store data and fit in memory. More...
struct	JRECONSTRUCTION::JStringEditor
	Auxiliary class to edit time offset of data per string. More...
struct	JRECONSTRUCTION::JPerth
	Thread pool for fits to data. More...
struct	JRECONSTRUCTION::JPerth_t
	Auxiliary data structure for chi2 function object. More...

Namespaces
	JRECONSTRUCTION
	Model fits to data.

Typedefs
typedef std::vector< JHitW0 >	JRECONSTRUCTION::buffer_type
	hits More...
typedef std::map< int, buffer_type >	JRECONSTRUCTION::map_type
	string -> hits More...
typedef std::vector< event_type >	JRECONSTRUCTION::data_type
typedef JLANG::JSTDObjectReader< const event_type >	JRECONSTRUCTION::input_type
typedef JFIT::JRegressorStorage < JFIT::JLine3Z, JFIT::JGandalf >	JRECONSTRUCTION::JRegressorStorage_t
typedef JFIT::JRegressor < JFIT::JLine3Z, JFIT::JGandalf >	JRECONSTRUCTION::JRegressor_t

Functions
int	main (int argc, char **argv)

Detailed Description

Program to determine inter-string time calibration.

Author

mdejong

Definition in file JPerth.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 277 of file JPerth.cc.

{
  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 JTYPELIST<JCOMPASS::JOrientation,
                    JACOUSTICS::JEvt>::typelist                    calibration_types;
  typedef JMultipleFileScanner<calibration_types>                  JCalibration_t;

  JMultipleFileScanner<>   inputFile;
  JLimit_t&                numberOfEvents = inputFile.getLimit();
  string                   detectorFile;
  JCalibration_t           calibrationFile;
  double                   Tmax_s;
  string                   pdfFile;
  JMuonGandalfParameters_t parameters;
  bool                     overwriteDetector;
  set<int>                 ids;                                  // string identifiers
  int                      number_of_iterations    =  1000;
  int                      number_of_extra_steps   =     0;
  double                   epsilon                 =   1.0e-4;
  double                   T_ns                    =   2.5;      // step size
  size_t                   threads;                              // number of threads
  int                      debug;

  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 inter-string time calibration.");
    
    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['P'] = make_field(pdfFile);
    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(ids,               "string identifiers (none = all)") = JPARSER::initialised();
    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);
  }


  JDetector detector;

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

  unique_ptr<JDynamics> dynamics;

  try {

    dynamics.reset(new JDynamics(detector, Tmax_s));

    dynamics->load(calibrationFile);
  }
  catch(const exception& error) {
    if (!calibrationFile.empty()) {
      FATAL(error.what());
    }
  }

  const JModuleRouter router(dynamics ? dynamics->getDetector() : detector);

  if (ids.empty()) {
    ids = getStringIDs(detector);
  }

  NOTICE("Reading PDFs... " << flush);
  
  const JRegressorStorage_t storage(pdfFile, parameters.TTS_ns);

  NOTICE("OK" << endl);

  JRegressor_t::debug              = debug;
  JRegressor_t::T_ns.setRange(parameters.TMin_ns, parameters.TMax_ns);
  JRegressor_t::Vmax_npe           = parameters.VMax_npe;
  JRegressor_t::MAXIMUM_ITERATIONS = parameters.NMax;


  // preprocess input data
  
  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, parameters.R_Hz);

    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;
      const JFIT::JEvt*  evt = ps;

      summary.update(*tev);

      if (dynamics) {
        dynamics->update(*tev);
      }

      if (!evt->empty()) {

        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.hasW(JSTART_LENGTH_METRES)       &&
            fit.getW(JSTART_LENGTH_METRES) > 0.0) {
          Z_m = JZRange(parameters.ZMin_m, parameters.ZMax_m + fit.getW(JSTART_LENGTH_METRES));
        }

        const JFIT::JModel<JLine1Z> match(tz, parameters.roadWidth_m, JRegressor_t::T_ns, Z_m);

        // hit selection based on start value

        buffer_type buffer;

        for (vector<JHitL0>::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {

          JHitW0 hit(*i, summary.getRate(i->getPMTIdentifier()));

          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());
          
          map[module.getString()].push_back(*hit);
        }
        
        data.push_back({map, tz});
      }
    }
  }
  STATUS(endl);


  // fit

  JGradient fit(number_of_iterations, number_of_extra_steps, epsilon, 3);

  for (int id : ids) {
    fit.push_back(JModifier_t(MAKE_STRING("string: " << FILL(4,'0') << id << FILL()), new JStringEditor(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 JStringEditor* p = dynamic_cast<const JStringEditor*>(fit[i].get());

    if (p != NULL) {
      STATUS("string: " << FILL(4,'0') << p->id << FILL() << ' ' << FIXED(9,3) << p->t0 << " [ns]" << endl);
    }
  }


  if (overwriteDetector) {

    detector.comment.add(JMeta(argc, argv));

    map<int, double> calibration;

    double t0 = 0.0;

    for (size_t i = 0; i != fit.size(); ++i) {

      const JStringEditor* p = dynamic_cast<const JStringEditor*>(fit[i].get());

      if (p != NULL) {

        calibration[p->id] = p->t0;

        t0 += p->t0;
      }
    }

    t0 /= calibration.size(); 

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

      if (!module->empty()) {

        map<int, double>::const_iterator p = calibration.find(module->getString());

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