Jpp master_rocky-44-g75b7c4f75
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 * \return number of modules
284 */
286 {
287 std::set<int> buffer;
288
289 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
290 buffer.insert(module->getID());
291 }
292
293 return buffer.size();
294 }
295
296
297 /**
298 * Get rotation over X axis in Geant4 coordinate system
299 *
300 * \param dir direction
301 * \return X-rotation [deg]
302 */
303 inline double GetXrotationG4(const JVersor3D& dir)
304 {
305 using namespace JPP;
306
307 const double phi = atan2(dir.getDY(), dir.getDZ())*(180.0/PI);
308
309 if (phi < 0.0){
310 return phi + 360.0;
311 }
312 else{
313 return phi;
314 }
315 }
316
317
318 /**
319 * Get rotation over Y axis in Geant4 coordinate system
320 *
321 * \param dir direction
322 * \return Y-rotation [deg]
323 */
324 inline double GetYrotationG4(const JVersor3D& dir)
325 {
326 using namespace JPP;
327
328 return asin(-dir.getDX())*(180.0/PI);
329 }
330
331
332 inline void read_gdml(std::istream&, JDetector&)
333 {
334 THROW(JException, "Operation not defined.");
335 }
336
337
338 /**
339 * Writes KM3Sim GDML input file from detector
340 *
341 * \param out output stream
342 * \param detector detector
343 */
344 inline void write_gdml(std::ostream& out, const JDetector& detector)
345 {
346 using namespace std;
347 using namespace JPP;
348
349 const double DEFAULT_CAN_MARGIN_M = 350.0; // default can margin [m]
350 const double DEFAULT_WORLD_MARGIN_M = 50.0; // default world margin [m]
351
352 const JCylinder3D cylinder(detector.begin(), detector.end());
353
354 double can_margin;
355
356 if (CAN_MARGIN_M) {
357 can_margin = atof(CAN_MARGIN_M);
358 } else {
359 cerr << "CAN_MARGIN_M not defined! Setting standard CAN_MARGIN_M " << DEFAULT_CAN_MARGIN_M << " m." << endl;
360 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)
361 }
362
363 const double WorldCylinderHeight = 2*(cylinder.getZmax() - cylinder.getZmin() + can_margin + DEFAULT_WORLD_MARGIN_M);
364 const double WorldRadius = cylinder.getLength() + cylinder.getRadius() + can_margin + DEFAULT_WORLD_MARGIN_M;
365
366 const double Crust_Z_size = WorldCylinderHeight/2;
367 const double Crust_Z_position = -WorldCylinderHeight/4;
368
369 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=\"";
370 if (!GDML_SCHEMA) {
371 cerr << "GDML_SCHEMA_DIR NOT DEFINED! Setting default path." << endl;
372 out << G4GDML_DEFAULT_SCHEMALOCATION << "\">\n\n\n";
373 } else {
374 out << GDML_SCHEMA << "gdml.xsd\">\n\n\n";
375 }
376 out << "<!-- Jpp version: "<< getGITVersion() << " -->\n";
377 out << "<define>\n";
378 out << "<rotation name=\"identity\"/>\n<position name=\"zero\"/>\n";
379
380 int PMTs_NO = 0;
381
382 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
383
384 if(module->empty()) continue;
385
386 const JVector3D center = module->getCenter();
387
388 out << "<position name=\"PosString" << module->getString() << "_Dom" << module->getID() << "\" unit=\"m\" x=\"" << module->getX() << "\" y=\"" << module->getY() << "\" z=\"" << module->getZ() << "\"/>\n";
389
390 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
391
392 const JVector3D vec = static_cast<JVector3D>(*pmt).sub(center);
393 out << "<position name=\"CathodPosition" << pmt->getID() << "_" << module->getID() << "\" unit=\"m\" x=\"" << vec.getX() << "\" y=\"" << vec.getY() << "\" z=\"" << vec.getZ() << "\"/>\n";
394 out << "<rotation name=\"CathodRotation" << pmt->getID() << "_" << module->getID() << "\" unit=\"deg\" x=\"" << GetXrotationG4(*pmt) << "\" y=\"" << GetYrotationG4(*pmt) << "\" z=\"0.000000\"/>\n";
395 out << "<constant name=\"CathodID_" << PMTs_NO << "\" value=\"" << pmt->getID() << "\"/>\n";
396 PMTs_NO++;
397 }
398
399 }
400
401 out << "<position name=\"OMShift\" unit=\"m\" x=\"0\" y=\"0\" z=\"0.0392\"/>\n";
402 out << "\n\n\n";
403 out << "<!-- end of DU position definitions -->\n<position name=\"CrustPosition\" unit=\"m\" x=\"0\" y=\"0\" z=\"" << Crust_Z_position << "\"/>\n";
404
405 out << "</define>\n";
406 out << "<materials>\n";
407 out << "</materials>\n";
408 out << "<solids>\n";
409 out << " <box name=\"CrustBox\" lunit=\"m\" x=\"2200\" y=\"2200\" z=\"" << Crust_Z_size << "\"/>\n";
410 out << " <box name=\"StoreyBox\" lunit=\"m\" x=\"0.6\" y=\"0.6\" z=\"0.6\"/>\n";
411 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";
412 out << " <tube name=\"CathodTube\" lunit=\"cm\" aunit=\"deg\" rmin=\"0.0\" rmax=\"4.7462\" z=\"0.5\" startphi=\"0.0\" deltaphi=\"360.0\"/>\n";
413 out << " <tube name=\"WorldTube\" lunit=\"m\" aunit=\"deg\" rmin=\"0.0\" rmax=\"" << WorldRadius << "\" z=\"" << WorldCylinderHeight << "\" startphi=\"0.0\" deltaphi=\"360.0\"/>\n";
414 out << "</solids>\n\n\n";
415
416 out << "<structure>\n";
417 out << " <volume name=\"CathodVolume\">\n";
418 out << " <materialref ref=\"Cathod\"/>\n";
419 out << " <solidref ref=\"CathodTube\"/>\n";
420 out << " </volume>\n";
421
422 out << "<!-- OMVolume(s) construction -->\n";
423
424 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
425
426 if(module->empty()) continue;
427 out << " <volume name=\"OMVolume" << module->getID() << "\">\n";
428 out << " <materialref ref=\"Water\"/>\n";
429 out << " <solidref ref=\"OMSphere\"/>\n";
430
431 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
432 out << " <physvol>\n";
433 out << " <volumeref ref=\"CathodVolume\"/>\n";
434 out << " <positionref ref=\"CathodPosition" << pmt->getID() << "_" << module->getID() << "\"/>\n";
435 out << " <rotationref ref=\"CathodRotation" << pmt->getID() << "_" << module->getID() << "\"/>\n";
436 out << " </physvol>\n";
437 }
438
439 out << " </volume>\n";
440 }
441
442 out << "<!-- StoreyVolume(s) construction -->\n";
443
444 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
445 if(module->empty()) continue;
446 out << " <volume name=\"StoreyVolume" << module->getID() << "\">\n";
447 out << " <materialref ref=\"Water\"/>\n";
448 out << " <solidref ref=\"StoreyBox\"/>\n";
449 out << " <physvol>\n";
450 out << " <volumeref ref=\"OMVolume" << module->getID() << "\"/>\n";
451 out << " <positionref ref=\"OMShift\"/>\n";
452 out << " <rotationref ref=\"identity\"/>\n";
453 out << " </physvol>\n";
454 out << " </volume>\n";
455 }
456
457 out << "<!-- Crust Volume construction-->\n";
458 out << "<volume name=\"CrustVolume\">\n";
459 out << " <materialref ref=\"Crust\"/>\n";
460 out << " <solidref ref=\"CrustBox\"/>\n";
461 out << "</volume>\n";
462
463 out << "<!-- World Volume construction-->\n";
464 out << "<volume name=\"WorldVolume\">\n";
465 out << " <materialref ref=\"Water\"/>\n";
466 out << " <solidref ref=\"WorldTube\"/>\n";
467
468 out << " <physvol>\n";
469 out << " <volumeref ref=\"CrustVolume\"/>\n";
470 out << " <positionref ref=\"CrustPosition\"/>\n";
471 out << " <rotationref ref=\"identity\"/>\n";
472 out << " </physvol>\n";
473
474 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
475 if(module->empty()) continue;
476 out << " <physvol>\n";
477 out << " <volumeref ref=\"StoreyVolume" << module->getID() << "\"/>\n";
478 out << " <positionref ref=\"PosString" << module->getString() << "_Dom" << module->getID() << "\"/>\n";
479 out << " <rotationref ref=\"identity\"/>\n";
480 out << " </physvol>\n";
481 }
482
483 out << "</volume>\n";
484
485 out << "</structure>\n";
486 out << "<setup name=\"Default\" version=\"1.0\">\n";
487 out << "<world ref=\"WorldVolume\"/>\n";
488 out << "</setup>\n";
489 out << "</gdml>\n";
490 }
491
492
493 /**
494 * Load detector from input file.
495 *
496 * Supported file formats:
497 * - ASCII file, extension JDETECTOR::GENDET_DETECTOR_FILE_FORMAT, gendet format;
498 * - binary file, extension JDETECTOR::BINARY_DETECTOR_FILE_FORMAT, Jpp internal format;
499 * - ASCII file, extension JDETECTOR::KM3NET_DETECTOR_FILE_FORMAT, %KM3NeT standard format;
500 * - gzipped file, extension JDETECTOR::ZIPPED_DETECTOR_FILE_FORMAT, %KM3NeT standard format;
501 *
502 * \param file_name file name
503 * \param detector detector
504 */
505 inline void load(const std::string& file_name, JDetector& detector)
506 {
507 using namespace std;
508 using namespace JPP;
509
510 if (getFilenameExtension(file_name) == GENDET_DETECTOR_FILE_FORMAT) {
511
512 JMonteCarloDetector buffer(true);
513
514 ifstream in(file_name.c_str());
515
516 if (!in) {
517 THROW(JFileOpenException, "File not opened for reading: " << file_name);
518 }
519
520 in >> buffer;
521
522 in.close();
523
524 detector.swap(buffer);
525
526 } else if (getFilenameExtension(file_name) == BINARY_DETECTOR_FILE_FORMAT[0] ||
527 getFilenameExtension(file_name) == BINARY_DETECTOR_FILE_FORMAT[1]) {
528
529 JFileStreamReader in(file_name.c_str());
530
531 if (!in) {
532 THROW(JFileOpenException, "File not opened for reading: " << file_name);
533 }
534
535 try {
536
537 detector.read(in);
538
539 } catch(const exception& error) {
540
541 // recovery procedure for old format of comments
542
544
545 in.clear();
546 in.rewind();
547
548 detector.read(in);
549
551 }
552
553 in.close();
554
555 } else if (getFilenameExtension(file_name) == KM3NET_DETECTOR_FILE_FORMAT) {
556
557 ifstream in(file_name.c_str());
558
559 if (!in) {
560 THROW(JFileOpenException, "File not opened for reading: " << file_name);
561 }
562
563 in >> detector;
564
565 in.close();
566
567 } else if (getFilenameExtension(file_name) == ZIPPED_DETECTOR_FILE_FORMAT) {
568
569 igzstream in(file_name.c_str());
570
571 if (!in) {
572 THROW(JFileOpenException, "File not opened for reading: " << file_name);
573 }
574
575 in >> detector;
576
577 in.close();
578
579 } else {
580
581 THROW(JProtocolException, "Protocol not defined: " << file_name);
582 }
583
584 if (!detector.hasVersion()) {
585 THROW(JValueOutOfRange, "Invalid version <" << detector.getVersion() << ">");
586 }
587 }
588
589
590 /**
591 * Store detector to output file.
592 *
593 * Supported file formats:
594 * - binary file, extension JDETECTOR::BINARY_DETECTOR_FILE_FORMAT, Jpp internal format;
595 * - ASCII file, extension JDETECTOR::KM3NET_DETECTOR_FILE_FORMAT, %KM3NeT standard format;
596 * - gzipped file, extension JDETECTOR::ZIPPED_DETECTOR_FILE_FORMAT, %KM3NeT standard format;
597 * - gdml file, extension JDETECTOR::GDML_DETECTOR_FILE_FORMAT, KM3Sim input format;
598 *
599 * \param file_name file name
600 * \param detector detector
601 */
602 inline void store(const std::string& file_name, const JDetector& detector)
603 {
604 using namespace std;
605 using namespace JPP;
606
607 if (getFilenameExtension(file_name) == BINARY_DETECTOR_FILE_FORMAT[1]) {
608
609 JFileStreamWriter out(file_name.c_str());
610
611 if (!out) {
612 THROW(JFileOpenException, "File not opened for writing: " << file_name);
613 }
614
615 detector.write(out);
616
617 out.close();
618
619 } else if (getFilenameExtension(file_name) == KM3NET_DETECTOR_FILE_FORMAT) {
620
621 std::ofstream out(file_name.c_str());
622
623 if (!out) {
624 THROW(JFileOpenException, "File not opened for writing: " << file_name);
625 }
626
627 out << detector;
628
629 out.close();
630
631 } else if (getFilenameExtension(file_name) == ZIPPED_DETECTOR_FILE_FORMAT) {
632
633 ogzstream out(file_name.c_str());
634
635 if (!out) {
636 THROW(JFileOpenException, "File not opened for writing: " << file_name);
637 }
638
639 out << detector;
640
641 out.close();
642
643 } else if (getFilenameExtension(file_name) == GDML_DETECTOR_FILE_FORMAT) {
644
645 std::ofstream out(file_name.c_str());
646
647 if (!out) {
648 THROW(JFileOpenException, "File not opened for writing: " << file_name);
649 }
650
651 write_gdml(out, detector);
652
653 out.close();
654
655 } else {
656
657 THROW(JProtocolException, "Protocol not defined: " << file_name);
658 }
659 }
660
661
662 /**
663 * Auxiliary class to get rotation matrix between two optical modules.
664 */
665 struct JRotation :
666 public JRotation3D
667 {
668
669 static const size_t NUMBER_OF_DIMENSIONS = 3; //!< Number of dimensions
670
671
672 /**
673 * Get rotation matrix to go from first to second module.
674 *
675 * \param first first module
676 * \param second second module
677 * \return rotation matrix
678 */
679 const JRotation3D& operator()(const JModule& first, const JModule& second)
680 {
681 this->setIdentity();
682
683 if (first.size() == second.size()) {
684
685 const size_t N = first.size();
686
687 if (N >= NUMBER_OF_DIMENSIONS) {
688
689 in .resize(N);
690 out.resize(N);
691
692 for (size_t i = 0; i != N; ++i) {
693 in [i] = first .getPMT(i).getDirection();
694 out[i] = second.getPMT(i).getDirection();
695 }
696
697 for (size_t i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
698 if (!orthonormalise(i)) {
699 THROW(JException, "Failure to orthonormalise direction " << i);
700 }
701 }
702
703 this->a00 = out[0].getX() * in[0].getX() + out[1].getX() * in[1].getX() + out[2].getX() * in[2].getX();
704 this->a01 = out[0].getX() * in[0].getY() + out[1].getX() * in[1].getY() + out[2].getX() * in[2].getY();
705 this->a02 = out[0].getX() * in[0].getZ() + out[1].getX() * in[1].getZ() + out[2].getX() * in[2].getZ();
706
707 this->a10 = out[0].getY() * in[0].getX() + out[1].getY() * in[1].getX() + out[2].getY() * in[2].getX();
708 this->a11 = out[0].getY() * in[0].getY() + out[1].getY() * in[1].getY() + out[2].getY() * in[2].getY();
709 this->a12 = out[0].getY() * in[0].getZ() + out[1].getY() * in[1].getZ() + out[2].getY() * in[2].getZ();
710
711 this->a20 = out[0].getZ() * in[0].getX() + out[1].getZ() * in[1].getX() + out[2].getZ() * in[2].getX();
712 this->a21 = out[0].getZ() * in[0].getY() + out[1].getZ() * in[1].getY() + out[2].getZ() * in[2].getY();
713 this->a22 = out[0].getZ() * in[0].getZ() + out[1].getZ() * in[1].getZ() + out[2].getZ() * in[2].getZ();
714
715 } else {
716
717 THROW(JException, "Module " << first.getID() << " size " << N << " < " << NUMBER_OF_DIMENSIONS);
718 }
719
720 } else {
721
722 THROW(JException, "Module " << first.getID() << " size " << first.size() << " != " << second.size());
723 }
724
725 return *this;
726 }
727
728 private:
729 /**
730 * Put normalised primary direction at specified index and orthoganilise following directions.\n
731 * This procedure follows Gram-Schmidt process.
732 *
733 * \param index index
734 * \param precision precision
735 * \return true if primary direction exists; else false
736 */
737 bool orthonormalise(const size_t index, const double precision = std::numeric_limits<double>::epsilon())
738 {
739 using namespace std;
740
741 size_t pos = index;
742
743 for (size_t i = index + 1; i != in.size(); ++i) {
744 if (in[i].getLengthSquared() > in[pos].getLengthSquared()) {
745 pos = i;
746 }
747 }
748
749 const double u = in[pos].getLength();
750
751 if (u > precision) {
752
753 in [pos] /= u;
754 out[pos] /= u;
755
756 if (pos != index) {
757 swap(in [pos], in [index]);
758 swap(out[pos], out[index]);
759 }
760
761 for (size_t i = index + 1; i != in.size(); ++i) {
762
763 const double dot = in[index].getDot(in[i]);
764
765 in [i] -= dot * in [index];
766 out[i] -= dot * out[index];
767 }
768
769 return true;
770
771 } else {
772
773 return false;
774 }
775 }
776
777
780 };
781
782
783 /**
784 * Function object to get rotation matrix to go from first to second module.
785 */
787
788
789 /**
790 * Get position to go from first to second module.
791 *
792 * \param first first module
793 * \param second second module
794 * \return position
795 */
796 inline JPosition3D getPosition(const JModule& first, const JModule& second)
797 {
798 return second.getPosition() - first.getPosition();
799 }
800
801
802 /**
803 * Get calibration to go from first to second calibration.
804 *
805 * \param first first calibration
806 * \param second second calibration
807 * \return calibration
808 */
809 inline JCalibration getCalibration(const JCalibration& first, const JCalibration& second)
810 {
811 return JCalibration(second.getT0() - first.getT0());
812 }
813}
814
815#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
int getNumberOfModules(const JDetector &detector)
Get number of modules.
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.
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.