Jpp
 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 <limits>
5 #include <ostream>
6 
7 #include "JDetector/JDetector.hh"
9 
10 #include "JCompass/JEvt.hh"
11 #include "JCompass/JEvtToolkit.hh"
12 
13 #include "JAcoustics/JModel.hh"
14 #include "JAcoustics/JGeometry.hh"
15 #include "JAcoustics/JEvt.hh"
16 
17 #include "JLang/JObjectIterator.hh"
18 #include "JLang/JException.hh"
19 
21 #include "JGeometry3D/JEigen3D.hh"
22 
23 #include "JTools/JElement.hh"
24 #include "JTools/JCollection.hh"
25 #include "JTools/JPolfit.hh"
26 #include "JTools/JHashMap.hh"
27 
28 
29 /**
30  * \file
31  *
32  * Dynamic detector calibration.
33  * \author mdejong
34  */
35 
36 namespace JDYNAMICS {}
37 namespace JPP { using namespace JDYNAMICS; }
38 
39 /**
40  * Main namespace for dynamic position and orientation calibration.
41  */
42 namespace JDYNAMICS {
43 
47 
48 
49  /**
50  * Dynamic detector calibration.
51  */
52  struct JDynamics {
53 
54  /**
55  * Get minimal abscissa of mapped function objects.
56  *
57  * \param input input
58  * \return minimal abscissa
59  */
60  template<class JKey_t, class JValue_t, class JEvaluator_t>
61  static inline double getXmin(const JTOOLS::JHashMap<JKey_t, JValue_t, JEvaluator_t>& input)
62  {
63  double xmin = std::numeric_limits<double>::max();
64 
65  for (typename JTOOLS::JHashMap<JKey_t, JValue_t, JEvaluator_t>::const_iterator i = input.begin(); i != input.end(); ++i) {
66  if (!i->second.empty() && i->second.getXmin() < xmin) {
67  xmin = i->second.getXmin();
68  }
69  }
70 
71  return xmin;
72  }
73 
74 
75  /**
76  * Get maximal abscissa of mapped function objects.
77  *
78  * \param input input
79  * \return maximal abscissa
80  */
81  template<class JKey_t, class JValue_t, class JEvaluator_t>
82  static inline double getXmax(const JTOOLS::JHashMap<JKey_t, JValue_t, JEvaluator_t>& input)
83  {
84  double xmax = std::numeric_limits<double>::lowest();
85 
86  for (typename JTOOLS::JHashMap<JKey_t, JValue_t, JEvaluator_t>::const_iterator i = input.begin(); i != input.end(); ++i) {
87  if (!i->second.empty() && i->second.getXmax() > xmax) {
88  xmax = i->second.getXmax();
89  }
90  }
91 
92  return xmax;
93  }
94 
95 
96  /**
97  * Dynamic orientation calibration.
98  */
99  struct JOrientation {
100 
101  enum {
102  NUMBER_OF_POINTS = 20, //!< number of points for interpolation
103  NUMBER_OF_DEGREES = 1 //!< number of degrees for interpolation
104  };
105 
110  typedef typename function_type::collection_type::container_type container_type;
111 
114 
117 
118 
119  /**
120  * Constructor.
121  *
122  * \param detector detector
123  * \param Tmax_s applicability period of calibration [s]
124  */
126  const double Tmax_s) :
127  Tmax_s(Tmax_s),
128  t0_s(std::numeric_limits<double>::lowest())
129  {
130  using namespace JPP;
131 
132  if (hasDetectorAddressMap(detector.getID())) {
133 
134  const JDetectorAddressMap& demo = getDetectorAddressMap(detector.getID());
135 
136  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
137  if (module->getFloor() != 0) {
138  buffer[module->getID()] = getRotation(getModule(demo.get(module->getID())), *module);
139  }
140  }
141 
142  } else {
143 
144  THROW(JValueOutOfRange, "No detector address map for detector identier " << detector.getID());
145  }
146  }
147 
148 
149  /**
150  * Load calibration data.
151  *
152  * \param input calibration data
153  */
155  {
156  t0_s = std::numeric_limits<double>::lowest();
157 
158  while (input.hasNext()) {
159 
160  const JCOMPASS::JOrientation* orientation = input.next();
161 
162  static_cast<container_type&>(calibration[orientation->id]).push_back(element_type(orientation->t, getQuaternion(*orientation)));
163  }
164 
165  for (data_type::iterator i = calibration.begin(); i != calibration.end(); ++i) {
166  i->second.sort();
167  i->second.compile();
168  }
169  }
170 
171 
172  bool empty() const { return calibration.empty(); } //!< empty
173  const_iterator begin() const { return calibration.begin(); } //!< begin of calibration data
174  const_iterator end() const { return calibration.end(); } //!< end of calibration data
175  const_reverse_iterator rbegin() const { return calibration.rbegin(); } //!< begin of reverse of calibration data
176  const_reverse_iterator rend() const { return calibration.rend(); } //!< begin of reverse of calibration data
177 
178 
179  /**
180  * Get applicability period.
181  *
182  * \return applicability period [s]
183  */
184  double getTmax() const
185  {
186  return Tmax_s;
187  }
188 
189 
190  /**
191  * Set applicability period.
192  *
193  * \param Tmax_s applicability period [s]
194  */
195  void setTmax(const double Tmax_s)
196  {
197  this->Tmax_s = Tmax_s;
198  }
199 
200 
201  /**
202  * Get coverage.
203  *
204  * \param detector detector
205  * \param t1_s time [s]
206  * \return coverage
207  */
208  double getCoverage(const JDetector& detector, const double t1_s) const
209  {
210  int n0 = 0;
211  int n1 = 0;
212 
213  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
214 
215  if (module->getFloor() != 0) {
216 
217  ++n0;
218 
219  if (calibration.has(module->getID()) && !calibration.get(module->getID()).empty()) {
220  if (calibration[module->getID()].getXmin() <= t1_s &&
221  calibration[module->getID()].getXmax() >= t1_s) {
222  ++n1;
223  }
224  }
225  }
226  }
227 
228  return (double) n1 / (double) n0;
229  }
230 
231 
232  /**
233  * Calibrate given detector at given time.
234  *
235  * \param detector detector (I/O)
236  * \param t1_s time [s]
237  */
238  void update(JDetector& detector, const double t1_s)
239  {
240  using namespace std;
241  using namespace JPP;
242 
243  if (!calibration.empty()) {
244 
245  if (fabs(t1_s - t0_s) > Tmax_s) {
246 
247  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
248 
249  if (module->getFloor() != 0 && calibration.has(module->getID()) && buffer.has(module->getID())) {
250 
251  const function_type& f1 = calibration.get(module->getID());
252 
253  if (!f1.empty()) {
254 
255  if (t1_s >= f1.getXmin() && t1_s <= f1.getXmax()) {
256 
257  JQuaternion3D Q0 = buffer[module->getID()];
258  JQuaternion3D Q1 = f1(t1_s);
259 
260  const JPosition3D center = module->getPosition();
261 
262  module->sub(center);
263 
264  module->rotate(Q1 * Q0.conjugate());
265 
266  module->add(center);
267 
268  buffer[module->getID()] = Q1;
269  }
270  }
271  }
272  }
273 
274  t0_s = t1_s;
275  }
276  }
277  }
278 
279 
280  /**
281  * Get minimal abscissa.
282  *
283  * \return minimal abscissa
284  */
285  double getXmin() const
286  {
287  return JDynamics::getXmin(this->calibration);
288  }
289 
290 
291  /**
292  * Get maximal abscissa.
293  *
294  * \return maximal abscissa
295  */
296  double getXmax() const
297  {
298  return JDynamics::getXmax(this->calibration);
299  }
300 
301  protected:
304  double Tmax_s;
305 
306  private:
307  double t0_s;
308  };
309 
310 
311  /**
312  * Dynamic position calibration.
313  */
314  struct JPosition {
315 
316  enum {
317  NUMBER_OF_POINTS = 7, //!< number of points for interpolation
318  NUMBER_OF_DEGREES = 2 //!< number of degrees for interpolation
319  };
320 
322 
327 
329 
332 
333 
334  /**
335  * Constructor.
336  *
337  * \param detector detector
338  * \param Tmax_s applicability period of calibration [s]
339  */
341  const double Tmax_s) :
342  geometry(detector),
343  Tmax_s(Tmax_s),
344  t0_s(std::numeric_limits<double>::lowest())
345  {}
346 
347 
348  /**
349  * Load calibration data.
350  *
351  * \param input calibration data
352  */
354  {
355  using namespace JPP;
356 
357  t0_s = std::numeric_limits<double>::lowest();
358 
359  while (input.hasNext()) {
360 
361  const JACOUSTICS::JEvt* evt = input.next();
362  const double t1_s = 0.5 * (evt->UNIXTimeStart + evt->UNIXTimeStop);
363 
364  for (JACOUSTICS::JEvt::const_iterator i = evt->begin(); i != evt->end(); ++i) {
365  calibration[i->id][t1_s] = JMODEL::JString(i->tx, i->ty);
366  }
367  }
368 
369  for (data_type::iterator i = calibration.begin(); i != calibration.end(); ++i) {
370  i->second.compile();
371  }
372  }
373 
374 
375  bool empty() const { return calibration.empty(); } //!< empty
376  const_iterator begin() const { return calibration.begin(); } //!< begin of calibration data
377  const_iterator end() const { return calibration.end(); } //!< end of calibration data
378  const_reverse_iterator rbegin() const { return calibration.rbegin(); } //!< begin of reverse of calibration data
379  const_reverse_iterator rend() const { return calibration.rend(); } //!< begin of reverse of calibration data
380 
381 
382  /**
383  * Get applicability period.
384  *
385  * \return applicability period [s]
386  */
387  double getTmax() const
388  {
389  return Tmax_s;
390  }
391 
392 
393  /**
394  * Set applicability period.
395  *
396  * \param Tmax_s applicability period [s]
397  */
398  void setTmax(const double Tmax_s)
399  {
400  this->Tmax_s = Tmax_s;
401  }
402 
403 
404  /**
405  * Get coverage.
406  *
407  * \param detector detector
408  * \param t1_s time [s]
409  * \return coverage
410  */
411  double getCoverage(const JDetector& detector, const double t1_s) const
412  {
413  int n0 = 0;
414  int n1 = 0;
415 
416  std::set<int> string;
417 
418  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
419  if (module->getFloor() != 0) {
420  string.insert(module->getString());
421  }
422  }
423 
424  for (std::set<int>::const_iterator i = string.begin(); i != string.end(); ++i) {
425 
426  ++n0;
427 
428  if (calibration.has(*i) && !calibration.get(*i).empty()) {
429  if (calibration[*i].getXmin() <= t1_s &&
430  calibration[*i].getXmax() >= t1_s) {
431  ++n1;
432  }
433  }
434  }
435 
436  return (double) n1 / (double) n0;
437  }
438 
439 
440  /**
441  * Get detector geometry.
442  *
443  * \return detector geometry
444  */
445  const JGeometry& getGeometry() const
446  {
447  return geometry;
448  }
449 
450 
451  /**
452  * Calibrate given detector at given time.
453  *
454  * \param detector detector (I/O)
455  * \param t1_s time [s]
456  */
457  void update(JDetector& detector, const double t1_s)
458  {
459  using namespace std;
460  using namespace JPP;
461 
462  if (!calibration.empty()) {
463 
464  if (fabs(t1_s - t0_s) > Tmax_s) {
465 
467 
468  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
469 
470  if (module->getFloor() != 0) {
471 
472  if (!buffer.has(module->getString())) {
473  if (calibration.has(module->getString())) {
474  buffer[module->getString()] = calibration[module->getString()](t1_s);
475  }
476  }
477 
478  if (buffer.has(module->getString())) {
479  module->set(geometry[module->getString()].getPosition(buffer[module->getString()], module->getFloor()) - getPiezoPosition());
480  }
481  }
482  }
483 
484  t0_s = t1_s;
485  }
486  }
487  }
488 
489 
490  /**
491  * Get minimal abscissa.
492  *
493  * \return minimal abscissa
494  */
495  double getXmin() const
496  {
497  return JDynamics::getXmin(this->calibration);
498  }
499 
500 
501  /**
502  * Get maximal abscissa.
503  *
504  * \return maximal abscissa
505  */
506  double getXmax() const
507  {
508  return JDynamics::getXmax(this->calibration);
509  }
510 
511  protected:
514  double Tmax_s;
515 
516  private:
517  double t0_s;
518  };
519 
520 
521  /**
522  * Constructor.
523  *
524  * \param detector detector
525  * \param Tmax_s applicability period of calibration [s]
526  */
528  const double Tmax_s) :
529  detector (detector),
530  orientation(detector, Tmax_s),
531  position (detector, Tmax_s)
532  {}
533 
534 
535  /**
536  * Load calibration data.
537  *
538  * \param input detector calibration data
539  */
540  template<class JObjectIterator_t>
541  void load(JObjectIterator_t& input)
542  {
543  try { orientation.load(dynamic_cast<JObjectIterator<JCOMPASS::JOrientation>&>(input)); } catch(const std::exception& error) {}
544  try { position .load(dynamic_cast<JObjectIterator<JACOUSTICS::JEvt>&> (input)); } catch(const std::exception& error) {}
545  }
546 
547 
548  /**
549  * Get detector calibrated at given time.
550  *
551  * \param t1_s time [s]
552  * \return detector
553  */
554  const JDetector& operator()(const double t1_s)
555  {
556  orientation.update(this->detector, t1_s);
557  position .update(this->detector, t1_s);
558 
559  return this->detector;
560  }
561 
562  protected:
564 
565  public:
568  };
569 }
570 
571 #endif
void setTmax(const double Tmax_s)
Set applicability period.
Definition: JDynamics.hh:195
JTOOLS::JElement2D< double, JACOUSTICS::JMODEL::JString > element_type
Definition: JDynamics.hh:323
number of points for interpolation
Definition: JDynamics.hh:102
Exceptions.
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:379
The elements in a collection are sorted according to their abscissa values and a given distance opera...
void update(JDetector &detector, const double t1_s)
Calibrate given detector at given time.
Definition: JDynamics.hh:457
Dynamic orientation calibration.
Definition: JDynamics.hh:99
Dynamic position calibration.
Definition: JDynamics.hh:314
container_type::const_reverse_iterator const_reverse_iterator
Definition: JHashMap.hh:86
JDynamics(const JDetector &detector, const double Tmax_s)
Constructor.
Definition: JDynamics.hh:527
General purpose class for hash map of unique keys.
Definition: JHashMap.hh:71
JTOOLS::JHashMap< int, function_type > data_type
Definition: JDynamics.hh:328
const_reverse_iterator rend() const
begin of reverse of calibration data
Definition: JDynamics.hh:176
JQuaternion3D getQuaternion(const JQuaternion &Q)
Get quaternion.
JTOOLS::JPolfitFunction1D< NUMBER_OF_POINTS, NUMBER_OF_DEGREES, element_type, JTOOLS::JCollection > function_type
Definition: JDynamics.hh:109
function_type::collection_type::container_type container_type
Definition: JDynamics.hh:110
General purpose class for hash map of unique elements.
Detector data structure.
Definition: JDetector.hh:80
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:670
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:378
JOrientation orientation
Definition: JDynamics.hh:566
const_iterator begin() const
begin of calibration data
Definition: JDynamics.hh:173
static JRotation getRotation
Function object to get rotation matrix to go from first to second module.
Lookup table for PMT addresses in detector.
Interface of object iteration for a single data type.
data_type::const_reverse_iterator const_reverse_iterator
Definition: JDynamics.hh:331
Data structure for detector geometry and calibration.
double getXmin() const
Get minimal abscissa.
Definition: JDynamics.hh:495
static double getXmax(const JTOOLS::JHashMap< JKey_t, JValue_t, JEvaluator_t > &input)
Get maximal abscissa of mapped function objects.
Definition: JDynamics.hh:82
JTOOLS::JHashMap< int, JGEOMETRY3D::JQuaternion3D > buffer_type
Definition: JDynamics.hh:112
static double getXmin(const JTOOLS::JHashMap< JKey_t, JValue_t, JEvaluator_t > &input)
Get minimal abscissa of mapped function objects.
Definition: JDynamics.hh:61
const JModuleAddressMap & get(const int id) const
Get module address map.
const JGeometry & getGeometry() const
Get detector geometry.
Definition: JDynamics.hh:445
JACOUSTICS::JGeometry JGeometry
Definition: JDynamics.hh:321
void load(JObjectIterator_t &input)
Load calibration data.
Definition: JDynamics.hh:541
const_reverse_iterator rbegin() const
begin of reverse of calibration data
Definition: JDynamics.hh:175
data_type::const_iterator const_iterator
Definition: JDynamics.hh:330
virtual mapped_type & get(typename JClass< key_type >::argument_type key) override
Get mapped value.
Definition: JHashMap.hh:118
double UNIXTimeStop
stop time
void load(JObjectIterator< JACOUSTICS::JEvt > &input)
Load calibration data.
Definition: JDynamics.hh:353
JVector3D & sub(const JVector3D &vector)
Subtract vector.
Definition: JVector3D.hh:158
const_iterator end() const
end of calibration data
Definition: JDynamics.hh:377
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:196
Acoustic event fit.
double getXmin() const
Get minimal abscissa.
Definition: JDynamics.hh:285
JPosition(const JDetector &detector, const double Tmax_s)
Constructor.
Definition: JDynamics.hh:340
int getID() const
Get identifier.
Definition: JObjectID.hh:50
void setTmax(const double Tmax_s)
Set applicability period.
Definition: JDynamics.hh:398
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
void update(JDetector &detector, const double t1_s)
Calibrate given detector at given time.
Definition: JDynamics.hh:238
double getCoverage(const JDetector &detector, const double t1_s) const
Get coverage.
Definition: JDynamics.hh:208
void load(JObjectIterator< JCOMPASS::JOrientation > &input)
Load calibration data.
Definition: JDynamics.hh:154
JDetectorAddressMap & getDetectorAddressMap()
Get detector address map.
double getXmax() const
Get maximal abscissa.
Definition: JDynamics.hh:506
JTOOLS::JElement2D< double, JGEOMETRY3D::JQuaternion3D > element_type
Definition: JDynamics.hh:106
data_type::const_iterator const_iterator
Definition: JDynamics.hh:115
const_iterator end() const
end of calibration data
Definition: JDynamics.hh:174
data_type::const_reverse_iterator const_reverse_iterator
Definition: JDynamics.hh:116
Data structure for unit quaternion in three dimensions.
bool empty() const
empty
Definition: JDynamics.hh:375
double getTmax() const
Get applicability period.
Definition: JDynamics.hh:184
Dynamic detector calibration.
Definition: JDynamics.hh:52
double UNIXTimeStart
start time
double getTmax() const
Get applicability period.
Definition: JDynamics.hh:387
JPosition3D getPiezoPosition()
Get relative position of piezo in optical module.
JTOOLS::JHashMap< int, function_type > data_type
Definition: JDynamics.hh:113
const_iterator begin() const
begin of calibration data
Definition: JDynamics.hh:376
JOrientation(const JDetector &detector, const double Tmax_s)
Constructor.
Definition: JDynamics.hh:125
Template class for polynomial interpolation in 1D.
Definition: JPolfit.hh:144
const JDetector & operator()(const double t1_s)
Get detector calibrated at given time.
Definition: JDynamics.hh:554
2D Element.
Definition: JElement.hh:46
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:85
double getCoverage(const JDetector &detector, const double t1_s) const
Get coverage.
Definition: JDynamics.hh:411
Model for fit to acoutsics data.
Acoustic event fit.
double getXmax() const
Get maximal abscissa.
Definition: JDynamics.hh:296
bool has(const T &value) const
Test whether given value is present.
number of degrees for interpolation
Definition: JDynamics.hh:103
JVector3D & add(const JVector3D &vector)
Add vector.
Definition: JVector3D.hh:142
const JModule & getModule(const JDetector &detector, const JModuleLocation &location)
find module with a given string and floor number
JTOOLS::JPolfitFunction1D< NUMBER_OF_POINTS, NUMBER_OF_DEGREES, element_type, JTOOLS::JCollection > function_type
Definition: JDynamics.hh:326
JQuaternion3D & conjugate()
Conjugate quaternion.
container_type::iterator iterator
Definition: JHashMap.hh:87