Jpp  18.0.0-rc.3
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JGeometry.hh
Go to the documentation of this file.
1 #ifndef __JACOUSTICS__JGEOMETRY__
2 #define __JACOUSTICS__JGEOMETRY__
3 
4 #include <ostream>
5 #include <iomanip>
6 #include <vector>
7 #include <map>
8 #include <algorithm>
9 #include <cmath>
10 
11 #include "JLang/JException.hh"
12 #include "JLang/JPredicate.hh"
13 #include "JLang/JComparator.hh"
14 #include "JLang/JManip.hh"
15 
17 
18 #include "JTools/JHashMap.hh"
19 
20 #include "JDetector/JLocation.hh"
21 #include "JDetector/JHydrophone.hh"
22 #include "JDetector/JDetector.hh"
24 
27 #include "JAcoustics/JMechanics.hh"
28 #include "JAcoustics/JModel.hh"
29 
30 
31 /**
32  * \file
33  *
34  * Acoustic geometries.
35  * \author mdejong
36  */
37 namespace JACOUSTICS {}
38 namespace JPP { using namespace JACOUSTICS; }
39 
40 namespace JACOUSTICS {
41 
43  using JLANG::JNoValue;
46  using JTOOLS::JHashMap;
48 
49 
50  /**
51  * Auxiliary namespace to encapsulate different geometries.\n
52  */
53  namespace JGEOMETRY {
54 
55  /**
56  * Floor geometry.
57  */
58  struct JFloor {
59  /**
60  * Default constructor.
61  */
62  JFloor() :
63  height(0.0)
64  {}
65 
66 
67  /**
68  * Constructor.\n
69  * The height refers to the maximal distance from the top of the so-called T-bar of the parent string to the actual sensor.
70  *
71  * \param height height
72  */
73  JFloor(const double height) :
74  height(height)
75  {}
76 
77 
78  /**
79  * Get height of this floor.\n
80  * The height refers to the maximal distance from
81  * the fixed reference point of the parent string.
82  *
83  * \return height
84  */
85  double getHeight() const
86  {
87  return height;
88  }
89 
90 
91  /**
92  * Write floor parameters to output stream.
93  *
94  * \param out output stream
95  * \param floor floor
96  * \return output stream
97  */
98  friend inline std::ostream& operator<<(std::ostream& out, const JFloor& floor)
99  {
100  return out << FIXED(7,2) << floor.getHeight();
101  }
102 
103 
104  protected:
105  double height;
106  };
107 
108 
109  /**
110  * String geometry.
111  *
112  * This data structure provides for the implementation of the dynamical geometry of a detector string.\n
113  * In this,
114  * the position of floor > 0 corresponds to the piezo sensor which is mounted in the optical module and \n
115  * the position of floor = 0 to the hydrophone which is optionally mounted on the anchor.\n
116  * The reference position of a string corresponds to the top of the T-bar located on the anchor.\n
117  * The position of a piezo sensor depends on the tilt of the string and the mechanical model.\n
118  * The position of the hydrophone is fixed.
119  * Its value is relative to the reference position of the string.
120  */
121  struct JString :
122  public JPosition3D,
123  public std::vector<JFloor>
124  {
128 
129  static constexpr double PRECISION_M = 1.0e-4; //!< precision of height evaluation [m]
130 
131  /**
132  * Get approximate length of string.
133  *
134  * \param parameters parameters
135  * \param mechanics mechanics
136  * \param height height
137  * \return length
138  */
139  static inline double getLength(const JMODEL::JString& parameters,
140  const JMechanics& mechanics,
141  const double height)
142  {
143  const double t2 = parameters.getLengthSquared();
144 
145  return 0.5*t2*mechanics.b * log(1.0 - mechanics.a*height) + sqrt(1.0 + t2) * height;
146  }
147 
148 
149  /**
150  * Get approximate height of string.
151  *
152  * \param parameters parameters
153  * \param mechanics mechanics
154  * \param length length
155  * \param precision precision
156  * \return height
157  */
158  static inline double getHeight(const JMODEL::JString& parameters,
159  const JMechanics& mechanics,
160  const double length,
161  const double precision = PRECISION_M)
162  {
163  const size_t MAXIMUM_NUMBER_OF_ITERATIONS = 10;
164 
165  const double ts = 0.5 * parameters.getLengthSquared();
166 
167  double z = length; // start value
168 
169  for (size_t i = 0; i != MAXIMUM_NUMBER_OF_ITERATIONS; ++i) {
170 
171  const double ls = getLength(parameters, mechanics, z);
172 
173  if (fabs(ls - length) <= precision) {
174  break;
175  }
176 
177  z -= (ls - length) / (1.0 + ts * (1.0 - mechanics.a*mechanics.b / (1.0 - mechanics.a*z)));
178  }
179 
180  return z;
181  }
182 
183  /**
184  * Get position at given height according to given string model parameters and mechanics.
185  *
186  * \param parameters parameters
187  * \param mechanics mechanics
188  * \param height height
189  * \return position
190  */
192  const JMechanics& mechanics,
193  const double height)
194  {
195  const double h1 = height * (1.0 + parameters.vs);
196  const double z1 = mechanics.getHeight(h1);
197 
198  return JPosition3D(parameters.tx * z1 + parameters.tx2 * h1*h1,
199  parameters.ty * z1 + parameters.ty2 * h1*h1,
200  getHeight(parameters, mechanics, h1));
201  }
202 
203 
204  /**
205  * Default constructor.
206  */
208  has_hydrophone(false)
209  {}
210 
211 
212  /**
213  * Constructor.
214  *
215  * The given position corresponds to the reference point of the string from
216  * which the positions of the piezo sensors and hydrophone are calculated.
217  *
218  * \param position position
219  * \param mechanics mechanical model parameters
220  */
221  JString(const JVector3D& position,
222  const JMechanics& mechanics) :
223  JPosition3D(position),
224  has_hydrophone(false),
225  mechanics(mechanics)
226  {}
227 
228 
229  /**
230  * Constructor.
231  *
232  * The given position corresponds to the reference position of the string from
233  * which the positions of the piezo sensors and hydrophone are calculated.
234  *
235  * The template parameter should correspond to a data type which implements the following policy methods.
236  * <pre>
237  * int %getFloor();
238  * JGEOMETRY3d::JPosition3D %getPosition();
239  * </pre>
240  * In this, the position should correspond to the center of the optical module.
241  *
242  * Note that the position of the piezo is offset by JDETECTOR::getPiezoPosition with respect to the center of the optical module.\n
243  * The position of the hydrophone should separately be set.
244  *
245  * \param position position
246  * \param mechanics mechanical model parameters
247  * \param __begin begin of optical modules
248  * \param __end end of optical modules
249  */
250  template<class T>
251  JString(const JVector3D& position,
252  const JMechanics& mechanics,
253  T __begin,
254  T __end) :
255  JPosition3D(position),
256  has_hydrophone(false),
257  mechanics(mechanics)
258  {
259  for (T i = __begin; i != __end; ++i) {
260  (*this)[i->getFloor()] = this->getDistance(i->getPosition() + JDETECTOR::getPiezoPosition());
261  }
262  }
263 
264 
265  /**
266  * Get mechanical model parameters.
267  *
268  * \return mechanical model parameters
269  */
270  const JMechanics& getMechanics() const
271  {
272  return mechanics;
273  }
274 
275 
276  /**
277  * Get floor data.
278  *
279  * \param floor floor number
280  * \return floor data
281  */
282  JFloor& operator[](size_t floor)
283  {
284  if (floor >= this->size()) {
285  this->resize(floor + 1);
286  }
287 
288  return static_cast<std::vector<JFloor>&>(*this)[floor];
289  }
290 
291 
292  /**
293  * Check if this string has given floor.
294  *
295  * \param floor floor
296  * \return true if floor present; else false
297  */
298  bool hasFloor(size_t floor) const
299  {
300  if (floor == 0)
301  return has_hydrophone;
302  else
303  return (floor < this->size());
304  }
305 
306 
307  /**
308  * Get height of given floor.
309  *
310  * \param floor floor
311  * \return height
312  */
313  double getHeight(const size_t floor) const
314  {
315  if (floor == 0) {
316 
317  if (has_hydrophone) {
318  return hydrophone.getZ();
319  }
320 
321  } else if (floor < this->size()) {
322 
323  return (*this)[floor].getHeight();
324  }
325 
326  THROW(JValueOutOfRange, "Invalid floor " << floor);
327  }
328 
329 
330  /**
331  * Get position of given floor according to given string model parameters.
332  *
333  * \param parameters parameters
334  * \param floor floor
335  * \return position
336  */
338  const size_t floor) const
339  {
340  if (floor == 0) {
341 
342  if (has_hydrophone) {
343  return this->getPosition() + hydrophone.getPosition();
344  }
345 
346  } else if (floor < this->size()) {
347 
348  return this->getPosition() + getPosition(parameters, mechanics, (*this)[floor].getHeight());
349  }
350 
351  THROW(JValueOutOfRange, "Invalid floor " << floor);
352  }
353 
354 
355  /**
356  * Get position of given floor according to default string model parameters.
357  *
358  * \param floor floor
359  * \return position
360  */
361  JPosition3D getPosition(const size_t floor) const
362  {
363  return getPosition(JMODEL::JString(), floor);
364  }
365 
366 
367  /**
368  * Get distance between given position and floor according to given string model parameters.
369  *
370  * \param parameters parameters
371  * \param position position
372  * \param floor floor
373  * \return distance
374  */
376  const JVector3D& position,
377  const size_t floor) const
378  {
379  return this->getPosition(parameters, floor).getDistance(position);
380  }
381 
382 
383  /**
384  * Get model gradient of distance between given position and floor according to given string model parameters.
385  *
386  * \param parameters parameters
387  * \param position position
388  * \param floor floor
389  * \return gradient
390  */
392  const JVector3D& position,
393  const size_t floor) const
394  {
395  if (floor == 0) {
396 
397  return JMODEL::JString();
398 
399  } else if (floor < this->size()) {
400 
401  const JPosition3D pos = this->getPosition(parameters, floor);
402  const double height = (*this)[floor].getHeight();
403  const double h1 = height * (1.0 + parameters.vs);
404  const double z1 = this->mechanics.getHeight(h1);
405 
406  const double tx = parameters.tx;
407  const double ty = parameters.ty;
408  const double tz = sqrt(1.0 - tx*tx - ty*ty);
409 
410  const double dx = pos.getX() - position.getX();
411  const double dy = pos.getY() - position.getY();
412  const double dz = pos.getZ() - position.getZ();
413 
414  const double D = sqrt(dx*dx + dy*dy + dz*dz);
415  const double vw = 1.0 - mechanics.a * mechanics.b / (1.0 - mechanics.a*h1);
416  const double vs = 1.0 + 0.5 * parameters.getLengthSquared() * vw;
417 
418  return JMODEL::JString(z1 * dx / D - height * (tx / tz) * dz / D,
419  z1 * dy / D - height * (ty / tz) * dz / D,
420  h1*h1 * dx / D,
421  h1*h1 * dy / D,
422  height * vw * (tx * dx + ty * dy) / D + h1 * (dz / vs) / D);
423 
424  } else {
425 
426  THROW(JValueOutOfRange, "Invalid floor " << floor);
427  }
428  }
429 
430 
431  /**
432  * Write string parameters to output stream.
433  *
434  * \param out output stream
435  * \param string string
436  * \return output stream
437  */
438  friend inline std::ostream& operator<<(std::ostream& out, const JString& string)
439  {
440  using namespace std;
441 
442  for (size_t i = 0; i != string.size(); ++i) {
443  if (string.hasFloor(i)) {
444  out << setw(2) << i << ' '
445  << FIXED(7,3) << string[i] << " | "
446  << string.getPosition(i) << ' '
447  << string.mechanics << endl;
448  }
449  }
450 
451  return out;
452  }
453 
454 
455  /**
456  * Hydrophone.
457  *
458  * The position of the hydrophone is relative to the reference position of the string.
459  */
462 
463  protected:
464  /**
465  * Mechanical data.
466  */
468  };
469 
470 
471  /**
472  * Detector geometry.
473  */
474  struct JDetector :
475  public JHashMap<int, JString>
476  {
477  /**
478  * Default constructor.
479  */
481  {}
482 
483 
484  /**
485  * Constructor.
486  *
487  * Note that the positions of the base modules correspond to the reference position of the string.\n
488  * As a consequence, a base module (i.e.\ floor = 0) is required for each string.\n
489  * Missing base modules should therefore be added beforehand (e.g.\ using application JDetectorDB.cc).
490  *
491  * Note that if the position of a hydrophone is not available,
492  * it is assumed that there is no hydrophone on that string.\n
493  * If the position of the hydrophone is manually set,
494  * the corresponding parameter JString::has_hydrophone should be set to <tt>true</tt>.
495  *
496  * \param detector detector
497  * \param hydrophones container with data of hydrophones
498  */
501  {
502  using namespace std;
503  using namespace JPP;
504 
505  map<int, vector<module_type> > buffer; // string -> modules
506 
507  for (JDETECTOR::JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
508  buffer[module->getString()].push_back(module_type(module->getLocation(), module->getPosition()));
509  }
510 
511  for (map<int, vector<module_type> >::iterator i = buffer.begin(); i != buffer.end(); ++i) {
512 
513  sort(i->second.begin(), i->second.end(), make_comparator(&module_type::getFloor));
514 
515  if (i->second[0].getFloor() == 0) {
516 
517  vector<module_type>::iterator p = i->second.begin();
518 
519  ++p;
520 
521  (*this)[i->first] = JGEOMETRY::JString(i->second[0].getPosition(), getMechanics(i->first), p, i->second.end());
522 
523  try {
524 
525  (*this)[i->first].hydrophone = getPosition(hydrophones.begin(),
526  hydrophones.end(),
528 
529  (*this)[i->first].has_hydrophone = true;
530  }
531  catch(const exception&) {
532  (*this)[i->first].has_hydrophone = false;
533  }
534 
535  } else {
536 
537  THROW(JNoValue, "No floor 0 in string " << i->first << "; use e.g. JDetectorDB -W.");
538  }
539  }
540  }
541 
542 
543  /**
544  * Check if this detector has given string.
545  *
546  * \param string string
547  * \return true if string present; else false
548  */
549  bool hasString(int string) const
550  {
551  return this->has(string);
552  }
553 
554 
555  /**
556  * Check if this detector has given location.
557  *
558  * \param location location
559  * \return true if location present; else false
560  */
561  bool hasLocation(const JLocation& location) const
562  {
563  return this->hasString(location.getString()) && (*this)[location.getString()].hasFloor(location.getFloor());
564  }
565 
566 
567  /**
568  * Write detector parameters to output stream.
569  *
570  * \param out output stream
571  * \param detector detector
572  * \return output stream
573  */
574  friend inline std::ostream& operator<<(std::ostream& out, const JDetector& detector)
575  {
576  using namespace std;
577 
578  for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
579  out << setw(4) << i->first << endl << i->second;
580  }
581 
582  return out;
583  }
584 
585 
586  /**
587  * Auxiliary data structure for module location and position.
588  */
589  struct module_type :
590  public JLocation,
591  public JPosition3D
592  {
593  /**
594  * Constructor.
595  *
596  * \param location module location
597  * \param position module position
598  */
599  module_type(const JLocation& location,
600  const JPosition3D& position) :
601  JLocation (location),
602  JPosition3D(position)
603  {}
604 
605 
606  /**
607  * Less-than operator.
608  *
609  * \param first first module
610  * \param second second module
611  * \return true if floor of first module less than that of second; else false
612  */
613  friend inline bool operator<(const module_type& first, const module_type& second)
614  {
615  return first.getFloor() < second.getFloor();
616  }
617  };
618  };
619  }
620 
621 
622  /**
623  * Type definition of detector geometry.
624  */
626 }
627 
628 #endif
static double getHeight(const JMODEL::JString &parameters, const JMechanics &mechanics, const double length, const double precision=PRECISION_M)
Get approximate height of string.
Definition: JGeometry.hh:158
Mechanical modelling of string.
Exceptions.
JPredicate< JResult_t T::*, JComparison::eq > make_predicate(JResult_t T::*member, const JResult_t value)
Helper method to create predicate for data member.
Definition: JPredicate.hh:128
const JMechanics & getMechanics() const
Get mechanical model parameters.
Definition: JGeometry.hh:270
JComparator< JResult_t T::*, JComparison::lt > make_comparator(JResult_t T::*member)
Helper method to create comparator between values of data member.
int getFloor() const
Get floor number.
Definition: JLocation.hh:145
static JPosition3D getPosition(const JMODEL::JString &parameters, const JMechanics &mechanics, const double height)
Get position at given height according to given string model parameters and mechanics.
Definition: JGeometry.hh:191
General purpose class for hash map of unique keys.
Definition: JHashMap.hh:72
friend std::ostream & operator<<(std::ostream &out, const JFloor &floor)
Write floor parameters to output stream.
Definition: JGeometry.hh:98
General purpose class for hash map of unique elements.
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
double getHeight(const size_t floor) const
Get height of given floor.
Definition: JGeometry.hh:313
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
Definition: diff-Tuna.sh:38
JString()
Default constructor.
Definition: JGeometry.hh:207
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
double getDistance(const JVector3D &pos) const
Get distance to point.
Definition: JVector3D.hh:270
JFloor()
Default constructor.
Definition: JGeometry.hh:62
Data structure for detector geometry and calibration.
JPosition3D getPosition(const size_t floor) const
Get position of given floor according to default string model parameters.
Definition: JGeometry.hh:361
Auxiliary data structure for module location and position.
Definition: JGeometry.hh:589
static constexpr double PRECISION_M
precision of height evaluation [m]
Definition: JGeometry.hh:129
then echo The file $DIR KM3NeT_00000001_00000000 root already please rename or remove it first
Data structure for hydrophone.
JFloor(const double height)
Constructor.
Definition: JGeometry.hh:73
friend std::ostream & operator<<(std::ostream &out, const JString &string)
Write string parameters to output stream.
Definition: JGeometry.hh:438
bool hasString(int string) const
Check if this detector has given string.
Definition: JGeometry.hh:549
JGEOMETRY::JDetector JGeometry
Type definition of detector geometry.
Definition: JGeometry.hh:625
Exception for missing value.
Definition: JException.hh:198
module_type(const JLocation &location, const JPosition3D &position)
Constructor.
Definition: JGeometry.hh:599
JString(const JVector3D &position, const JMechanics &mechanics)
Constructor.
Definition: JGeometry.hh:221
Detector file.
Definition: JHead.hh:226
Acoustics support kit.
std::vector< double > vs
Detector support kit.
Data structure for vector in three dimensions.
Definition: JVector3D.hh:34
Logical location of module.
Definition: JLocation.hh:37
double getDistance(const JMODEL::JString &parameters, const JVector3D &position, const size_t floor) const
Get distance between given position and floor according to given string model parameters.
Definition: JGeometry.hh:375
Acoustics toolkit.
JPosition3D getPosition(const JMODEL::JString &parameters, const size_t floor) const
Get position of given floor according to given string model parameters.
Definition: JGeometry.hh:337
JMODEL::JString getGradient(const JMODEL::JString &parameters, const JVector3D &position, const size_t floor) const
Get model gradient of distance between given position and floor according to given string model param...
Definition: JGeometry.hh:391
do set_variable OUTPUT_DIRECTORY $WORKDIR T
static double getLength(const JMODEL::JString &parameters, const JMechanics &mechanics, const double height)
Get approximate length of string.
Definition: JGeometry.hh:139
double a
0 &lt;= a &lt; (maximal height)⁻1; [m^-1]
Definition: JMechanics.hh:100
double getHeight(const double height) const
Get effective height for given actual height.
Definition: JMechanics.hh:68
double getLengthSquared() const
Get length squared.
double getY() const
Get y position.
Definition: JVector3D.hh:104
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
double getLength() const
Get length.
Definition: JVector3D.hh:246
Logical location of module.
I/O manipulators.
JString(const JVector3D &position, const JMechanics &mechanics, T __begin, T __end)
Constructor.
Definition: JGeometry.hh:251
double b
0 &lt;= b; [m]
Definition: JMechanics.hh:101
JDetector(const JDETECTOR::JDetector &detector, const std::vector< JDETECTOR::JHydrophone > &hydrophones=std::vector< JDETECTOR::JHydrophone >())
Constructor.
Definition: JGeometry.hh:499
JFloor & operator[](size_t floor)
Get floor data.
Definition: JGeometry.hh:282
bool hasFloor(size_t floor) const
Check if this string has given floor.
Definition: JGeometry.hh:298
int getString() const
Get string number.
Definition: JLocation.hh:134
friend std::ostream & operator<<(std::ostream &out, const JDetector &detector)
Write detector parameters to output stream.
Definition: JGeometry.hh:574
JPosition3D hydrophone
Hydrophone.
Definition: JGeometry.hh:460
JMechanics mechanics
Mechanical data.
Definition: JGeometry.hh:467
JPosition3D getPiezoPosition()
Get relative position of piezo in optical module.
JPosition3D()
Default constructor.
Definition: JPosition3D.hh:48
friend bool operator<(const module_type &first, const module_type &second)
Less-than operator.
Definition: JGeometry.hh:613
JDetector()
Default constructor.
Definition: JGeometry.hh:480
double getX() const
Get x position.
Definition: JVector3D.hh:94
bool hasLocation(const JLocation &location) const
Check if this detector has given location.
Definition: JGeometry.hh:561
Auxiliary data structure to list files in directory.
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
double getHeight() const
Get height of this floor.
Definition: JGeometry.hh:85
Exception for accessing a value in a collection that is outside of its range.
Definition: JException.hh:162
container_type::const_iterator const_iterator
Definition: JHashMap.hh:86
then cat $TRIPOD_INITIAL<< EOF1 256877.5 4743716.7-2438.42 256815.5 4743395.0-2435.53 257096.2 4743636.0-2439.5EOFfiif[[!-f $DETECTOR]];then JEditDetector-a $DETECTOR_INITIAL-s"-1 addz -6.9"-o $DETECTOR--!eval`JPrintDetector-a $DETECTOR-O SUMMARY`for STRING in ${STRINGS[*]};do set_variable MODULE`getModule-a $DETECTOR-L"$STRING 0"`JEditDetector-a $DETECTOR-M"$MODULE setz -2.9"-o $DETECTOR--!donefiif[[!-f $TRIPOD]];then cp-p $TRIPOD_INITIAL $TRIPODfiJAcoustics.sh $DETECTOR_IDcat > acoustics_trigger_parameters txt<< EOFQ=0.0;TMax_s=0.020;numberOfHits=90;EOFJAcousticsEventBuilder.sh $DETECTOR $RUNS[*]INPUT_FILES=(`ls KM3NeT_ ${(l:8::0::0:) DETECTOR_ID}_0 *${^RUNS}_event.root`) cd $WORKDIRsource $JPP_DIR/examples/JAcoustics/acoustics-fit-toolkit.shif[!$HOMEDIR-ef $WORKDIR];then cp-p $HOMEDIR/$DETECTOR $WORKDIR cp-p $HOMEDIR/${^ACOUSTICS_KEYS}.txt $WORKDIR cp-p $HOMEDIR/${^INPUT_FILES}$WORKDIRfitimer_startinitialise stage_b 1 0 100.0e-6 10.0 0.002 0.1 0 > &stage log
Model for fit to acoutsics data.
do echo Generating $dir eval D
Definition: JDrawLED.sh:53
double getZ() const
Get z position.
Definition: JVector3D.hh:115
Auxiliary data structure for parameters of mechanical model.
Definition: JMechanics.hh:39
container_type::iterator iterator
Definition: JHashMap.hh:88