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());
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;
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;
376 out <<
"<!-- Jpp version: "<< getGITVersion() <<
" -->\n";
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());
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());
683 if (first.size() == second.size()) {
685 const size_t N = first.size();
692 for (
size_t i = 0; i != N; ++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();
722 THROW(
JException,
"Module " << first.
getID() <<
" size " << first.size() <<
" != " << second.size());
737 bool orthonormalise(
const size_t index,
const double precision = std::numeric_limits<double>::epsilon())
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.
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
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).
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.