314{
318
320 typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type;
322
324
325
327 JLimit_t& numberOfEvents = inputFile.getLimit();
328 string detectorFile;
329 JCalibration_t calibrationFile;
330 double Tmax_s;
331 string pdfFile;
333 bool overwriteDetector;
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;
343 size_t threads;
345
346 parameters.
ZMin_m = -1.0e4;
347 parameters.
ZMax_m = +1.0e4;
348
349 const int DEFAULT_ID = -1;
350
351 try {
352
354
359
360 JParser<> zap(
"Program to determine string or optical module time calibrations.");
361
370 zap[
'A'] =
make_field(overwriteDetector,
"overwrite detector file provided through '-a' with fitted time offsets.");
376 zap[
'N'] =
make_field(threads,
"number of threads") = 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
394
396
398
399 try {
401 }
404 }
405
406 unique_ptr<JDynamics> dynamics;
407
408 if (!calibrationFile.empty()) {
409
410 try {
411
413
414 dynamics->load(calibrationFile);
415 }
416 catch(const exception& error) {
418 }
419 }
420
422
423 NOTICE(
"Reading PDFs... " << flush);
424
426
428
429 JRegressor_t::debug =
debug;
430 JRegressor_t::Vmax_npe = parameters.
VMax_npe;
431 JRegressor_t::MAXIMUM_ITERATIONS = parameters.
NMax;
432
433
434
435 NOTICE(
"Reading data" << endl);
436
438
440
443
445
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
456
457 summary.update(*tev);
458
459 if (dynamics) {
460 dynamics->update(*tev);
461 }
462
464
465 if (evt->begin() != __end) {
466
468
470
471 buildL0(*tev, router, true, back_inserter(dataL0));
472
474
478
482 }
483
485
486
487
489
490 for (vector<JHitL0>::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
491
493
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
510
512
513 buffer_type::const_iterator __end = unique(buffer.begin(), buffer.end(), equal_to<JDAQPMTIdentifier>());
514
515
516
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 }
535
536
537
538
539 JGradient fit(number_of_iterations, number_of_extra_steps, epsilon, 3);
540
544
548
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 {
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
567
569 double Y[2] = { 0.0, 0.0 };
570
571 for (size_t i = 0; i != fit.size(); ++i) {
572
574
575 if (p != NULL) {
576
577 const double x = p->
id;
578 const double y = p->
t0;
579
580 if (X(x)) {
581
584
587 }
588 }
589 }
590
591 try {
592
594
595 const double a = V.
a00 * Y[0] + V.
a01 * Y[1];
596
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
604
605 cout <<
"this->push_back(" <<
FIXED(6,2) << getTimeSlewing[i] -
y <<
");" << endl;
606 }
607 }
608 catch(const exception& error) {
610 }
611 }
612
613 if (overwriteDetector) {
614
616
618
619 for (size_t i = 0; i != fit.size(); ++i) {
620
622
623 if (p != NULL) {
625 }
626 }
627
629
630 for (JDetector::iterator module =
detector.begin(); module !=
detector.end(); ++module) {
631
632 if (!module->empty()) {
633
639 module->getPMT(pmt).addT0(p->second - t0);
640 }
641 }
642 }
643 }
644
645 NOTICE(
"Store calibration data on file " << detectorFile << endl);
646
648 }
649}
#define DEBUG(A)
Message macros.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
#define MAKE_STRING(A)
Make string.
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
int getString() const
Get string number.
Router for direct addressing of module data in detector data structure.
Data structure for a composite optical module.
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.
int getID() const
Get identifier.
void invert()
Invert matrix.
Utility class to parse command line options.
Auxiliary class for a hit with background rate value.
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.
unsigned char JTOT_t
time over threshold [ns]
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.
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.
Long64_t counter_type
Type definition for counter.
KM3NeT DAQ data structures and auxiliaries.
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Auxiliary data structure for sequence of same character.
Auxiliary data structure for floating point format specification.
Data structure for measured coincidence rates of all pairs of PMTs in optical module.
Dynamic detector calibration.
Auxiliary class to test history.
Auxiliary class to match data points with given model.
Auxiliary data structure for editable parameter.
Template data structure for storage of internal data.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Auxiliary class for editing time offset.
double t0
time offset [ns]
Data structure for fit parameters.
double TTS_ns
transition-time spread [ns]
double TMin_ns
minimal time w.r.t. Cherenkov hypothesis [ns]
double roadWidth_m
road width [m]
double TMax_ns
maximal time w.r.t. Cherenkov hypothesis [ns]
double VMax_npe
maximum number of of photo-electrons
double ZMax_m
maximal z-positon [m]
double ZMin_m
minimal z-positon [m]
int NMax
maximum number of iterations
double R_Hz
default rate [Hz]
Auxiliary data structure for chi2 function object.
Auxiliary class for defining the range of iterations of objects.
static counter_type max()
Get maximum counter value.
Auxiliary data structure for sorting of hits.