1#ifndef __JDETECTOR__JDETECTORTOOLKIT__
2#define __JDETECTOR__JDETECTORTOOLKIT__
65 static const char*
const GDML_SCHEMA = getenv(
"GDML_SCHEMA_DIR");
84 for (JDetector::const_iterator i1 =
detector.begin(); i1 !=
detector.end(); ++i1) {
85 for (JDetector::const_iterator i2 =
detector.begin(); i2 != i1; ++i2) {
87 if (option || (i1->getFloor() != 0 && i2->getFloor() != 0)) {
89 const double ds = getDistance(i1->getPosition(), i2->getPosition());
130 return (sqrt((Dmax_m + roadWidth_m*getSinThetaC()) *
131 (Dmax_m - roadWidth_m*getSinThetaC())) +
132 roadWidth_m * getSinThetaC() * getTanThetaC()) * getInverseSpeedOfLight();
147 if (!module.empty()) {
151 for (JModule::const_iterator pmt = module.begin(); pmt != module.end(); ++pmt) {
176 return module.size();
188 int number_of_pmts = 0;
190 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
194 return number_of_pmts;
208 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
209 buffer.insert(module->getString());
227 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
228 if (option || !module->empty()) {
229 buffer.insert(module->getID());
247 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
248 buffer.insert(module->getFloor());
251 return buffer.size();
271 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
272 buffer.
include(module->getFloor());
290 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
291 if (option || !module->empty()) {
292 buffer.insert(module->getID());
296 return buffer.size();
310 const double phi = atan2(dir.
getDY(), dir.
getDZ())*(180.0/PI);
331 return asin(-dir.
getDX())*(180.0/PI);
352 const double DEFAULT_CAN_MARGIN_M = 350.0;
353 const double DEFAULT_WORLD_MARGIN_M = 50.0;
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;
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;
369 const double Crust_Z_size = WorldCylinderHeight/2;
370 const double Crust_Z_position = -WorldCylinderHeight/4;
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=\"";
374 cerr <<
"GDML_SCHEMA_DIR NOT DEFINED! Setting default path." << endl;
379 out <<
"<!-- Jpp version: "<< getGITVersion() <<
" -->\n";
381 out <<
"<rotation name=\"identity\"/>\n<position name=\"zero\"/>\n";
385 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
387 if(module->empty())
continue;
389 const JVector3D center =
module->getCenter();
391 out <<
"<position name=\"PosString" << module->getString() <<
"_Dom" <<
module->getID() << "\" unit=\"m\" x=\"" << module->getX() << "\" y=\"" << module->getY() << "\" z=\"" << module->getZ() << "\"/>\n";
393 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
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";
404 out <<
"<position name=\"OMShift\" unit=\"m\" x=\"0\" y=\"0\" z=\"0.0392\"/>\n";
406 out <<
"<!-- end of DU position definitions -->\n<position name=\"CrustPosition\" unit=\"m\" x=\"0\" y=\"0\" z=\"" << Crust_Z_position <<
"\"/>\n";
408 out <<
"</define>\n";
409 out <<
"<materials>\n";
410 out <<
"</materials>\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";
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";
425 out <<
"<!-- OMVolume(s) construction -->\n";
427 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
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";
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";
442 out <<
" </volume>\n";
445 out <<
"<!-- StoreyVolume(s) construction -->\n";
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";
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";
466 out <<
"<!-- World Volume construction-->\n";
467 out <<
"<volume name=\"WorldVolume\">\n";
468 out <<
" <materialref ref=\"Water\"/>\n";
469 out <<
" <solidref ref=\"WorldTube\"/>\n";
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";
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";
486 out <<
"</volume>\n";
488 out <<
"</structure>\n";
489 out <<
"<setup name=\"Default\" version=\"1.0\">\n";
490 out <<
"<world ref=\"WorldVolume\"/>\n";
517 ifstream in(file_name.c_str());
542 }
catch(
const exception& error) {
560 ifstream in(file_name.c_str());
624 std::ofstream out(file_name.c_str());
648 std::ofstream out(file_name.c_str());
686 if (first.size() == second.size()) {
688 const size_t N = first.size();
695 for (
size_t i = 0; i != N; ++i) {
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();
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();
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();
725 THROW(
JException,
"Module " << first.
getID() <<
" size " << first.size() <<
" != " << second.size());
740 bool orthonormalise(
const size_t index,
const double precision = std::numeric_limits<double>::epsilon())
746 for (
size_t i = index + 1; i !=
in.size(); ++i) {
747 if (
in[i].getLengthSquared() >
in[pos].getLengthSquared()) {
752 const double u =
in[pos].getLength();
760 swap(
in [pos],
in [index]);
761 swap(
out[pos],
out[index]);
764 for (
size_t i = index + 1; i !=
in.size(); ++i) {
766 const double dot =
in[index].getDot(
in[i]);
768 in [i] -= dot *
in [index];
769 out[i] -= dot *
out[index];
Data structure for detector geometry and calibration.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Logical location of module.
Auxiliary class to define a range between two values.
Jpp environment information.
Data structure for time calibration.
double getT0() const
Get time offset.
static void setStartOfComment(const start_of_comment_type option)
Set option for short start of comment in binary I/O.
Data structure for a composite optical module.
const JPMT & getPMT(const int index) const
Get PMT.
Monte Carlo detector (i.e. ".det" file).
double getRadius() const
Get radius.
double getLength() const
Get length.
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.
JVector3D & sub(const JVector3D &vector)
Subtract vector.
Data structure for normalised vector in three dimensions.
double getDY() const
Get y direction.
double getDX() const
Get x direction.
double getDZ() const
Get z direction.
Binary buffered file input.
virtual void clear() override
Clear status of reader.
Binary buffered file output.
Exception for opening of file.
int getID() const
Get identifier.
Exception for accessing a value in a collection that is outside of its range.
JMatrix3D & setIdentity()
Set to identity matrix.
file Auxiliary data structures and methods for detector calibration.
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).
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.