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());
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());
289 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
290 buffer.insert(module->getID());
293 return buffer.size();
307 const double phi = atan2(dir.
getDY(), dir.
getDZ())*(180.0/
PI);
328 return asin(-dir.
getDX())*(180.0/
PI);
349 const double DEFAULT_CAN_MARGIN_M = 350.0;
350 const double DEFAULT_WORLD_MARGIN_M = 50.0;
352 const JCylinder3D cylinder(detector.begin(), detector.end());
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;
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;
366 const double Crust_Z_size = WorldCylinderHeight/2;
367 const double Crust_Z_position = -WorldCylinderHeight/4;
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=\"";
371 cerr <<
"GDML_SCHEMA_DIR NOT DEFINED! Setting default path." << endl;
378 out <<
"<rotation name=\"identity\"/>\n<position name=\"zero\"/>\n";
382 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
384 if(module->empty())
continue;
386 const JVector3D center = module->getCenter();
388 out <<
"<position name=\"PosString" << module->getString() <<
"_Dom" << module->getID() <<
"\" unit=\"m\" x=\"" << module->
getX() <<
"\" y=\"" << module->getY() <<
"\" z=\"" << module->getZ() <<
"\"/>\n";
390 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
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";
401 out <<
"<position name=\"OMShift\" unit=\"m\" x=\"0\" y=\"0\" z=\"0.0392\"/>\n";
403 out <<
"<!-- end of DU position definitions -->\n<position name=\"CrustPosition\" unit=\"m\" x=\"0\" y=\"0\" z=\"" << Crust_Z_position <<
"\"/>\n";
405 out <<
"</define>\n";
406 out <<
"<materials>\n";
407 out <<
"</materials>\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";
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";
422 out <<
"<!-- OMVolume(s) construction -->\n";
424 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
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";
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";
439 out <<
" </volume>\n";
442 out <<
"<!-- StoreyVolume(s) construction -->\n";
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";
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";
463 out <<
"<!-- World Volume construction-->\n";
464 out <<
"<volume name=\"WorldVolume\">\n";
465 out <<
" <materialref ref=\"Water\"/>\n";
466 out <<
" <solidref ref=\"WorldTube\"/>\n";
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";
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";
483 out <<
"</volume>\n";
485 out <<
"</structure>\n";
486 out <<
"<setup name=\"Default\" version=\"1.0\">\n";
487 out <<
"<world ref=\"WorldVolume\"/>\n";
514 ifstream in(file_name.c_str());
524 detector.swap(buffer);
539 }
catch(
const exception& error) {
557 ifstream in(file_name.c_str());
621 std::ofstream out(file_name.c_str());
645 std::ofstream out(file_name.c_str());
669 static const size_t NUMBER_OF_DIMENSIONS = 3;
683 if (first.size() == second.size()) {
685 const size_t N = first.size();
687 if (N >= NUMBER_OF_DIMENSIONS) {
692 for (
size_t i = 0; i != N; ++i) {
697 for (
size_t i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
698 if (!orthonormalise(i)) {
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();
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();
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();
717 THROW(
JException,
"Module " << first.
getID() <<
" size " << N <<
" < " << NUMBER_OF_DIMENSIONS);
722 THROW(
JException,
"Module " << first.
getID() <<
" size " << first.size() <<
" != " << second.size());
743 for (
size_t i = index + 1; i != in.size(); ++i) {
744 if (in[i].getLengthSquared() > in[pos].getLengthSquared()) {
749 const double u = in[pos].getLength();
757 swap(in [pos], in [index]);
758 swap(out[pos], out[index]);
761 for (
size_t i = index + 1; i != in.size(); ++i) {
763 const double dot = in[index].getDot(in[i]);
765 in [i] -= dot * in [index];
766 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.
std::set< int > getModuleIDs(const JDetector &detector, const bool option=false)
Get list of modules identifiers.
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.