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  * Calibrate given detector at given time.
203  *
204  * \param detector detector (I/O)
205  * \param t1_s time [s]
206  */
207  void update(JDetector& detector, const double t1_s)
208  {
209  using namespace std;
210  using namespace JPP;
211 
212  if (!calibration.empty()) {
213 
214  if (fabs(t1_s - t0_s) > Tmax_s) {
215 
216  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
217 
218  if (module->getFloor() != 0 && calibration.has(module->getID()) && buffer.has(module->getID())) {
219 
220  const function_type& f1 = calibration.get(module->getID());
221 
222  if (!f1.empty()) {
223 
224  if (t1_s >= f1.getXmin() && t1_s <= f1.getXmax()) {
225 
226  JQuaternion3D Q0 = buffer[module->getID()];
227  JQuaternion3D Q1 = f1(t1_s);
228 
229  const JPosition3D center = module->getPosition();
230 
231  module->sub(center);
232 
233  module->rotate(Q1 * Q0.conjugate());
234 
235  module->add(center);
236 
237  buffer[module->getID()] = Q1;
238  }
239  }
240  }
241  }
242 
243  t0_s = t1_s;
244  }
245  }
246  }
247 
248 
249  /**
250  * Get minimal abscissa.
251  *
252  * \return minimal abscissa
253  */
254  double getXmin() const
255  {
256  return JDynamics::getXmin(this->calibration);
257  }
258 
259 
260  /**
261  * Get maximal abscissa.
262  *
263  * \return maximal abscissa
264  */
265  double getXmax() const
266  {
267  return JDynamics::getXmax(this->calibration);
268  }
269 
270 
271  /**
272  * Write calibration to output stream.
273  *
274  * \param out output stream
275  * \param calibration
276  * \return output stream
277  */
278  friend inline std::ostream& operator<<(std::ostream& out, const JOrientation& calibration)
279  {
280  if (!calibration.empty())
281  return out << "[" << FIXED(15,0) << calibration.getXmin()
282  << "," << FIXED(15,0) << calibration.getXmax()
283  << "]";
284  else
285  return out << "[,]";
286  }
287 
288  protected:
291  double Tmax_s;
292 
293  private:
294  double t0_s;
295  };
296 
297 
298  /**
299  * Dynamic position calibration.
300  */
301  struct JPosition {
302 
303  enum {
304  NUMBER_OF_POINTS = 7, //!< number of points for interpolation
305  NUMBER_OF_DEGREES = 2 //!< number of degrees for interpolation
306  };
307 
309 
314 
316 
319 
320 
321  /**
322  * Constructor.
323  *
324  * \param detector detector
325  * \param Tmax_s applicability period of calibration [s]
326  */
328  const double Tmax_s) :
329  geometry(detector),
330  Tmax_s(Tmax_s),
331  t0_s(std::numeric_limits<double>::lowest())
332  {}
333 
334 
335  /**
336  * Load calibration data.
337  *
338  * \param input calibration data
339  */
341  {
342  using namespace JPP;
343 
344  t0_s = std::numeric_limits<double>::lowest();
345 
346  while (input.hasNext()) {
347 
348  const JACOUSTICS::JEvt* evt = input.next();
349  const double t1_s = 0.5 * (evt->UNIXTimeStart + evt->UNIXTimeStop);
350 
351  for (JACOUSTICS::JEvt::const_iterator i = evt->begin(); i != evt->end(); ++i) {
352  calibration[i->id][t1_s] = JMODEL::JString(i->tx, i->ty);
353  }
354  }
355 
356  for (data_type::iterator i = calibration.begin(); i != calibration.end(); ++i) {
357  i->second.compile();
358  }
359  }
360 
361 
362  bool empty() const { return calibration.empty(); } //!< empty
363  const_iterator begin() const { return calibration.begin(); } //!< begin of calibration data
364  const_iterator end() const { return calibration.end(); } //!< end of calibration data
365  const_reverse_iterator rbegin() const { return calibration.rbegin(); } //!< begin of reverse of calibration data
366  const_reverse_iterator rend() const { return calibration.rend(); } //!< begin of reverse of calibration data
367 
368 
369  /**
370  * Get applicability period.
371  *
372  * \return applicability period [s]
373  */
374  double getTmax() const
375  {
376  return Tmax_s;
377  }
378 
379 
380  /**
381  * Set applicability period.
382  *
383  * \param Tmax_s applicability period [s]
384  */
385  void setTmax(const double Tmax_s)
386  {
387  this->Tmax_s = Tmax_s;
388  }
389 
390 
391  /**
392  * Get detector geometry.
393  *
394  * \return detector geometry
395  */
396  const JGeometry& getGeometry() const
397  {
398  return geometry;
399  }
400 
401 
402  /**
403  * Calibrate given detector at given time.
404  *
405  * \param detector detector (I/O)
406  * \param t1_s time [s]
407  */
408  void update(JDetector& detector, const double t1_s)
409  {
410  using namespace std;
411  using namespace JPP;
412 
413  if (!calibration.empty()) {
414 
415  if (fabs(t1_s - t0_s) > Tmax_s) {
416 
418 
419  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
420 
421  if (module->getFloor() != 0) {
422 
423  if (!buffer.has(module->getString())) {
424  if (calibration.has(module->getString())) {
425  buffer[module->getString()] = calibration[module->getString()](t1_s);
426  }
427  }
428 
429  if (buffer.has(module->getString())) {
430  module->set(geometry[module->getString()].getPosition(buffer[module->getString()], module->getFloor()) - getPiezoPosition());
431  }
432  }
433  }
434 
435  t0_s = t1_s;
436  }
437  }
438  }
439 
440 
441  /**
442  * Get minimal abscissa.
443  *
444  * \return minimal abscissa
445  */
446  double getXmin() const
447  {
448  return JDynamics::getXmin(this->calibration);
449  }
450 
451 
452  /**
453  * Get maximal abscissa.
454  *
455  * \return maximal abscissa
456  */
457  double getXmax() const
458  {
459  return JDynamics::getXmax(this->calibration);
460  }
461 
462 
463  /**
464  * Write calibration to output stream.
465  *
466  * \param out output stream
467  * \param calibration
468  * \return output stream
469  */
470  friend inline std::ostream& operator<<(std::ostream& out, const JPosition& calibration)
471  {
472  if (!calibration.empty())
473  return out << "[" << FIXED(15,0) << calibration.getXmin()
474  << "," << FIXED(15,0) << calibration.getXmax()
475  << "]";
476  else
477  return out << "[,]";
478  }
479 
480  protected:
483  double Tmax_s;
484 
485  private:
486  double t0_s;
487  };
488 
489 
490  /**
491  * Constructor.
492  *
493  * \param detector detector
494  * \param Tmax_s applicability period of calibration [s]
495  */
497  const double Tmax_s) :
498  detector (detector),
499  orientation(detector, Tmax_s),
500  position (detector, Tmax_s)
501  {}
502 
503 
504  /**
505  * Load calibration data.
506  *
507  * \param input detector calibration data
508  */
509  template<class JObjectIterator_t>
510  void load(JObjectIterator_t& input)
511  {
512  try { orientation.load(dynamic_cast<JObjectIterator<JCOMPASS::JOrientation>&>(input)); } catch(const std::exception& error) {}
513  try { position .load(dynamic_cast<JObjectIterator<JACOUSTICS::JEvt>&> (input)); } catch(const std::exception& error) {}
514  }
515 
516 
517  /**
518  * Get detector calibrated at given time.
519  *
520  * \param t1_s time [s]
521  * \return detector
522  */
523  const JDetector& operator()(const double t1_s)
524  {
525  orientation.update(this->detector, t1_s);
526  position .update(this->detector, t1_s);
527 
528  return this->detector;
529  }
530 
531  protected:
533 
534  public:
537  };
538 }
539 
540 #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:310
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:366
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:408
Dynamic orientation calibration.
Definition: JDynamics.hh:99
Dynamic position calibration.
Definition: JDynamics.hh:301
container_type::const_reverse_iterator const_reverse_iterator
Definition: JHashMap.hh:86
JDynamics(const JDetector &detector, const double Tmax_s)
Constructor.
Definition: JDynamics.hh:496
General purpose class for hash map of unique keys.
Definition: JHashMap.hh:71
friend std::ostream & operator<<(std::ostream &out, const JOrientation &calibration)
Write calibration to output stream.
Definition: JDynamics.hh:278
JTOOLS::JHashMap< int, function_type > data_type
Definition: JDynamics.hh:315
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:365
JOrientation orientation
Definition: JDynamics.hh:535
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.
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:445
Interface of object iteration for a single data type.
data_type::const_reverse_iterator const_reverse_iterator
Definition: JDynamics.hh:318
Data structure for detector geometry and calibration.
double getXmin() const
Get minimal abscissa.
Definition: JDynamics.hh:446
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:396
JACOUSTICS::JGeometry JGeometry
Definition: JDynamics.hh:308
void load(JObjectIterator_t &input)
Load calibration data.
Definition: JDynamics.hh:510
const_reverse_iterator rbegin() const
begin of reverse of calibration data
Definition: JDynamics.hh:175
data_type::const_iterator const_iterator
Definition: JDynamics.hh:317
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:340
JVector3D & sub(const JVector3D &vector)
Subtract vector.
Definition: JVector3D.hh:158
const_iterator end() const
end of calibration data
Definition: JDynamics.hh:364
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:254
JPosition(const JDetector &detector, const double Tmax_s)
Constructor.
Definition: JDynamics.hh:327
int getID() const
Get identifier.
Definition: JObjectID.hh:50
void setTmax(const double Tmax_s)
Set applicability period.
Definition: JDynamics.hh:385
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:207
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:457
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:362
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:374
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:363
JOrientation(const JDetector &detector, const double Tmax_s)
Constructor.
Definition: JDynamics.hh:125
Template class for polynomial interpolation in 1D.
Definition: JPolfit.hh:144
friend std::ostream & operator<<(std::ostream &out, const JPosition &calibration)
Write calibration to output stream.
Definition: JDynamics.hh:470
const JDetector & operator()(const double t1_s)
Get detector calibrated at given time.
Definition: JDynamics.hh:523
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
Model for fit to acoutsics data.
Acoustic event fit.
double getXmax() const
Get maximal abscissa.
Definition: JDynamics.hh:265
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:313
JQuaternion3D & conjugate()
Conjugate quaternion.
container_type::iterator iterator
Definition: JHashMap.hh:87