Jpp 19.3.0-rc.2
the software that should make you happy
Loading...
Searching...
No Matches
JDetectorToolkit.hh
Go to the documentation of this file.
1#ifndef __JDETECTOR__JDETECTORTOOLKIT__
2#define __JDETECTOR__JDETECTORTOOLKIT__
3
4#include <limits>
5#include <istream>
6#include <ostream>
7#include <iostream>
8#include <fstream>
9#include <cmath>
10#include <vector>
11#include <set>
12#include <algorithm>
13#include <cstdlib>
14
26#include "JMath/JConstants.hh"
27#include "JMath/JMathToolkit.hh"
28#include "JTools/JRange.hh"
29#include "JIO/JFileStreamIO.hh"
30#include "Jeep/JeepToolkit.hh"
31#include "JLang/JException.hh"
32#include "JLang/gzstream.h"
33#include "JLang/JManip.hh"
34#include "JLang/Jpp.hh"
35
36
37/**
38 * \author mdejong
39 */
40
41namespace JDETECTOR {}
42namespace JPP { using namespace JDETECTOR; }
43
44/**
45 * Auxiliary classes and methods for detector calibration and simulation.
46 */
47namespace JDETECTOR {
48
54
55
56 /**
57 * File name extensions.
58 */
59 static const char* const GENDET_DETECTOR_FILE_FORMAT = "det"; //!< file format used by gendet
60 static const char* const BINARY_DETECTOR_FILE_FORMAT[] = { "dat", "datx" }; //!< JIO binary file format
61 static const char* const KM3NET_DETECTOR_FILE_FORMAT = "detx"; //!< %KM3NeT standard ASCII format
62 static const char* const ZIPPED_DETECTOR_FILE_FORMAT = "gz"; //!< zipped %KM3NeT standard ASCII format
63 static const char* const GDML_DETECTOR_FILE_FORMAT = "gdml"; //!< KM3Sim input format
64
65 static const char* const GDML_SCHEMA = getenv("GDML_SCHEMA_DIR"); //!< directory necessary for correct GDML header output
66 static const char* const CAN_MARGIN_M = getenv("CAN_MARGIN_M"); //!< extension of the detector size to comply with the can definition
67 static const char* const G4GDML_DEFAULT_SCHEMALOCATION = "http://service-spi.web.cern.ch/service-spi/app/releases/GDML/schema/gdml.xsd";
68
69
70 /**
71 * Get maximal distance between modules in detector.
72 * The option can be used to include base modules, if any.
73 *
74 * \param detector detector
75 * \param option option
76 * \return maximal distance [m]
77 */
78 inline double getMaximalDistance(const JDetector& detector, const bool option = false)
79 {
80 using namespace JPP;
81
82 double dmax = 0.0;
83
84 for (JDetector::const_iterator i1 = detector.begin(); i1 != detector.end(); ++i1) {
85 for (JDetector::const_iterator i2 = detector.begin(); i2 != i1; ++i2) {
86
87 if (option || (i1->getFloor() != 0 && i2->getFloor() != 0)) {
88
89 const double ds = getDistance(i1->getPosition(), i2->getPosition());
90
91 if (ds > dmax) {
92 dmax = ds;
93 }
94 }
95 }
96 }
97
98 return dmax;
99 }
100
101
102 /**
103 * Get maximal time between optical modules in detector following causality.
104 *
105 * \param detector detector
106 * \return maximal time [ns]
107 */
108 inline double getMaximalTime(const JDetector& detector)
109 {
110 using namespace JPP;
111
112 return getMaximalDistance(detector) * getIndexOfRefraction() * getInverseSpeedOfLight();
113 }
114
115
116 /**
117 * Get maximal time between optical modules in detector following causality.
118 * The road width corresponds to the maximal distance traveled by the light.
119 *
120 * \param detector detector
121 * \param roadWidth_m road width [m]
122 * \return maximal time [ns]
123 */
124 inline double getMaximalTime(const JDetector& detector, const double roadWidth_m)
125 {
126 using namespace JPP;
127
128 const double Dmax_m = getMaximalDistance(detector);
129
130 return (sqrt((Dmax_m + roadWidth_m*getSinThetaC()) *
131 (Dmax_m - roadWidth_m*getSinThetaC())) +
132 roadWidth_m * getSinThetaC() * getTanThetaC()) * getInverseSpeedOfLight();
133 }
134
135
136 /**
137 * Get de-calibrated time range.
138 *
139 * The de-calibrated time range is corrected for minimal and maximal time offset of PMTs in given module.
140 *
141 * \param timeRange time range [ns]
142 * \param module module
143 * \return time range [ns]
144 */
145 inline JTimeRange getTimeRange(const JTimeRange& timeRange, const JModule& module)
146 {
147 if (!module.empty()) {
148
150
151 for (JModule::const_iterator pmt = module.begin(); pmt != module.end(); ++pmt) {
152
153 const JCalibration& calibration = pmt->getCalibration();
154
155 time_range.include(putTime(timeRange.getLowerLimit(), calibration));
156 time_range.include(putTime(timeRange.getUpperLimit(), calibration));
157 }
158
159 return time_range;
160
161 } else {
162
163 return timeRange;
164 }
165 }
166
167
168 /**
169 * Get number of PMTs.
170 *
171 * \param module module
172 * \return number of PMTs
173 */
174 inline int getNumberOfPMTs(const JModule& module)
175 {
176 return module.size();
177 }
178
179
180 /**
181 * Get number of PMTs.
182 *
183 * \param detector detector
184 * \return number of PMTs
185 */
187 {
188 int number_of_pmts = 0;
189
190 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
191 number_of_pmts += getNumberOfPMTs(*module);
192 }
193
194 return number_of_pmts;
195 }
196
197
198 /**
199 * Get list of strings identifiers.
200 *
201 * \param detector detector
202 * \return list of string identifiers
203 */
205 {
206 std::set<int> buffer;
207
208 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
209 buffer.insert(module->getString());
210 }
211
212 return buffer;
213 }
214
215
216 /**
217 * Get list of modules identifiers.
218 *
219 * \param detector detector
220 * \param option option to include base modules
221 * \return list of module identifiers
222 */
223 inline std::set<int> getModuleIDs(const JDetector& detector, const bool option = false)
224 {
225 std::set<int> buffer;
226
227 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
228 if (option || !module->empty()) {
229 buffer.insert(module->getID());
230 }
231 }
232
233 return buffer;
234 }
235
236
237 /**
238 * Get number of floors.
239 *
240 * \param detector detector
241 * \return number of floors
242 */
244 {
245 std::set<int> buffer;
246
247 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
248 buffer.insert(module->getFloor());
249 }
250
251 return buffer.size();
252 }
253
254
255 /**
256 * Type definition for range of floors.
257 */
259
260
261 /**
262 * Get range of floors.
263 *
264 * \param detector detector
265 * \return range of floors
266 */
268 {
270
271 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
272 buffer.include(module->getFloor());
273 }
274
275 return buffer;
276 }
277
278
279 /**
280 * Get number of modules.
281 *
282 * \param detector detector
283 * \param option option to include base modules
284 * \return number of modules
285 */
286 inline int getNumberOfModules(const JDetector& detector, const bool option = false)
287 {
288 std::set<int> buffer;
289
290 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
291 if (option || !module->empty()) {
292 buffer.insert(module->getID());
293 }
294 }
295
296 return buffer.size();
297 }
298
299
300 /**
301 * Get rotation over X axis in Geant4 coordinate system
302 *
303 * \param dir direction
304 * \return X-rotation [deg]
305 */
306 inline double GetXrotationG4(const JVersor3D& dir)
307 {
308 using namespace JPP;
309
310 const double phi = atan2(dir.getDY(), dir.getDZ())*(180.0/PI);
311
312 if (phi < 0.0){
313 return phi + 360.0;
314 }
315 else{
316 return phi;
317 }
318 }
319
320
321 /**
322 * Get rotation over Y axis in Geant4 coordinate system
323 *
324 * \param dir direction
325 * \return Y-rotation [deg]
326 */
327 inline double GetYrotationG4(const JVersor3D& dir)
328 {
329 using namespace JPP;
330
331 return asin(-dir.getDX())*(180.0/PI);
332 }
333
334
335 inline void read_gdml(std::istream&, JDetector&)
336 {
337 THROW(JException, "Operation not defined.");
338 }
339
340
341 /**
342 * Writes KM3Sim GDML input file from detector
343 *
344 * \param out output stream
345 * \param detector detector
346 */
347 inline void write_gdml(std::ostream& out, const JDetector& detector)
348 {
349 using namespace std;
350 using namespace JPP;
351
352 const double DEFAULT_CAN_MARGIN_M = 350.0; // default can margin [m]
353 const double DEFAULT_WORLD_MARGIN_M = 50.0; // default world margin [m]
354
355 const JCylinder3D cylinder(detector.begin(), detector.end());
356
357 double can_margin;
358
359 if (CAN_MARGIN_M) {
360 can_margin = atof(CAN_MARGIN_M);
361 } else {
362 cerr << "CAN_MARGIN_M not defined! Setting standard CAN_MARGIN_M " << DEFAULT_CAN_MARGIN_M << " m." << endl;
363 can_margin = DEFAULT_CAN_MARGIN_M; //this is given in meters like all the dimensions in the GDML file (look at the unit field of every position and solid definition)
364 }
365
366 const double WorldCylinderHeight = 2*(cylinder.getZmax() - cylinder.getZmin() + can_margin + DEFAULT_WORLD_MARGIN_M);
367 const double WorldRadius = cylinder.getLength() + cylinder.getRadius() + can_margin + DEFAULT_WORLD_MARGIN_M;
368
369 const double Crust_Z_size = WorldCylinderHeight/2;
370 const double Crust_Z_position = -WorldCylinderHeight/4;
371
372 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=\"";
373 if (!GDML_SCHEMA) {
374 cerr << "GDML_SCHEMA_DIR NOT DEFINED! Setting default path." << endl;
375 out << G4GDML_DEFAULT_SCHEMALOCATION << "\">\n\n\n";
376 } else {
377 out << GDML_SCHEMA << "gdml.xsd\">\n\n\n";
378 }
379 out << "<!-- Jpp version: "<< getGITVersion() << " -->\n";
380 out << "<define>\n";
381 out << "<rotation name=\"identity\"/>\n<position name=\"zero\"/>\n";
382
383 int PMTs_NO = 0;
384
385 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
386
387 if(module->empty()) continue;
388
389 const JVector3D center = module->getCenter();
390
391 out << "<position name=\"PosString" << module->getString() << "_Dom" << module->getID() << "\" unit=\"m\" x=\"" << module->getX() << "\" y=\"" << module->getY() << "\" z=\"" << module->getZ() << "\"/>\n";
392
393 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
394
395 const JVector3D vec = static_cast<JVector3D>(*pmt).sub(center);
396 out << "<position name=\"CathodPosition" << pmt->getID() << "_" << module->getID() << "\" unit=\"m\" x=\"" << vec.getX() << "\" y=\"" << vec.getY() << "\" z=\"" << vec.getZ() << "\"/>\n";
397 out << "<rotation name=\"CathodRotation" << pmt->getID() << "_" << module->getID() << "\" unit=\"deg\" x=\"" << GetXrotationG4(*pmt) << "\" y=\"" << GetYrotationG4(*pmt) << "\" z=\"0.000000\"/>\n";
398 out << "<constant name=\"CathodID_" << PMTs_NO << "\" value=\"" << pmt->getID() << "\"/>\n";
399 PMTs_NO++;
400 }
401
402 }
403
404 out << "<position name=\"OMShift\" unit=\"m\" x=\"0\" y=\"0\" z=\"0.0392\"/>\n";
405 out << "\n\n\n";
406 out << "<!-- end of DU position definitions -->\n<position name=\"CrustPosition\" unit=\"m\" x=\"0\" y=\"0\" z=\"" << Crust_Z_position << "\"/>\n";
407
408 out << "</define>\n";
409 out << "<materials>\n";
410 out << "</materials>\n";
411 out << "<solids>\n";
412 out << " <box name=\"CrustBox\" lunit=\"m\" x=\"2200\" y=\"2200\" z=\"" << Crust_Z_size << "\"/>\n";
413 out << " <box name=\"StoreyBox\" lunit=\"m\" x=\"0.6\" y=\"0.6\" z=\"0.6\"/>\n";
414 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";
415 out << " <tube name=\"CathodTube\" lunit=\"cm\" aunit=\"deg\" rmin=\"0.0\" rmax=\"4.7462\" z=\"0.5\" startphi=\"0.0\" deltaphi=\"360.0\"/>\n";
416 out << " <tube name=\"WorldTube\" lunit=\"m\" aunit=\"deg\" rmin=\"0.0\" rmax=\"" << WorldRadius << "\" z=\"" << WorldCylinderHeight << "\" startphi=\"0.0\" deltaphi=\"360.0\"/>\n";
417 out << "</solids>\n\n\n";
418
419 out << "<structure>\n";
420 out << " <volume name=\"CathodVolume\">\n";
421 out << " <materialref ref=\"Cathod\"/>\n";
422 out << " <solidref ref=\"CathodTube\"/>\n";
423 out << " </volume>\n";
424
425 out << "<!-- OMVolume(s) construction -->\n";
426
427 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
428
429 if(module->empty()) continue;
430 out << " <volume name=\"OMVolume" << module->getID() << "\">\n";
431 out << " <materialref ref=\"Water\"/>\n";
432 out << " <solidref ref=\"OMSphere\"/>\n";
433
434 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
435 out << " <physvol>\n";
436 out << " <volumeref ref=\"CathodVolume\"/>\n";
437 out << " <positionref ref=\"CathodPosition" << pmt->getID() << "_" << module->getID() << "\"/>\n";
438 out << " <rotationref ref=\"CathodRotation" << pmt->getID() << "_" << module->getID() << "\"/>\n";
439 out << " </physvol>\n";
440 }
441
442 out << " </volume>\n";
443 }
444
445 out << "<!-- StoreyVolume(s) construction -->\n";
446
447 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
448 if(module->empty()) continue;
449 out << " <volume name=\"StoreyVolume" << module->getID() << "\">\n";
450 out << " <materialref ref=\"Water\"/>\n";
451 out << " <solidref ref=\"StoreyBox\"/>\n";
452 out << " <physvol>\n";
453 out << " <volumeref ref=\"OMVolume" << module->getID() << "\"/>\n";
454 out << " <positionref ref=\"OMShift\"/>\n";
455 out << " <rotationref ref=\"identity\"/>\n";
456 out << " </physvol>\n";
457 out << " </volume>\n";
458 }
459
460 out << "<!-- Crust Volume construction-->\n";
461 out << "<volume name=\"CrustVolume\">\n";
462 out << " <materialref ref=\"Crust\"/>\n";
463 out << " <solidref ref=\"CrustBox\"/>\n";
464 out << "</volume>\n";
465
466 out << "<!-- World Volume construction-->\n";
467 out << "<volume name=\"WorldVolume\">\n";
468 out << " <materialref ref=\"Water\"/>\n";
469 out << " <solidref ref=\"WorldTube\"/>\n";
470
471 out << " <physvol>\n";
472 out << " <volumeref ref=\"CrustVolume\"/>\n";
473 out << " <positionref ref=\"CrustPosition\"/>\n";
474 out << " <rotationref ref=\"identity\"/>\n";
475 out << " </physvol>\n";
476
477 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
478 if(module->empty()) continue;
479 out << " <physvol>\n";
480 out << " <volumeref ref=\"StoreyVolume" << module->getID() << "\"/>\n";
481 out << " <positionref ref=\"PosString" << module->getString() << "_Dom" << module->getID() << "\"/>\n";
482 out << " <rotationref ref=\"identity\"/>\n";
483 out << " </physvol>\n";
484 }
485
486 out << "</volume>\n";
487
488 out << "</structure>\n";
489 out << "<setup name=\"Default\" version=\"1.0\">\n";
490 out << "<world ref=\"WorldVolume\"/>\n";
491 out << "</setup>\n";
492 out << "</gdml>\n";
493 }
494
495
496 /**
497 * Load detector from input file.
498 *
499 * Supported file formats:
500 * - ASCII file, extension JDETECTOR::GENDET_DETECTOR_FILE_FORMAT, gendet format;
501 * - binary file, extension JDETECTOR::BINARY_DETECTOR_FILE_FORMAT, Jpp internal format;
502 * - ASCII file, extension JDETECTOR::KM3NET_DETECTOR_FILE_FORMAT, %KM3NeT standard format;
503 * - gzipped file, extension JDETECTOR::ZIPPED_DETECTOR_FILE_FORMAT, %KM3NeT standard format;
504 *
505 * \param file_name file name
506 * \param detector detector
507 */
508 inline void load(const std::string& file_name, JDetector& detector)
509 {
510 using namespace std;
511 using namespace JPP;
512
513 if (getFilenameExtension(file_name) == GENDET_DETECTOR_FILE_FORMAT) {
514
515 JMonteCarloDetector buffer(true);
516
517 ifstream in(file_name.c_str());
518
519 if (!in) {
520 THROW(JFileOpenException, "File not opened for reading: " << file_name);
521 }
522
523 in >> buffer;
524
525 in.close();
526
527 detector.swap(buffer);
528
529 } else if (getFilenameExtension(file_name) == BINARY_DETECTOR_FILE_FORMAT[0] ||
530 getFilenameExtension(file_name) == BINARY_DETECTOR_FILE_FORMAT[1]) {
531
532 JFileStreamReader in(file_name.c_str());
533
534 if (!in) {
535 THROW(JFileOpenException, "File not opened for reading: " << file_name);
536 }
537
538 try {
539
540 detector.read(in);
541
542 } catch(const exception& error) {
543
544 // recovery procedure for old format of comments
545
547
548 in.clear();
549 in.rewind();
550
551 detector.read(in);
552
554 }
555
556 in.close();
557
558 } else if (getFilenameExtension(file_name) == KM3NET_DETECTOR_FILE_FORMAT) {
559
560 ifstream in(file_name.c_str());
561
562 if (!in) {
563 THROW(JFileOpenException, "File not opened for reading: " << file_name);
564 }
565
566 in >> detector;
567
568 in.close();
569
570 } else if (getFilenameExtension(file_name) == ZIPPED_DETECTOR_FILE_FORMAT) {
571
572 igzstream in(file_name.c_str());
573
574 if (!in) {
575 THROW(JFileOpenException, "File not opened for reading: " << file_name);
576 }
577
578 in >> detector;
579
580 in.close();
581
582 } else {
583
584 THROW(JProtocolException, "Protocol not defined: " << file_name);
585 }
586
587 if (!detector.hasVersion()) {
588 THROW(JValueOutOfRange, "Invalid version <" << detector.getVersion() << ">");
589 }
590 }
591
592
593 /**
594 * Store detector to output file.
595 *
596 * Supported file formats:
597 * - binary file, extension JDETECTOR::BINARY_DETECTOR_FILE_FORMAT, Jpp internal format;
598 * - ASCII file, extension JDETECTOR::KM3NET_DETECTOR_FILE_FORMAT, %KM3NeT standard format;
599 * - gzipped file, extension JDETECTOR::ZIPPED_DETECTOR_FILE_FORMAT, %KM3NeT standard format;
600 * - gdml file, extension JDETECTOR::GDML_DETECTOR_FILE_FORMAT, KM3Sim input format;
601 *
602 * \param file_name file name
603 * \param detector detector
604 */
605 inline void store(const std::string& file_name, const JDetector& detector)
606 {
607 using namespace std;
608 using namespace JPP;
609
610 if (getFilenameExtension(file_name) == BINARY_DETECTOR_FILE_FORMAT[1]) {
611
612 JFileStreamWriter out(file_name.c_str());
613
614 if (!out) {
615 THROW(JFileOpenException, "File not opened for writing: " << file_name);
616 }
617
618 detector.write(out);
619
620 out.close();
621
622 } else if (getFilenameExtension(file_name) == KM3NET_DETECTOR_FILE_FORMAT) {
623
624 std::ofstream out(file_name.c_str());
625
626 if (!out) {
627 THROW(JFileOpenException, "File not opened for writing: " << file_name);
628 }
629
630 out << detector;
631
632 out.close();
633
634 } else if (getFilenameExtension(file_name) == ZIPPED_DETECTOR_FILE_FORMAT) {
635
636 ogzstream out(file_name.c_str());
637
638 if (!out) {
639 THROW(JFileOpenException, "File not opened for writing: " << file_name);
640 }
641
642 out << detector;
643
644 out.close();
645
646 } else if (getFilenameExtension(file_name) == GDML_DETECTOR_FILE_FORMAT) {
647
648 std::ofstream out(file_name.c_str());
649
650 if (!out) {
651 THROW(JFileOpenException, "File not opened for writing: " << file_name);
652 }
653
654 write_gdml(out, detector);
655
656 out.close();
657
658 } else {
659
660 THROW(JProtocolException, "Protocol not defined: " << file_name);
661 }
662 }
663
664
665 /**
666 * Auxiliary class to get rotation matrix between two optical modules.
667 */
668 struct JRotation :
669 public JRotation3D
670 {
671
672 static const size_t NUMBER_OF_DIMENSIONS = 3; //!< Number of dimensions
673
674
675 /**
676 * Get rotation matrix to go from first to second module.
677 *
678 * \param first first module
679 * \param second second module
680 * \return rotation matrix
681 */
682 const JRotation3D& operator()(const JModule& first, const JModule& second)
683 {
684 this->setIdentity();
685
686 if (first.size() == second.size()) {
687
688 const size_t N = first.size();
689
690 if (N >= NUMBER_OF_DIMENSIONS) {
691
692 in .resize(N);
693 out.resize(N);
694
695 for (size_t i = 0; i != N; ++i) {
696 in [i] = first .getPMT(i).getDirection();
697 out[i] = second.getPMT(i).getDirection();
698 }
699
700 for (size_t i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
701 if (!orthonormalise(i)) {
702 THROW(JException, "Failure to orthonormalise direction " << i);
703 }
704 }
705
706 this->a00 = out[0].getX() * in[0].getX() + out[1].getX() * in[1].getX() + out[2].getX() * in[2].getX();
707 this->a01 = out[0].getX() * in[0].getY() + out[1].getX() * in[1].getY() + out[2].getX() * in[2].getY();
708 this->a02 = out[0].getX() * in[0].getZ() + out[1].getX() * in[1].getZ() + out[2].getX() * in[2].getZ();
709
710 this->a10 = out[0].getY() * in[0].getX() + out[1].getY() * in[1].getX() + out[2].getY() * in[2].getX();
711 this->a11 = out[0].getY() * in[0].getY() + out[1].getY() * in[1].getY() + out[2].getY() * in[2].getY();
712 this->a12 = out[0].getY() * in[0].getZ() + out[1].getY() * in[1].getZ() + out[2].getY() * in[2].getZ();
713
714 this->a20 = out[0].getZ() * in[0].getX() + out[1].getZ() * in[1].getX() + out[2].getZ() * in[2].getX();
715 this->a21 = out[0].getZ() * in[0].getY() + out[1].getZ() * in[1].getY() + out[2].getZ() * in[2].getY();
716 this->a22 = out[0].getZ() * in[0].getZ() + out[1].getZ() * in[1].getZ() + out[2].getZ() * in[2].getZ();
717
718 } else {
719
720 THROW(JException, "Module " << first.getID() << " size " << N << " < " << NUMBER_OF_DIMENSIONS);
721 }
722
723 } else {
724
725 THROW(JException, "Module " << first.getID() << " size " << first.size() << " != " << second.size());
726 }
727
728 return *this;
729 }
730
731 private:
732 /**
733 * Put normalised primary direction at specified index and orthoganilise following directions.\n
734 * This procedure follows Gram-Schmidt process.
735 *
736 * \param index index
737 * \param precision precision
738 * \return true if primary direction exists; else false
739 */
740 bool orthonormalise(const size_t index, const double precision = std::numeric_limits<double>::epsilon())
741 {
742 using namespace std;
743
744 size_t pos = index;
745
746 for (size_t i = index + 1; i != in.size(); ++i) {
747 if (in[i].getLengthSquared() > in[pos].getLengthSquared()) {
748 pos = i;
749 }
750 }
751
752 const double u = in[pos].getLength();
753
754 if (u > precision) {
755
756 in [pos] /= u;
757 out[pos] /= u;
758
759 if (pos != index) {
760 swap(in [pos], in [index]);
761 swap(out[pos], out[index]);
762 }
763
764 for (size_t i = index + 1; i != in.size(); ++i) {
765
766 const double dot = in[index].getDot(in[i]);
767
768 in [i] -= dot * in [index];
769 out[i] -= dot * out[index];
770 }
771
772 return true;
773
774 } else {
775
776 return false;
777 }
778 }
779
780
783 };
784
785
786 /**
787 * Function object to get rotation matrix to go from first to second module.
788 */
790
791
792 /**
793 * Get position to go from first to second module.
794 *
795 * \param first first module
796 * \param second second module
797 * \return position
798 */
799 inline JPosition3D getPosition(const JModule& first, const JModule& second)
800 {
801 return second.getPosition() - first.getPosition();
802 }
803
804
805 /**
806 * Get calibration to go from first to second calibration.
807 *
808 * \param first first calibration
809 * \param second second calibration
810 * \return calibration
811 */
812 inline JCalibration getCalibration(const JCalibration& first, const JCalibration& second)
813 {
814 return JCalibration(second.getT0() - first.getT0());
815 }
816}
817
818#endif
Data structure for detector geometry and calibration.
Exceptions.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Logical location of module.
I/O manipulators.
Auxiliary methods for geometrical methods.
Mathematical constants.
Physics constants.
Auxiliary class to define a range between two values.
Auxiliary methods for handling file names, type names and environment.
Jpp environment information.
Data structure for time calibration.
double getT0() const
Get time offset.
Detector data structure.
Definition JDetector.hh:96
static void setStartOfComment(const start_of_comment_type option)
Set option for short start of comment in binary I/O.
Definition JDetector.hh:589
Data structure for a composite optical module.
Definition JModule.hh:75
const JPMT & getPMT(const int index) const
Get PMT.
Definition JModule.hh:172
Monte Carlo detector (i.e. ".det" file).
double getRadius() const
Get radius.
Definition JCircle2D.hh:144
double getLength() const
Get length.
Definition JVector2D.hh:199
double getZmin() const
Get minimal z position.
double getZmax() const
Get maximal z position.
const JDirection3D & getDirection() const
Get direction.
Data structure for position in three dimensions.
const JPosition3D & getPosition() const
Get position.
Data structure for vector in three dimensions.
Definition JVector3D.hh:36
JVector3D & sub(const JVector3D &vector)
Subtract vector.
Definition JVector3D.hh:158
Data structure for normalised vector in three dimensions.
Definition JVersor3D.hh:28
double getDY() const
Get y direction.
Definition JVersor3D.hh:106
double getDX() const
Get x direction.
Definition JVersor3D.hh:95
double getDZ() const
Get z direction.
Definition JVersor3D.hh:117
Binary buffered file input.
virtual void clear() override
Clear status of reader.
Binary buffered file output.
void close()
Close file.
General exception.
Definition JException.hh:24
Exception for opening of file.
int getID() const
Get identifier.
Definition JObjectID.hh:50
Protocol exception.
Exception for accessing a value in a collection that is outside of its range.
JMatrix3D & setIdentity()
Set to identity matrix.
static JRange< double, std::less< double > > DEFAULT_RANGE()
Definition JRange.hh:555
range_type & include(argument_type x)
Include given value to range.
Definition JRange.hh:397
T getLowerLimit() const
Get lower limit.
Definition JRange.hh:202
T getUpperLimit() const
Get upper limit.
Definition JRange.hh:213
void close()
Definition gzstream.h:185
file Auxiliary data structures and methods for detector calibration.
Definition JAnchor.hh:12
int getNumberOfPMTs(const JModule &module)
Get number of PMTs.
static JRotation getRotation
Function object to get rotation matrix to go from first to second module.
int getNumberOfFloors(const JDetector &detector)
Get number of floors.
floor_range getRangeOfFloors(const JDetector &detector)
Get range of floors.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
static const char *const GDML_DETECTOR_FILE_FORMAT
KM3Sim input format.
JCalibration getCalibration(const JCalibration &first, const JCalibration &second)
Get calibration to go from first to second calibration.
void write_gdml(std::ostream &out, const JDetector &detector)
Writes KM3Sim GDML input file from detector.
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
static const char *const ZIPPED_DETECTOR_FILE_FORMAT
zipped KM3NeT standard ASCII format
std::set< int > getStringIDs(const JDetector &detector)
Get list of strings identifiers.
static const char *const G4GDML_DEFAULT_SCHEMALOCATION
JPosition3D getPosition(const JModule &first, const JModule &second)
Get position to go from first to second module.
JTOOLS::JRange< int > floor_range
Type definition for range of floors.
double getMaximalTime(const JDetector &detector)
Get maximal time between optical modules in detector following causality.
static const char *const CAN_MARGIN_M
extension of the detector size to comply with the can definition
static const char *const KM3NET_DETECTOR_FILE_FORMAT
KM3NeT standard ASCII format
static const char *const GDML_SCHEMA
directory necessary for correct GDML header output
static const char *const GENDET_DETECTOR_FILE_FORMAT
File name extensions.
double GetYrotationG4(const JVersor3D &dir)
Get rotation over Y axis in Geant4 coordinate system.
int getNumberOfModules(const JDetector &detector, const bool option=false)
Get number of modules.
double getMaximalDistance(const JDetector &detector, const bool option=false)
Get maximal distance between modules in detector.
double putTime(const T &t1, const JCalibration &cal)
Get de-calibrated time.
JTimeRange getTimeRange(const JTimeRange &timeRange, const JModule &module)
Get de-calibrated time range.
std::set< int > getModuleIDs(const JDetector &detector, const bool option=false)
Get list of modules identifiers.
double GetXrotationG4(const JVersor3D &dir)
Get rotation over X axis in Geant4 coordinate system.
void read_gdml(std::istream &, JDetector &)
static const char *const BINARY_DETECTOR_FILE_FORMAT[]
JIO binary file format.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Calibration.
Definition JHead.hh:330
Detector file.
Definition JHead.hh:227
Auxiliary class to get rotation matrix between two optical modules.
std::vector< JVector3D > out
const JRotation3D & operator()(const JModule &first, const JModule &second)
Get rotation matrix to go from first to second module.
std::vector< JVector3D > in
bool orthonormalise(const size_t index, const double precision=std::numeric_limits< double >::epsilon())
Put normalised primary direction at specified index and orthoganilise following directions.
static const size_t NUMBER_OF_DIMENSIONS
Number of dimensions.