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