Jpp 20.0.0-195-g190c9e876
the software that should make you happy
Loading...
Searching...
No Matches
JPerth.cc File Reference

Program to determine string or optical module time calibrations. 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 "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"

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::JEditor
 Auxiliary class for editing time offset. More...
 
struct  JRECONSTRUCTION::JPerth
 Thread pool for fits to data. More...
 
struct  JRECONSTRUCTION::JPerth_t
 Auxiliary data structure for chi2 function object. More...
 

Namespaces

namespace  JRECONSTRUCTION
 Support.
 

Typedefs

typedef std::vector< JHitW0JRECONSTRUCTION::buffer_type
 hits
 
typedef std::map< int, buffer_typeJRECONSTRUCTION::map_type
 identifier -> hits
 
typedef std::vector< event_typeJRECONSTRUCTION::data_type
 
typedef JLANG::JSTDObjectReader< const event_typeJRECONSTRUCTION::input_type
 
typedef JFIT::JRegressorStorage< JFIT::JLine3Z, JFIT::JGandalfJRECONSTRUCTION::JRegressorStorage_t
 
typedef JFIT::JRegressor< JFIT::JLine3Z, JFIT::JGandalfJRECONSTRUCTION::JRegressor_t
 

Functions

std::set< int > JRECONSTRUCTION::getTimeOverThresholdIDs ()
 Get list of time-over-threshold identifiers.
 
int main (int argc, char **argv)
 

Detailed Description

Program to determine string or optical module time calibrations.

Author
mdejong

Definition in file JPerth.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 313 of file JPerth.cc.

