Jpp  17.3.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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"
20 #include "JDetector/JLocation.hh"
24 #include "JGeometry3D/JVector3D.hh"
25 #include "JGeometry3D/JVersor3D.hh"
27 #include "JPhysics/JConstants.hh"
28 #include "JMath/JConstants.hh"
29 #include "JMath/JMathToolkit.hh"
30 #include "JTools/JRange.hh"
31 #include "JIO/JFileStreamIO.hh"
32 #include "Jeep/JeepToolkit.hh"
33 #include "JLang/JException.hh"
34 #include "JLang/gzstream.h"
35 #include "JLang/JManip.hh"
36 #include "JLang/Jpp.hh"
37 #include "JLang/JType.hh"
38 
39 
40 /**
41  * \author mdejong
42  */
43 
44 namespace JDETECTOR {}
45 namespace JPP { using namespace JDETECTOR; }
46 
47 /**
48  * Auxiliary classes and methods for detector calibration and simulation.
49  */
50 namespace JDETECTOR {
51 
52  using JLANG::JException;
55  using JLANG::JType;
56 
57 
58  /**
59  * File name extensions.
60  */
61  static const char* const GENDET_DETECTOR_FILE_FORMAT = "det"; //!< file format used by gendet
62  static const char* const BINARY_DETECTOR_FILE_FORMAT[] = { "dat", "datx" }; //!< JIO binary file format
63  static const char* const KM3NET_DETECTOR_FILE_FORMAT = "detx"; //!< %KM3NeT standard ASCII format
64  static const char* const ZIPPED_DETECTOR_FILE_FORMAT = "gz"; //!< zipped %KM3NeT standard ASCII format
65  static const char* const GDML_DETECTOR_FILE_FORMAT = "gdml"; //!< KM3Sim input format
66 
67  static const char* const GDML_SCHEMA = getenv("GDML_SCHEMA_DIR"); //!< directory necessary for correct GDML header output
68  static const char* const CAN_MARGIN_M = getenv("CAN_MARGIN_M"); //!< extension of the detector size to comply with the can definition
69  static const char* const G4GDML_DEFAULT_SCHEMALOCATION = "http://service-spi.web.cern.ch/service-spi/app/releases/GDML/schema/gdml.xsd";
70 
71 
72  /**
73  * Get maximal distance between modules in detector.
74  * The option can be used to include base modules, if any.
75  *
76  * \param detector detector
77  * \param option option
78  * \return maximal distance [m]
79  */
80  inline double getMaximalDistance(const JDetector& detector, const bool option = false)
81  {
82  using namespace JPP;
83 
84  double dmax = 0.0;
85 
86  for (JDetector::const_iterator i1 = detector.begin(); i1 != detector.end(); ++i1) {
87  for (JDetector::const_iterator i2 = detector.begin(); i2 != i1; ++i2) {
88 
89  if (option || (i1->getFloor() != 0 && i2->getFloor() != 0)) {
90 
91  const double ds = getDistance(i1->getPosition(), i2->getPosition());
92 
93  if (ds > dmax) {
94  dmax = ds;
95  }
96  }
97  }
98  }
99 
100  return dmax;
101  }
102 
103 
104  /**
105  * Get rotation over X axis in Geant4 coordinate system
106  *
107  * \param dir direction
108  * \return X-rotation [deg]
109  */
110  inline double GetXrotationG4(const JVersor3D dir)
111  {
112  using namespace JPP;
113 
114  const double phi = atan2(dir.getDY(), dir.getDZ())*(180.0/PI);
115 
116  if (phi < 0.0){
117  return phi + 360.0;
118  }
119  else{
120  return phi;
121  }
122  }
123 
124 
125  /**
126  * Get rotation over Y axis in Geant4 coordinate system
127  *
128  * \param dir direction
129  * \return Y-rotation [deg]
130  */
131  inline double GetYrotationG4(const JVersor3D dir)
132  {
133  using namespace JPP;
134 
135  return asin(-dir.getDX())*(180.0/PI);
136  }
137 
138 
139  inline void read_gdml(std::istream&, JDetector&)
140  {
141  THROW(JException, "Operation not defined.");
142  }
143 
144 
145  /**
146  * Writes KM3Sim GDML input file from detector
147  *
148  * \param out output stream
149  * \param detector detector
150  */
151  inline void write_gdml(std::ostream& out, const JDetector& detector)
152  {
153  using namespace std;
154  using namespace JPP;
155 
156  const double DEFAULT_CAN_MARGIN_M = 350.0; // default can margin [m]
157  const double DEFAULT_WORLD_MARGIN_M = 50.0; // default world margin [m]
158 
159  const JCylinder3D cylinder(detector.begin(), detector.end());
160 
161  double can_margin;
162 
163  if (CAN_MARGIN_M) {
164  can_margin = atof(CAN_MARGIN_M);
165  } else {
166  cerr << "CAN_MARGIN_M not defined! Setting standard CAN_MARGIN_M " << DEFAULT_CAN_MARGIN_M << " m." << endl;
167  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)
168  }
169 
170  const double WorldCylinderHeight = 2*(cylinder.getZmax() - cylinder.getZmin() + can_margin + DEFAULT_WORLD_MARGIN_M);
171  const double WorldRadius = cylinder.getLength() + cylinder.getRadius() + can_margin + DEFAULT_WORLD_MARGIN_M;
172 
173  const double Crust_Z_size = WorldCylinderHeight/2;
174  const double Crust_Z_position = -WorldCylinderHeight/4;
175 
176  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=\"";
177  if (!GDML_SCHEMA) {
178  cerr << "GDML_SCHEMA_DIR NOT DEFINED! Setting default path." << endl;
179  out << G4GDML_DEFAULT_SCHEMALOCATION << "\">\n\n\n";
180  } else {
181  out << GDML_SCHEMA << "gdml.xsd\">\n\n\n";
182  }
183  out << "<!-- Jpp version: "<< getGITVersion() << " -->\n";
184  out << "<define>\n";
185  out << "<rotation name=\"identity\"/>\n<position name=\"zero\"/>\n";
186 
187  int PMTs_NO = 0;
188 
189  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
190 
191  if(module->empty()) continue;
192 
193  const JVector3D center = module->getCenter();
194 
195  out << "<position name=\"PosString" << module->getString() << "_Dom" << module->getID() << "\" unit=\"m\" x=\"" << module->getX() << "\" y=\"" << module->getY() << "\" z=\"" << module->getZ() << "\"/>\n";
196 
197  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
198 
199  const JVector3D vec = static_cast<JVector3D>(*pmt).sub(center);
200  out << "<position name=\"CathodPosition" << pmt->getID() << "_" << module->getID() << "\" unit=\"m\" x=\"" << vec.getX() << "\" y=\"" << vec.getY() << "\" z=\"" << vec.getZ() << "\"/>\n";
201  out << "<rotation name=\"CathodRotation" << pmt->getID() << "_" << module->getID() << "\" unit=\"deg\" x=\"" << GetXrotationG4(*pmt) << "\" y=\"" << GetYrotationG4(*pmt) << "\" z=\"0.000000\"/>\n";
202  out << "<constant name=\"CathodID_" << PMTs_NO << "\" value=\"" << pmt->getID() << "\"/>\n";
203  PMTs_NO++;
204  }
205 
206  }
207 
208  out << "<position name=\"OMShift\" unit=\"m\" x=\"0\" y=\"0\" z=\"0.0392\"/>\n";
209  out << "\n\n\n";
210  out << "<!-- end of DU position definitions -->\n<position name=\"CrustPosition\" unit=\"m\" x=\"0\" y=\"0\" z=\"" << Crust_Z_position << "\"/>\n";
211 
212  out << "</define>\n";
213  out << "<materials>\n";
214  out << "</materials>\n";
215  out << "<solids>\n";
216  out << " <box name=\"CrustBox\" lunit=\"m\" x=\"2200\" y=\"2200\" z=\"" << Crust_Z_size << "\"/>\n";
217  out << " <box name=\"StoreyBox\" lunit=\"m\" x=\"0.6\" y=\"0.6\" z=\"0.6\"/>\n";
218  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";
219  out << " <tube name=\"CathodTube\" lunit=\"cm\" aunit=\"deg\" rmin=\"0.0\" rmax=\"4.7462\" z=\"0.5\" startphi=\"0.0\" deltaphi=\"360.0\"/>\n";
220  out << " <tube name=\"WorldTube\" lunit=\"m\" aunit=\"deg\" rmin=\"0.0\" rmax=\"" << WorldRadius << "\" z=\"" << WorldCylinderHeight << "\" startphi=\"0.0\" deltaphi=\"360.0\"/>\n";
221  out << "</solids>\n\n\n";
222 
223  out << "<structure>\n";
224  out << " <volume name=\"CathodVolume\">\n";
225  out << " <materialref ref=\"Cathod\"/>\n";
226  out << " <solidref ref=\"CathodTube\"/>\n";
227  out << " </volume>\n";
228 
229  out << "<!-- OMVolume(s) construction -->\n";
230 
231  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
232 
233  if(module->empty()) continue;
234  out << " <volume name=\"OMVolume" << module->getID() << "\">\n";
235  out << " <materialref ref=\"Water\"/>\n";
236  out << " <solidref ref=\"OMSphere\"/>\n";
237 
238  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
239  out << " <physvol>\n";
240  out << " <volumeref ref=\"CathodVolume\"/>\n";
241  out << " <positionref ref=\"CathodPosition" << pmt->getID() << "_" << module->getID() << "\"/>\n";
242  out << " <rotationref ref=\"CathodRotation" << pmt->getID() << "_" << module->getID() << "\"/>\n";
243  out << " </physvol>\n";
244  }
245 
246  out << " </volume>\n";
247  }
248 
249  out << "<!-- StoreyVolume(s) construction -->\n";
250 
251  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
252  if(module->empty()) continue;
253  out << " <volume name=\"StoreyVolume" << module->getID() << "\">\n";
254  out << " <materialref ref=\"Water\"/>\n";
255  out << " <solidref ref=\"StoreyBox\"/>\n";
256  out << " <physvol>\n";
257  out << " <volumeref ref=\"OMVolume" << module->getID() << "\"/>\n";
258  out << " <positionref ref=\"OMShift\"/>\n";
259  out << " <rotationref ref=\"identity\"/>\n";
260  out << " </physvol>\n";
261  out << " </volume>\n";
262  }
263 
264  out << "<!-- Crust Volume construction-->\n";
265  out << "<volume name=\"CrustVolume\">\n";
266  out << " <materialref ref=\"Crust\"/>\n";
267  out << " <solidref ref=\"CrustBox\"/>\n";
268  out << "</volume>\n";
269 
270  out << "<!-- World Volume construction-->\n";
271  out << "<volume name=\"WorldVolume\">\n";
272  out << " <materialref ref=\"Water\"/>\n";
273  out << " <solidref ref=\"WorldTube\"/>\n";
274 
275  out << " <physvol>\n";
276  out << " <volumeref ref=\"CrustVolume\"/>\n";
277  out << " <positionref ref=\"CrustPosition\"/>\n";
278  out << " <rotationref ref=\"identity\"/>\n";
279  out << " </physvol>\n";
280 
281  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
282  if(module->empty()) continue;
283  out << " <physvol>\n";
284  out << " <volumeref ref=\"StoreyVolume" << module->getID() << "\"/>\n";
285  out << " <positionref ref=\"PosString" << module->getString() << "_Dom" << module->getID() << "\"/>\n";
286  out << " <rotationref ref=\"identity\"/>\n";
287  out << " </physvol>\n";
288  }
289 
290  out << "</volume>\n";
291 
292  out << "</structure>\n";
293  out << "<setup name=\"Default\" version=\"1.0\">\n";
294  out << "<world ref=\"WorldVolume\"/>\n";
295  out << "</setup>\n";
296  out << "</gdml>\n";
297  }
298 
299 
300  /**
301  * Get maximal time between optical modules in detector following causality.
302  *
303  * \param detector detector
304  * \return maximal time [ns]
305  */
306  inline double getMaximalTime(const JDetector& detector)
307  {
308  using namespace JPP;
309 
311  }
312 
313 
314  /**
315  * Get maximal time between optical modules in detector following causality.
316  * The road width corresponds to the maximal distance traveled by the light.
317  *
318  * \param detector detector
319  * \param roadWidth_m road width [m]
320  * \return maximal time [ns]
321  */
322  inline double getMaximalTime(const JDetector& detector, const double roadWidth_m)
323  {
324  using namespace JPP;
325 
326  const double Dmax_m = getMaximalDistance(detector);
327 
328  return (sqrt((Dmax_m + roadWidth_m*getSinThetaC()) *
329  (Dmax_m - roadWidth_m*getSinThetaC())) +
330  roadWidth_m * getSinThetaC() * getTanThetaC()) * getInverseSpeedOfLight();
331  }
332 
333 
334  /**
335  * Get de-calibrated time range.
336  *
337  * The de-calibrated time range is corrected for minimal and maximal time offset of PMTs in given module.
338  *
339  * \param timeRange time range [ns]
340  * \param module module
341  * \return time range [ns]
342  */
343  inline JTimeRange getTimeRange(const JTimeRange& timeRange, const JModule& module)
344  {
345  if (!module.empty()) {
346 
348 
349  for (JModule::const_iterator pmt = module.begin(); pmt != module.end(); ++pmt) {
350 
351  const JCalibration& calibration = pmt->getCalibration();
352 
353  time_range.include(putTime(timeRange.getLowerLimit(), calibration));
354  time_range.include(putTime(timeRange.getUpperLimit(), calibration));
355  }
356 
357  return time_range;
358 
359  } else {
360 
361  return timeRange;
362  }
363  }
364 
365 
366  /**
367  * Get number of PMTs.
368  *
369  * \param module module
370  * \return number of PMTs
371  */
372  inline int getNumberOfPMTs(const JModule& module)
373  {
374  return module.size();
375  }
376 
377 
378  /**
379  * Get number of PMTs.
380  *
381  * \param detector detector
382  * \return number of PMTs
383  */
384  inline int getNumberOfPMTs(const JDetector& detector)
385  {
386  int number_of_pmts = 0;
387 
388  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
389  number_of_pmts += getNumberOfPMTs(*module);
390  }
391 
392  return number_of_pmts;
393  }
394 
395 
396  /**
397  * Get list of strings identifiers.
398  *
399  * \param detector detector
400  * \return list of string identifiers
401  */
403  {
404  std::set<int> buffer;
405 
406  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
407  buffer.insert(module->getString());
408  }
409 
410  return buffer;
411  }
412 
413 
414  /**
415  * Get number of floors.
416  *
417  * \param detector detector
418  * \return number of floors
419  */
421  {
422  std::set<int> buffer;
423 
424  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
425  buffer.insert(module->getFloor());
426  }
427 
428  return buffer.size();
429  }
430 
431 
432  /**
433  * Type definition for range of floors.
434  */
436 
437 
438  /**
439  * Get range of floors.
440  *
441  * \param detector detector
442  * \return range of floors
443  */
444  inline floor_range getRangeOfFloors(const JDetector& detector)
445  {
446  floor_range buffer = floor_range::DEFAULT_RANGE;
447 
448  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
449  buffer.include(module->getFloor());
450  }
451 
452  return buffer;
453  }
454 
455 
456  /**
457  * Get number of modules.
458  *
459  * \param detector detector
460  * \return number of modules
461  */
463  {
464  std::set<int> buffer;
465 
466  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
467  buffer.insert(module->getID());
468  }
469 
470  return buffer.size();
471  }
472 
473 
474  /**
475  * Load detector from input file.
476  *
477  * Supported file formats:
478  * - ASCII file, extension JDETECTOR::GENDET_DETECTOR_FILE_FORMAT, gendet format;
479  * - binary file, extension JDETECTOR::BINARY_DETECTOR_FILE_FORMAT, Jpp internal format;
480  * - ASCII file, extension JDETECTOR::KM3NET_DETECTOR_FILE_FORMAT, %KM3NeT standard format;
481  * - gzipped file, extension JDETECTOR::ZIPPED_DETECTOR_FILE_FORMAT, %KM3NeT standard format;
482  *
483  * \param file_name file name
484  * \param detector detector
485  */
486  inline void load(const std::string& file_name, JDetector& detector)
487  {
488  using namespace std;
489  using namespace JPP;
490 
492 
493  JMonteCarloDetector buffer(true);
494 
495  ifstream in(file_name.c_str());
496 
497  if (!in) {
498  THROW(JFileOpenException, "File not opened for reading: " << file_name);
499  }
500 
501  in >> buffer;
502 
503  in.close();
504 
505  detector.swap(buffer);
506 
507  } else if (getFilenameExtension(file_name) == BINARY_DETECTOR_FILE_FORMAT[0] ||
509 
510  JFileStreamReader in(file_name.c_str());
511 
512  if (!in) {
513  THROW(JFileOpenException, "File not opened for reading: " << file_name);
514  }
515 
516  try {
517 
518  detector.read(in);
519 
520  } catch(const exception& error) {
521 
522  // recovery procedure for old format of comments
523 
525 
526  in.clear();
527  in.rewind();
528 
529  detector.read(in);
530 
532  }
533 
534  in.close();
535 
536  } else if (getFilenameExtension(file_name) == KM3NET_DETECTOR_FILE_FORMAT) {
537 
538  ifstream in(file_name.c_str());
539 
540  if (!in) {
541  THROW(JFileOpenException, "File not opened for reading: " << file_name);
542  }
543 
544  in >> detector;
545 
546  in.close();
547 
548  } else if (getFilenameExtension(file_name) == ZIPPED_DETECTOR_FILE_FORMAT) {
549 
550  igzstream in(file_name.c_str());
551 
552  if (!in) {
553  THROW(JFileOpenException, "File not opened for reading: " << file_name);
554  }
555 
556  in >> detector;
557 
558  in.close();
559 
560  } else {
561 
562  THROW(JProtocolException, "Protocol not defined: " << file_name);
563  }
564  }
565 
566 
567  /**
568  * Store detector to output file.
569  *
570  * Supported file formats:
571  * - binary file, extension JDETECTOR::BINARY_DETECTOR_FILE_FORMAT, Jpp internal format;
572  * - ASCII file, extension JDETECTOR::KM3NET_DETECTOR_FILE_FORMAT, %KM3NeT standard format;
573  * - gzipped file, extension JDETECTOR::ZIPPED_DETECTOR_FILE_FORMAT, %KM3NeT standard format;
574  * - gdml file, extension JDETECTOR::GDML_DETECTOR_FILE_FORMAT, KM3Sim input format;
575  *
576  * \param file_name file name
577  * \param detector detector
578  */
579  inline void store(const std::string& file_name, const JDetector& detector)
580  {
581  using namespace std;
582  using namespace JPP;
583 
584  if (getFilenameExtension(file_name) == BINARY_DETECTOR_FILE_FORMAT[1]) {
585 
586  JFileStreamWriter out(file_name.c_str());
587 
588  if (!out) {
589  THROW(JFileOpenException, "File not opened for writing: " << file_name);
590  }
591 
592  detector.write(out);
593 
594  out.close();
595 
596  } else if (getFilenameExtension(file_name) == KM3NET_DETECTOR_FILE_FORMAT) {
597 
598  std::ofstream out(file_name.c_str());
599 
600  if (!out) {
601  THROW(JFileOpenException, "File not opened for writing: " << file_name);
602  }
603 
604  out << detector;
605 
606  out.close();
607 
608  } else if (getFilenameExtension(file_name) == ZIPPED_DETECTOR_FILE_FORMAT) {
609 
610  ogzstream out(file_name.c_str());
611 
612  if (!out) {
613  THROW(JFileOpenException, "File not opened for writing: " << file_name);
614  }
615 
616  out << detector;
617 
618  out.close();
619 
620  } else if (getFilenameExtension(file_name) == GDML_DETECTOR_FILE_FORMAT) {
621 
622  std::ofstream out(file_name.c_str());
623 
624  if (!out) {
625  THROW(JFileOpenException, "File not opened for writing: " << file_name);
626  }
627 
628  write_gdml(out, detector);
629 
630  out.close();
631 
632  } else {
633 
634  THROW(JProtocolException, "Protocol not defined: " << file_name);
635  }
636  }
637 
638 
639  /**
640  * Get module according module address map.
641  *
642  * \param memo module address map
643  * \param id module identifier
644  * \param location module location
645  * \return module
646  */
647  inline const JModule& getModule(const JModuleAddressMap& memo,
648  const int id = -1,
649  const JLocation& location = JLocation())
650  {
651  static JModule module;
652 
653 
654  module.setID(id);
655 
656  module.setLocation(location);
657 
658  module.resize(memo.size());
659 
660  if (memo.has( 0)) { module[memo[ 0].tdc] = JPMT( 1, JAxis3D(JVector3D(+0.000, +0.000, -0.200), JVersor3D(+0.000, +0.000, -1.000))); }
661 
662  if (memo.has( 1)) { module[memo[ 1].tdc] = JPMT( 2, JAxis3D(JVector3D(+0.000, +0.105, -0.170), JVersor3D(+0.000, +0.527, -0.850))); }
663  if (memo.has( 2)) { module[memo[ 2].tdc] = JPMT( 3, JAxis3D(JVector3D(+0.091, +0.053, -0.170), JVersor3D(+0.456, +0.263, -0.850))); }
664  if (memo.has( 3)) { module[memo[ 3].tdc] = JPMT( 4, JAxis3D(JVector3D(+0.091, -0.053, -0.170), JVersor3D(+0.456, -0.263, -0.850))); }
665  if (memo.has( 4)) { module[memo[ 4].tdc] = JPMT( 5, JAxis3D(JVector3D(+0.000, -0.105, -0.170), JVersor3D(+0.000, -0.527, -0.850))); }
666  if (memo.has( 5)) { module[memo[ 5].tdc] = JPMT( 6, JAxis3D(JVector3D(-0.091, -0.053, -0.170), JVersor3D(-0.456, -0.263, -0.850))); }
667  if (memo.has( 6)) { module[memo[ 6].tdc] = JPMT( 7, JAxis3D(JVector3D(-0.091, +0.053, -0.170), JVersor3D(-0.456, +0.263, -0.850))); }
668 
669  if (memo.has( 7)) { module[memo[ 7].tdc] = JPMT( 8, JAxis3D(JVector3D(+0.083, +0.144, -0.111), JVersor3D(+0.416, +0.720, -0.555))); }
670  if (memo.has( 8)) { module[memo[ 8].tdc] = JPMT( 9, JAxis3D(JVector3D(+0.166, +0.000, -0.111), JVersor3D(+0.832, +0.000, -0.555))); }
671  if (memo.has( 9)) { module[memo[ 9].tdc] = JPMT(10, JAxis3D(JVector3D(+0.083, -0.144, -0.111), JVersor3D(+0.416, -0.720, -0.555))); }
672  if (memo.has(10)) { module[memo[10].tdc] = JPMT(11, JAxis3D(JVector3D(-0.083, -0.144, -0.111), JVersor3D(-0.416, -0.720, -0.555))); }
673  if (memo.has(11)) { module[memo[11].tdc] = JPMT(12, JAxis3D(JVector3D(-0.166, +0.000, -0.111), JVersor3D(-0.832, +0.000, -0.555))); }
674  if (memo.has(12)) { module[memo[12].tdc] = JPMT(13, JAxis3D(JVector3D(-0.083, +0.144, -0.111), JVersor3D(-0.416, +0.720, -0.555))); }
675 
676  if (memo.has(13)) { module[memo[13].tdc] = JPMT(14, JAxis3D(JVector3D(+0.000, +0.191, -0.059), JVersor3D(+0.000, +0.955, -0.295))); }
677  if (memo.has(14)) { module[memo[14].tdc] = JPMT(15, JAxis3D(JVector3D(+0.165, +0.096, -0.059), JVersor3D(+0.827, +0.478, -0.295))); }
678  if (memo.has(15)) { module[memo[15].tdc] = JPMT(16, JAxis3D(JVector3D(+0.165, -0.096, -0.059), JVersor3D(+0.827, -0.478, -0.295))); }
679  if (memo.has(16)) { module[memo[16].tdc] = JPMT(17, JAxis3D(JVector3D(+0.000, -0.191, -0.059), JVersor3D(+0.000, -0.955, -0.295))); }
680  if (memo.has(17)) { module[memo[17].tdc] = JPMT(18, JAxis3D(JVector3D(-0.165, -0.096, -0.059), JVersor3D(-0.827, -0.478, -0.295))); }
681  if (memo.has(18)) { module[memo[18].tdc] = JPMT(19, JAxis3D(JVector3D(-0.165, +0.096, -0.059), JVersor3D(-0.827, +0.478, -0.295))); }
682 
683  if (memo.has(19)) { module[memo[19].tdc] = JPMT(20, JAxis3D(JVector3D(+0.096, +0.165, +0.059), JVersor3D(+0.478, +0.827, +0.295))); }
684  if (memo.has(20)) { module[memo[20].tdc] = JPMT(21, JAxis3D(JVector3D(+0.191, +0.000, +0.059), JVersor3D(+0.955, +0.000, +0.295))); }
685  if (memo.has(21)) { module[memo[21].tdc] = JPMT(22, JAxis3D(JVector3D(+0.096, -0.165, +0.059), JVersor3D(+0.478, -0.827, +0.295))); }
686  if (memo.has(22)) { module[memo[22].tdc] = JPMT(23, JAxis3D(JVector3D(-0.096, -0.165, +0.059), JVersor3D(-0.478, -0.827, +0.295))); }
687  if (memo.has(23)) { module[memo[23].tdc] = JPMT(24, JAxis3D(JVector3D(-0.191, +0.000, +0.059), JVersor3D(-0.955, +0.000, +0.295))); }
688  if (memo.has(24)) { module[memo[24].tdc] = JPMT(25, JAxis3D(JVector3D(-0.096, +0.165, +0.059), JVersor3D(-0.478, +0.827, +0.295))); }
689 
690  if (memo.has(25)) { module[memo[25].tdc] = JPMT(26, JAxis3D(JVector3D(+0.000, +0.166, +0.111), JVersor3D(+0.000, +0.832, +0.555))); }
691  if (memo.has(26)) { module[memo[26].tdc] = JPMT(27, JAxis3D(JVector3D(+0.144, +0.083, +0.111), JVersor3D(+0.720, +0.416, +0.555))); }
692  if (memo.has(27)) { module[memo[27].tdc] = JPMT(28, JAxis3D(JVector3D(+0.144, -0.083, +0.111), JVersor3D(+0.720, -0.416, +0.555))); }
693  if (memo.has(28)) { module[memo[28].tdc] = JPMT(29, JAxis3D(JVector3D(+0.000, -0.166, +0.111), JVersor3D(+0.000, -0.832, +0.555))); }
694  if (memo.has(29)) { module[memo[29].tdc] = JPMT(30, JAxis3D(JVector3D(-0.144, -0.083, +0.111), JVersor3D(-0.720, -0.416, +0.555))); }
695  if (memo.has(30)) { module[memo[30].tdc] = JPMT(31, JAxis3D(JVector3D(-0.144, +0.083, +0.111), JVersor3D(-0.720, +0.416, +0.555))); }
696 
697  module.compile();
698 
699  return module;
700  }
701 
702 
703  /**
704  * Get module according given detector type.
705  *
706  * \param type detector type
707  * \param id module identifier
708  * \param location module location
709  * \return module
710  */
711  template<class JDetector_t>
712  inline const JModule& getModule(const JType<JDetector_t> type,
713  const int id,
714  const JLocation& location = JLocation())
715  {
716  return getModule(getModuleAddressMap<JDetector_t>(id), id, location);
717  }
718 
719 
720  /**
721  * Get module according Antares detector type.
722  *
723  * \param type Antares detector type
724  * \param id module identifier
725  * \param location module location
726  * \return module
727  */
728  inline const JModule& getModule(const JType<JAntares_t> type,
729  const int id,
730  const JLocation& location = JLocation())
731  {
732  static JModule module;
733 
734  module.setID(id);
735 
736  module.setLocation(location);
737 
738  if (module.empty()) {
739 
740  module.resize(3);
741 
742  const double R = 0.5; // [m]
743 
744  const double st = sin(0.75*PI);
745  const double ct = cos(0.75*PI);
746 
747  for (int i = 0; i != 3; ++i) {
748 
749  const double phi = (2.0*PI*i) / 3.0;
750  const double cp = cos(phi);
751  const double sp = sin(phi);
752 
753  module[i] = JPMT(i, JAxis3D(JVector3D(R*st*cp, R*st*sp, R*ct), JVersor3D(st*cp, st*sp, ct)));
754  }
755  }
756 
757  return module;
758  }
759 
760 
761  /**
762  * Get module according detector template.
763  *
764  * \param id module identifier
765  * \param location module location
766  * \return module
767  */
768  template<class JDetector_t>
769  inline const JModule& getModule(const int id,
770  const JLocation& location = JLocation())
771  {
772  return getModule(JType<JDetector_t>(), id, location);
773  }
774 
775 
776  /**
777  * Auxiliary class to get rotation matrix between two optical modules.
778  */
779  struct JRotation :
780  public JRotation3D
781  {
782 
783  static const size_t NUMBER_OF_DIMENSIONS = 3; //!< Number of dimensions
784 
785 
786  /**
787  * Get rotation matrix to go from first to second module.
788  *
789  * \param first first module
790  * \param second second module
791  * \return rotation matrix
792  */
793  const JRotation3D& operator()(const JModule& first, const JModule& second)
794  {
795  this->setIdentity();
796 
797  if (first.size() == second.size()) {
798 
799  const size_t N = first.size();
800 
801  if (N >= NUMBER_OF_DIMENSIONS) {
802 
803  in .resize(N);
804  out.resize(N);
805 
806  for (size_t i = 0; i != N; ++i) {
807  in [i] = first .getPMT(i).getDirection();
808  out[i] = second.getPMT(i).getDirection();
809  }
810 
811  for (size_t i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
812  if (!orthonormalise(i)) {
813  THROW(JException, "Failure to orthonormalise direction " << i);
814  }
815  }
816 
817  this->a00 = out[0].getX() * in[0].getX() + out[1].getX() * in[1].getX() + out[2].getX() * in[2].getX();
818  this->a01 = out[0].getX() * in[0].getY() + out[1].getX() * in[1].getY() + out[2].getX() * in[2].getY();
819  this->a02 = out[0].getX() * in[0].getZ() + out[1].getX() * in[1].getZ() + out[2].getX() * in[2].getZ();
820 
821  this->a10 = out[0].getY() * in[0].getX() + out[1].getY() * in[1].getX() + out[2].getY() * in[2].getX();
822  this->a11 = out[0].getY() * in[0].getY() + out[1].getY() * in[1].getY() + out[2].getY() * in[2].getY();
823  this->a12 = out[0].getY() * in[0].getZ() + out[1].getY() * in[1].getZ() + out[2].getY() * in[2].getZ();
824 
825  this->a20 = out[0].getZ() * in[0].getX() + out[1].getZ() * in[1].getX() + out[2].getZ() * in[2].getX();
826  this->a21 = out[0].getZ() * in[0].getY() + out[1].getZ() * in[1].getY() + out[2].getZ() * in[2].getY();
827  this->a22 = out[0].getZ() * in[0].getZ() + out[1].getZ() * in[1].getZ() + out[2].getZ() * in[2].getZ();
828 
829  } else {
830 
831  THROW(JException, "Module " << first.getID() << " size " << N << " < " << NUMBER_OF_DIMENSIONS);
832  }
833 
834  } else {
835 
836  THROW(JException, "Module " << first.getID() << " size " << first.size() << " != " << second.size());
837  }
838 
839  return *this;
840  }
841 
842  private:
843  /**
844  * Put normalised primary direction at specified index and orthoganilise following directions.\n
845  * This procedure follows Gram-Schmidt process.
846  *
847  * \param index index
848  * \param precision precision
849  * \return true if primary direction exists; else false
850  */
851  bool orthonormalise(const size_t index, const double precision = std::numeric_limits<double>::epsilon())
852  {
853  using namespace std;
854 
855  size_t pos = index;
856 
857  for (size_t i = index + 1; i != in.size(); ++i) {
858  if (in[i].getLengthSquared() > in[pos].getLengthSquared()) {
859  pos = i;
860  }
861  }
862 
863  const double u = in[pos].getLength();
864 
865  if (u > precision) {
866 
867  in [pos] /= u;
868  out[pos] /= u;
869 
870  if (pos != index) {
871  swap(in [pos], in [index]);
872  swap(out[pos], out[index]);
873  }
874 
875  for (size_t i = index + 1; i != in.size(); ++i) {
876 
877  const double dot = in[index].getDot(in[i]);
878 
879  in [i] -= dot * in [index];
880  out[i] -= dot * out[index];
881  }
882 
883  return true;
884 
885  } else {
886 
887  return false;
888  }
889  }
890 
891 
894  };
895 
896 
897  /**
898  * Function object to get rotation matrix to go from first to second module.
899  */
901 
902 
903  /**
904  * Get position to go from first to second module.
905  *
906  * \param first first module
907  * \param second second module
908  * \return position
909  */
910  inline JPosition3D getPosition(const JModule& first, const JModule& second)
911  {
912  return second.getPosition() - first.getPosition();
913  }
914 
915 
916  /**
917  * Get calibration to go from first to second calibration.
918  *
919  * \param first first calibration
920  * \param second second calibration
921  * \return calibration
922  */
924  {
925  return JCalibration(second.getT0() - first.getT0());
926  }
927 }
928 
929 #endif
Exception for opening of file.
Definition: JException.hh:342
static const JRange< double, std::less< double > > DEFAULT_RANGE
Default range.
Definition: JRange.hh:555
double GetYrotationG4(const JVersor3D dir)
Get rotation over Y axis in Geant4 coordinate system.
General exception.
Definition: JException.hh:23
Exceptions.
T getLowerLimit() const
Get lower limit.
Definition: JRange.hh:202
Auxiliary methods for geometrical methods.
Data structure for a composite optical module.
Definition: JModule.hh:68
range_type & include(argument_type x)
Include given value to range.
Definition: JRange.hh:397
void setLocation(const JLocation &location)
Set location.
Definition: JLocation.hh:91
const JCalibration & getCalibration() const
Get calibration.
double getIndexOfRefraction()
Get average index of refraction of water corresponding to group velocity.
Detector data structure.
Definition: JDetector.hh:89
const JDirection3D & getDirection() const
Get direction.
Rotation matrix.
Definition: JRotation3D.hh:111
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:696
std::set< int > getStringIDs(const JDetector &detector)
Get list of strings identifiers.
JTOOLS::JRange< int > floor_range
Type definition for range of floors.
then JShowerPostfit f $INPUT_FILE o $OUTPUT_FILE N
int getNumberOfPMTs(const JModule &module)
Get number of PMTs.
static void setStartOfComment(const start_of_comment_type option)
Set option for short start of comment in binary I/O.
Definition: JDetector.hh:573
static JRotation getRotation
Function object to get rotation matrix to go from first to second module.
JCalibration getCalibration(const JCalibration &first, const JCalibration &second)
Get calibration to go from first to second calibration.
Jpp environment information.
static const char *const BINARY_DETECTOR_FILE_FORMAT[]
JIO binary file format.
Data structure for time calibration.
Auxiliary class for a type holder.
Definition: JType.hh:19
void read_gdml(std::istream &, JDetector &)
static const char *const GDML_SCHEMA
directory necessary for correct GDML header output
static const JModuleCounter getNumberOfModules
Function object to count unique modules.
static const char *const CAN_MARGIN_M
extension of the detector size to comply with the can definition
floor_range getRangeOfFloors(const JDetector &detector)
Get range of floors.
Data structure for detector geometry and calibration.
then echo The file $DIR KM3NeT_00000001_00000000 root already please rename or remove it first
double getDistance(const JFirst_t &first, const JSecond_t &second)
Get distance between objects.
JTimeRange getTimeRange(const Evt &event)
Get time range (i.e. time between earliest and latest hit) of Monte Carlo event.
Axis object.
Definition: JAxis3D.hh:38
then usage set_variable ACOUSTICS_WORKDIR $WORKDIR set_variable FORMULA sin([0]+2 *$PI *([1]+[2]*x)*x)" set_variable DY 1.0e-8 mkdir $WORKDIR for DETECTOR in $DETECTORS[*]
Monte Carlo detector (i.e. so-called .det file).
Lookup table for PMT addresses in optical module.
double putTime(const T &t1, const JCalibration &cal)
Get de-calibrated time.
JVector3D & sub(const JVector3D &vector)
Subtract vector.
Definition: JVector3D.hh:158
Detector specific mapping between logical positions and readout channels of PMTs in optical modules...
Cylinder object.
Definition: JCylinder3D.hh:39
virtual JReader & read(JReader &in) override
Read from input.
Definition: JDetector.hh:455
Detector file.
Definition: JHead.hh:226
std::string getGITVersion(const std::string &tag)
Get GIT version for given GIT tag.
Data structure for vector in three dimensions.
Definition: JVector3D.hh:34
virtual JWriter & write(JWriter &out) const override
Write to output.
Definition: JDetector.hh:507
Logical location of module.
Definition: JLocation.hh:37
Mathematical constants.
double getDY() const
Get y direction.
Definition: JVersor3D.hh:106
Auxiliary class to get rotation matrix between two optical modules.
double getDX() const
Get x direction.
Definition: JVersor3D.hh:95
Auxiliary methods for handling file names, type names and environment.
int getID() const
Get identifier.
Definition: JObjectID.hh:50
Data structure for PMT geometry, calibration and status.
Definition: JPMT.hh:43
void compile()
Compile module data.
Definition: JModule.hh:282
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
JPosition3D getPosition(const Vec &pos)
Get position.
T getUpperLimit() const
Get upper limit.
Definition: JRange.hh:213
int getNumberOfFloors(const JDetector &detector)
Get number of floors.
then awk string
Physics constants.
static const double PI
Mathematical constants.
double getY() const
Get y position.
Definition: JVector3D.hh:104
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:130
double getMaximalTime(const double R_Hz)
Get maximal time for given rate.
const JPMT & getPMT(const int index) const
Get PMT.
Definition: JModule.hh:173
Logical location of module.
Protocol exception.
Definition: JException.hh:648
const JRotation3D & operator()(const JModule &first, const JModule &second)
Get rotation matrix to go from first to second module.
std::vector< JVector3D > in
I/O manipulators.
static const char *const KM3NET_DETECTOR_FILE_FORMAT
KM3NeT standard ASCII format
static const char *const ZIPPED_DETECTOR_FILE_FORMAT
zipped KM3NeT standard ASCII format
then JCookie sh JDataQuality D $DETECTOR_ID R
Definition: JDataQuality.sh:41
static const char *const GENDET_DETECTOR_FILE_FORMAT
File name extensions.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Auxiliary class to define a range between two values.
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
const double getInverseSpeedOfLight()
Get inverse speed of light.
void setID(const int id)
Set identifier.
Definition: JObjectID.hh:72
double getX() const
Get x position.
Definition: JVector3D.hh:94
double getMaximalDistance(const JDetector &detector, const bool option=false)
Get maximal distance between modules in detector.
static const char *const G4GDML_DEFAULT_SCHEMALOCATION
then cp
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
do set_variable DETECTOR_TXT $WORKDIR detector
Data structure for normalised vector in three dimensions.
Definition: JVersor3D.hh:26
double u[N+1]
Definition: JPolint.hh:776
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
Definition: JCanberra.sh:46
Binary buffered file output.
double GetXrotationG4(const JVersor3D dir)
Get rotation over X axis in Geant4 coordinate system.
static const char *const GDML_DETECTOR_FILE_FORMAT
KM3Sim input format.
double getZ() const
Get z position.
Definition: JVector3D.hh:115
Binary buffered file input.
bool has(const int index) const
Test whether index is available.
double getDZ() const
Get z direction.
Definition: JVersor3D.hh:117
std::vector< JVector3D > out
const double epsilon
Definition: JQuadrature.cc:21
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 JModule & getModule(const JDetector &detector, const JModuleLocation &location)
find module with a given string and floor number
double getSinThetaC()
Get average sine of Cherenkov angle of water corresponding to group velocity.
void write_gdml(std::ostream &out, const JDetector &detector)
Writes KM3Sim GDML input file from detector.
double getT0() const
Get time offset.