Jpp
 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 <fstream>
8 #include <cmath>
9 #include <set>
10 
11 #include "JDetector/JDetector.hh"
13 #include "JDetector/JTimeRange.hh"
17 #include "JGeometry3D/JVector3D.hh"
18 #include "JGeometry3D/JVersor3D.hh"
20 #include "JTools/JConstants.hh"
21 #include "JIO/JFileStreamIO.hh"
22 #include "JLang/JString.hh"
23 #include "JLang/JException.hh"
24 #include "JLang/gzstream.h"
26 #include "JFit/JPoint3D.hh"
27 
28 
29 /**
30  * \author mdejong
31  */
32 
33 namespace JDETECTOR {}
34 namespace JPP { using namespace JDETECTOR; }
35 
36 /**
37  * Auxiliary classes and methods for detector calibration and simulation.
38  */
39 namespace JDETECTOR {
40 
47  using JTOOLS::PI;
48  using JLANG::JString;
49  using JLANG::JException;
53  using JFIT::JPoint3D;
54 
55 
56 
57  /**
58  * File name extensions.
59  */
60  static const char* const GENDET_DETECTOR_FILE_FORMAT = ".det"; //!< file format used by gendet
61  static const char* const BINARY_DETECTOR_FILE_FORMAT = ".dat"; //!< JIO binary file format
62  static const char* const KM3NET_DETECTOR_FILE_FORMAT = ".detx"; //!< KM3NeT standard ASCII format
63  static const char* const ZIPPED_DETECTOR_FILE_FORMAT = ".gz"; //!< zipped KM3NeT standard ASCII format
64  static const char* const GDML_DETECTOR_FILE_FORMAT = ".xml"; //!< KM3Sim input format
65 
66 
67  /**
68  * Get maximal distance between modules in detector.
69  *
70  * \param detector detector
71  * \return maximal distance [m]
72  */
73  inline double getMaximalDistance(const JDetector& detector)
74  {
75  double dmax = 0.0;
76 
77  for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
78  for (JDetector::const_iterator j = detector.begin(); j != i; ++j) {
79  if (i->getDistance(*j) > dmax) {
80  dmax = i->getDistance(*j);
81  }
82  }
83  }
84 
85  return dmax;
86  }
87 
88 
89  /**
90  * Get rotation over X axis in Geant4 coordinate system
91  *
92  * \param dir direction
93  * \return X-rotation [deg]
94  */
95  inline double GetXrotationG4(const JVersor3D dir)
96  {
97  const double phi = atan2(dir.getDY(), dir.getDZ())*(180.0/PI);
98 
99  if (phi < 0.0){
100  return phi + 360.0;
101  }
102  else{
103  return phi;
104  }
105  }
106 
107 
108  /**
109  * Get rotation over Y axis in Geant4 coordinate system
110  *
111  * \param dir direction
112  * \return Y-rotation [deg]
113  */
114  inline double GetYrotationG4(const JVersor3D dir)
115  {
116  return asin(-dir.getDX())*(180.0/PI);
117  }
118 
119 
120  inline void read_gdml(std::istream&, const JDetector&)
121  {}
122 
123 
124  /**
125  * Writes KM3Sim GDML input file from detector
126  *
127  * \param out output stream
128  * \param detector detector
129  */
130  inline void write_gdml(std::ostream& out, const JDetector& detector)
131  {
132  const JCylinder3D cylinder(detector.begin(), detector.end());
133 
134  // Detector height is bascially just a distance between the highest and the lowest DOM
135  //double Detector_Height = cylinder.getZmax() - cylinder.getZmin();
136  // Distance from the center of the lowest DOM to the seabed
137  //double Seabed_Distance = cylinder.getZmin();
138 
139  // To translate Z coordinate from detx format which corresponds to the distance from the seabed one has to perform a following calculation
140  // Dom_Z_GDML = Dom_Z_Detx - Distance to the seabed from the lowest DOM - Detector_height/2
141 
142  //double DistanceOfSeabedFromLowestOM = Seabed_Distance;
143  // For now, the size of the WorldBox is fixed to 2200 meters
144  const double WorldBoxHeight = 2200.;
145  // Below formula comes directly from the note by Apostolos
146  //double Crust_Z_size = WorldBoxHeight/2-Detector_Height/2-DistanceOfSeabedFromLowestOM; Formula for obsolete coordinate covention, where Z_seabed != 0
147  const double Crust_Z_size = WorldBoxHeight/2;
148  //double Crust_Z_position = -WorldBoxHeight/4-Detector_Height/4-DistanceOfSeabedFromLowestOM/2; Formula for obsolete coordinate covention, where Z_seabed != 0
149  const double Crust_Z_position= -WorldBoxHeight/4;
150 
151 
152  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=\"/afs/cern.ch/sw/geant4/releases/share/geant4.9.5/source/persistency/gdml/schema/gdml.xsd\">\n\n\n";
153  out<<"<define>\n";
154  out<<"<rotation name=\"identity\"/>\n<position name=\"zero\"/>\n";
155 
156  int PMTs_NO = 0;
157  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
158 
159  const JFIT::JEstimator<JPoint3D> center(module->begin(), module->end());
160 
161  //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
162  out<<"<position name=\"PosString"<<module->getString()<<"_Dom"<<module->getID()<<"\" unit=\"m\" x=\""<<module->getX()<<"\" y=\""<<module->getY()<<"\" z=\""<<module->getZ()<<"\"/>\n";
163 
164  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
165  //JVector3D vec;
166  //static_cast<JVector3D> (vec) = module->sub(*pmt);
167 
168  const JVector3D vec = static_cast<JVector3D>(*pmt).sub(center);
169  out<<"<position name=\"CathodPosition"<<pmt->getID()<<"_"<<module->getID()<<"\" unit=\"m\" x=\""<<vec.getX()<<"\" y=\""<<vec.getY()<<"\" z=\""<<vec.getZ()<<"\"/>\n";
170  out<<"<rotation name=\"CathodRotation"<<pmt->getID()<<"_"<<module->getID()<<"\" unit=\"deg\" x=\""<<GetXrotationG4(*pmt)<<"\" y=\""<<GetYrotationG4(*pmt)<<"\" z=\"0.000000\"/>\n";
171  out<<"<constant name=\"CathodID_"<<PMTs_NO<<"\" value=\""<<pmt->getID()<<"\"/>\n";
172  PMTs_NO++;
173  }
174 
175  }
176 
177  out<<"<position name=\"OMShift\" unit=\"m\" x=\"0\" y=\"0\" z=\"0.0392\"/>\n";
178  out<<"\n\n\n";
179  out<<"<!-- end of DU position definitions -->\n<position name=\"CrustPosition\" unit=\"m\" x=\"0\" y=\"0\" z=\""<<Crust_Z_position<<"\"/>\n";
180 
181  out<<"</define>\n";
182  out<<"<materials>\n";
183  out<<"</materials>\n";
184  out<<"<solids>\n";
185  out<<" <box name=\"WorldBox\" lunit=\"m\" x=\"2200\" y=\"2200\" z=\"2200\"/>\n";
186  out<<" <box name=\"CrustBox\" lunit=\"m\" x=\"2200\" y=\"2200\" z=\""<<Crust_Z_size<<"\"/>\n";
187  out<<" <box name=\"StoreyBox\" lunit=\"m\" x=\"0.6\" y=\"0.6\" z=\"0.6\"/>\n";
188  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";
189  out<<" <tube name=\"CathodTube\" lunit=\"cm\" aunit=\"deg\" rmin=\"0.0\" rmax=\"4.7462\" z=\"0.5\" startphi=\"0.0\" deltaphi=\"360.0\"/>\n";
190  out<<"</solids>\n\n\n";
191 
192  out<<"<structure>\n";
193  out<<" <volume name=\"CathodVolume\">\n";
194  out<<" <materialref ref=\"Cathod\"/>\n";
195  out<<" <solidref ref=\"CathodTube\"/>\n";
196  out<<" </volume>\n";
197 
198  out<<"<!-- OMVolume(s) construction -->\n";
199 
200  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
201 
202  out<<" <volume name=\"OMVolume"<<module->getID()<<"\">\n";
203  out<<" <materialref ref=\"Water\"/>\n";
204  out<<" <solidref ref=\"OMSphere\"/>\n";
205 
206  for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
207  out<<" <physvol>\n";
208  out<<" <volumeref ref=\"CathodVolume\"/>\n";
209  out<<" <positionref ref=\"CathodPosition"<<pmt->getID()<<"_"<<module->getID()<<"\"/>\n";
210  out<<" <rotationref ref=\"CathodRotation"<<pmt->getID()<<"_"<<module->getID()<<"\"/>\n";
211  out<<" </physvol>\n";
212  }
213 
214  out<<" </volume>\n";
215  }
216 
217  out<<"<!-- StoreyVolume(s) construction -->\n";
218 
219  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
220  out<<" <volume name=\"StoreyVolume"<<module->getID()<<"\">\n";
221  out<<" <materialref ref=\"Water\"/>\n";
222  out<<" <solidref ref=\"StoreyBox\"/>\n";
223  out<<" <physvol>\n";
224  out<<" <volumeref ref=\"OMVolume"<<module->getID()<<"\"/>\n";
225  out<<" <positionref ref=\"OMShift\"/>\n";
226  //out<<" <positionref ref=\"zero\"/>\n";
227  out<<" <rotationref ref=\"identity\"/>\n";
228  out<<" </physvol>\n";
229  out<<" </volume>\n";
230  }
231 
232  out<<"<!-- Crust Volume construction-->\n";
233  out<<"<volume name=\"CrustVolume\">\n";
234  out<<" <materialref ref=\"Crust\"/>\n";
235  out<<" <solidref ref=\"CrustBox\"/>\n";
236  out<<"</volume>\n";
237 
238  out<<"<!-- World Volume construction-->\n";
239  out<<"<volume name=\"WorldVolume\">\n";
240  out<<" <materialref ref=\"Water\"/>\n";
241  out<<" <solidref ref=\"WorldBox\"/>\n";
242 
243  out<<" <physvol>\n";
244  out<<" <volumeref ref=\"CrustVolume\"/>\n";
245  out<<" <positionref ref=\"CrustPosition\"/>\n";
246  out<<" <rotationref ref=\"identity\"/>\n";
247  out<<" </physvol>\n";
248 
249  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
250  out<<" <physvol>\n";
251  out<<" <volumeref ref=\"StoreyVolume"<<module->getID()<<"\"/>\n";
252  out<<" <positionref ref=\"PosString"<< module->getString() <<"_Dom"<< module->getID() <<"\"/>\n";
253  //out<<" <rotation x="{}" y="{}" z="0.0"/>\n'.format(Dom_Phi[i], Dom_Theta[i]))
254  out<<" <rotationref ref=\"identity\"/>\n";
255  out<<" </physvol>\n";
256  }
257 
258  out<<"</volume>\n";
259 
260  out<<"</structure>\n";
261  out<<"<setup name=\"Default\" version=\"1.0\">\n";
262  out<<"<world ref=\"WorldVolume\"/>\n";
263  out<<"</setup>\n";
264  out<<"</gdml>\n";
265  }
266 
267 
268  /**
269  * Get maximal time between modules in detector following causality.
270  *
271  * \param detector detector
272  * \return maximal time [ns]
273  */
274  inline double getMaximalTime(const JDetector& detector)
275  {
277  }
278 
279 
280  /**
281  * Get maximal time between modules in detector following causality.
282  * The road width corresponds to the maximal distance traveled by the light.
283  *
284  * \param detector detector
285  * \param roadWidth_m road width [m]
286  * \return maximal time [ns]
287  */
288  inline double getMaximalTime(const JDetector& detector, const double roadWidth_m)
289  {
290  const double Dmax_m = getMaximalDistance(detector);
291 
292  return (sqrt((Dmax_m + roadWidth_m*getSinThetaC()) *
293  (Dmax_m - roadWidth_m*getSinThetaC())) +
294  roadWidth_m * getSinThetaC() * getTanThetaC()) * getInverseSpeedOfLight();
295  }
296 
297 
298  /**
299  * Get de-calibrated time range.
300  *
301  * The de-calibrated time range is corrected for minimal and maximal time offset of PMTs in given module.
302  *
303  * \param timeRange time range [ns]
304  * \param module module
305  * \return time range [ns]
306  */
307  inline JTimeRange getTimeRange(const JTimeRange& timeRange, const JModule& module)
308  {
309  if (!module.empty()) {
310 
312 
313  for (JModule::const_iterator pmt = module.begin(); pmt != module.end(); ++pmt) {
314 
315  const JCalibration& calibration = pmt->getCalibration();
316 
317  time_range.include(putTime(timeRange.getLowerLimit(), calibration));
318  time_range.include(putTime(timeRange.getUpperLimit(), calibration));
319  }
320 
321  return time_range;
322 
323  } else {
324 
325  return timeRange;
326  }
327  }
328 
329 
330  /**
331  * Get number of PMTs.
332  *
333  * \param module module
334  * \return number of PMTs
335  */
336  inline int getNumberOfPMTs(const JModule& module)
337  {
338  return module.size();
339  }
340 
341 
342  /**
343  * Get number of PMTs.
344  *
345  * \param detector detector
346  * \return number of PMTs
347  */
348  inline int getNumberOfPMTs(const JDetector& detector)
349  {
350  int number_of_pmts = 0;
351 
352  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
353  number_of_pmts += getNumberOfPMTs(*module);
354  }
355 
356  return number_of_pmts;
357  }
358 
359 
360  /**
361  * Get number of strings.
362  *
363  * \param detector detector
364  * \return number of strings
365  */
366  inline int getNumberOfStrings(const JDetector& detector)
367  {
368  std::set<int> buffer;
369 
370  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
371  buffer.insert(module->getString());
372  }
373 
374  return buffer.size();
375  }
376 
377 
378  /**
379  * Get list of strings IDs.
380  *
381  * \param detector detector
382  * \return list of string IDs
383  */
384  inline std::set<int> getStringIDs(const JDetector& detector)
385  {
386  std::set<int> buffer;
387 
388  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
389  buffer.insert(module->getString());
390  }
391 
392  return buffer;
393  }
394 
395 
396  /**
397  * Get number of floors.
398  *
399  * \param detector detector
400  * \return number of floors
401  */
402  inline int getNumberOfFloors(const JDetector& detector)
403  {
404  std::set<int> buffer;
405 
406  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
407  buffer.insert(module->getFloor());
408  }
409 
410  return buffer.size();
411  }
412 
413 
414  /**
415  * Get number of modules.
416  *
417  * \param detector detector
418  * \return number of modules
419  */
420  inline int getNumberOfModules(const JDetector& detector)
421  {
422  std::set<int> buffer;
423 
424  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
425  buffer.insert(module->getID());
426  }
427 
428  return buffer.size();
429  }
430 
431 
432  /**
433  * Load detector from input file.
434  *
435  * Supported file formats:
436  * - ASCII file, extension JDETECTOR::GENDET_DETECTOR_FILE_FORMAT, gendet format;
437  * - binary file, extension JDETECTOR::BINARY_DETECTOR_FILE_FORMAT, Jpp internal format;
438  * - ASCII file, extension JDETECTOR::KM3NET_DETECTOR_FILE_FORMAT, KM3NeT standard format;
439  * - gzipped file, extension JDETECTOR::ZIPPED_DETECTOR_FILE_FORMAT, KM3NeT standard format;
440  *
441  * \param file_name file name
442  * \param detector detector
443  */
444  inline void load(const JString& file_name, JDetector& detector)
445  {
446  using namespace std;
447  using namespace JIO;
448 
449  if (file_name.endsWith(GENDET_DETECTOR_FILE_FORMAT)) {
450 
451  JMonteCarloDetector buffer(true);
452 
453  ifstream in(file_name.c_str());
454 
455  if (!in) {
456  THROW(JFileOpenException, "File not opened: " << file_name);
457  }
458 
459  in >> buffer;
460 
461  in.close();
462 
463  detector.swap(buffer);
464 
465  } else if (file_name.endsWith(BINARY_DETECTOR_FILE_FORMAT)) {
466 
467  JFileStreamReader in(file_name.c_str());
468 
469  if (!in) {
470  THROW(JFileOpenException, "File not opened: " << file_name);
471  }
472 
473  detector.read(in);
474 
475  in.close();
476 
477  } else if (file_name.endsWith(KM3NET_DETECTOR_FILE_FORMAT)) {
478 
479  ifstream in(file_name.c_str());
480 
481  if (!in) {
482  THROW(JFileOpenException, "File not opened: " << file_name);
483  }
484 
485  in >> detector;
486 
487  in.close();
488 
489  } else if (file_name.endsWith(ZIPPED_DETECTOR_FILE_FORMAT)) {
490 
491  igzstream 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 {
502 
503  THROW(JProtocolException, "Protocol not defined: " << file_name);
504  }
505  }
506 
507 
508  /**
509  * Store detector to output file.
510  *
511  * Supported file formats:
512  * - binary file, extension JDETECTOR::BINARY_DETECTOR_FILE_FORMAT, Jpp internal format;
513  * - ASCII file, extension JDETECTOR::KM3NET_DETECTOR_FILE_FORMAT, KM3NeT standard format;
514  * - gzipped file, extension JDETECTOR::ZIPPED_DETECTOR_FILE_FORMAT, KM3NeT standard format;
515  * - gdml file, extension JDETECTOR::GDML_DETECTOR_FILE_FORMAT, KM3Sim input format;
516  *
517  * \param file_name file name
518  * \param detector detector
519  */
520  inline void store(const JString& file_name, const JDetector& detector)
521  {
522  using namespace std;
523  using namespace JIO;
524 
525  if (file_name.endsWith(BINARY_DETECTOR_FILE_FORMAT)) {
526 
527  JFileStreamWriter out(file_name.c_str());
528 
529  detector.write(out);
530 
531  out.close();
532 
533  } else if (file_name.endsWith(KM3NET_DETECTOR_FILE_FORMAT)) {
534 
535  std::ofstream out(file_name.c_str());
536 
537  out << detector;
538 
539  out.close();
540 
541  } else if (file_name.endsWith(ZIPPED_DETECTOR_FILE_FORMAT)) {
542 
543  ogzstream out(file_name.c_str());
544 
545  out << detector;
546 
547  out.close();
548 
549  } else if (file_name.endsWith(GDML_DETECTOR_FILE_FORMAT)) {
550 
551  std::ofstream out(file_name.c_str());
552 
553  write_gdml(out, detector);
554 
555  out.close();
556 
557  } else {
558 
559  throw JProtocolException("Protocol not defined.");
560  }
561  }
562 
563 
564  /**
565  * Get module according module address map.
566  *
567  * \param memo module address map
568  * \param id module identifier
569  * \param location module location
570  * \return module
571  */
572  inline const JModule& getModule(const JModuleAddressMap& memo,
573  const int id,
574  const JModuleLocation& location = JModuleLocation())
575  {
576  static JModule module;
577 
578 
579  module.setID(id);
580 
581  module.setLocation(location);
582 
583  module.resize(memo.size());
584 
585  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))); }
586 
587  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))); }
588  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))); }
589  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))); }
590  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))); }
591  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))); }
592  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))); }
593 
594  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))); }
595  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))); }
596  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))); }
597  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))); }
598  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))); }
599  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))); }
600 
601  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))); }
602  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))); }
603  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))); }
604  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))); }
605  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))); }
606  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))); }
607 
608  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))); }
609  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))); }
610  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))); }
611  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))); }
612  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))); }
613  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))); }
614 
615  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))); }
616  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))); }
617  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))); }
618  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))); }
619  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))); }
620  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))); }
621 
622  return module;
623  }
624 
625 
626  /**
627  * Get module corresponding to Antares storey.
628  *
629  * \param id module identifier
630  * \param location module location
631  * \return module
632  */
633  inline const JModule& getModule(const int id,
634  const JModuleLocation& location = JModuleLocation())
635  {
636  using JTOOLS::PI;
637 
638  static JModule module;
639 
640  module.setID(id);
641 
642  module.setLocation(location);
643 
644  if (module.empty()) {
645 
646  module.resize(3);
647 
648  const double R = 0.5; // [m]
649 
650  const double st = sin(0.75*PI);
651  const double ct = cos(0.75*PI);
652 
653  for (int i = 0; i != 3; ++i) {
654 
655  const double phi = (2.0*PI*i) / 3.0;
656  const double cp = cos(phi);
657  const double sp = sin(phi);
658 
659  module[i] = JPMT(i, JAxis3D(JVector3D(R*st*cp, R*st*sp, R*ct), JVersor3D(st*cp, st*sp, ct)));
660  }
661  }
662 
663  return module;
664  }
665 
666 
667  /**
668  * Get module label (DU-floor) for JMonitor applications
669  *
670  * \param location module location
671  */
672  inline std::string getModuleLabel(const JModuleLocation& location) {
673  using namespace std;
674  ostringstream os;
675  os << "DU" << setfill('0') << setw(3) << location.getString();
676  os << "F" << setfill('0') << setw(2) << location.getFloor();
677  return os.str();
678  }
679 
680 
681 
682 
683 
684 }
685 
686 #endif
Exception for opening of file.
Definition: JException.hh:324
static const JRange< double, std::less< double > > DEFAULT_RANGE
Default range.
Definition: JRange.hh:506
double GetYrotationG4(const JVersor3D dir)
Get rotation over Y axis in Geant4 coordinate system.
General exception.
Definition: JException.hh:40
Exceptions.
Wrapper class around STL string class.
Definition: JString.hh:28
Data structure for a composite optical module.
Definition: JModule.hh:47
Logical location of module.
int getNumberOfStrings(const JDetector &detector)
Get number of strings.
double getIndexOfRefraction()
Get average index of refraction of water.
Definition: JConstants.hh:111
const JCalibration & getCalibration() const
Get calibration.
Detector data structure.
Definition: JDetector.hh:77
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:633
std::set< int > getStringIDs(const JDetector &detector)
Get list of strings IDs.
virtual JWriter & write(JWriter &out) const
Write to output.
Definition: JDetector.hh:397
static const double PI
Constants.
Definition: JConstants.hh:20
int getNumberOfPMTs(const JModule &module)
Get number of PMTs.
virtual JReader & read(JReader &in)
Read from input.
Definition: JDetector.hh:367
Data structure for PMT calibration.
Data structure for detector geometry and calibration.
Data structure for position fit.
Definition: JPoint3D.hh:22
bool endsWith(const std::string &suffix) const
Test if this string ends with the specified suffix.
Definition: JString.hh:175
JRange< double > JTimeRange
Type definition for time range.
JTimeRange getTimeRange(const Evt &event)
Get time range (i.e.
Axis object.
Definition: JAxis3D.hh:37
Monte Carlo detector (i.e.
double getSinThetaC()
Get average sine of Cherenkov angle of water.
Definition: JConstants.hh:144
const double getInverseSpeedOfLight()
Get inverse speed of light.
Definition: JConstants.hh:100
int getFloor() const
Get floor number.
int getString() const
Get string number.
Lookup table for PMT addresses in optical module.
double putTime(const T &t1, const JCalibration &cal)
Get de-calibrated time.
Constants.
JVector3D & sub(const JVector3D &vector)
Subtract vector.
Definition: JVector3D.hh:156
double getMaximalDistance(const JDetector &detector)
Get maximal distance between modules in detector.
Cylinder object.
Definition: JCylinder3D.hh:37
Data structure for vector in three dimensions.
Definition: JVector3D.hh:32
double getDY() const
Get y direction.
Definition: JVersor3D.hh:101
double getDX() const
Get x direction.
Definition: JVersor3D.hh:90
void setLocation(const JModuleLocation &location)
Set location.
Data structure for PMT geometry and calibration.
Definition: JPMT.hh:52
int getNumberOfFloors(const JDetector &detector)
Get number of floors.
double getY() const
Get y position.
Definition: JVector3D.hh:102
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
double getMaximalTime(const JDetector &detector)
Get maximal time between modules in detector following causality.
Protocol exception.
Definition: JException.hh:612
static const char *const KM3NET_DETECTOR_FILE_FORMAT
KM3NeT standard ASCII format.
std::string getModuleLabel(const JModuleLocation &location)
Get module label (DU-floor) for JMonitor applications.
static const char *const ZIPPED_DETECTOR_FILE_FORMAT
zipped KM3NeT standard ASCII format
Linear fit of JFIT::JPoint3D.
static const char *const GENDET_DETECTOR_FILE_FORMAT
File name extensions.
static const char *const BINARY_DETECTOR_FILE_FORMAT
JIO binary file format.
void setID(const int id)
Set identifier.
Definition: JObjectID.hh:65
double getX() const
Get x position.
Definition: JVector3D.hh:92
void store(const JString &file_name, const JDetector &detector)
Store detector to output file.
Data structure for normalised vector in three dimensions.
Definition: JVersor3D.hh:23
double getTanThetaC()
Get average tangent of Cherenkov angle of water.
Definition: JConstants.hh:122
void read_gdml(std::istream &, const JDetector &)
Binary buffered file output.
int getNumberOfModules(const JDetector &detector)
Get number of modules.
Linear fit of crossing point (position) between axes (objects with position and direction).
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:113
Binary buffered file input.
bool has(const int index) const
Test whether index is available.
double getDZ() const
Get z direction.
Definition: JVersor3D.hh:112
const JModule & getModule(const JDetector &detector, const JModuleLocation &location)
find module with a given string and floor number
Logical location of module.
void write_gdml(std::ostream &out, const JDetector &detector)
Writes KM3Sim GDML input file from detector.