1 #ifndef __JDETECTOR__JDETECTORTOOLKIT__
2 #define __JDETECTOR__JDETECTORTOOLKIT__
41 namespace JDETECTOR {}
42 namespace JPP {
using namespace JDETECTOR; }
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);
508 JFileStreamReader
in(file_name.c_str());
518 }
catch(
const exception& error) {
536 ifstream
in(file_name.c_str());
588 JFileStreamWriter out(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];
Exception for opening of file.
Data structure for a composite optical module.
const JCalibration & getCalibration() const
Get calibration.
double getIndexOfRefraction()
Get average index of refraction of water corresponding to group velocity.
const JDirection3D & getDirection() const
Get direction.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
std::set< int > getStringIDs(const JDetector &detector)
Get list of strings identifiers.
JTOOLS::JRange< int > floor_range
Type definition for range of floors.
int getNumberOfPMTs(const JModule &module)
Get number of PMTs.
static void setStartOfComment(const start_of_comment_type option)
Set option for short start of comment in binary I/O.
static JRotation getRotation
Function object to get rotation matrix to go from first to second module.
JCalibration getCalibration(const JCalibration &first, const JCalibration &second)
Get calibration to go from first to second calibration.
Jpp environment information.
const std::string & getVersion() const
Get version.
static const char *const BINARY_DETECTOR_FILE_FORMAT[]
JIO binary file format.
Data structure for time calibration.
void read_gdml(std::istream &, JDetector &)
static const char *const GDML_SCHEMA
directory necessary for correct GDML header output
static const JModuleCounter getNumberOfModules
Function object to count unique modules.
double GetYrotationG4(const JVersor3D &dir)
Get rotation over Y axis in Geant4 coordinate system.
static const char *const CAN_MARGIN_M
extension of the detector size to comply with the can definition
floor_range getRangeOfFloors(const JDetector &detector)
Get range of floors.
Data structure for detector geometry and calibration.
then echo The file $DIR KM3NeT_00000001_00000000 root already please rename or remove it first
double getDistance(const JFirst_t &first, const JSecond_t &second)
Get distance between objects.
JTimeRange getTimeRange(const Evt &event)
Get time range (i.e. time between earliest and latest hit) of Monte Carlo event.
Monte Carlo detector (i.e. ".det" file).
double putTime(const T &t1, const JCalibration &cal)
Get de-calibrated time.
JVector3D & sub(const JVector3D &vector)
Subtract vector.
virtual JReader & read(JReader &in) override
Read from input.
std::string getGITVersion(const std::string &tag)
Get GIT version for given GIT tag.
Data structure for vector in three dimensions.
virtual JWriter & write(JWriter &out) const override
Write to output.
double getDY() const
Get y direction.
Auxiliary class to get rotation matrix between two optical modules.
double getDX() const
Get x direction.
int getID() const
Get identifier.
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
JPosition3D getPosition(const Vec &pos)
Get position.
int getNumberOfFloors(const JDetector &detector)
Get number of floors.
static const double PI
Mathematical constants.
double getY() const
Get y position.
const JPosition3D & getPosition() const
Get position.
double getMaximalTime(const double R_Hz)
Get maximal time for given rate.
const JPMT & getPMT(const int index) const
Get PMT.
Logical location of module.
const JRotation3D & operator()(const JModule &first, const JModule &second)
Get rotation matrix to go from first to second module.
std::vector< JVector3D > in
double GetXrotationG4(const JVersor3D &dir)
Get rotation over X axis in Geant4 coordinate system.
static const char *const KM3NET_DETECTOR_FILE_FORMAT
KM3NeT standard ASCII format
static const char *const ZIPPED_DETECTOR_FILE_FORMAT
zipped KM3NeT standard ASCII format
static const char *const GENDET_DETECTOR_FILE_FORMAT
File name extensions.
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Auxiliary class to define a range between two values.
then fatal The output file must have the wildcard in the e g root fi eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
std::string getFilenameExtension(const std::string &file_name)
Get file name extension, i.e. part after last JEEP::FILENAME_SEPARATOR if any.
const double getInverseSpeedOfLight()
Get inverse speed of light.
double getX() const
Get x position.
double getMaximalDistance(const JDetector &detector, const bool option=false)
Get maximal distance between modules in detector.
static const char *const G4GDML_DEFAULT_SCHEMALOCATION
double getTanThetaC()
Get average tangent of Cherenkov angle of water corresponding to group velocity.
Data structure for position in three dimensions.
Exception for accessing a value in a collection that is outside of its range.
do set_variable DETECTOR_TXT $WORKDIR detector
Data structure for normalised vector in three dimensions.
bool hasVersion() const
Check validity of version.
static const char *const GDML_DETECTOR_FILE_FORMAT
KM3Sim input format.
double getZ() const
Get z position.
double getDZ() const
Get z direction.
std::vector< JVector3D > out
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.
double getSinThetaC()
Get average sine of Cherenkov angle of water corresponding to group velocity.
void write_gdml(std::ostream &out, const JDetector &detector)
Writes KM3Sim GDML input file from detector.
double getT0() const
Get time offset.