Jpp  19.1.0-rc.1
the software that should make you happy
JDetectorToolkit.hh
Go to the documentation of this file.
1 #ifndef __JDETECTOR__JDETECTORTOOLKIT__
2 #define __JDETECTOR__JDETECTORTOOLKIT__
3 
4 #include <limits>
5 #include <istream>
6 #include <ostream>
7 #include <iostream>
8 #include <fstream>
9 #include <cmath>
10 #include <vector>
11 #include <set>
12 #include <algorithm>
13 #include <cstdlib>
14 
15 #include "JDetector/JDetector.hh"
17 #include "JDetector/JTimeRange.hh"
18 #include "JDetector/JLocation.hh"
21 #include "JGeometry3D/JVector3D.hh"
22 #include "JGeometry3D/JVersor3D.hh"
25 #include "JPhysics/JConstants.hh"
26 #include "JMath/JConstants.hh"
27 #include "JMath/JMathToolkit.hh"
28 #include "JTools/JRange.hh"
29 #include "JIO/JFileStreamIO.hh"
30 #include "Jeep/JeepToolkit.hh"
31 #include "JLang/JException.hh"
32 #include "JLang/gzstream.h"
33 #include "JLang/JManip.hh"
34 #include "JLang/Jpp.hh"
35 
36 
37 /**
38  * \author mdejong
39  */
40 
41 namespace JDETECTOR {}
42 namespace JPP { using namespace JDETECTOR; }
43 
44 /**
45  * Auxiliary classes and methods for detector calibration and simulation.
46  */
47 namespace JDETECTOR {
48 
49  using JLANG::JException;
54 
55 
56  /**
57  * File name extensions.
58  */
59  static const char* const GENDET_DETECTOR_FILE_FORMAT = "det"; //!< file format used by gendet
60  static const char* const BINARY_DETECTOR_FILE_FORMAT[] = { "dat", "datx" }; //!< JIO binary file format
61  static const char* const KM3NET_DETECTOR_FILE_FORMAT = "detx"; //!< %KM3NeT standard ASCII format
62  static const char* const ZIPPED_DETECTOR_FILE_FORMAT = "gz"; //!< zipped %KM3NeT standard ASCII format
63  static const char* const GDML_DETECTOR_FILE_FORMAT = "gdml"; //!< KM3Sim input format
64 
65  static const char* const GDML_SCHEMA = getenv("GDML_SCHEMA_DIR"); //!< directory necessary for correct GDML header output
66  static const char* const CAN_MARGIN_M = getenv("CAN_MARGIN_M"); //!< extension of the detector size to comply with the can definition
67  static const char* const G4GDML_DEFAULT_SCHEMALOCATION = "http://service-spi.web.cern.ch/service-spi/app/releases/GDML/schema/gdml.xsd";
68 
69 
70  /**
71  * Get maximal distance between modules in detector.
72  * The option can be used to include base modules, if any.
73  *
74  * \param detector detector
75  * \param option option
76  * \return maximal distance [m]
77  */
78  inline double getMaximalDistance(const JDetector& detector, const bool option = false)
79  {
80  using namespace JPP;
81 
82  double dmax = 0.0;
83 
84  for (JDetector::const_iterator i1 = detector.begin(); i1 != detector.end(); ++i1) {
85  for (JDetector::const_iterator i2 = detector.begin(); i2 != i1; ++i2) {
86 
87  if (option || (i1->getFloor() != 0 && i2->getFloor() != 0)) {
88 
89  const double ds = getDistance(i1->getPosition(), i2->getPosition());
90 
91  if (ds > dmax) {
92  dmax = ds;
93  }
94  }
95  }
96  }
97 
98  return dmax;
99  }
100 
101 
102  /**
103  * Get maximal time between optical modules in detector following causality.
104  *
105  * \param detector detector
106  * \return maximal time [ns]
107  */
108  inline double getMaximalTime(const JDetector& detector)
109  {
110  using namespace JPP;
111 
113  }
114 
115 
116  /**
117  * Get maximal time between optical modules in detector following causality.
118  * The road width corresponds to the maximal distance traveled by the light.
119  *
120  * \param detector detector
121  * \param roadWidth_m road width [m]
122  * \return maximal time [ns]
123  */
124  inline double getMaximalTime(const JDetector& detector, const double roadWidth_m)
125  {
126  using namespace JPP;
127 
128  const double Dmax_m = getMaximalDistance(detector);
129 
130  return (sqrt((Dmax_m + roadWidth_m*getSinThetaC()) *
131  (Dmax_m - roadWidth_m*getSinThetaC())) +
132  roadWidth_m * getSinThetaC() * getTanThetaC()) * getInverseSpeedOfLight();
133  }
134 
135 
136  /**
137  * Get de-calibrated time range.
138  *
139  * The de-calibrated time range is corrected for minimal and maximal time offset of PMTs in given module.
140  *
141  * \param timeRange time range [ns]
142  * \param module module
143  * \return time range [ns]
144  */
145  inline JTimeRange getTimeRange(const JTimeRange& timeRange, const JModule& module)
146  {
147  if (!module.empty()) {
148 
150 
151  for (JModule::const_iterator pmt = module.begin(); pmt != module.end(); ++pmt) {
152 
153  const JCalibration& calibration = pmt->getCalibration();
154 
155  time_range.include(putTime(timeRange.getLowerLimit(), calibration));
156  time_range.include(putTime(timeRange.getUpperLimit(), calibration));
157  }
158 
159  return time_range;
160 
161  } else {
162 
163  return timeRange;
164  }
165  }
166 
167 
168  /**
169  * Get number of PMTs.
170  *
171  * \param module module
172  * \return number of PMTs
173  */
174  inline int getNumberOfPMTs(const JModule& module)
175  {
176  return module.size();
177  }
178 
179 
180  /**
181  * Get number of PMTs.
182  *
183  * \param detector detector
184  * \return number of PMTs
185  */
186  inline int getNumberOfPMTs(const JDetector& detector)
187  {
188  int number_of_pmts = 0;
189 
190  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
191  number_of_pmts += getNumberOfPMTs(*module);
192  }
193 
194  return number_of_pmts;
195  }
196 
197 
198  /**
199  * Get list of strings identifiers.
200  *
201  * \param detector detector
202  * \return list of string identifiers
203  */
204  inline std::set<int> getStringIDs(const JDetector& detector)
205  {
206  std::set<int> buffer;
207 
208  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
209  buffer.insert(module->getString());
210  }
211 
212  return buffer;
213  }
214 
215 
216  /**
217  * Get number of floors.
218  *
219  * \param detector detector
220  * \return number of floors
221  */
222  inline int getNumberOfFloors(const JDetector& detector)
223  {
224  std::set<int> buffer;
225 
226  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
227  buffer.insert(module->getFloor());
228  }
229 
230  return buffer.size();
231  }
232 
233 
234  /**
235  * Type definition for range of floors.
236  */
238 
239 
240  /**
241  * Get range of floors.
242  *
243  * \param detector detector
244  * \return range of floors
245  */
246  inline floor_range getRangeOfFloors(const JDetector& detector)
247  {
249 
250  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
251  buffer.include(module->getFloor());
252  }
253 
254  return buffer;
255  }
256 
257 
258  /**
259  * Get number of modules.
260  *
261  * \param detector detector
262  * \return number of modules
263  */
264  inline int getNumberOfModules(const JDetector& detector)
265  {
266  std::set<int> buffer;
267 
268  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
269  buffer.insert(module->getID());
270  }
271 
272  return buffer.size();
273  }
274 
275 
276  /**
277  * Get rotation over X axis in Geant4 coordinate system
278  *
279  * \param dir direction
280  * \return X-rotation [deg]
281  */
282  inline double GetXrotationG4(const JVersor3D& dir)
283  {
284  using namespace JPP;
285 
286  const double phi = atan2(dir.getDY(), dir.getDZ())*(180.0/PI);
287 
288  if (phi < 0.0){
289  return phi + 360.0;
290  }
291  else{
292  return phi;
293  }
294  }
295 
296 
297  /**
298  * Get rotation over Y axis in Geant4 coordinate system
299  *
300  * \param dir direction
301  * \return Y-rotation [deg]
302  */
303  inline double GetYrotationG4(const JVersor3D& dir)
304  {
305  using namespace JPP;
306 
307  return asin(-dir.getDX())*(180.0/PI);
308  }
309 
310 
311  inline void read_gdml(std::istream&, JDetector&)
312  {
313  THROW(JException, "Operation not defined.");
314  }
315 
316 
317  /**
318  * Writes KM3Sim GDML input file from detector
319  *
320  * \param out output stream
321  * \param detector detector
322  */
323  inline void write_gdml(std::ostream& out, const JDetector& detector)
324  {
325  using namespace std;
326  using namespace JPP;
327 
328  const double DEFAULT_CAN_MARGIN_M = 350.0; // default can margin [m]
329  const double DEFAULT_WORLD_MARGIN_M = 50.0; // default world margin [m]
330 
331  const JCylinder3D cylinder(detector.begin(), detector.end());
332 
333  double can_margin;
334 
335  if (CAN_MARGIN_M) {
336  can_margin = atof(CAN_MARGIN_M);
337  } else {
338  cerr << "CAN_MARGIN_M not defined! Setting standard CAN_MARGIN_M " << DEFAULT_CAN_MARGIN_M << " m." << endl;
339  can_margin = DEFAULT_CAN_MARGIN_M; //this is given in meters like all the dimensions in the GDML file (look at the unit field of every position and solid definition)
340  }
341 
342  const double WorldCylinderHeight = 2*(cylinder.getZmax() - cylinder.getZmin() + can_margin + DEFAULT_WORLD_MARGIN_M);
343  const double WorldRadius = cylinder.getLength() + cylinder.getRadius() + can_margin + DEFAULT_WORLD_MARGIN_M;
344 
345  const double Crust_Z_size = WorldCylinderHeight/2;
346  const double Crust_Z_position = -WorldCylinderHeight/4;
347 
348  out << "<?xml version=\"1.0\" encoding=\"UTF-8\" ?>\n<gdml xmlns:gdml=\"http://cern.ch/2001/Schemas/GDML\" xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\" xsi:noNamespaceSchemaLocation=\"";
349  if (!GDML_SCHEMA) {
350  cerr << "GDML_SCHEMA_DIR NOT DEFINED! Setting default path." << endl;
351  out << G4GDML_DEFAULT_SCHEMALOCATION << "\">\n\n\n";
352  } else {
353  out << GDML_SCHEMA << "gdml.xsd\">\n\n\n";
354  }
355  out << "<!-- Jpp version: "<< getGITVersion() << " -->\n";
356  out << "<define>\n";
357  out << "<rotation name=\"identity\"/>\n<position name=\"zero\"/>\n";
358 
359  int PMTs_NO = 0;
360 
361  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
362 
363  if(module->empty()) continue;
364 
365  const JVector3D center = module->getCenter();
366 
367  out << "<position name=\"PosString" << module->getString() << "_Dom" << module->getID() << "\" unit=\"m\" x=\"" << module->getX() << "\" y=\"" << module->getY() << "\" z=\"" << module->getZ() << "\"/>\n";
368 
369  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
370 
371  const JVector3D vec = static_cast<JVector3D>(*pmt).sub(center);
372  out << "<position name=\"CathodPosition" << pmt->getID() << "_" << module->getID() << "\" unit=\"m\" x=\"" << vec.getX() << "\" y=\"" << vec.getY() << "\" z=\"" << vec.getZ() << "\"/>\n";
373  out << "<rotation name=\"CathodRotation" << pmt->getID() << "_" << module->getID() << "\" unit=\"deg\" x=\"" << GetXrotationG4(*pmt) << "\" y=\"" << GetYrotationG4(*pmt) << "\" z=\"0.000000\"/>\n";
374  out << "<constant name=\"CathodID_" << PMTs_NO << "\" value=\"" << pmt->getID() << "\"/>\n";
375  PMTs_NO++;
376  }
377 
378  }
379 
380  out << "<position name=\"OMShift\" unit=\"m\" x=\"0\" y=\"0\" z=\"0.0392\"/>\n";
381  out << "\n\n\n";
382  out << "<!-- end of DU position definitions -->\n<position name=\"CrustPosition\" unit=\"m\" x=\"0\" y=\"0\" z=\"" << Crust_Z_position << "\"/>\n";
383 
384  out << "</define>\n";
385  out << "<materials>\n";
386  out << "</materials>\n";
387  out << "<solids>\n";
388  out << " <box name=\"CrustBox\" lunit=\"m\" x=\"2200\" y=\"2200\" z=\"" << Crust_Z_size << "\"/>\n";
389  out << " <box name=\"StoreyBox\" lunit=\"m\" x=\"0.6\" y=\"0.6\" z=\"0.6\"/>\n";
390  out << " <sphere name=\"OMSphere\" lunit=\"cm\" aunit=\"deg\" rmin=\"0.0\" rmax=\"21.6\" startphi=\"0.0\" deltaphi=\"360.0\" starttheta=\"0.0\" deltatheta=\"180.0\"/>\n";
391  out << " <tube name=\"CathodTube\" lunit=\"cm\" aunit=\"deg\" rmin=\"0.0\" rmax=\"4.7462\" z=\"0.5\" startphi=\"0.0\" deltaphi=\"360.0\"/>\n";
392  out << " <tube name=\"WorldTube\" lunit=\"m\" aunit=\"deg\" rmin=\"0.0\" rmax=\"" << WorldRadius << "\" z=\"" << WorldCylinderHeight << "\" startphi=\"0.0\" deltaphi=\"360.0\"/>\n";
393  out << "</solids>\n\n\n";
394 
395  out << "<structure>\n";
396  out << " <volume name=\"CathodVolume\">\n";
397  out << " <materialref ref=\"Cathod\"/>\n";
398  out << " <solidref ref=\"CathodTube\"/>\n";
399  out << " </volume>\n";
400 
401  out << "<!-- OMVolume(s) construction -->\n";
402 
403  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
404 
405  if(module->empty()) continue;
406  out << " <volume name=\"OMVolume" << module->getID() << "\">\n";
407  out << " <materialref ref=\"Water\"/>\n";
408  out << " <solidref ref=\"OMSphere\"/>\n";
409 
410  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
411  out << " <physvol>\n";
412  out << " <volumeref ref=\"CathodVolume\"/>\n";
413  out << " <positionref ref=\"CathodPosition" << pmt->getID() << "_" << module->getID() << "\"/>\n";
414  out << " <rotationref ref=\"CathodRotation" << pmt->getID() << "_" << module->getID() << "\"/>\n";
415  out << " </physvol>\n";
416  }
417 
418  out << " </volume>\n";
419  }
420 
421  out << "<!-- StoreyVolume(s) construction -->\n";
422 
423  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
424  if(module->empty()) continue;
425  out << " <volume name=\"StoreyVolume" << module->getID() << "\">\n";
426  out << " <materialref ref=\"Water\"/>\n";
427  out << " <solidref ref=\"StoreyBox\"/>\n";
428  out << " <physvol>\n";
429  out << " <volumeref ref=\"OMVolume" << module->getID() << "\"/>\n";
430  out << " <positionref ref=\"OMShift\"/>\n";
431  out << " <rotationref ref=\"identity\"/>\n";
432  out << " </physvol>\n";
433  out << " </volume>\n";
434  }
435 
436  out << "<!-- Crust Volume construction-->\n";
437  out << "<volume name=\"CrustVolume\">\n";
438  out << " <materialref ref=\"Crust\"/>\n";
439  out << " <solidref ref=\"CrustBox\"/>\n";
440  out << "</volume>\n";
441 
442  out << "<!-- World Volume construction-->\n";
443  out << "<volume name=\"WorldVolume\">\n";
444  out << " <materialref ref=\"Water\"/>\n";
445  out << " <solidref ref=\"WorldTube\"/>\n";
446 
447  out << " <physvol>\n";
448  out << " <volumeref ref=\"CrustVolume\"/>\n";
449  out << " <positionref ref=\"CrustPosition\"/>\n";
450  out << " <rotationref ref=\"identity\"/>\n";
451  out << " </physvol>\n";
452 
453  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
454  if(module->empty()) continue;
455  out << " <physvol>\n";
456  out << " <volumeref ref=\"StoreyVolume" << module->getID() << "\"/>\n";
457  out << " <positionref ref=\"PosString" << module->getString() << "_Dom" << module->getID() << "\"/>\n";
458  out << " <rotationref ref=\"identity\"/>\n";
459  out << " </physvol>\n";
460  }
461 
462  out << "</volume>\n";
463 
464  out << "</structure>\n";
465  out << "<setup name=\"Default\" version=\"1.0\">\n";
466  out << "<world ref=\"WorldVolume\"/>\n";
467  out << "</setup>\n";
468  out << "</gdml>\n";
469  }
470 
471 
472  /**
473  * Load detector from input file.
474  *
475  * Supported file formats:
476  * - ASCII file, extension JDETECTOR::GENDET_DETECTOR_FILE_FORMAT, gendet format;
477  * - binary file, extension JDETECTOR::BINARY_DETECTOR_FILE_FORMAT, Jpp internal format;
478  * - ASCII file, extension JDETECTOR::KM3NET_DETECTOR_FILE_FORMAT, %KM3NeT standard format;
479  * - gzipped file, extension JDETECTOR::ZIPPED_DETECTOR_FILE_FORMAT, %KM3NeT standard format;
480  *
481  * \param file_name file name
482  * \param detector detector
483  */
484  inline void load(const std::string& file_name, JDetector& detector)
485  {
486  using namespace std;
487  using namespace JPP;
488 
490 
491  JMonteCarloDetector buffer(true);
492 
493  ifstream in(file_name.c_str());
494 
495  if (!in) {
496  THROW(JFileOpenException, "File not opened for reading: " << file_name);
497  }
498 
499  in >> buffer;
500 
501  in.close();
502 
503  detector.swap(buffer);
504 
505  } else if (getFilenameExtension(file_name) == BINARY_DETECTOR_FILE_FORMAT[0] ||
507 
508  JFileStreamReader in(file_name.c_str());
509 
510  if (!in) {
511  THROW(JFileOpenException, "File not opened for reading: " << file_name);
512  }
513 
514  try {
515 
516  detector.read(in);
517 
518  } catch(const exception& error) {
519 
520  // recovery procedure for old format of comments
521 
523 
524  in.clear();
525  in.rewind();
526 
527  detector.read(in);
528 
530  }
531 
532  in.close();
533 
534  } else if (getFilenameExtension(file_name) == KM3NET_DETECTOR_FILE_FORMAT) {
535 
536  ifstream in(file_name.c_str());
537 
538  if (!in) {
539  THROW(JFileOpenException, "File not opened for reading: " << file_name);
540  }
541 
542  in >> detector;
543 
544  in.close();
545 
546  } else if (getFilenameExtension(file_name) == ZIPPED_DETECTOR_FILE_FORMAT) {
547 
548  igzstream in(file_name.c_str());
549 
550  if (!in) {
551  THROW(JFileOpenException, "File not opened for reading: " << file_name);
552  }
553 
554  in >> detector;
555 
556  in.close();
557 
558  } else {
559 
560  THROW(JProtocolException, "Protocol not defined: " << file_name);
561  }
562 
563  if (!detector.hasVersion()) {
564  THROW(JValueOutOfRange, "Invalid version <" << detector.getVersion() << ">");
565  }
566  }
567 
568 
569  /**
570  * Store detector to output file.
571  *
572  * Supported file formats:
573  * - binary file, extension JDETECTOR::BINARY_DETECTOR_FILE_FORMAT, Jpp internal format;
574  * - ASCII file, extension JDETECTOR::KM3NET_DETECTOR_FILE_FORMAT, %KM3NeT standard format;
575  * - gzipped file, extension JDETECTOR::ZIPPED_DETECTOR_FILE_FORMAT, %KM3NeT standard format;
576  * - gdml file, extension JDETECTOR::GDML_DETECTOR_FILE_FORMAT, KM3Sim input format;
577  *
578  * \param file_name file name
579  * \param detector detector
580  */
581  inline void store(const std::string& file_name, const JDetector& detector)
582  {
583  using namespace std;
584  using namespace JPP;
585 
586  if (getFilenameExtension(file_name) == BINARY_DETECTOR_FILE_FORMAT[1]) {
587 
588  JFileStreamWriter out(file_name.c_str());
589 
590  if (!out) {
591  THROW(JFileOpenException, "File not opened for writing: " << file_name);
592  }
593 
594  detector.write(out);
595 
596  out.close();
597 
598  } else if (getFilenameExtension(file_name) == KM3NET_DETECTOR_FILE_FORMAT) {
599 
600  std::ofstream out(file_name.c_str());
601 
602  if (!out) {
603  THROW(JFileOpenException, "File not opened for writing: " << file_name);
604  }
605 
606  out << detector;
607 
608  out.close();
609 
610  } else if (getFilenameExtension(file_name) == ZIPPED_DETECTOR_FILE_FORMAT) {
611 
612  ogzstream out(file_name.c_str());
613 
614  if (!out) {
615  THROW(JFileOpenException, "File not opened for writing: " << file_name);
616  }
617 
618  out << detector;
619 
620  out.close();
621 
622  } else if (getFilenameExtension(file_name) == GDML_DETECTOR_FILE_FORMAT) {
623 
624  std::ofstream out(file_name.c_str());
625 
626  if (!out) {
627  THROW(JFileOpenException, "File not opened for writing: " << file_name);
628  }
629 
630  write_gdml(out, detector);
631 
632  out.close();
633 
634  } else {
635 
636  THROW(JProtocolException, "Protocol not defined: " << file_name);
637  }
638  }
639 
640 
641  /**
642  * Auxiliary class to get rotation matrix between two optical modules.
643  */
644  struct JRotation :
645  public JRotation3D
646  {
647 
648  static const size_t NUMBER_OF_DIMENSIONS = 3; //!< Number of dimensions
649 
650 
651  /**
652  * Get rotation matrix to go from first to second module.
653  *
654  * \param first first module
655  * \param second second module
656  * \return rotation matrix
657  */
658  const JRotation3D& operator()(const JModule& first, const JModule& second)
659  {
660  this->setIdentity();
661 
662  if (first.size() == second.size()) {
663 
664  const size_t N = first.size();
665 
666  if (N >= NUMBER_OF_DIMENSIONS) {
667 
668  in .resize(N);
669  out.resize(N);
670 
671  for (size_t i = 0; i != N; ++i) {
672  in [i] = first .getPMT(i).getDirection();
673  out[i] = second.getPMT(i).getDirection();
674  }
675 
676  for (size_t i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
677  if (!orthonormalise(i)) {
678  THROW(JException, "Failure to orthonormalise direction " << i);
679  }
680  }
681 
682  this->a00 = out[0].getX() * in[0].getX() + out[1].getX() * in[1].getX() + out[2].getX() * in[2].getX();
683  this->a01 = out[0].getX() * in[0].getY() + out[1].getX() * in[1].getY() + out[2].getX() * in[2].getY();
684  this->a02 = out[0].getX() * in[0].getZ() + out[1].getX() * in[1].getZ() + out[2].getX() * in[2].getZ();
685 
686  this->a10 = out[0].getY() * in[0].getX() + out[1].getY() * in[1].getX() + out[2].getY() * in[2].getX();
687  this->a11 = out[0].getY() * in[0].getY() + out[1].getY() * in[1].getY() + out[2].getY() * in[2].getY();
688  this->a12 = out[0].getY() * in[0].getZ() + out[1].getY() * in[1].getZ() + out[2].getY() * in[2].getZ();
689 
690  this->a20 = out[0].getZ() * in[0].getX() + out[1].getZ() * in[1].getX() + out[2].getZ() * in[2].getX();
691  this->a21 = out[0].getZ() * in[0].getY() + out[1].getZ() * in[1].getY() + out[2].getZ() * in[2].getY();
692  this->a22 = out[0].getZ() * in[0].getZ() + out[1].getZ() * in[1].getZ() + out[2].getZ() * in[2].getZ();
693 
694  } else {
695 
696  THROW(JException, "Module " << first.getID() << " size " << N << " < " << NUMBER_OF_DIMENSIONS);
697  }
698 
699  } else {
700 
701  THROW(JException, "Module " << first.getID() << " size " << first.size() << " != " << second.size());
702  }
703 
704  return *this;
705  }
706 
707  private:
708  /**
709  * Put normalised primary direction at specified index and orthoganilise following directions.\n
710  * This procedure follows Gram-Schmidt process.
711  *
712  * \param index index
713  * \param precision precision
714  * \return true if primary direction exists; else false
715  */
716  bool orthonormalise(const size_t index, const double precision = std::numeric_limits<double>::epsilon())
717  {
718  using namespace std;
719 
720  size_t pos = index;
721 
722  for (size_t i = index + 1; i != in.size(); ++i) {
723  if (in[i].getLengthSquared() > in[pos].getLengthSquared()) {
724  pos = i;
725  }
726  }
727 
728  const double u = in[pos].getLength();
729 
730  if (u > precision) {
731 
732  in [pos] /= u;
733  out[pos] /= u;
734 
735  if (pos != index) {
736  swap(in [pos], in [index]);
737  swap(out[pos], out[index]);
738  }
739 
740  for (size_t i = index + 1; i != in.size(); ++i) {
741 
742  const double dot = in[index].getDot(in[i]);
743 
744  in [i] -= dot * in [index];
745  out[i] -= dot * out[index];
746  }
747 
748  return true;
749 
750  } else {
751 
752  return false;
753  }
754  }
755 
756 
759  };
760 
761 
762  /**
763  * Function object to get rotation matrix to go from first to second module.
764  */
766 
767 
768  /**
769  * Get position to go from first to second module.
770  *
771  * \param first first module
772  * \param second second module
773  * \return position
774  */
775  inline JPosition3D getPosition(const JModule& first, const JModule& second)
776  {
777  return second.getPosition() - first.getPosition();
778  }
779 
780 
781  /**
782  * Get calibration to go from first to second calibration.
783  *
784  * \param first first calibration
785  * \param second second calibration
786  * \return calibration
787  */
788  inline JCalibration getCalibration(const JCalibration& first, const JCalibration& second)
789  {
790  return JCalibration(second.getT0() - first.getT0());
791  }
792 }
793 
794 #endif
Data structure for detector geometry and calibration.
Exceptions.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:712
Logical location of module.
I/O manipulators.
Auxiliary methods for geometrical methods.
Mathematical constants.
Physics constants.
Auxiliary class to define a range between two values.
Auxiliary methods for handling file names, type names and environment.
Jpp environment information.
Data structure for time calibration.
double getT0() const
Get time offset.
Detector data structure.
Definition: JDetector.hh:96
static void setStartOfComment(const start_of_comment_type option)
Set option for short start of comment in binary I/O.
Definition: JDetector.hh:589
virtual JWriter & write(JWriter &out) const override
Write to output.
Definition: JDetector.hh:523
virtual JReader & read(JReader &in) override
Read from input.
Definition: JDetector.hh:473
Data structure for a composite optical module.
Definition: JModule.hh:75
const JPMT & getPMT(const int index) const
Get PMT.
Definition: JModule.hh:172
Monte Carlo detector (i.e. ".det" file).
double getRadius() const
Get radius.
Definition: JCircle2D.hh:144
double getLength() const
Get length.
Definition: JVector2D.hh:199
Cylinder object.
Definition: JCylinder3D.hh:41
double getZmin() const
Get minimal z position.
Definition: JCylinder3D.hh:111
double getZmax() const
Get maximal z position.
Definition: JCylinder3D.hh:122
const JDirection3D & getDirection() const
Get direction.
Data structure for position in three dimensions.
Definition: JPosition3D.hh:38
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
Rotation matrix.
Definition: JRotation3D.hh:114
Data structure for vector in three dimensions.
Definition: JVector3D.hh:36
double getY() const
Get y position.
Definition: JVector3D.hh:104
double getZ() const
Get z position.
Definition: JVector3D.hh:115
JVector3D & sub(const JVector3D &vector)
Subtract vector.
Definition: JVector3D.hh:158
double getX() const
Get x position.
Definition: JVector3D.hh:94
Data structure for normalised vector in three dimensions.
Definition: JVersor3D.hh:28
double getDY() const
Get y direction.
Definition: JVersor3D.hh:106
double getDX() const
Get x direction.
Definition: JVersor3D.hh:95
double getDZ() const
Get z direction.
Definition: JVersor3D.hh:117
Binary buffered file input.
virtual void clear() override
Clear status of reader.
Binary buffered file output.
void close()
Close file.
General exception.
Definition: JException.hh:24
Exception for opening of file.
Definition: JException.hh:360
int getID() const
Get identifier.
Definition: JObjectID.hh:50
Protocol exception.
Definition: JException.hh:666
Exception for accessing a value in a collection that is outside of its range.
Definition: JException.hh:180
range_type & include(argument_type x)
Include given value to range.
Definition: JRange.hh:397
T getLowerLimit() const
Get lower limit.
Definition: JRange.hh:202
static JRange< double, std::less< double > > DEFAULT_RANGE()
Default range.
Definition: JRange.hh:555
T getUpperLimit() const
Get upper limit.
Definition: JRange.hh:213
void close()
Definition: gzstream.h:185
const double epsilon
Definition: JQuadrature.cc:21
file Auxiliary data structures and methods for detector calibration.
Definition: JAnchor.hh:12
static JRotation getRotation
Function object to get rotation matrix to go from first to second module.
std::set< int > getStringIDs(const JDetector &detector)
Get list of strings identifiers.
int getNumberOfFloors(const JDetector &detector)
Get number of floors.
floor_range getRangeOfFloors(const JDetector &detector)
Get range of floors.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
static const char *const GDML_DETECTOR_FILE_FORMAT
KM3Sim input format.
JCalibration getCalibration(const JCalibration &first, const JCalibration &second)
Get calibration to go from first to second calibration.
void write_gdml(std::ostream &out, const JDetector &detector)
Writes KM3Sim GDML input file from detector.
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
static const char *const ZIPPED_DETECTOR_FILE_FORMAT
zipped KM3NeT standard ASCII format
static const char *const G4GDML_DEFAULT_SCHEMALOCATION
int getNumberOfModules(const JDetector &detector)
Get number of modules.
JPosition3D getPosition(const JModule &first, const JModule &second)
Get position to go from first to second module.
JTOOLS::JRange< int > floor_range
Type definition for range of floors.
int getNumberOfPMTs(const JDetector &detector)
Get number of PMTs.
static const char *const CAN_MARGIN_M
extension of the detector size to comply with the can definition
static const char *const KM3NET_DETECTOR_FILE_FORMAT
KM3NeT standard ASCII format
double getMaximalTime(const JDetector &detector, const double roadWidth_m)
Get maximal time between optical modules in detector following causality.
static const char *const GDML_SCHEMA
directory necessary for correct GDML header output
static const char *const GENDET_DETECTOR_FILE_FORMAT
File name extensions.
double GetYrotationG4(const JVersor3D &dir)
Get rotation over Y axis in Geant4 coordinate system.
double getMaximalDistance(const JDetector &detector, const bool option=false)
Get maximal distance between modules in detector.
double putTime(const T &t1, const JCalibration &cal)
Get de-calibrated time.
JTimeRange getTimeRange(const JTimeRange &timeRange, const JModule &module)
Get de-calibrated time range.
double GetXrotationG4(const JVersor3D &dir)
Get rotation over X axis in Geant4 coordinate system.
void read_gdml(std::istream &, JDetector &)
static const char *const BINARY_DETECTOR_FILE_FORMAT[]
JIO binary file format.
std::string getFilenameExtension(const std::string &file_name)
Get file name extension, i.e. part after last JEEP::FILENAME_SEPARATOR if any.
Definition: JeepToolkit.hh:109
std::string getGITVersion(const std::string &tag)
Get GIT version for given GIT tag.
static const double PI
Mathematical constants.
double getDistance(const JFirst_t &first, const JSecond_t &second)
Get distance between objects.
double getIndexOfRefraction()
Get average index of refraction of water corresponding to group velocity.
const double getInverseSpeedOfLight()
Get inverse speed of light.
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
double getSinThetaC()
Get average sine of Cherenkov angle of water corresponding to group velocity.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
double u[N+1]
Definition: JPolint.hh:865
Definition: JSTDTypes.hh:14
Calibration.
Definition: JHead.hh:330
bool hasVersion() const
Check validity of version.
Auxiliary class to get rotation matrix between two optical modules.
std::vector< JVector3D > out
std::vector< JVector3D > in
bool orthonormalise(const size_t index, const double precision=std::numeric_limits< double >::epsilon())
Put normalised primary direction at specified index and orthoganilise following directions.
const JRotation3D & operator()(const JModule &first, const JModule &second)
Get rotation matrix to go from first to second module.
const std::string & getVersion() const
Get version.