314{
315 using namespace std;
316 using namespace JPP;
317 using namespace KM3NETDAQ;
318
319 typedef JParallelFileScanner< JTypeList<JDAQEvent, JFIT::JEvt> > JParallelFileScanner_t;
320 typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
321 typedef JMultipleFileScanner<calibration_types> JCalibration_t;
322
324
325
326 JMultipleFileScanner<> inputFile;
327 JLimit_t& numberOfEvents = inputFile.getLimit();
328 string detectorFile;
329 JCalibration_t calibrationFile;
330 double Tmax_s;
331 string pdfFile;
332 JMuonGandalfParameters_t parameters;
333 bool overwriteDetector;
334 set<int> strings; // string identifiers
335 set<int> modules; // module identifiers
336 set<int> indices; // time-over-threshold identifiers
337 range_type X; // time-over-threshold values
338 int number_of_iterations = 1000;
339 int number_of_extra_steps = 0;
340 double epsilon = 1.0e-4;
341 double T_ns = 2.5; // step size
342 JPMTParametersMap pmtParameters;
343 size_t threads; // number of threads
344 int debug;
345
346 parameters.ZMin_m = -1.0e4; // set range hit selection
347 parameters.ZMax_m = +1.0e4;
348
349 const int DEFAULT_ID = -1;
350
351 try {
352
353 JProperties properties;
354
355 properties.insert(gmake_property(number_of_iterations));
356 properties.insert(gmake_property(number_of_extra_steps));
357 properties.insert(gmake_property(epsilon));
358 properties.insert(gmake_property(T_ns));
359
360 JParser<> zap("Program to determine string or optical module time calibrations.");
361
362 zap['f'] = make_field(inputFile);
363 zap['a'] = make_field(detectorFile);
364 zap['+'] = make_field(calibrationFile) = JPARSER::initialised();
365 zap['T'] = make_field(Tmax_s) = 100.0;
366 zap['n'] = make_field(numberOfEvents) = JLimit::max();
367 zap['F'] = make_field(pdfFile);
368 zap['P'] = make_field(pmtParameters, "PMT simulation data (or corresponding file name)") = JPARSER::initialised();
369 zap['@'] = make_field(parameters) = JPARSER::initialised();
370 zap['A'] = make_field(overwriteDetector, "overwrite detector file provided through '-a' with fitted time offsets.");
371 zap['S'] = make_field(strings, "string identifiers (" << DEFAULT_ID << " => all)") = JPARSER::initialised();
372 zap['M'] = make_field(modules, "module identifiers (" << DEFAULT_ID << " => all)") = JPARSER::initialised();
373 zap['I'] = make_field(indices, "time-over-threshold indices (" << DEFAULT_ID << " => all)") = JPARSER::initialised();
374 zap['x'] = make_field(X, "fit range for time slewing correction") = range_type(4.0, 175.0);
375 zap['%'] = make_field(properties, "fit parameters") = JPARSER::initialised();
376 zap['N'] = make_field(threads, "number of threads") = 1;
377 zap['d'] = make_field(debug) = 1;
378
379 zap(argc, argv);
380 }
381 catch(const exception& error) {
382 FATAL(error.what() << endl);
383 }
384
385
386 if ((strings.empty() ? 0 : 1) +
387 (modules.empty() ? 0 : 1) +
388 (indices.empty() ? 0 : 1) != 1) {
389 FATAL("Set either strings (option -S) or modules (option -M) or indices (option -I)." << endl);
390 }
391
392
393 typedef JDAQHit::JTOT_t JTOT_t;
394
395 const JRange<double> range(0.0, numeric_limits<JTOT_t>::max());
396
398
399 try {
400 load(detectorFile, detector);
401 }
402 catch(const JException& error) {
403 FATAL(error);
404 }
405
406 unique_ptr<JDynamics> dynamics;
407
408 if (!calibrationFile.empty()) {
409
410 try {
411
412 dynamics.reset(new JDynamics(detector, Tmax_s));
413
414 dynamics->load(calibrationFile);
415 }
416 catch(const exception& error) {
417 FATAL(error.what());
418 }
419 }
420
421 const JModuleRouter router(dynamics ? dynamics->getDetector() : detector);
422
423 NOTICE("Reading PDFs... " << flush);
424
425 const JRegressorStorage_t storage(pdfFile, JTimeRange(parameters.TMin_ns, parameters.TMax_ns), parameters.TTS_ns);
426
427 NOTICE("OK" << endl);
428
429 JRegressor_t::debug = debug;
430 JRegressor_t::Vmax_npe = parameters.VMax_npe;
431 JRegressor_t::MAXIMUM_ITERATIONS = parameters.NMax;
432
433 // preprocess input data
434
435 NOTICE("Reading data" << endl);
436
437 const JBuildL0<JHitL0> buildL0;
438
440
441 counter_type skip = inputFile.getLowerLimit();
442 counter_type counter = inputFile.getLowerLimit();
443
444 for (JMultipleFileScanner<>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
445
446 JSummaryFileRouter summary(*i);
447
448 for (JParallelFileScanner_t in(*i); (skip -= in.skip(skip)) == 0 && in.hasNext() && counter != numberOfEvents; ++counter) {
449
450 STATUS("event: " << setw(10) << counter << '\r'); DEBUG(endl);
451
452 multi_pointer_type ps = in.next();
453
454 const JDAQEvent* tev = ps;
455 JFIT::JEvt* evt = ps;
456
457 summary.update(*tev);
458
459 if (dynamics) {
460 dynamics->update(*tev);
461 }
462
463 JFIT::JEvt::iterator __end = partition(evt->begin(), evt->end(), JHistory::is_application(JMUONGANDALF));
464
465 if (evt->begin() != __end) {
466
467 sort(evt->begin(), __end, qualitySorter);
468
469 vector<JHitL0> dataL0;
470
471 buildL0(*tev, router, true, back_inserter(dataL0));
472
473 const JFIT::JFit& fit = (*evt)[0];
474
475 const JRotation3D R (getDirection(fit));
476 const JLine1Z tz(getPosition (fit).rotate(R), fit.getT());
477 JRange<double> Z_m;
478
479 if (fit.getW(JSTART_LENGTH_METRES, 0.0) > 0.0) {
480 Z_m = JZRange(fit.getW(JSTART_ZMIN_M) + parameters.ZMin_m,
481 fit.getW(JSTART_ZMAX_M) + parameters.ZMax_m);
482 }
483
484 const JFIT::JModel<JLine1Z> match(tz, parameters.roadWidth_m, JTimeRange(parameters.TMin_ns, parameters.TMax_ns), Z_m);
485
486 // hit selection based on start value
487
488 buffer_type buffer;
489
490 for (vector<JHitL0>::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
491
492 const JPMTIdentifier id(i->getModuleID(), i->getPMTAddress());
493
494 const JPMTParameters& wip = pmtParameters.getPMTParameters(id);
495
496 const int type = wip.getType();
497 const double QE = wip.QE;
498 const double R_Hz = summary.getRate(i->getPMTIdentifier(), parameters.R_Hz);
499
500 JHitW0 hit(*i, type, QE, R_Hz);
501
502 hit.rotate(R);
503
504 if (match(hit)) {
505 buffer.push_back(hit);
506 }
507 }
508
509 // select first hit
510
511 sort(buffer.begin(), buffer.end(), JHitL0::compare);
512
513 buffer_type::const_iterator __end = unique(buffer.begin(), buffer.end(), equal_to<JDAQPMTIdentifier>());
514
515 // re-organise data
516
517 map_type map;
518
519 for (buffer_type::const_iterator hit = buffer.begin(); hit != __end; ++hit) {
520
521 const JModule& module = router.getModule(hit->getModuleID());
522 const JTOT_t index = (JTOT_t) range.constrain(hit->getToT());
523
524 if (!strings.empty()) { map[module.getString()].push_back(*hit); }
525 if (!modules.empty()) { map[module.getID()] .push_back(*hit); }
526 if (!indices.empty()) { map[index] .push_back(*hit); }
527 }
528
529 data.push_back({map, tz, true});
530 }
531 }
532 }
533 STATUS(endl);
534 NOTICE("OK" << endl);
535
536
537 // fit
538
539 JGradient fit(number_of_iterations, number_of_extra_steps, epsilon, 3);
540
541 if (strings.count(DEFAULT_ID) != 0) { strings = getStringIDs(detector); }
542 if (modules.count(DEFAULT_ID) != 0) { modules = getModuleIDs(detector); }
543 if (indices.count(DEFAULT_ID) != 0) { indices = getTimeOverThresholdIDs(); }
544
545 for (int id : strings) { fit.push_back(JModifier_t(MAKE_STRING("string: " << FILL(4,'0') << id << FILL()), new JEditor(data, id), T_ns)); }
546 for (int id : modules) { fit.push_back(JModifier_t(MAKE_STRING("module: " << setw(8) << id), new JEditor(data, id), T_ns)); }
547 for (int id : indices) { fit.push_back(JModifier_t(MAKE_STRING("index: " << setw(3) << id), new JEditor(data, id), T_ns)); }
548
549 const JPerth_t perth = {storage, data, threads};
550
551 const double chi2 = fit(perth);
552
553 STATUS("result: " << FIXED(12,6) << chi2 << ' ' << setw(6) << fit.numberOfIterations << endl);
554
555 for (size_t i = 0; i != fit.size(); ++i) {
556 {
557 const JEditor* p = dynamic_cast<const JEditor*>(fit[i].get());
558
559 if (p != NULL) { STATUS(fit[i].name << ' ' << FIXED(9,3) << p->t0 << " [ns]" << endl); }
560 }
561 }
562
563
564 if (!indices.empty()) {
565
566 // solve y = ax + b
567
568 JMatrix2S V;
569 double Y[2] = { 0.0, 0.0 };
570
571 for (size_t i = 0; i != fit.size(); ++i) {
572
573 const JEditor* p = dynamic_cast<const JEditor*>(fit[i].get());
574
575 if (p != NULL) {
576
577 const double x = p->id;
578 const double y = p->t0;
579
580 if (X(x)) {
581
582 V.a00 += x*x; V.a01 += x;
583 V.a10 += x; V.a11 += 1;
584
585 Y[0] += x*y;
586 Y[1] += y;
587 }
588 }
589 }
590
591 try {
592
593 V.invert();
594
595 const double a = V.a00 * Y[0] + V.a01 * Y[1];
596 //const double b = V.a10 * Y[0] + V.a11 * Y[1];
597
598 cout << "// Produced by JPerth.cc; to be included in JTimeSlewing.hh" << endl;
599
600 for (size_t i = 0; i != getTimeSlewing.size(); ++i) {
601
602 const double x = i;
603 const double y = a * (X.constrain(x) - TIME_OVER_THRESHOLD_NS);
604
605 cout << "this->push_back(" << FIXED(6,2) << getTimeSlewing[i] - y << ");" << endl;
606 }
607 }
608 catch(const exception& error) {
609 FATAL(error.what());
610 }
611 }
612
613 if (overwriteDetector) {
614
615 detector.comment.add(JMeta(argc, argv));
616
618
619 for (size_t i = 0; i != fit.size(); ++i) {
620
621 const JEditor* p = dynamic_cast<const JEditor*>(fit[i].get());
622
623 if (p != NULL) {
624 calibration[p->id] = p->t0;
625 }
626 }
627
628 const double t0 = getAverage(get_values(calibration), 0.0);
629
630 for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
631
632 if (!module->empty()) {
633
634 map<int, double>::const_iterator p = (!strings.empty() ? calibration.find(module->getString()) :
635 !modules.empty() ? calibration.find(module->getID()) :
636 calibration.end());
637 if (p != calibration.end()) {
638 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
639 module->getPMT(pmt).addT0(p->second - t0);
640 }
641 }
642 }
643 }
644
645 NOTICE("Store calibration data on file " << detectorFile << endl);
646
647 store(detectorFile, detector);
648 }
649}
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
#define STATUS(A)
Definition JMessage.hh:63
#define NOTICE(A)
Definition JMessage.hh:64
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:74
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
#define MAKE_STRING(A)
Make string.
Definition JPrint.hh:63
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
Detector data structure.
Definition JDetector.hh:96
int getString() const
Get string number.
Definition JLocation.hh:135
Router for direct addressing of module data in detector data structure.
Data structure for a composite optical module.
Definition JModule.hh:76
Auxiliary class for map of PMT parameters.
const JPMTParameters & getPMTParameters(const JPMTIdentifier &id) const
Get PMT parameters.
Data structure for PMT parameters.
double QE
relative quantum efficiency
int getType() const
Get type for for time-slewing correction.
Utility class to parse parameter values.
Data structure for set of track fit results.
Data structure for track fit results with history and optional associated values.
const std::vector< double > & getW() const
Get associated values.
double getT() const
Get time.
Data structure for fit of straight line paralel to z-axis.
Definition JLine1Z.hh:29
General exception.
Definition JException.hh:24
int getID() const
Get identifier.
Definition JObjectID.hh:50
2 x 2 symmetric matrix
Definition JMatrix2S.hh:28
void invert()
Invert matrix.
Definition JMatrix2S.hh:66
Utility class to parse command line options.
Definition JParser.hh:1698
Auxiliary class for a hit with background rate value.
Definition JHitW0.hh:25
General purpose class for object reading from a list of file names.
General purpose class for parallel reading of objects from a single file or multiple files.
File router for fast addressing of summary data.
Range of values.
Definition JRange.hh:42
T constrain(argument_type x) const
Constrain value to range.
Definition JRange.hh:350
Template L0 hit builder.
Definition JBuildL0.hh:38
unsigned char JTOT_t
time over threshold [ns]
Definition JDAQHit.hh:40
const char * map
Definition elog.cc:87
static const int JSTART_ZMAX_M
end position of track see JRECONSTRUCTION::JMuonStart
static const int JSTART_LENGTH_METRES
distance between projected positions on the track of optical modules for which the response does not ...
static const int JSTART_ZMIN_M
start position of track see JRECONSTRUCTION::JMuonStart
JDirection3D getDirection(const Vec &dir)
Get direction.
JPosition3D getPosition(const Vec &pos)
Get position.
JTOOLS::JRange< double > JTimeRange
Type definition for time range (unit [s]).
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
std::set< int > getStringIDs(const JDetector &detector)
Get list of strings identifiers.
const double TIME_OVER_THRESHOLD_NS
Specification for time-over-threshold corresponding to a one photo-electron pulse.
std::set< int > getModuleIDs(const JDetector &detector, const bool option=false)
Get list of modules identifiers.
JTOOLS::JRange< double > JZRange
const array_type< JValue_t > & get_values(const std::map< JKey_t, JValue_t, JComparator_t, JAllocator_t > &data)
Method to create array of values of map.
std::iterator_traits< T >::value_type getAverage(T __begin, T __end)
Get average.
Definition JMath.hh:494
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
bool qualitySorter(const JFit &first, const JFit &second)
Comparison of fit results.
std::set< int > getTimeOverThresholdIDs()
Get list of time-over-threshold identifiers.
Definition JPerth.cc:84
Long64_t counter_type
Type definition for counter.
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
JRange< int > range_type
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition JDAQ.hh:26
Auxiliary data structure for sequence of same character.
Definition JManip.hh:330
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
Calibration.
Definition JHead.hh:330
Detector file.
Definition JHead.hh:227
Data structure for measured coincidence rates of all pairs of PMTs in optical module.
Definition JFitK40.hh:103
Dynamic detector calibration.
Definition JDynamics.hh:81
Conjugate gradient fit.
Definition JGradient.hh:75
Auxiliary class to test history.
Definition JHistory.hh:188
Auxiliary class to match data points with given model.
Auxiliary data structure for editable parameter.
Definition JGradient.hh:49
Template data structure for storage of internal data.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
Auxiliary class for editing time offset.
Definition JPerth.cc:114
double t0
time offset [ns]
Definition JPerth.cc:151
double TMin_ns
minimal time w.r.t. Cherenkov hypothesis [ns]
double TMax_ns
maximal time w.r.t. Cherenkov hypothesis [ns]
double VMax_npe
maximum number of of photo-electrons
Auxiliary data structure for chi2 function object.
Definition JPerth.cc:285
Auxiliary class for defining the range of iterations of objects.
Definition JLimit.hh:45
static counter_type max()
Get maximum counter value.
Definition JLimit.hh:128
Auxiliary class for ROOT I/O of application specific meta data.
Definition JMeta.hh:72
Auxiliary data structure for sorting of hits.
Definition JHitL0.hh:85