Jpp  master_rocky
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 list of modules identifiers.
218  *
219  * \param detector detector
220  * \param option option to include base modules
221  * \return list of module identifiers
222  */
223  inline std::set<int> getModuleIDs(const JDetector& detector, const bool option = false)
224  {
225  std::set<int> buffer;
226 
227  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
228  if (option || !module->empty()) {
229  buffer.insert(module->getID());
230  }
231  }
232 
233  return buffer;
234  }
235 
236 
237  /**
238  * Get number of floors.
239  *
240  * \param detector detector
241  * \return number of floors
242  */
243  inline int getNumberOfFloors(const JDetector& detector)
244  {
245  std::set<int> buffer;
246 
247  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
248  buffer.insert(module->getFloor());
249  }
250 
251  return buffer.size();
252  }
253 
254 
255  /**
256  * Type definition for range of floors.
257  */
259 
260 
261  /**
262  * Get range of floors.
263  *
264  * \param detector detector
265  * \return range of floors
266  */
267  inline floor_range getRangeOfFloors(const JDetector& detector)
268  {
270 
271  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
272  buffer.include(module->getFloor());
273  }
274 
275  return buffer;
276  }
277 
278 
279  /**
280  * Get number of modules.
281  *
282  * \param detector detector
283  * \return number of modules
284  */
285  inline int getNumberOfModules(const JDetector& detector)
286  {
287  std::set<int> buffer;
288 
289  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
290  buffer.insert(module->getID());
291  }
292 
293  return buffer.size();
294  }
295 
296 
297  /**
298  * Get rotation over X axis in Geant4 coordinate system
299  *
300  * \param dir direction
301  * \return X-rotation [deg]
302  */
303  inline double GetXrotationG4(const JVersor3D& dir)
304  {
305  using namespace JPP;
306 
307  const double phi = atan2(dir.getDY(), dir.getDZ())*(180.0/PI);
308 
309  if (phi < 0.0){
310  return phi + 360.0;
311  }
312  else{
313  return phi;
314  }
315  }
316 
317 
318  /**
319  * Get rotation over Y axis in Geant4 coordinate system
320  *
321  * \param dir direction
322  * \return Y-rotation [deg]
323  */
324  inline double GetYrotationG4(const JVersor3D& dir)
325  {
326  using namespace JPP;
327 
328  return asin(-dir.getDX())*(180.0/PI);
329  }
330 
331 
332  inline void read_gdml(std::istream&, JDetector&)
333  {
334  THROW(JException, "Operation not defined.");
335  }
336 
337 
338  /**
339  * Writes KM3Sim GDML input file from detector
340  *
341  * \param out output stream
342  * \param detector detector
343  */
344  inline void write_gdml(std::ostream& out, const JDetector& detector)
345  {
346  using namespace std;
347  using namespace JPP;
348 
349  const double DEFAULT_CAN_MARGIN_M = 350.0; // default can margin [m]
350  const double DEFAULT_WORLD_MARGIN_M = 50.0; // default world margin [m]
351 
352  const JCylinder3D cylinder(detector.begin(), detector.end());
353 
354  double can_margin;
355 
356  if (CAN_MARGIN_M) {
357  can_margin = atof(CAN_MARGIN_M);
358  } else {
359  cerr << "CAN_MARGIN_M not defined! Setting standard CAN_MARGIN_M " << DEFAULT_CAN_MARGIN_M << " m." << endl;
360  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)
361  }
362 
363  const double WorldCylinderHeight = 2*(cylinder.getZmax() - cylinder.getZmin() + can_margin + DEFAULT_WORLD_MARGIN_M);
364  const double WorldRadius = cylinder.getLength() + cylinder.getRadius() + can_margin + DEFAULT_WORLD_MARGIN_M;
365 
366  const double Crust_Z_size = WorldCylinderHeight/2;
367  const double Crust_Z_position = -WorldCylinderHeight/4;
368 
369  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=\"";
370  if (!GDML_SCHEMA) {
371  cerr << "GDML_SCHEMA_DIR NOT DEFINED! Setting default path." << endl;
372  out << G4GDML_DEFAULT_SCHEMALOCATION << "\">\n\n\n";
373  } else {
374  out << GDML_SCHEMA << "gdml.xsd\">\n\n\n";
375  }
376  out << "<!-- Jpp version: "<< getGITVersion() << " -->\n";
377  out << "<define>\n";
378  out << "<rotation name=\"identity\"/>\n<position name=\"zero\"/>\n";
379 
380  int PMTs_NO = 0;
381 
382  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
383 
384  if(module->empty()) continue;
385 
386  const JVector3D center = module->getCenter();
387 
388  out << "<position name=\"PosString" << module->getString() << "_Dom" << module->getID() << "\" unit=\"m\" x=\"" << module->getX() << "\" y=\"" << module->getY() << "\" z=\"" << module->getZ() << "\"/>\n";
389 
390  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
391 
392  const JVector3D vec = static_cast<JVector3D>(*pmt).sub(center);
393  out << "<position name=\"CathodPosition" << pmt->getID() << "_" << module->getID() << "\" unit=\"m\" x=\"" << vec.getX() << "\" y=\"" << vec.getY() << "\" z=\"" << vec.getZ() << "\"/>\n";
394  out << "<rotation name=\"CathodRotation" << pmt->getID() << "_" << module->getID() << "\" unit=\"deg\" x=\"" << GetXrotationG4(*pmt) << "\" y=\"" << GetYrotationG4(*pmt) << "\" z=\"0.000000\"/>\n";
395  out << "<constant name=\"CathodID_" << PMTs_NO << "\" value=\"" << pmt->getID() << "\"/>\n";
396  PMTs_NO++;
397  }
398 
399  }
400 
401  out << "<position name=\"OMShift\" unit=\"m\" x=\"0\" y=\"0\" z=\"0.0392\"/>\n";
402  out << "\n\n\n";
403  out << "<!-- end of DU position definitions -->\n<position name=\"CrustPosition\" unit=\"m\" x=\"0\" y=\"0\" z=\"" << Crust_Z_position << "\"/>\n";
404 
405  out << "</define>\n";
406  out << "<materials>\n";
407  out << "</materials>\n";
408  out << "<solids>\n";
409  out << " <box name=\"CrustBox\" lunit=\"m\" x=\"2200\" y=\"2200\" z=\"" << Crust_Z_size << "\"/>\n";
410  out << " <box name=\"StoreyBox\" lunit=\"m\" x=\"0.6\" y=\"0.6\" z=\"0.6\"/>\n";
411  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";
412  out << " <tube name=\"CathodTube\" lunit=\"cm\" aunit=\"deg\" rmin=\"0.0\" rmax=\"4.7462\" z=\"0.5\" startphi=\"0.0\" deltaphi=\"360.0\"/>\n";
413  out << " <tube name=\"WorldTube\" lunit=\"m\" aunit=\"deg\" rmin=\"0.0\" rmax=\"" << WorldRadius << "\" z=\"" << WorldCylinderHeight << "\" startphi=\"0.0\" deltaphi=\"360.0\"/>\n";
414  out << "</solids>\n\n\n";
415 
416  out << "<structure>\n";
417  out << " <volume name=\"CathodVolume\">\n";
418  out << " <materialref ref=\"Cathod\"/>\n";
419  out << " <solidref ref=\"CathodTube\"/>\n";
420  out << " </volume>\n";
421 
422  out << "<!-- OMVolume(s) construction -->\n";
423 
424  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
425 
426  if(module->empty()) continue;
427  out << " <volume name=\"OMVolume" << module->getID() << "\">\n";
428  out << " <materialref ref=\"Water\"/>\n";
429  out << " <solidref ref=\"OMSphere\"/>\n";
430 
431  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
432  out << " <physvol>\n";
433  out << " <volumeref ref=\"CathodVolume\"/>\n";
434  out << " <positionref ref=\"CathodPosition" << pmt->getID() << "_" << module->getID() << "\"/>\n";
435  out << " <rotationref ref=\"CathodRotation" << pmt->getID() << "_" << module->getID() << "\"/>\n";
436  out << " </physvol>\n";
437  }
438 
439  out << " </volume>\n";
440  }
441 
442  out << "<!-- StoreyVolume(s) construction -->\n";
443 
444  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
445  if(module->empty()) continue;
446  out << " <volume name=\"StoreyVolume" << module->getID() << "\">\n";
447  out << " <materialref ref=\"Water\"/>\n";
448  out << " <solidref ref=\"StoreyBox\"/>\n";
449  out << " <physvol>\n";
450  out << " <volumeref ref=\"OMVolume" << module->getID() << "\"/>\n";
451  out << " <positionref ref=\"OMShift\"/>\n";
452  out << " <rotationref ref=\"identity\"/>\n";
453  out << " </physvol>\n";
454  out << " </volume>\n";
455  }
456 
457  out << "<!-- Crust Volume construction-->\n";
458  out << "<volume name=\"CrustVolume\">\n";
459  out << " <materialref ref=\"Crust\"/>\n";
460  out << " <solidref ref=\"CrustBox\"/>\n";
461  out << "</volume>\n";
462 
463  out << "<!-- World Volume construction-->\n";
464  out << "<volume name=\"WorldVolume\">\n";
465  out << " <materialref ref=\"Water\"/>\n";
466  out << " <solidref ref=\"WorldTube\"/>\n";
467 
468  out << " <physvol>\n";
469  out << " <volumeref ref=\"CrustVolume\"/>\n";
470  out << " <positionref ref=\"CrustPosition\"/>\n";
471  out << " <rotationref ref=\"identity\"/>\n";
472  out << " </physvol>\n";
473 
474  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
475  if(module->empty()) continue;
476  out << " <physvol>\n";
477  out << " <volumeref ref=\"StoreyVolume" << module->getID() << "\"/>\n";
478  out << " <positionref ref=\"PosString" << module->getString() << "_Dom" << module->getID() << "\"/>\n";
479  out << " <rotationref ref=\"identity\"/>\n";
480  out << " </physvol>\n";
481  }
482 
483  out << "</volume>\n";
484 
485  out << "</structure>\n";
486  out << "<setup name=\"Default\" version=\"1.0\">\n";
487  out << "<world ref=\"WorldVolume\"/>\n";
488  out << "</setup>\n";
489  out << "</gdml>\n";
490  }
491 
492 
493  /**
494  * Load detector from input file.
495  *
496  * Supported file formats:
497  * - ASCII file, extension JDETECTOR::GENDET_DETECTOR_FILE_FORMAT, gendet format;
498  * - binary file, extension JDETECTOR::BINARY_DETECTOR_FILE_FORMAT, Jpp internal format;
499  * - ASCII file, extension JDETECTOR::KM3NET_DETECTOR_FILE_FORMAT, %KM3NeT standard format;
500  * - gzipped file, extension JDETECTOR::ZIPPED_DETECTOR_FILE_FORMAT, %KM3NeT standard format;
501  *
502  * \param file_name file name
503  * \param detector detector
504  */
505  inline void load(const std::string& file_name, JDetector& detector)
506  {
507  using namespace std;
508  using namespace JPP;
509 
511 
512  JMonteCarloDetector buffer(true);
513 
514  ifstream in(file_name.c_str());
515 
516  if (!in) {
517  THROW(JFileOpenException, "File not opened for reading: " << file_name);
518  }
519 
520  in >> buffer;
521 
522  in.close();
523 
524  detector.swap(buffer);
525 
526  } else if (getFilenameExtension(file_name) == BINARY_DETECTOR_FILE_FORMAT[0] ||
528 
529  JFileStreamReader in(file_name.c_str());
530 
531  if (!in) {
532  THROW(JFileOpenException, "File not opened for reading: " << file_name);
533  }
534 
535  try {
536 
537  detector.read(in);
538 
539  } catch(const exception& error) {
540 
541  // recovery procedure for old format of comments
542 
544 
545  in.clear();
546  in.rewind();
547 
548  detector.read(in);
549 
551  }
552 
553  in.close();
554 
555  } else if (getFilenameExtension(file_name) == KM3NET_DETECTOR_FILE_FORMAT) {
556 
557  ifstream in(file_name.c_str());
558 
559  if (!in) {
560  THROW(JFileOpenException, "File not opened for reading: " << file_name);
561  }
562 
563  in >> detector;
564 
565  in.close();
566 
567  } else if (getFilenameExtension(file_name) == ZIPPED_DETECTOR_FILE_FORMAT) {
568 
569  igzstream in(file_name.c_str());
570 
571  if (!in) {
572  THROW(JFileOpenException, "File not opened for reading: " << file_name);
573  }
574 
575  in >> detector;
576 
577  in.close();
578 
579  } else {
580 
581  THROW(JProtocolException, "Protocol not defined: " << file_name);
582  }
583 
584  if (!detector.hasVersion()) {
585  THROW(JValueOutOfRange, "Invalid version <" << detector.getVersion() << ">");
586  }
587  }
588 
589 
590  /**
591  * Store detector to output file.
592  *
593  * Supported file formats:
594  * - binary file, extension JDETECTOR::BINARY_DETECTOR_FILE_FORMAT, Jpp internal format;
595  * - ASCII file, extension JDETECTOR::KM3NET_DETECTOR_FILE_FORMAT, %KM3NeT standard format;
596  * - gzipped file, extension JDETECTOR::ZIPPED_DETECTOR_FILE_FORMAT, %KM3NeT standard format;
597  * - gdml file, extension JDETECTOR::GDML_DETECTOR_FILE_FORMAT, KM3Sim input format;
598  *
599  * \param file_name file name
600  * \param detector detector
601  */
602  inline void store(const std::string& file_name, const JDetector& detector)
603  {
604  using namespace std;
605  using namespace JPP;
606 
607  if (getFilenameExtension(file_name) == BINARY_DETECTOR_FILE_FORMAT[1]) {
608 
609  JFileStreamWriter out(file_name.c_str());
610 
611  if (!out) {
612  THROW(JFileOpenException, "File not opened for writing: " << file_name);
613  }
614 
615  detector.write(out);
616 
617  out.close();
618 
619  } else if (getFilenameExtension(file_name) == KM3NET_DETECTOR_FILE_FORMAT) {
620 
621  std::ofstream out(file_name.c_str());
622 
623  if (!out) {
624  THROW(JFileOpenException, "File not opened for writing: " << file_name);
625  }
626 
627  out << detector;
628 
629  out.close();
630 
631  } else if (getFilenameExtension(file_name) == ZIPPED_DETECTOR_FILE_FORMAT) {
632 
633  ogzstream out(file_name.c_str());
634 
635  if (!out) {
636  THROW(JFileOpenException, "File not opened for writing: " << file_name);
637  }
638 
639  out << detector;
640 
641  out.close();
642 
643  } else if (getFilenameExtension(file_name) == GDML_DETECTOR_FILE_FORMAT) {
644 
645  std::ofstream out(file_name.c_str());
646 
647  if (!out) {
648  THROW(JFileOpenException, "File not opened for writing: " << file_name);
649  }
650 
651  write_gdml(out, detector);
652 
653  out.close();
654 
655  } else {
656 
657  THROW(JProtocolException, "Protocol not defined: " << file_name);
658  }
659  }
660 
661 
662  /**
663  * Auxiliary class to get rotation matrix between two optical modules.
664  */
665  struct JRotation :
666  public JRotation3D
667  {
668 
669  static const size_t NUMBER_OF_DIMENSIONS = 3; //!< Number of dimensions
670 
671 
672  /**
673  * Get rotation matrix to go from first to second module.
674  *
675  * \param first first module
676  * \param second second module
677  * \return rotation matrix
678  */
679  const JRotation3D& operator()(const JModule& first, const JModule& second)
680  {
681  this->setIdentity();
682 
683  if (first.size() == second.size()) {
684 
685  const size_t N = first.size();
686 
687  if (N >= NUMBER_OF_DIMENSIONS) {
688 
689  in .resize(N);
690  out.resize(N);
691 
692  for (size_t i = 0; i != N; ++i) {
693  in [i] = first .getPMT(i).getDirection();
694  out[i] = second.getPMT(i).getDirection();
695  }
696 
697  for (size_t i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
698  if (!orthonormalise(i)) {
699  THROW(JException, "Failure to orthonormalise direction " << i);
700  }
701  }
702 
703  this->a00 = out[0].getX() * in[0].getX() + out[1].getX() * in[1].getX() + out[2].getX() * in[2].getX();
704  this->a01 = out[0].getX() * in[0].getY() + out[1].getX() * in[1].getY() + out[2].getX() * in[2].getY();
705  this->a02 = out[0].getX() * in[0].getZ() + out[1].getX() * in[1].getZ() + out[2].getX() * in[2].getZ();
706 
707  this->a10 = out[0].getY() * in[0].getX() + out[1].getY() * in[1].getX() + out[2].getY() * in[2].getX();
708  this->a11 = out[0].getY() * in[0].getY() + out[1].getY() * in[1].getY() + out[2].getY() * in[2].getY();
709  this->a12 = out[0].getY() * in[0].getZ() + out[1].getY() * in[1].getZ() + out[2].getY() * in[2].getZ();
710 
711  this->a20 = out[0].getZ() * in[0].getX() + out[1].getZ() * in[1].getX() + out[2].getZ() * in[2].getX();
712  this->a21 = out[0].getZ() * in[0].getY() + out[1].getZ() * in[1].getY() + out[2].getZ() * in[2].getY();
713  this->a22 = out[0].getZ() * in[0].getZ() + out[1].getZ() * in[1].getZ() + out[2].getZ() * in[2].getZ();
714 
715  } else {
716 
717  THROW(JException, "Module " << first.getID() << " size " << N << " < " << NUMBER_OF_DIMENSIONS);
718  }
719 
720  } else {
721 
722  THROW(JException, "Module " << first.getID() << " size " << first.size() << " != " << second.size());
723  }
724 
725  return *this;
726  }
727 
728  private:
729  /**
730  * Put normalised primary direction at specified index and orthoganilise following directions.\n
731  * This procedure follows Gram-Schmidt process.
732  *
733  * \param index index
734  * \param precision precision
735  * \return true if primary direction exists; else false
736  */
737  bool orthonormalise(const size_t index, const double precision = std::numeric_limits<double>::epsilon())
738  {
739  using namespace std;
740 
741  size_t pos = index;
742 
743  for (size_t i = index + 1; i != in.size(); ++i) {
744  if (in[i].getLengthSquared() > in[pos].getLengthSquared()) {
745  pos = i;
746  }
747  }
748 
749  const double u = in[pos].getLength();
750 
751  if (u > precision) {
752 
753  in [pos] /= u;
754  out[pos] /= u;
755 
756  if (pos != index) {
757  swap(in [pos], in [index]);
758  swap(out[pos], out[index]);
759  }
760 
761  for (size_t i = index + 1; i != in.size(); ++i) {
762 
763  const double dot = in[index].getDot(in[i]);
764 
765  in [i] -= dot * in [index];
766  out[i] -= dot * out[index];
767  }
768 
769  return true;
770 
771  } else {
772 
773  return false;
774  }
775  }
776 
777 
780  };
781 
782 
783  /**
784  * Function object to get rotation matrix to go from first to second module.
785  */
787 
788 
789  /**
790  * Get position to go from first to second module.
791  *
792  * \param first first module
793  * \param second second module
794  * \return position
795  */
796  inline JPosition3D getPosition(const JModule& first, const JModule& second)
797  {
798  return second.getPosition() - first.getPosition();
799  }
800 
801 
802  /**
803  * Get calibration to go from first to second calibration.
804  *
805  * \param first first calibration
806  * \param second second calibration
807  * \return calibration
808  */
809  inline JCalibration getCalibration(const JCalibration& first, const JCalibration& second)
810  {
811  return JCalibration(second.getT0() - first.getT0());
812  }
813 }
814 
815 #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.
std::set< int > getModuleIDs(const JDetector &detector, const bool option=false)
Get list of modules identifiers.
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.