Jpp  test_elongated_shower_pde
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JDynamics.hh
Go to the documentation of this file.
1 #ifndef __JDYNAMICS__JDYNAMICS__
2 #define __JDYNAMICS__JDYNAMICS__
3 
4 #include <ostream>
5 #include <set>
6 
7 #include "JDetector/JDetector.hh"
10 #include "JUTC/JUTCTimeRange.hh"
11 
12 #include "JCompass/JEvt.hh"
13 #include "JCompass/JEvtToolkit.hh"
14 #include "JCompass/JSupport.hh"
15 
16 #include "JAcoustics/JModel.hh"
17 #include "JAcoustics/JGeometry.hh"
18 #include "JAcoustics/JEvt.hh"
19 #include "JAcoustics/JSupport.hh"
20 
21 #include "JLang/JObjectIterator.hh"
22 #include "JLang/JException.hh"
23 
25 #include "JGeometry3D/JEigen3D.hh"
26 
27 #include "JTools/JElement.hh"
28 #include "JTools/JCollection.hh"
29 #include "JTools/JPolfit.hh"
30 #include "JTools/JHashMap.hh"
31 #include "JTools/JRange.hh"
32 
35 
36 
37 /**
38  * \file
39  *
40  * Dynamic detector calibration.
41  * \author mdejong
42  */
43 
44 namespace JDYNAMICS {}
45 namespace JPP { using namespace JDYNAMICS; }
46 
47 /**
48  * Main namespace for dynamic position and orientation calibration.
49  */
50 namespace JDYNAMICS {
51 
56  using JUTC::JUTCTimeRange;
57 
58 
59  /**
60  * Dynamic detector calibration.
61  */
62  struct JDynamics :
63  public JDetector
64  {
65  /**
66  * Auxiliary data structure to track applicability period of calibration data.
67  */
68  struct JUTCTracker :
69  public JUTCTimeRange
70  {
71  /**
72  * Constructor.
73  *
74  * \param Tmax_s applicability period [s]
75  */
76  JUTCTracker(const double Tmax_s) :
77  JUTCTimeRange(-Tmax_s, +Tmax_s)
78  {}
79 
80 
81  /**
82  * Reset.
83  */
84  void reset()
85  {
86  set(0.0);
87  }
88 
89 
90  /**
91  * Set.
92  *
93  * \param t0_s UTC time [s]
94  */
95  void set(const double t0_s)
96  {
97  const double Tmax_s = 0.5 * (getUpperLimit() - getLowerLimit());
98 
99  setRange(t0_s - Tmax_s, t0_s + Tmax_s);
100  }
101  };
102 
103 
104  /**
105  * Dynamic orientation calibration.
106  */
107  struct JOrientation :
108  public JUTCTracker
109  {
110  enum {
111  NUMBER_OF_POINTS = 20, //!< number of points for interpolation
112  NUMBER_OF_DEGREES = 1 //!< number of degrees for interpolation
113  };
114 
119  typedef typename function_type::collection_type::container_type container_type;
120 
123 
126 
127 
128  /**
129  * Constructor.
130  *
131  * \param detector detector
132  * \param Tmax_s applicability period of calibration [s]
133  */
135  const double Tmax_s) :
136  JUTCTracker(Tmax_s)
137  {
138  using namespace JPP;
139 
141  THROW(JValueOutOfRange, "Detector version " << detector.getVersion() << " < " << JDetectorVersion::V4);
142  }
143 
144  if (hasDetectorAddressMap(detector.getID())) {
145 
146  const JDetectorAddressMap& demo = getDetectorAddressMap(detector.getID());
147 
148  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
149  if (module->getFloor() != 0) {
150  buffer[module->getID()] = getRotation(JDETECTOR::getModule(demo.get(module->getID())), *module);
151  }
152  }
153 
154  } else {
155 
156  THROW(JValueOutOfRange, "No detector address map for detector identier " << detector.getID());
157  }
158  }
159 
160 
161  /**
162  * Load calibration data.
163  *
164  * \param input calibration data
165  */
167  {
168  this->reset();
169 
170  while (input.hasNext()) {
171 
172  const JCOMPASS::JOrientation* orientation = input.next();
173 
174  static_cast<container_type&>(calibration[orientation->id]).push_back(element_type(orientation->t, getQuaternion(*orientation)));
175  }
176 
177  for (data_type::iterator i = calibration.begin(); i != calibration.end(); ++i) {
178  i->second.sort();
179  i->second.compile();
180  }
181  }
182 
183 
184  bool empty() const { return calibration.empty(); } //!< empty
185  size_t size() const { return calibration.size(); } //!< size
186  const_iterator begin() const { return calibration.begin(); } //!< begin of calibration data
187  const_iterator end() const { return calibration.end(); } //!< end of calibration data
188  const_reverse_iterator rbegin() const { return calibration.rbegin(); } //!< begin of reverse of calibration data
189  const_reverse_iterator rend() const { return calibration.rend(); } //!< begin of reverse of calibration data
190 
191 
192  /**
193  * Get coverage.
194  *
195  * \param detector detector
196  * \param t1_s UTC time [s]
197  * \return coverage [0,1]
198  */
199  double getCoverage(const JDetector& detector, const double t1_s) const
200  {
201  int n0 = 0;
202  int n1 = 0;
203 
204  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
205 
206  if (module->getFloor() != 0 && !module->has(COMPASS_DISABLE)) {
207 
208  ++n0;
209 
210  if (calibration.has(module->getID()) && !calibration.get(module->getID()).empty()) {
211  if (calibration[module->getID()].getXmin() <= t1_s &&
212  calibration[module->getID()].getXmax() >= t1_s) {
213  ++n1;
214  }
215  }
216  }
217  }
218 
219  if (n0 != 0)
220  return (double) n1 / (double) n0;
221  else
222  return 0.0;
223  }
224 
225 
226  /**
227  * Calibrate given detector at given UTC time.
228  *
229  * \param detector detector (I/O)
230  * \param t1_s UTC time [s]
231  * \return true if updated; else false
232  */
233  bool update(JDetector& detector, const double t1_s)
234  {
235  using namespace std;
236  using namespace JPP;
237 
238  if (!calibration.empty()) {
239 
240  if (!in_range(t1_s)) {
241 
242  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
243 
244  if (module->getFloor() != 0 && !module->has(COMPASS_DISABLE) && calibration.has(module->getID()) && buffer.has(module->getID())) {
245 
246  const function_type& f1 = calibration.get(module->getID());
247 
248  if (!f1.empty()) {
249 
250  if (t1_s >= f1.getXmin() && t1_s <= f1.getXmax()) {
251 
252  JQuaternion3D Q0 = buffer[module->getID()];
253  JQuaternion3D Q1 = module->getQuaternion() * f1(t1_s);
254 
255  const JPosition3D center = module->getPosition();
256 
257  module->sub(center);
258 
259  module->rotate(Q1 * Q0.conjugate());
260 
261  module->add(center);
262 
263  buffer[module->getID()] = Q1;
264  }
265  }
266  }
267  }
268 
269  set(t1_s);
270 
271  return true;
272  }
273  }
274 
275  return false;
276  }
277 
278  protected:
281  };
282 
283 
284  /**
285  * Dynamic position calibration.
286  */
287  struct JPosition :
288  public JUTCTracker
289  {
290  enum {
291  NUMBER_OF_POINTS = 7, //!< number of points for interpolation
292  NUMBER_OF_DEGREES = 2 //!< number of degrees for interpolation
293  };
294 
296 
301 
303 
306 
307 
308  /**
309  * Constructor.
310  *
311  * \param detector detector
312  * \param Tmax_s applicability period of calibration [s]
313  */
315  const double Tmax_s) :
316  JUTCTracker(Tmax_s),
317  geometry(detector)
318  {
319  using namespace JPP;
320 
322  THROW(JValueOutOfRange, "Detector version " << detector.getVersion() << " < " << JDetectorVersion::V4);
323  }
324  }
325 
326 
327  /**
328  * Load calibration data.
329  *
330  * \param input calibration data
331  */
333  {
334  using namespace JPP;
335 
336  this->reset();
337 
338  while (input.hasNext()) {
339 
340  const JACOUSTICS::JEvt* evt = input.next();
341  const double t1_s = 0.5 * (evt->UNIXTimeStart + evt->UNIXTimeStop);
342 
343  for (JACOUSTICS::JEvt::const_iterator i = evt->begin(); i != evt->end(); ++i) {
344  calibration[i->id][t1_s] = JMODEL::JString(i->tx, i->ty);
345  }
346  }
347 
348  for (data_type::iterator i = calibration.begin(); i != calibration.end(); ++i) {
349  i->second.compile();
350  }
351  }
352 
353 
354  bool empty() const { return calibration.empty(); } //!< empty
355  size_t size() const { return calibration.size(); } //!< size
356  const_iterator begin() const { return calibration.begin(); } //!< begin of calibration data
357  const_iterator end() const { return calibration.end(); } //!< end of calibration data
358  const_reverse_iterator rbegin() const { return calibration.rbegin(); } //!< begin of reverse of calibration data
359  const_reverse_iterator rend() const { return calibration.rend(); } //!< begin of reverse of calibration data
360 
361 
362  /**
363  * Get coverage.
364  *
365  * \param detector detector
366  * \param t1_s UTC time [s]
367  * \return coverage [0,1]
368  */
369  double getCoverage(const JDetector& detector, const double t1_s) const
370  {
371  int n0 = 0;
372  int n1 = 0;
373 
374  std::set<int> string;
375 
376  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
377  if (module->getFloor() != 0) {
378  string.insert(module->getString());
379  }
380  }
381 
382  for (std::set<int>::const_iterator i = string.begin(); i != string.end(); ++i) {
383 
384  ++n0;
385 
386  if (calibration.has(*i) && !calibration.get(*i).empty()) {
387  if (calibration[*i].getXmin() <= t1_s &&
388  calibration[*i].getXmax() >= t1_s) {
389  ++n1;
390  }
391  }
392  }
393 
394  if (n0 != 0)
395  return (double) n1 / (double) n0;
396  else
397  return 0.0;
398  }
399 
400 
401  /**
402  * Get detector geometry.
403  *
404  * \return detector geometry
405  */
406  const JGeometry& getGeometry() const
407  {
408  return geometry;
409  }
410 
411 
412  /**
413  * Calibrate given detector at given UTC time.
414  *
415  * \param detector detector (I/O)
416  * \param t1_s UTC time [s]
417  * \return true if updated; else false
418  */
419  bool update(JDetector& detector, const double t1_s)
420  {
421  using namespace std;
422  using namespace JPP;
423 
424  if (!calibration.empty()) {
425 
426  if (!in_range(t1_s)) {
427 
429 
430  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
431 
432  if (module->getFloor() != 0) {
433 
434  if (!buffer.has(module->getString())) {
435 
436  if (calibration.has(module->getString())) {
437 
438  const function_type& f1 = calibration.get(module->getString());
439 
440  if (!f1.empty()) {
441  if (t1_s >= f1.getXmin() && t1_s <= f1.getXmax()) {
442  buffer[module->getString()] = f1(t1_s);
443  }
444  }
445  }
446  }
447 
448  if (buffer.has(module->getString())) {
449  module->set(geometry[module->getString()].getPosition(buffer[module->getString()], module->getFloor()) - getPiezoPosition());
450  }
451  }
452  }
453 
454  set(t1_s);
455 
456  return true;
457  }
458  }
459 
460  return false;
461  }
462 
463  protected:
466  };
467 
468 
469  /**
470  * Constructor.
471  *
472  * \param detector detector
473  * \param Tmax_s applicability period of calibration [s]
474  */
476  const double Tmax_s) :
477  JDetector (detector),
478  orientation(detector, Tmax_s),
479  position (detector, Tmax_s)
480  {
482  }
483 
484 
485  /**
486  * Get actual detector.
487  *
488  * \return detector
489  */
490  const JDetector& getDetector() const
491  {
492  return static_cast<const JDetector&>(*this);
493  }
494 
495 
496  /**
497  * Load calibration data.
498  *
499  * \param input calibration data
500  */
501  template<class JObjectIterator_t>
502  void load(JObjectIterator_t& input)
503  {
505 
506  try { orientation.load(dynamic_cast<JObjectIterator<JCOMPASS::JOrientation>&>(input)); } catch(const std::exception&) {}
507  try { position .load(dynamic_cast<JObjectIterator<JACOUSTICS::JEvt>&> (input)); } catch(const std::exception&) {}
508  }
509 
510 
511  /**
512  * Get detector calibrated at given time.
513  *
514  * \param t1_s UTC time [s]
515  * \return true if updated; else false
516  */
517  bool update(const double t1_s)
518  {
519  bool is_updated = false;
520 
521  if (!in_range(t1_s)) {
522 
524 
525  if (orientation.update(static_cast<JDetector&>(*this), t1_s)) { is_updated = true; }
526  if (position .update(static_cast<JDetector&>(*this), t1_s)) { is_updated = true; }
527 
528  range.join(orientation.getUTCTimeRange());
529  range.join(position .getUTCTimeRange());
530 
531  setUTCTimeRange(range);
532  }
533 
534  return is_updated;
535  }
536 
537 
538  /**
539  * Get detector calibrated at given time.
540  *
541  * \param chronometer chronometer
542  * \return true if updated; else false
543  */
544  bool update(const JDAQChronometer& chronometer)
545  {
546  return update(chronometer.getTimesliceStart().getUTCseconds());
547  }
548 
549 
550  /**
551  * Data structure for coverage of dynamic calibrations.
552  */
553  struct coverage_type {
554  double orientation; //!< coverage of detector by available orientation calibration [0,1]
555  double position; //!< coverage of detector by available position calibration [0,1]
556  };
557 
558 
559  /**
560  * Get coverage at given time.
561  *
562  * \param t1_s UTC time [s]
563  * \return coverage
564  */
565  coverage_type getCoverage(const double t1_s) const
566  {
567  return { orientation.getCoverage(*this, t1_s), position.getCoverage(*this, t1_s) };
568  }
569 
570 
571  /**
572  * Get coverage at given time.
573  *
574  * \param chronometer chronometer
575  * \return coverage
576  */
577  coverage_type getCoverage(const JDAQChronometer& chronometer) const
578  {
579  return getCoverage(chronometer.getTimesliceStart().getUTCseconds());
580  }
581 
582 
583  /**
584  * Get actual coverage.
585  *
586  * \return coverage
587  */
589  {
590  return getCoverage(0.5 * (getLowerLimit() + getUpperLimit()));
591  }
592 
593 
594  JOrientation orientation; //!< orientation calibration
595  JPosition position; //!< position calibration
596  };
597 }
598 
599 #endif
JTOOLS::JElement2D< double, JACOUSTICS::JMODEL::JString > element_type
Definition: JDynamics.hh:297
static const JRange< T, JComparator_t > DEFAULT_RANGE
Default range.
Definition: JRange.hh:555
void set(const double t0_s)
Set.
Definition: JDynamics.hh:95
Exceptions.
T getLowerLimit() const
Get lower limit.
Definition: JRange.hh:202
JDAQUTCExtended getTimesliceStart() const
Get start of timeslice.
int id
module identifier
General purpose class for collection of elements, see: &lt;a href=&quot;JTools.PDF&quot;;&gt;Collection of elements...
Definition: JCollection.hh:73
const_reverse_iterator rend() const
begin of reverse of calibration data
Definition: JDynamics.hh:359
The elements in a collection are sorted according to their abscissa values and a given distance opera...
coverage_type getCoverage(const JDAQChronometer &chronometer) const
Get coverage at given time.
Definition: JDynamics.hh:577
coverage_type getCoverage(const double t1_s) const
Get coverage at given time.
Definition: JDynamics.hh:565
Dynamic orientation calibration.
Definition: JDynamics.hh:107
Dynamic position calibration.
Definition: JDynamics.hh:287
double position
coverage of detector by available position calibration [0,1]
Definition: JDynamics.hh:555
const JDetector & getDetector() const
Get actual detector.
Definition: JDynamics.hh:490
container_type::const_reverse_iterator const_reverse_iterator
Definition: JHashMap.hh:87
JDynamics(const JDetector &detector, const double Tmax_s)
Constructor.
Definition: JDynamics.hh:475
JTOOLS::JHashMap< int, function_type > data_type
Definition: JDynamics.hh:302
const_reverse_iterator rend() const
begin of reverse of calibration data
Definition: JDynamics.hh:189
JQuaternion3D getQuaternion(const JQuaternion &Q)
Get quaternion.
JTOOLS::JPolfitFunction1D< NUMBER_OF_POINTS, NUMBER_OF_DEGREES, element_type, JTOOLS::JCollection > function_type
Definition: JDynamics.hh:118
function_type::collection_type::container_type container_type
Definition: JDynamics.hh:119
General purpose class for hash map of unique elements.
ROOT TTree parameter settings.
JPosition position
position calibration
Definition: JDynamics.hh:595
Detector data structure.
Definition: JDetector.hh:89
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:696
Acoustic geometries.
bool hasDetectorAddressMap(const int id)
Check if detector address map is available.
const_reverse_iterator rbegin() const
begin of reverse of calibration data
Definition: JDynamics.hh:358
JOrientation orientation
orientation calibration
Definition: JDynamics.hh:594
const_iterator begin() const
begin of calibration data
Definition: JDynamics.hh:186
static JRotation getRotation
Function object to get rotation matrix to go from first to second module.
Lookup table for PMT addresses in detector.
const std::string & getVersion() const
Get version.
Interface of object iteration for a single data type.
data_type::const_reverse_iterator const_reverse_iterator
Definition: JDynamics.hh:305
Data structure for detector geometry and calibration.
JTOOLS::JHashMap< int, JGEOMETRY3D::JQuaternion3D > buffer_type
Definition: JDynamics.hh:121
const JModuleAddressMap & get(const int id) const
Get module address map.
ROOT TTree parameter settings.
const JGeometry & getGeometry() const
Get detector geometry.
Definition: JDynamics.hh:406
number of points for interpolation
Definition: JDynamics.hh:111
const JQuaternion3D & getQuaternion() const
Get quaternion.
JACOUSTICS::JGeometry JGeometry
Definition: JDynamics.hh:295
static const int COMPASS_DISABLE
Enable (disable) use of compass if this status bit is 0 (1);.
double orientation
coverage of detector by available orientation calibration [0,1]
Definition: JDynamics.hh:554
void load(JObjectIterator_t &input)
Load calibration data.
Definition: JDynamics.hh:502
const_reverse_iterator rbegin() const
begin of reverse of calibration data
Definition: JDynamics.hh:188
JUTCTracker(const double Tmax_s)
Constructor.
Definition: JDynamics.hh:76
void setUTCTimeRange(const JRange< double > &timerange)
Set UTC time range.
data_type::const_iterator const_iterator
Definition: JDynamics.hh:304
virtual mapped_type & get(typename JClass< key_type >::argument_type key) override
Get mapped value.
Definition: JHashMap.hh:119
double UNIXTimeStop
stop time
void load(JObjectIterator< JACOUSTICS::JEvt > &input)
Load calibration data.
Definition: JDynamics.hh:332
JVector3D & sub(const JVector3D &vector)
Subtract vector.
Definition: JVector3D.hh:158
const_iterator end() const
end of calibration data
Definition: JDynamics.hh:357
Detector specific mapping between logical positions and readout channels of PMTs in optical modules...
virtual const pointer_type & next()=0
Get next element.
Detector file.
Definition: JHead.hh:224
Acoustic event fit.
Version with quaternion and time offset per module.
JPosition(const JDetector &detector, const double Tmax_s)
Constructor.
Definition: JDynamics.hh:314
void setRange(const range_type &range)
Set range.
Definition: JRange.hh:146
int getID() const
Get identifier.
Definition: JObjectID.hh:50
coverage_type getCoverage() const
Get actual coverage.
Definition: JDynamics.hh:588
T getUpperLimit() const
Get upper limit.
Definition: JRange.hh:213
virtual bool hasNext()=0
Check availability of next element.
General purpose class for a collection of sorted elements.
Compass event fit.
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
double getCoverage(const JDetector &detector, const double t1_s) const
Get coverage.
Definition: JDynamics.hh:199
void load(JObjectIterator< JCOMPASS::JOrientation > &input)
Load calibration data.
Definition: JDynamics.hh:166
JUINT32_t getUTCseconds() const
Get time.
JDetectorAddressMap & getDetectorAddressMap()
Get detector address map.
JTOOLS::JElement2D< double, JGEOMETRY3D::JQuaternion3D > element_type
Definition: JDynamics.hh:115
bool in_range(argument_type x) const
Test whether value is inside range.
Definition: JRange.hh:323
Data structure for coverage of dynamic calibrations.
Definition: JDynamics.hh:553
static const JGetDetectorVersion getDetectorVersion
Function object to map detector version to numerical value.
data_type::const_iterator const_iterator
Definition: JDynamics.hh:124
bool update(JDetector &detector, const double t1_s)
Calibrate given detector at given UTC time.
Definition: JDynamics.hh:233
const_iterator end() const
end of calibration data
Definition: JDynamics.hh:187
data_type::const_reverse_iterator const_reverse_iterator
Definition: JDynamics.hh:125
bool update(const double t1_s)
Get detector calibrated at given time.
Definition: JDynamics.hh:517
Data structure for unit quaternion in three dimensions.
bool empty() const
empty
Definition: JDynamics.hh:354
bool update(const JDAQChronometer &chronometer)
Get detector calibrated at given time.
Definition: JDynamics.hh:544
z range($ZMAX-$ZMIN)< $MINIMAL_DZ." fi fi typeset -Z 4 STRING typeset -Z 2 FLOOR JPlot1D -f $
UTC time range [s].
size_t size() const
size
Definition: JDynamics.hh:355
Dynamic detector calibration.
Definition: JDynamics.hh:62
double UNIXTimeStart
start time
Auxiliary class to define a range between two values.
JPosition3D getPiezoPosition()
Get relative position of piezo in optical module.
JTOOLS::JHashMap< int, function_type > data_type
Definition: JDynamics.hh:122
const_iterator begin() const
begin of calibration data
Definition: JDynamics.hh:356
JOrientation(const JDetector &detector, const double Tmax_s)
Constructor.
Definition: JDynamics.hh:134
Template class for polynomial interpolation in 1D.
Definition: JPolfit.hh:149
bool update(JDetector &detector, const double t1_s)
Calibrate given detector at given UTC time.
Definition: JDynamics.hh:419
2D Element.
Definition: JElement.hh:46
number of degrees for interpolation
Definition: JDynamics.hh:112
Orientation of module.
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
Exception for accessing a value in a collection that is outside of its range.
Definition: JException.hh:162
Compass event data types.
container_type::const_iterator const_iterator
Definition: JHashMap.hh:86
double getCoverage(const JDetector &detector, const double t1_s) const
Get coverage.
Definition: JDynamics.hh:369
Model for fit to acoutsics data.
const JUTCTimeRange & getUTCTimeRange() const
Get UTC time range.
Acoustic event fit.
bool has(const T &value) const
Test whether given value is present.
JVector3D & add(const JVector3D &vector)
Add vector.
Definition: JVector3D.hh:142
const JModule & getModule(const JModuleAddressMap &memo, const int id=-1, const JLocation &location=JLocation())
Get module according module address map.
JTOOLS::JPolfitFunction1D< NUMBER_OF_POINTS, NUMBER_OF_DEGREES, element_type, JTOOLS::JCollection > function_type
Definition: JDynamics.hh:300
Auxiliary data structure to track applicability period of calibration data.
Definition: JDynamics.hh:68
then usage $script< detector file >< tripodfile >< stage > input file nInput files correspond to the output of JAcousticsEventBuilder[.sh] nFirst stage eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY eval JPrintDetector a $DETECTOR O CAN source JAcoustics sh $DETECTOR_ID typeset A CONFIGURATION for key in Tmax_s
range_type & join(const range_type &range)
Join ranges.
Definition: JRange.hh:415
JQuaternion3D & conjugate()
Conjugate quaternion.
container_type::iterator iterator
Definition: JHashMap.hh:88