Jpp
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 <fstream>
8 #include <cmath>
9 #include <vector>
10 #include <set>
11 #include <algorithm>
12 
13 #include "JDetector/JDetector.hh"
14 #include "JDetector/JMonteCarloDetector.hh"
15 #include "JDetector/JTimeRange.hh"
16 #include "JDetector/JModuleAddressMap.hh"
17 #include "JDetector/JModuleLocation.hh"
18 #include "JDetector/JModuleGeometry.hh"
19 #include "JGeometry3D/JVector3D.hh"
20 #include "JGeometry3D/JVersor3D.hh"
21 #include "JGeometry3D/JCylinder3D.hh"
22 #include "JTools/JConstants.hh"
23 #include "JIO/JFileStreamIO.hh"
24 #include "JLang/JString.hh"
25 #include "JLang/JException.hh"
26 #include "JLang/gzstream.h"
27 #include "JFit/JPoint3DEstimator.hh"
28 #include "JFit/JPoint3D.hh"
29 
30 
31 /**
32  * \author mdejong
33  */
34 
35 namespace JDETECTOR {}
36 namespace JPP { using namespace JDETECTOR; }
37 
38 /**
39  * Auxiliary classes and methods for detector calibration and simulation.
40  */
41 namespace JDETECTOR {
42 
43  using JGEOMETRY3D::JVector3D;
44  using JGEOMETRY3D::JVersor3D;
45  using JTOOLS::getInverseSpeedOfLight;
46  using JTOOLS::getIndexOfRefraction;
47  using JTOOLS::getTanThetaC;
48  using JTOOLS::getSinThetaC;
49  using JTOOLS::PI;
50  using JLANG::JString;
51  using JLANG::JException;
52  using JLANG::JFileOpenException;
53  using JLANG::JProtocolException;
54  using JGEOMETRY3D::JCylinder3D;
55  using JFIT::JPoint3D;
56 
57 
58 
59  /**
60  * File name extensions.
61  */
62  static const char* const GENDET_DETECTOR_FILE_FORMAT = ".det"; //!< file format used by gendet
63  static const char* const BINARY_DETECTOR_FILE_FORMAT = ".dat"; //!< JIO binary file format
64  static const char* const KM3NET_DETECTOR_FILE_FORMAT = ".detx"; //!< KM3NeT standard ASCII format
65  static const char* const ZIPPED_DETECTOR_FILE_FORMAT = ".gz"; //!< zipped KM3NeT standard ASCII format
66  static const char* const GDML_DETECTOR_FILE_FORMAT = ".xml"; //!< KM3Sim input format
67  static const char* const GDML_SCHEMA = getenv("GDML_SCHEMA_DIR"); //!< directory necessary for correct GDML header output
68 
69 
70  /**
71  * Get maximal distance between modules in detector.
72  *
73  * \param detector detector
74  * \return maximal distance [m]
75  */
76  inline double getMaximalDistance(const JDetector& detector)
77  {
78  double dmax = 0.0;
79 
80  for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
81  for (JDetector::const_iterator j = detector.begin(); j != i; ++j) {
82  if (i->getDistance(*j) > dmax) {
83  dmax = i->getDistance(*j);
84  }
85  }
86  }
87 
88  return dmax;
89  }
90 
91 
92  /**
93  * Get rotation over X axis in Geant4 coordinate system
94  *
95  * \param dir direction
96  * \return X-rotation [deg]
97  */
98  inline double GetXrotationG4(const JVersor3D dir)
99  {
100  const double phi = atan2(dir.getDY(), dir.getDZ())*(180.0/PI);
101 
102  if (phi < 0.0){
103  return phi + 360.0;
104  }
105  else{
106  return phi;
107  }
108  }
109 
110 
111  /**
112  * Get rotation over Y axis in Geant4 coordinate system
113  *
114  * \param dir direction
115  * \return Y-rotation [deg]
116  */
117  inline double GetYrotationG4(const JVersor3D dir)
118  {
119  return asin(-dir.getDX())*(180.0/PI);
120  }
121 
122 
123  inline void read_gdml(std::istream&, const JDetector&)
124  {}
125 
126 
127  /**
128  * Writes KM3Sim GDML input file from detector
129  *
130  * \param out output stream
131  * \param detector detector
132  */
133  inline void write_gdml(std::ostream& out, const JDetector& detector)
134  {
135 
136 
137  const JCylinder3D cylinder(detector.begin(), detector.end());
138 
139  // Detector height is bascially just a distance between the highest and the lowest DOM
140  //double Detector_Height = cylinder.getZmax() - cylinder.getZmin();
141  // Distance from the center of the lowest DOM to the seabed
142  //double Seabed_Distance = cylinder.getZmin();
143 
144  // To translate Z coordinate from detx format which corresponds to the distance from the seabed one has to perform a following calculation
145  // Dom_Z_GDML = Dom_Z_Detx - Distance to the seabed from the lowest DOM - Detector_height/2
146 
147  //double DistanceOfSeabedFromLowestOM = Seabed_Distance;
148  // For now, the size of the WorldBox is fixed to 2200 meters
149  const double WorldBoxHeight = 2200.;
150  // Below formula comes directly from the note by Apostolos
151  //double Crust_Z_size = WorldBoxHeight/2-Detector_Height/2-DistanceOfSeabedFromLowestOM; Formula for obsolete coordinate covention, where Z_seabed != 0
152  const double Crust_Z_size = WorldBoxHeight/2;
153  //double Crust_Z_position = -WorldBoxHeight/4-Detector_Height/4-DistanceOfSeabedFromLowestOM/2; Formula for obsolete coordinate covention, where Z_seabed != 0
154  const double Crust_Z_position= -WorldBoxHeight/4;
155 
156 
157  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=\"";
158  if(!GDML_SCHEMA){
159  std::cout<<"GDML_SCHEMA_DIR NOT DEFINED! Setting default path"<<std::endl;
160  out<<"/sps/km3net/users/jmanczak/km3sim/geant4.10.05/source/persistency/gdml/schema/gdml.xsd\">\n\n\n";
161  }
162  else{
163  out<<GDML_SCHEMA<<"gdml.xsd\">\n\n\n";
164  }
165  out<<"<define>\n";
166  out<<"<rotation name=\"identity\"/>\n<position name=\"zero\"/>\n";
167 
168  int PMTs_NO = 0;
169  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
170 
171  const JFIT::JEstimator<JPoint3D> center(module->begin(), module->end());
172 
173  //out<<"<position name=\"PosString"<<module->getString()<<"_Dom"<<module->getID()<<"\" unit=\"m\" x=\""<<module->getX()<<"\" y=\""<<module->getY()<<"\" z=\""<<module->getZ() - Seabed_Distance-Detector_Height/2<<"\"/>\n"; Formula for obsolete coordinate covention, where Z_seabed != 0
174  out<<"<position name=\"PosString"<<module->getString()<<"_Dom"<<module->getID()<<"\" unit=\"m\" x=\""<<module->getX()<<"\" y=\""<<module->getY()<<"\" z=\""<<module->getZ()<<"\"/>\n";
175 
176  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
177  //JVector3D vec;
178  //static_cast<JVector3D> (vec) = module->sub(*pmt);
179 
180  const JVector3D vec = static_cast<JVector3D>(*pmt).sub(center);
181  out<<"<position name=\"CathodPosition"<<pmt->getID()<<"_"<<module->getID()<<"\" unit=\"m\" x=\""<<vec.getX()<<"\" y=\""<<vec.getY()<<"\" z=\""<<vec.getZ()<<"\"/>\n";
182  out<<"<rotation name=\"CathodRotation"<<pmt->getID()<<"_"<<module->getID()<<"\" unit=\"deg\" x=\""<<GetXrotationG4(*pmt)<<"\" y=\""<<GetYrotationG4(*pmt)<<"\" z=\"0.000000\"/>\n";
183  out<<"<constant name=\"CathodID_"<<PMTs_NO<<"\" value=\""<<pmt->getID()<<"\"/>\n";
184  PMTs_NO++;
185  }
186 
187  }
188 
189  out<<"<position name=\"OMShift\" unit=\"m\" x=\"0\" y=\"0\" z=\"0.0392\"/>\n";
190  out<<"\n\n\n";
191  out<<"<!-- end of DU position definitions -->\n<position name=\"CrustPosition\" unit=\"m\" x=\"0\" y=\"0\" z=\""<<Crust_Z_position<<"\"/>\n";
192 
193  out<<"</define>\n";
194  out<<"<materials>\n";
195  out<<"</materials>\n";
196  out<<"<solids>\n";
197  out<<" <box name=\"WorldBox\" lunit=\"m\" x=\"2200\" y=\"2200\" z=\"2200\"/>\n";
198  out<<" <box name=\"CrustBox\" lunit=\"m\" x=\"2200\" y=\"2200\" z=\""<<Crust_Z_size<<"\"/>\n";
199  out<<" <box name=\"StoreyBox\" lunit=\"m\" x=\"0.6\" y=\"0.6\" z=\"0.6\"/>\n";
200  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";
201  out<<" <tube name=\"CathodTube\" lunit=\"cm\" aunit=\"deg\" rmin=\"0.0\" rmax=\"4.7462\" z=\"0.5\" startphi=\"0.0\" deltaphi=\"360.0\"/>\n";
202  out<<"</solids>\n\n\n";
203 
204  out<<"<structure>\n";
205  out<<" <volume name=\"CathodVolume\">\n";
206  out<<" <materialref ref=\"Cathod\"/>\n";
207  out<<" <solidref ref=\"CathodTube\"/>\n";
208  out<<" </volume>\n";
209 
210  out<<"<!-- OMVolume(s) construction -->\n";
211 
212  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
213 
214  out<<" <volume name=\"OMVolume"<<module->getID()<<"\">\n";
215  out<<" <materialref ref=\"Water\"/>\n";
216  out<<" <solidref ref=\"OMSphere\"/>\n";
217 
218  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
219  out<<" <physvol>\n";
220  out<<" <volumeref ref=\"CathodVolume\"/>\n";
221  out<<" <positionref ref=\"CathodPosition"<<pmt->getID()<<"_"<<module->getID()<<"\"/>\n";
222  out<<" <rotationref ref=\"CathodRotation"<<pmt->getID()<<"_"<<module->getID()<<"\"/>\n";
223  out<<" </physvol>\n";
224  }
225 
226  out<<" </volume>\n";
227  }
228 
229  out<<"<!-- StoreyVolume(s) construction -->\n";
230 
231  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
232  out<<" <volume name=\"StoreyVolume"<<module->getID()<<"\">\n";
233  out<<" <materialref ref=\"Water\"/>\n";
234  out<<" <solidref ref=\"StoreyBox\"/>\n";
235  out<<" <physvol>\n";
236  out<<" <volumeref ref=\"OMVolume"<<module->getID()<<"\"/>\n";
237  out<<" <positionref ref=\"OMShift\"/>\n";
238  //out<<" <positionref ref=\"zero\"/>\n";
239  out<<" <rotationref ref=\"identity\"/>\n";
240  out<<" </physvol>\n";
241  out<<" </volume>\n";
242  }
243 
244  out<<"<!-- Crust Volume construction-->\n";
245  out<<"<volume name=\"CrustVolume\">\n";
246  out<<" <materialref ref=\"Crust\"/>\n";
247  out<<" <solidref ref=\"CrustBox\"/>\n";
248  out<<"</volume>\n";
249 
250  out<<"<!-- World Volume construction-->\n";
251  out<<"<volume name=\"WorldVolume\">\n";
252  out<<" <materialref ref=\"Water\"/>\n";
253  out<<" <solidref ref=\"WorldBox\"/>\n";
254 
255  out<<" <physvol>\n";
256  out<<" <volumeref ref=\"CrustVolume\"/>\n";
257  out<<" <positionref ref=\"CrustPosition\"/>\n";
258  out<<" <rotationref ref=\"identity\"/>\n";
259  out<<" </physvol>\n";
260 
261  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
262  out<<" <physvol>\n";
263  out<<" <volumeref ref=\"StoreyVolume"<<module->getID()<<"\"/>\n";
264  out<<" <positionref ref=\"PosString"<< module->getString() <<"_Dom"<< module->getID() <<"\"/>\n";
265  //out<<" <rotation x="{}" y="{}" z="0.0"/>\n'.format(Dom_Phi[i], Dom_Theta[i]))
266  out<<" <rotationref ref=\"identity\"/>\n";
267  out<<" </physvol>\n";
268  }
269 
270  out<<"</volume>\n";
271 
272  out<<"</structure>\n";
273  out<<"<setup name=\"Default\" version=\"1.0\">\n";
274  out<<"<world ref=\"WorldVolume\"/>\n";
275  out<<"</setup>\n";
276  out<<"</gdml>\n";
277  }
278 
279 
280  /**
281  * Get maximal time between modules in detector following causality.
282  *
283  * \param detector detector
284  * \return maximal time [ns]
285  */
286  inline double getMaximalTime(const JDetector& detector)
287  {
288  return getMaximalDistance(detector) * getIndexOfRefraction() * getInverseSpeedOfLight();
289  }
290 
291 
292  /**
293  * Get maximal time between modules in detector following causality.
294  * The road width corresponds to the maximal distance traveled by the light.
295  *
296  * \param detector detector
297  * \param roadWidth_m road width [m]
298  * \return maximal time [ns]
299  */
300  inline double getMaximalTime(const JDetector& detector, const double roadWidth_m)
301  {
302  const double Dmax_m = getMaximalDistance(detector);
303 
304  return (sqrt((Dmax_m + roadWidth_m*getSinThetaC()) *
305  (Dmax_m - roadWidth_m*getSinThetaC())) +
306  roadWidth_m * getSinThetaC() * getTanThetaC()) * getInverseSpeedOfLight();
307  }
308 
309 
310  /**
311  * Get de-calibrated time range.
312  *
313  * The de-calibrated time range is corrected for minimal and maximal time offset of PMTs in given module.
314  *
315  * \param timeRange time range [ns]
316  * \param module module
317  * \return time range [ns]
318  */
319  inline JTimeRange getTimeRange(const JTimeRange& timeRange, const JModule& module)
320  {
321  if (!module.empty()) {
322 
323  JTimeRange time_range(JTimeRange::DEFAULT_RANGE);
324 
325  for (JModule::const_iterator pmt = module.begin(); pmt != module.end(); ++pmt) {
326 
327  const JCalibration& calibration = pmt->getCalibration();
328 
329  time_range.include(putTime(timeRange.getLowerLimit(), calibration));
330  time_range.include(putTime(timeRange.getUpperLimit(), calibration));
331  }
332 
333  return time_range;
334 
335  } else {
336 
337  return timeRange;
338  }
339  }
340 
341 
342  /**
343  * Get number of PMTs.
344  *
345  * \param module module
346  * \return number of PMTs
347  */
348  inline int getNumberOfPMTs(const JModule& module)
349  {
350  return module.size();
351  }
352 
353 
354  /**
355  * Get number of PMTs.
356  *
357  * \param detector detector
358  * \return number of PMTs
359  */
360  inline int getNumberOfPMTs(const JDetector& detector)
361  {
362  int number_of_pmts = 0;
363 
364  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
365  number_of_pmts += getNumberOfPMTs(*module);
366  }
367 
368  return number_of_pmts;
369  }
370 
371 
372  /**
373  * Get number of strings.
374  *
375  * \param detector detector
376  * \return number of strings
377  */
378  inline int getNumberOfStrings(const JDetector& detector)
379  {
380  std::set<int> buffer;
381 
382  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
383  buffer.insert(module->getString());
384  }
385 
386  return buffer.size();
387  }
388 
389 
390  /**
391  * Get list of strings IDs.
392  *
393  * \param detector detector
394  * \return list of string IDs
395  */
396  inline std::set<int> getStringIDs(const JDetector& detector)
397  {
398  std::set<int> buffer;
399 
400  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
401  buffer.insert(module->getString());
402  }
403 
404  return buffer;
405  }
406 
407 
408  /**
409  * Get number of floors.
410  *
411  * \param detector detector
412  * \return number of floors
413  */
414  inline int getNumberOfFloors(const JDetector& detector)
415  {
416  std::set<int> buffer;
417 
418  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
419  buffer.insert(module->getFloor());
420  }
421 
422  return buffer.size();
423  }
424 
425 
426  /**
427  * Get number of modules.
428  *
429  * \param detector detector
430  * \return number of modules
431  */
432  inline int getNumberOfModules(const JDetector& detector)
433  {
434  std::set<int> buffer;
435 
436  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
437  buffer.insert(module->getID());
438  }
439 
440  return buffer.size();
441  }
442 
443 
444  /**
445  * Load detector from input file.
446  *
447  * Supported file formats:
448  * - ASCII file, extension JDETECTOR::GENDET_DETECTOR_FILE_FORMAT, gendet format;
449  * - binary file, extension JDETECTOR::BINARY_DETECTOR_FILE_FORMAT, Jpp internal format;
450  * - ASCII file, extension JDETECTOR::KM3NET_DETECTOR_FILE_FORMAT, KM3NeT standard format;
451  * - gzipped file, extension JDETECTOR::ZIPPED_DETECTOR_FILE_FORMAT, KM3NeT standard format;
452  *
453  * \param file_name file name
454  * \param detector detector
455  */
456  inline void load(const JString& file_name, JDetector& detector)
457  {
458  using namespace std;
459  using namespace JIO;
460 
461  if (file_name.endsWith(GENDET_DETECTOR_FILE_FORMAT)) {
462 
463  JMonteCarloDetector buffer(true);
464 
465  ifstream in(file_name.c_str());
466 
467  if (!in) {
468  THROW(JFileOpenException, "File not opened: " << file_name);
469  }
470 
471  in >> buffer;
472 
473  in.close();
474 
475  detector.swap(buffer);
476 
477  } else if (file_name.endsWith(BINARY_DETECTOR_FILE_FORMAT)) {
478 
479  JFileStreamReader in(file_name.c_str());
480 
481  if (!in) {
482  THROW(JFileOpenException, "File not opened: " << file_name);
483  }
484 
485  detector.read(in);
486 
487  in.close();
488 
489  } else if (file_name.endsWith(KM3NET_DETECTOR_FILE_FORMAT)) {
490 
491  ifstream in(file_name.c_str());
492 
493  if (!in) {
494  THROW(JFileOpenException, "File not opened: " << file_name);
495  }
496 
497  in >> detector;
498 
499  in.close();
500 
501  } else if (file_name.endsWith(ZIPPED_DETECTOR_FILE_FORMAT)) {
502 
503  igzstream in(file_name.c_str());
504 
505  if (!in) {
506  THROW(JFileOpenException, "File not opened: " << file_name);
507  }
508 
509  in >> detector;
510 
511  in.close();
512 
513  } else {
514 
515  THROW(JProtocolException, "Protocol not defined: " << file_name);
516  }
517  }
518 
519 
520  /**
521  * Store detector to output file.
522  *
523  * Supported file formats:
524  * - binary file, extension JDETECTOR::BINARY_DETECTOR_FILE_FORMAT, Jpp internal format;
525  * - ASCII file, extension JDETECTOR::KM3NET_DETECTOR_FILE_FORMAT, KM3NeT standard format;
526  * - gzipped file, extension JDETECTOR::ZIPPED_DETECTOR_FILE_FORMAT, KM3NeT standard format;
527  * - gdml file, extension JDETECTOR::GDML_DETECTOR_FILE_FORMAT, KM3Sim input format;
528  *
529  * \param file_name file name
530  * \param detector detector
531  */
532  inline void store(const JString& file_name, const JDetector& detector)
533  {
534  using namespace std;
535  using namespace JIO;
536 
537  if (file_name.endsWith(BINARY_DETECTOR_FILE_FORMAT)) {
538 
539  JFileStreamWriter out(file_name.c_str());
540 
541  detector.write(out);
542 
543  out.close();
544 
545  } else if (file_name.endsWith(KM3NET_DETECTOR_FILE_FORMAT)) {
546 
547  std::ofstream out(file_name.c_str());
548 
549  out << detector;
550 
551  out.close();
552 
553  } else if (file_name.endsWith(ZIPPED_DETECTOR_FILE_FORMAT)) {
554 
555  ogzstream out(file_name.c_str());
556 
557  out << detector;
558 
559  out.close();
560 
561  } else if (file_name.endsWith(GDML_DETECTOR_FILE_FORMAT)) {
562 
563  std::ofstream out(file_name.c_str());
564 
565  write_gdml(out, detector);
566 
567  out.close();
568 
569  } else {
570 
571  throw JProtocolException("Protocol not defined.");
572  }
573  }
574 
575 
576  /**
577  * Get module according module address map.
578  *
579  * \param memo module address map
580  * \param id module identifier
581  * \param location module location
582  * \return module
583  */
584  inline const JModule& getModule(const JModuleAddressMap& memo,
585  const int id,
586  const JModuleLocation& location = JModuleLocation())
587  {
588  static JModule module;
589 
590 
591  module.setID(id);
592 
593  module.setLocation(location);
594 
595  module.resize(memo.size());
596 
597  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))); }
598 
599  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))); }
600  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))); }
601  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))); }
602  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))); }
603  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))); }
604  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))); }
605 
606  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))); }
607  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))); }
608  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))); }
609  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))); }
610  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))); }
611  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))); }
612 
613  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))); }
614  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))); }
615  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))); }
616  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))); }
617  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))); }
618  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))); }
619 
620  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))); }
621  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))); }
622  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))); }
623  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))); }
624  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))); }
625  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))); }
626 
627  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))); }
628  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))); }
629  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))); }
630  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))); }
631  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))); }
632  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))); }
633 
634  module.compile();
635 
636  return module;
637  }
638 
639 
640  /**
641  * Get module corresponding to Antares storey.
642  *
643  * \param id module identifier
644  * \param location module location
645  * \return module
646  */
647  inline const JModule& getModule(const int id,
648  const JModuleLocation& location = JModuleLocation())
649  {
650  using JTOOLS::PI;
651 
652  static JModule module;
653 
654  module.setID(id);
655 
656  module.setLocation(location);
657 
658  if (module.empty()) {
659 
660  module.resize(3);
661 
662  const double R = 0.5; // [m]
663 
664  const double st = sin(0.75*PI);
665  const double ct = cos(0.75*PI);
666 
667  for (int i = 0; i != 3; ++i) {
668 
669  const double phi = (2.0*PI*i) / 3.0;
670  const double cp = cos(phi);
671  const double sp = sin(phi);
672 
673  module[i] = JPMT(i, JAxis3D(JVector3D(R*st*cp, R*st*sp, R*ct), JVersor3D(st*cp, st*sp, ct)));
674  }
675  }
676 
677  return module;
678  }
679 
680 
681  /**
682  * Get module label for monitoring and other applications.\
683  * The format is "(XXX,YY)", where XXX is the string number and YY the floor.
684  *
685  * \param location module location
686  * \retun label
687  */
688  inline std::string getModuleLabel(const JModuleLocation& location)
689  {
690  using namespace std;
691 
692  ostringstream os;
693 
694  os << "(" << setfill('0') << setw(3) << location.getString() << "," << setfill('0') << setw(2) << location.getFloor() << ")";
695 
696  return os.str();
697  }
698 
699 
700  /**
701  * Auxiliary class to get rotation matrix between two optical modules.
702  */
703  struct JRotation :
704  public JRotation3D
705  {
706 
707  static const size_t NUMBER_OF_DIMENSIONS = 3; //!< Number of dimensions
708 
709 
710  /**
711  * Get rotation matrix to go from first to second module.
712  *
713  * \param first first module
714  * \param second second module
715  * \return rotation matrix
716  */
717  const JRotation3D& operator()(const JModule& first, const JModule& second)
718  {
719  this->setIdentity();
720 
721  if (first.size() == second.size()) {
722 
723  const size_t N = first.size();
724 
725  if (N >= NUMBER_OF_DIMENSIONS) {
726 
727  in .resize(N);
728  out.resize(N);
729 
730  for (size_t i = 0; i != N; ++i) {
731  in [i] = first .getPMT(i).getDirection();
732  out[i] = second.getPMT(i).getDirection();
733  }
734 
735  for (size_t i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
736  if (!orthonormalise(i)) {
737  THROW(JException, "Failure to orthonormalise direction " << i);
738  }
739  }
740 
741  this->a00 = out[0].getX() * in[0].getX() + out[1].getX() * in[1].getX() + out[2].getX() * in[2].getX();
742  this->a01 = out[0].getX() * in[0].getY() + out[1].getX() * in[1].getY() + out[2].getX() * in[2].getY();
743  this->a02 = out[0].getX() * in[0].getZ() + out[1].getX() * in[1].getZ() + out[2].getX() * in[2].getZ();
744 
745  this->a10 = out[0].getY() * in[0].getX() + out[1].getY() * in[1].getX() + out[2].getY() * in[2].getX();
746  this->a11 = out[0].getY() * in[0].getY() + out[1].getY() * in[1].getY() + out[2].getY() * in[2].getY();
747  this->a12 = out[0].getY() * in[0].getZ() + out[1].getY() * in[1].getZ() + out[2].getY() * in[2].getZ();
748 
749  this->a20 = out[0].getZ() * in[0].getX() + out[1].getZ() * in[1].getX() + out[2].getZ() * in[2].getX();
750  this->a21 = out[0].getZ() * in[0].getY() + out[1].getZ() * in[1].getY() + out[2].getZ() * in[2].getY();
751  this->a22 = out[0].getZ() * in[0].getZ() + out[1].getZ() * in[1].getZ() + out[2].getZ() * in[2].getZ();
752 
753  } else {
754 
755  THROW(JException, "Module size" << N << " < " << NUMBER_OF_DIMENSIONS);
756  }
757 
758  } else {
759 
760  THROW(JException, "Module size" << first.size() << " != " << second.size());
761  }
762 
763  return *this;
764  }
765 
766  private:
767  /**
768  * Put normalised primary direction at specified index and orthoganilise others.\n
769  * This procedure follows Gram-Schmidt process.
770  *
771  * \param index index
772  * \return true if primary direction exists; else false
773  */
774  bool orthonormalise(const size_t index)
775  {
776  using namespace std;
777 
778  size_t pos = index;
779 
780  for (size_t i = index + 1; i != in.size(); ++i) {
781  if (in[i].getLengthSquared() > in[pos].getLengthSquared()) {
782  pos = i;
783  }
784  }
785 
786  const double u = in[pos].getLength();
787 
788  if (u > 0.0) {
789 
790  in [pos] /= u;
791  out[pos] /= u;
792 
793  if (pos != index) {
794  swap(in [pos], in [index]);
795  swap(out[pos], out[index]);
796  }
797 
798  for (size_t i = index + 1; i != in.size(); ++i) {
799 
800  const double dot = in[index].getDot(in[i]);
801 
802  in [i] -= dot * in [index];
803  out[i] -= dot * out[index];
804  }
805 
806  return true;
807 
808  } else {
809 
810  return false;
811  }
812  }
813 
814 
815  std::vector<JVector3D> in;
816  std::vector<JVector3D> out;
817  };
818 
819 
820  /**
821  * Function object to get rotation matrix to go from first to second module.
822  */
823  static JRotation getRotation;
824 
825 
826  /**
827  * Get position to go from first to second module.
828  *
829  * \param first first module
830  * \param second second module
831  * \return position
832  */
833  inline JPosition3D getPosition(const JModule& first, const JModule& second)
834  {
835  return second.getPosition() - first.getPosition();
836  }
837 
838 
839  /**
840  * Get calibration to go from first to second calibration.
841  *
842  * \param first first calibration
843  * \param second second calibration
844  * \return calibration
845  */
846  inline JCalibration getCalibration(const JCalibration& first, const JCalibration& second)
847  {
848  return JCalibration(second.getT0() - first.getT0());
849  }
850 }
851 
852 #endif
JException.hh
JGEOMETRY3D::JVersor3D::getDZ
double getDZ() const
Get z direction.
Definition: JVersor3D.hh:113
JDETECTOR::KM3NET_DETECTOR_FILE_FORMAT
static const char *const KM3NET_DETECTOR_FILE_FORMAT
KM3NeT standard ASCII format.
Definition: JDetectorToolkit.hh:64
JIO::JFileStreamWriter
Binary buffered file output.
Definition: JFileStreamIO.hh:72
JDETECTOR::JRotation::out
std::vector< JVector3D > out
Definition: JDetectorToolkit.hh:816
JTOOLS::getSinThetaC
double getSinThetaC()
Get average sine of Cherenkov angle of water.
Definition: JConstants.hh:155
JDETECTOR::JRotation::orthonormalise
bool orthonormalise(const size_t index)
Put normalised primary direction at specified index and orthoganilise others.
Definition: JDetectorToolkit.hh:774
JIO::JFileStreamReader
Binary buffered file input.
Definition: JFileStreamIO.hh:22
JDETECTOR::read_gdml
void read_gdml(std::istream &, const JDetector &)
Definition: JDetectorToolkit.hh:123
JDETECTOR::getModuleLabel
std::string getModuleLabel(const JModuleLocation &location)
Get module label for monitoring and other applications.
Definition: JDetectorToolkit.hh:688
JVector3D.hh
JDETECTOR::getCalibration
JCalibration getCalibration(const JCalibration &first, const JCalibration &second)
Get calibration to go from first to second calibration.
Definition: JDetectorToolkit.hh:846
JDETECTOR::getMaximalDistance
double getMaximalDistance(const JDetector &detector)
Get maximal distance between modules in detector.
Definition: JDetectorToolkit.hh:76
JDETECTOR::JRotation
Auxiliary class to get rotation matrix between two optical modules.
Definition: JDetectorToolkit.hh:703
igzstream
Definition: gzstream.h:203
JFileStreamIO.hh
JDETECTOR::JModuleAddressMap::has
bool has(const int index) const
Test whether index is available.
Definition: JModuleAddressMap.hh:103
JGEOMETRY3D::JDirection3D::getDirection
const JDirection3D & getDirection() const
Get direction.
Definition: JDirection3D.hh:106
JGEOMETRY3D::JVector3D::getZ
double getZ() const
Get z position.
Definition: JVector3D.hh:114
JGEOMETRY3D::JAxis3D
Axis object.
Definition: JAxis3D.hh:38
JLANG::JFileOpenException
Exception for opening of file.
Definition: JException.hh:342
JDETECTOR::JCalibration
Data structure for PMT calibration.
Definition: JDetector/JCalibration.hh:35
getModule
const JModule & getModule(const JDetector &detector, const JModuleLocation &location)
find module with a given string and floor number
Definition: JBeaconSimulator.cc:927
JDETECTOR::load
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
Definition: JDetectorToolkit.hh:456
JPoint3D.hh
JDETECTOR::GetYrotationG4
double GetYrotationG4(const JVersor3D dir)
Get rotation over Y axis in Geant4 coordinate system.
Definition: JDetectorToolkit.hh:117
JTOOLS::u
double u[N+1]
Definition: JPolint.hh:706
JDETECTOR::putTime
double putTime(const T &t1, const JCalibration &cal)
Get de-calibrated time.
Definition: JDetector/JCalibration.hh:216
std::vector
Definition: JSTDTypes.hh:12
JDETECTOR::getMaximalTime
double getMaximalTime(const JDetector &detector, const double roadWidth_m)
Get maximal time between modules in detector following causality.
Definition: JDetectorToolkit.hh:300
JPoint3DEstimator.hh
JDETECTOR::JCalibration::getT0
double getT0() const
Get time offset.
Definition: JDetector/JCalibration.hh:83
JTOOLS::j
int j
Definition: JPolint.hh:634
JTOOLS::JRange< double >::DEFAULT_RANGE
static const JRange< double, std::less< double > > DEFAULT_RANGE
Default range.
Definition: JRange.hh:558
JDETECTOR::GDML_DETECTOR_FILE_FORMAT
static const char *const GDML_DETECTOR_FILE_FORMAT
KM3Sim input format.
Definition: JDetectorToolkit.hh:66
JGEOMETRY3D::JVersor3D
Data structure for normalised vector in three dimensions.
Definition: JVersor3D.hh:23
ogzstream
Definition: gzstream.h:213
JDETECTOR::JDetector::read
virtual JReader & read(JReader &in)
Read from input.
Definition: JDetector.hh:370
JDETECTOR::getNumberOfStrings
int getNumberOfStrings(const JDetector &detector)
Get number of strings.
Definition: JDetectorToolkit.hh:378
JFIT::JEstimator< JPoint3D >
Linear fit of crossing point (position) between axes (objects with position and direction).
Definition: JPoint3DEstimator.hh:55
JModuleAddressMap.hh
JDETECTOR::JModuleAddressMap
Lookup table for PMT addresses in optical module.
Definition: JModuleAddressMap.hh:82
JString.hh
JDETECTOR::getNumberOfFloors
int getNumberOfFloors(const JDetector &detector)
Get number of floors.
Definition: JDetectorToolkit.hh:414
JTimeRange.hh
std::set< int >
JDETECTOR::JDetector::write
virtual JWriter & write(JWriter &out) const
Write to output.
Definition: JDetector.hh:400
JTOOLS::JTimeRange
JRange< double > JTimeRange
Type definition for time range.
Definition: JTools/JTimeRange.hh:19
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:37
JDETECTOR::JModuleLocation
Logical location of module.
Definition: JModuleLocation.hh:36
JGEOMETRY3D::JCylinder3D
Cylinder object.
Definition: JCylinder3D.hh:37
JLANG::JString::endsWith
bool endsWith(const std::string &suffix) const
Test if this string ends with the specified suffix.
Definition: JString.hh:175
JGEOMETRY3D::JVector3D
Data structure for vector in three dimensions.
Definition: JVector3D.hh:33
JDETECTOR::JModule::getPMT
const JPMT & getPMT(const int index) const
Get PMT.
Definition: JModule.hh:173
JTOOLS::getInverseSpeedOfLight
const double getInverseSpeedOfLight()
Get inverse speed of light.
Definition: JConstants.hh:100
JLANG::JString
Wrapper class around STL string class.
Definition: JString.hh:28
JDETECTOR::getNumberOfPMTs
int getNumberOfPMTs(const JDetector &detector)
Get number of PMTs.
Definition: JDetectorToolkit.hh:360
JDETECTOR::getNumberOfModules
int getNumberOfModules(const JDetector &detector)
Get number of modules.
Definition: JDetectorToolkit.hh:432
JConstants.hh
JGEOMETRY3D::JPosition3D
Data structure for position in three dimensions.
Definition: JPosition3D.hh:35
JLANG::JObjectID::setID
void setID(const int id)
Set identifier.
Definition: JObjectID.hh:66
JDETECTOR::GDML_SCHEMA
static const char *const GDML_SCHEMA
directory necessary for correct GDML header output
Definition: JDetectorToolkit.hh:67
JTOOLS::getIndexOfRefraction
double getIndexOfRefraction()
Get average index of refraction of water.
Definition: JConstants.hh:111
JDETECTOR::JCalibration::getCalibration
const JCalibration & getCalibration() const
Get calibration.
Definition: JDetector/JCalibration.hh:61
JDETECTOR::ZIPPED_DETECTOR_FILE_FORMAT
static const char *const ZIPPED_DETECTOR_FILE_FORMAT
zipped KM3NeT standard ASCII format
Definition: JDetectorToolkit.hh:65
JDETECTOR::GENDET_DETECTOR_FILE_FORMAT
static const char *const GENDET_DETECTOR_FILE_FORMAT
File name extensions.
Definition: JDetectorToolkit.hh:62
JMonteCarloDetector.hh
THROW
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:669
JGEOMETRY3D::JVersor3D::getDX
double getDX() const
Get x direction.
Definition: JVersor3D.hh:91
JDETECTOR::JModule
Data structure for a composite optical module.
Definition: JModule.hh:49
JLANG::JProtocolException
Protocol exception.
Definition: JException.hh:630
JDETECTOR::BINARY_DETECTOR_FILE_FORMAT
static const char *const BINARY_DETECTOR_FILE_FORMAT
JIO binary file format.
Definition: JDetectorToolkit.hh:63
JDETECTOR::JRotation::in
std::vector< JVector3D > in
Definition: JDetectorToolkit.hh:815
JGEOMETRY3D::JVersor3D::getDY
double getDY() const
Get y direction.
Definition: JVersor3D.hh:102
JDETECTOR::JPMT
Data structure for PMT geometry and calibration.
Definition: JPMT.hh:53
JDETECTOR::JModuleLocation::setLocation
void setLocation(const JModuleLocation &location)
Set location.
Definition: JModuleLocation.hh:90
JDETECTOR::store
void store(const JString &file_name, const JDetector &detector)
Store detector to output file.
Definition: JDetectorToolkit.hh:532
JModuleLocation.hh
JTOOLS::getTanThetaC
double getTanThetaC()
Get average tangent of Cherenkov angle of water.
Definition: JConstants.hh:133
JDETECTOR::JRotation::operator()
const JRotation3D & operator()(const JModule &first, const JModule &second)
Get rotation matrix to go from first to second module.
Definition: JDetectorToolkit.hh:717
gzstream.h
JDETECTOR::JDetector
Detector data structure.
Definition: JDetector.hh:80
JGEOMETRY3D::JPosition3D::getPosition
const JPosition3D & getPosition() const
Get position.
Definition: JPosition3D.hh:129
std
Definition: jaanetDictionary.h:36
JGEOMETRY3D::JVector3D::getY
double getY() const
Get y position.
Definition: JVector3D.hh:103
JCylinder3D.hh
JTOOLS::PI
static const double PI
Constants.
Definition: JConstants.hh:20
JFIT::JPoint3D
Data structure for position fit.
Definition: JPoint3D.hh:22
JIO
Auxiliary classes and methods for binary I/O.
Definition: JBinaryFileReader.hh:17
JDETECTOR::getStringIDs
std::set< int > getStringIDs(const JDetector &detector)
Get list of strings IDs.
Definition: JDetectorToolkit.hh:396
JDETECTOR::getPosition
JPosition3D getPosition(const JModule &first, const JModule &second)
Get position to go from first to second module.
Definition: JDetectorToolkit.hh:833
JModuleGeometry.hh
JDetector.hh
JGEOMETRY3D::JVector3D::getX
double getX() const
Get x position.
Definition: JVector3D.hh:93
JVersor3D.hh
JDETECTOR::JModuleLocation::getFloor
int getFloor() const
Get floor number.
Definition: JModuleLocation.hh:144
JDETECTOR::getRotation
static JRotation getRotation
Function object to get rotation matrix to go from first to second module.
Definition: JDetectorToolkit.hh:823
JDETECTOR::JModule::compile
void compile()
Compile position of module from the positions of the PMTs.
Definition: JModule.hh:210
JDETECTOR::getTimeRange
JTimeRange getTimeRange(const JTimeRange &timeRange, const JModule &module)
Get de-calibrated time range.
Definition: JDetectorToolkit.hh:319
JGEOMETRY3D::JVector3D::sub
JVector3D & sub(const JVector3D &vector)
Subtract vector.
Definition: JVector3D.hh:157
JDETECTOR::GetXrotationG4
double GetXrotationG4(const JVersor3D dir)
Get rotation over X axis in Geant4 coordinate system.
Definition: JDetectorToolkit.hh:98
JDETECTOR
Auxiliary classes and methods for detector calibration.
Definition: JAnchor.hh:12
JLANG::JException
General exception.
Definition: JException.hh:40
JDETECTOR::JMonteCarloDetector
Monte Carlo detector (i.e.
Definition: JMonteCarloDetector.hh:43
JGEOMETRY3D::JRotation3D
Rotation matrix.
Definition: JRotation3D.hh:111
JDETECTOR::JModuleLocation::getString
int getString() const
Get string number.
Definition: JModuleLocation.hh:133
JDETECTOR::write_gdml
void write_gdml(std::ostream &out, const JDetector &detector)
Writes KM3Sim GDML input file from detector.
Definition: JDetectorToolkit.hh:133
Generated by   doxygen 1.8.16