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());
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());
226 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
227 buffer.insert(module->getFloor());
230 return buffer.size();
250 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
251 buffer.
include(module->getFloor());
268 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
269 buffer.insert(module->getID());
272 return buffer.size();
286 const double phi = atan2(dir.
getDY(), dir.
getDZ())*(180.0/
PI);
307 return asin(-dir.
getDX())*(180.0/
PI);
328 const double DEFAULT_CAN_MARGIN_M = 350.0;
329 const double DEFAULT_WORLD_MARGIN_M = 50.0;
331 const JCylinder3D cylinder(detector.begin(), detector.end());
338 cerr <<
"CAN_MARGIN_M not defined! Setting standard CAN_MARGIN_M " << DEFAULT_CAN_MARGIN_M <<
" m." << endl;
339 can_margin = DEFAULT_CAN_MARGIN_M;
342 const double WorldCylinderHeight = 2*(cylinder.
getZmax() - cylinder.
getZmin() + can_margin + DEFAULT_WORLD_MARGIN_M);
343 const double WorldRadius = cylinder.
getLength() + cylinder.
getRadius() + can_margin + DEFAULT_WORLD_MARGIN_M;
345 const double Crust_Z_size = WorldCylinderHeight/2;
346 const double Crust_Z_position = -WorldCylinderHeight/4;
348 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=\"";
350 cerr <<
"GDML_SCHEMA_DIR NOT DEFINED! Setting default path." << endl;
357 out <<
"<rotation name=\"identity\"/>\n<position name=\"zero\"/>\n";
361 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
363 if(module->empty())
continue;
365 const JVector3D center = module->getCenter();
367 out <<
"<position name=\"PosString" << module->getString() <<
"_Dom" << module->getID() <<
"\" unit=\"m\" x=\"" << module->
getX() <<
"\" y=\"" << module->getY() <<
"\" z=\"" << module->getZ() <<
"\"/>\n";
369 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
372 out <<
"<position name=\"CathodPosition" << pmt->getID() <<
"_" << module->getID() <<
"\" unit=\"m\" x=\"" << vec.
getX() <<
"\" y=\"" << vec.
getY() <<
"\" z=\"" << vec.
getZ() <<
"\"/>\n";
373 out <<
"<rotation name=\"CathodRotation" << pmt->getID() <<
"_" << module->getID() <<
"\" unit=\"deg\" x=\"" <<
GetXrotationG4(*pmt) <<
"\" y=\"" <<
GetYrotationG4(*pmt) <<
"\" z=\"0.000000\"/>\n";
374 out <<
"<constant name=\"CathodID_" << PMTs_NO <<
"\" value=\"" << pmt->getID() <<
"\"/>\n";
380 out <<
"<position name=\"OMShift\" unit=\"m\" x=\"0\" y=\"0\" z=\"0.0392\"/>\n";
382 out <<
"<!-- end of DU position definitions -->\n<position name=\"CrustPosition\" unit=\"m\" x=\"0\" y=\"0\" z=\"" << Crust_Z_position <<
"\"/>\n";
384 out <<
"</define>\n";
385 out <<
"<materials>\n";
386 out <<
"</materials>\n";
388 out <<
" <box name=\"CrustBox\" lunit=\"m\" x=\"2200\" y=\"2200\" z=\"" << Crust_Z_size <<
"\"/>\n";
389 out <<
" <box name=\"StoreyBox\" lunit=\"m\" x=\"0.6\" y=\"0.6\" z=\"0.6\"/>\n";
390 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";
391 out <<
" <tube name=\"CathodTube\" lunit=\"cm\" aunit=\"deg\" rmin=\"0.0\" rmax=\"4.7462\" z=\"0.5\" startphi=\"0.0\" deltaphi=\"360.0\"/>\n";
392 out <<
" <tube name=\"WorldTube\" lunit=\"m\" aunit=\"deg\" rmin=\"0.0\" rmax=\"" << WorldRadius <<
"\" z=\"" << WorldCylinderHeight <<
"\" startphi=\"0.0\" deltaphi=\"360.0\"/>\n";
393 out <<
"</solids>\n\n\n";
395 out <<
"<structure>\n";
396 out <<
" <volume name=\"CathodVolume\">\n";
397 out <<
" <materialref ref=\"Cathod\"/>\n";
398 out <<
" <solidref ref=\"CathodTube\"/>\n";
399 out <<
" </volume>\n";
401 out <<
"<!-- OMVolume(s) construction -->\n";
403 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
405 if(module->empty())
continue;
406 out <<
" <volume name=\"OMVolume" << module->getID() <<
"\">\n";
407 out <<
" <materialref ref=\"Water\"/>\n";
408 out <<
" <solidref ref=\"OMSphere\"/>\n";
410 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
411 out <<
" <physvol>\n";
412 out <<
" <volumeref ref=\"CathodVolume\"/>\n";
413 out <<
" <positionref ref=\"CathodPosition" << pmt->getID() <<
"_" << module->getID() <<
"\"/>\n";
414 out <<
" <rotationref ref=\"CathodRotation" << pmt->getID() <<
"_" << module->getID() <<
"\"/>\n";
415 out <<
" </physvol>\n";
418 out <<
" </volume>\n";
421 out <<
"<!-- StoreyVolume(s) construction -->\n";
423 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
424 if(module->empty())
continue;
425 out <<
" <volume name=\"StoreyVolume" << module->getID() <<
"\">\n";
426 out <<
" <materialref ref=\"Water\"/>\n";
427 out <<
" <solidref ref=\"StoreyBox\"/>\n";
428 out <<
" <physvol>\n";
429 out <<
" <volumeref ref=\"OMVolume" << module->getID() <<
"\"/>\n";
430 out <<
" <positionref ref=\"OMShift\"/>\n";
431 out <<
" <rotationref ref=\"identity\"/>\n";
432 out <<
" </physvol>\n";
433 out <<
" </volume>\n";
436 out <<
"<!-- Crust Volume construction-->\n";
437 out <<
"<volume name=\"CrustVolume\">\n";
438 out <<
" <materialref ref=\"Crust\"/>\n";
439 out <<
" <solidref ref=\"CrustBox\"/>\n";
440 out <<
"</volume>\n";
442 out <<
"<!-- World Volume construction-->\n";
443 out <<
"<volume name=\"WorldVolume\">\n";
444 out <<
" <materialref ref=\"Water\"/>\n";
445 out <<
" <solidref ref=\"WorldTube\"/>\n";
447 out <<
" <physvol>\n";
448 out <<
" <volumeref ref=\"CrustVolume\"/>\n";
449 out <<
" <positionref ref=\"CrustPosition\"/>\n";
450 out <<
" <rotationref ref=\"identity\"/>\n";
451 out <<
" </physvol>\n";
453 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
454 if(module->empty())
continue;
455 out <<
" <physvol>\n";
456 out <<
" <volumeref ref=\"StoreyVolume" << module->getID() <<
"\"/>\n";
457 out <<
" <positionref ref=\"PosString" << module->getString() <<
"_Dom" << module->getID() <<
"\"/>\n";
458 out <<
" <rotationref ref=\"identity\"/>\n";
459 out <<
" </physvol>\n";
462 out <<
"</volume>\n";
464 out <<
"</structure>\n";
465 out <<
"<setup name=\"Default\" version=\"1.0\">\n";
466 out <<
"<world ref=\"WorldVolume\"/>\n";
493 ifstream in(file_name.c_str());
503 detector.swap(buffer);
518 }
catch(
const exception& error) {
536 ifstream in(file_name.c_str());
600 std::ofstream out(file_name.c_str());
624 std::ofstream out(file_name.c_str());
648 static const size_t NUMBER_OF_DIMENSIONS = 3;
662 if (first.size() == second.size()) {
664 const size_t N = first.size();
666 if (N >= NUMBER_OF_DIMENSIONS) {
671 for (
size_t i = 0; i != N; ++i) {
676 for (
size_t i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
677 if (!orthonormalise(i)) {
682 this->a00 = out[0].getX() * in[0].getX() + out[1].getX() * in[1].getX() + out[2].getX() * in[2].getX();
683 this->a01 = out[0].getX() * in[0].getY() + out[1].getX() * in[1].getY() + out[2].getX() * in[2].getY();
684 this->a02 = out[0].getX() * in[0].getZ() + out[1].getX() * in[1].getZ() + out[2].getX() * in[2].getZ();
686 this->a10 = out[0].getY() * in[0].getX() + out[1].getY() * in[1].getX() + out[2].getY() * in[2].getX();
687 this->a11 = out[0].getY() * in[0].getY() + out[1].getY() * in[1].getY() + out[2].getY() * in[2].getY();
688 this->a12 = out[0].getY() * in[0].getZ() + out[1].getY() * in[1].getZ() + out[2].getY() * in[2].getZ();
690 this->a20 = out[0].getZ() * in[0].getX() + out[1].getZ() * in[1].getX() + out[2].getZ() * in[2].getX();
691 this->a21 = out[0].getZ() * in[0].getY() + out[1].getZ() * in[1].getY() + out[2].getZ() * in[2].getY();
692 this->a22 = out[0].getZ() * in[0].getZ() + out[1].getZ() * in[1].getZ() + out[2].getZ() * in[2].getZ();
696 THROW(
JException,
"Module " << first.
getID() <<
" size " << N <<
" < " << NUMBER_OF_DIMENSIONS);
701 THROW(
JException,
"Module " << first.
getID() <<
" size " << first.size() <<
" != " << second.size());
722 for (
size_t i = index + 1; i != in.size(); ++i) {
723 if (in[i].getLengthSquared() > in[pos].getLengthSquared()) {
728 const double u = in[pos].getLength();
736 swap(in [pos], in [index]);
737 swap(out[pos], out[index]);
740 for (
size_t i = index + 1; i != in.size(); ++i) {
742 const double dot = in[index].getDot(in[i]);
744 in [i] -= dot * in [index];
745 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.
virtual JWriter & write(JWriter &out) const override
Write to output.
virtual JReader & read(JReader &in) override
Read from input.
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.
double getY() const
Get y position.
double getZ() const
Get z position.
JVector3D & sub(const JVector3D &vector)
Subtract vector.
double getX() const
Get x position.
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.
file Auxiliary data structures and methods for detector calibration.
static JRotation getRotation
Function object to get rotation matrix to go from first to second module.
std::set< int > getStringIDs(const JDetector &detector)
Get list of strings identifiers.
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
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.
int getNumberOfPMTs(const JDetector &detector)
Get number of PMTs.
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
double getMaximalTime(const JDetector &detector, const double roadWidth_m)
Get maximal time between optical modules in detector following causality.
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.
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.
std::string getFilenameExtension(const std::string &file_name)
Get file name extension, i.e. part after last JEEP::FILENAME_SEPARATOR if any.
std::string getGITVersion(const std::string &tag)
Get GIT version for given GIT tag.
static const double PI
Mathematical constants.
double getDistance(const JFirst_t &first, const JSecond_t &second)
Get distance between objects.
double getIndexOfRefraction()
Get average index of refraction of water corresponding to group velocity.
const double getInverseSpeedOfLight()
Get inverse speed of light.
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
double getSinThetaC()
Get average sine of Cherenkov angle of water corresponding to group velocity.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
bool hasVersion() const
Check validity of version.
Auxiliary class to get rotation matrix between two optical modules.
std::vector< JVector3D > out
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.
const JRotation3D & operator()(const JModule &first, const JModule &second)
Get rotation matrix to go from first to second module.
const std::string & getVersion() const
Get version.