1 #ifndef __JDETECTOR__JDETECTORTOOLKIT__
2 #define __JDETECTOR__JDETECTORTOOLKIT__
44 namespace JDETECTOR {}
45 namespace JPP {
using namespace JDETECTOR; }
67 static const char*
const GDML_SCHEMA = getenv(
"GDML_SCHEMA_DIR");
84 for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {
85 for (JDetector::const_iterator
j = detector.begin();
j != i; ++
j) {
86 if (
getDistance(i->getPosition(),
j->getPosition()) > dmax) {
106 const double phi = atan2(dir.
getDY(), dir.
getDZ())*(180.0/
PI);
127 return asin(-dir.
getDX())*(180.0/
PI);
148 const double DEFAULT_CAN_MARGIN_M = 350.0;
149 const double DEFAULT_WORLD_MARGIN_M = 50.0;
151 const JCylinder3D cylinder(detector.begin(), detector.end());
158 cerr <<
"CAN_MARGIN_M not defined! Setting standard CAN_MARGIN_M " << DEFAULT_CAN_MARGIN_M <<
" m." << endl;
159 can_margin = DEFAULT_CAN_MARGIN_M;
162 const double WorldCylinderHeight = 2*(cylinder.getZmax() - cylinder.getZmin() + can_margin + DEFAULT_WORLD_MARGIN_M);
163 const double WorldRadius = cylinder.getLength() + cylinder.getRadius() + can_margin + DEFAULT_WORLD_MARGIN_M;
165 const double Crust_Z_size = WorldCylinderHeight/2;
166 const double Crust_Z_position = -WorldCylinderHeight/4;
168 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=\"";
170 cerr <<
"GDML_SCHEMA_DIR NOT DEFINED! Setting default path." << endl;
177 out <<
"<rotation name=\"identity\"/>\n<position name=\"zero\"/>\n";
181 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
183 if(module->empty())
continue;
185 const JVector3D center = module->getCenter();
187 out <<
"<position name=\"PosString" << module->getString() <<
"_Dom" << module->getID() <<
"\" unit=\"m\" x=\"" << module->
getX() <<
"\" y=\"" << module->getY() <<
"\" z=\"" << module->getZ() <<
"\"/>\n";
189 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
192 out <<
"<position name=\"CathodPosition" << pmt->getID() <<
"_" << module->getID() <<
"\" unit=\"m\" x=\"" << vec.
getX() <<
"\" y=\"" << vec.
getY() <<
"\" z=\"" << vec.
getZ() <<
"\"/>\n";
193 out <<
"<rotation name=\"CathodRotation" << pmt->getID() <<
"_" << module->getID() <<
"\" unit=\"deg\" x=\"" <<
GetXrotationG4(*pmt) <<
"\" y=\"" <<
GetYrotationG4(*pmt) <<
"\" z=\"0.000000\"/>\n";
194 out <<
"<constant name=\"CathodID_" << PMTs_NO <<
"\" value=\"" << pmt->getID() <<
"\"/>\n";
200 out <<
"<position name=\"OMShift\" unit=\"m\" x=\"0\" y=\"0\" z=\"0.0392\"/>\n";
202 out <<
"<!-- end of DU position definitions -->\n<position name=\"CrustPosition\" unit=\"m\" x=\"0\" y=\"0\" z=\"" << Crust_Z_position <<
"\"/>\n";
204 out <<
"</define>\n";
205 out <<
"<materials>\n";
206 out <<
"</materials>\n";
208 out <<
" <box name=\"CrustBox\" lunit=\"m\" x=\"2200\" y=\"2200\" z=\"" << Crust_Z_size <<
"\"/>\n";
209 out <<
" <box name=\"StoreyBox\" lunit=\"m\" x=\"0.6\" y=\"0.6\" z=\"0.6\"/>\n";
210 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";
211 out <<
" <tube name=\"CathodTube\" lunit=\"cm\" aunit=\"deg\" rmin=\"0.0\" rmax=\"4.7462\" z=\"0.5\" startphi=\"0.0\" deltaphi=\"360.0\"/>\n";
212 out <<
" <tube name=\"WorldTube\" lunit=\"m\" aunit=\"deg\" rmin=\"0.0\" rmax=\"" << WorldRadius <<
"\" z=\"" << WorldCylinderHeight <<
"\" startphi=\"0.0\" deltaphi=\"360.0\"/>\n";
213 out <<
"</solids>\n\n\n";
215 out <<
"<structure>\n";
216 out <<
" <volume name=\"CathodVolume\">\n";
217 out <<
" <materialref ref=\"Cathod\"/>\n";
218 out <<
" <solidref ref=\"CathodTube\"/>\n";
219 out <<
" </volume>\n";
221 out <<
"<!-- OMVolume(s) construction -->\n";
223 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
225 if(module->empty())
continue;
226 out <<
" <volume name=\"OMVolume" << module->getID() <<
"\">\n";
227 out <<
" <materialref ref=\"Water\"/>\n";
228 out <<
" <solidref ref=\"OMSphere\"/>\n";
230 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
231 out <<
" <physvol>\n";
232 out <<
" <volumeref ref=\"CathodVolume\"/>\n";
233 out <<
" <positionref ref=\"CathodPosition" << pmt->getID() <<
"_" << module->getID() <<
"\"/>\n";
234 out <<
" <rotationref ref=\"CathodRotation" << pmt->getID() <<
"_" << module->getID() <<
"\"/>\n";
235 out <<
" </physvol>\n";
238 out <<
" </volume>\n";
241 out <<
"<!-- StoreyVolume(s) construction -->\n";
243 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
244 if(module->empty())
continue;
245 out <<
" <volume name=\"StoreyVolume" << module->getID() <<
"\">\n";
246 out <<
" <materialref ref=\"Water\"/>\n";
247 out <<
" <solidref ref=\"StoreyBox\"/>\n";
248 out <<
" <physvol>\n";
249 out <<
" <volumeref ref=\"OMVolume" << module->getID() <<
"\"/>\n";
250 out <<
" <positionref ref=\"OMShift\"/>\n";
251 out <<
" <rotationref ref=\"identity\"/>\n";
252 out <<
" </physvol>\n";
253 out <<
" </volume>\n";
256 out <<
"<!-- Crust Volume construction-->\n";
257 out <<
"<volume name=\"CrustVolume\">\n";
258 out <<
" <materialref ref=\"Crust\"/>\n";
259 out <<
" <solidref ref=\"CrustBox\"/>\n";
260 out <<
"</volume>\n";
262 out <<
"<!-- World Volume construction-->\n";
263 out <<
"<volume name=\"WorldVolume\">\n";
264 out <<
" <materialref ref=\"Water\"/>\n";
265 out <<
" <solidref ref=\"WorldTube\"/>\n";
267 out <<
" <physvol>\n";
268 out <<
" <volumeref ref=\"CrustVolume\"/>\n";
269 out <<
" <positionref ref=\"CrustPosition\"/>\n";
270 out <<
" <rotationref ref=\"identity\"/>\n";
271 out <<
" </physvol>\n";
273 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
274 if(module->empty())
continue;
275 out <<
" <physvol>\n";
276 out <<
" <volumeref ref=\"StoreyVolume" << module->getID() <<
"\"/>\n";
277 out <<
" <positionref ref=\"PosString" << module->getString() <<
"_Dom" << module->getID() <<
"\"/>\n";
278 out <<
" <rotationref ref=\"identity\"/>\n";
279 out <<
" </physvol>\n";
282 out <<
"</volume>\n";
284 out <<
"</structure>\n";
285 out <<
"<setup name=\"Default\" version=\"1.0\">\n";
286 out <<
"<world ref=\"WorldVolume\"/>\n";
337 if (!module.empty()) {
341 for (JModule::const_iterator pmt = module.begin(); pmt != module.end(); ++pmt) {
345 time_range.include(
putTime(timeRange.getLowerLimit(), calibration));
346 time_range.include(
putTime(timeRange.getUpperLimit(), calibration));
366 return module.size();
378 int number_of_pmts = 0;
380 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
384 return number_of_pmts;
398 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
399 buffer.insert(module->getString());
416 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
417 buffer.insert(module->getFloor());
420 return buffer.size();
440 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
441 buffer.
include(module->getFloor());
458 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
459 buffer.insert(module->getID());
462 return buffer.size();
487 ifstream
in(file_name.c_str());
497 detector.swap(buffer);
512 }
catch(
const exception& error) {
530 ifstream
in(file_name.c_str());
586 std::ofstream out(file_name.c_str());
602 std::ofstream out(file_name.c_str());
634 module.resize(memo.size());
687 template<
class JDetector_t>
692 return getModule(getModuleAddressMap<JDetector_t>(
id),
id, location);
714 if (module.empty()) {
718 const double R = 0.5;
720 const double st = sin(0.75*
PI);
721 const double ct = cos(0.75*
PI);
723 for (
int i = 0; i != 3; ++i) {
725 const double phi = (2.0*
PI*i) / 3.0;
726 const double cp = cos(phi);
727 const double sp = sin(phi);
744 template<
class JDetector_t>
759 static const size_t NUMBER_OF_DIMENSIONS = 3;
773 if (first.size() == second.size()) {
775 const size_t N = first.size();
777 if (N >= NUMBER_OF_DIMENSIONS) {
782 for (
size_t i = 0; i !=
N; ++i) {
787 for (
size_t i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
788 if (!orthonormalise(i)) {
793 this->a00 = out[0].getX() *
in[0].getX() + out[1].getX() *
in[1].getX() + out[2].getX() *
in[2].getX();
794 this->a01 = out[0].getX() *
in[0].getY() + out[1].getX() *
in[1].getY() + out[2].getX() *
in[2].getY();
795 this->a02 = out[0].getX() *
in[0].getZ() + out[1].getX() *
in[1].getZ() + out[2].getX() *
in[2].getZ();
797 this->a10 = out[0].getY() *
in[0].getX() + out[1].getY() *
in[1].getX() + out[2].getY() *
in[2].getX();
798 this->a11 = out[0].getY() *
in[0].getY() + out[1].getY() *
in[1].getY() + out[2].getY() *
in[2].getY();
799 this->a12 = out[0].getY() *
in[0].getZ() + out[1].getY() *
in[1].getZ() + out[2].getY() *
in[2].getZ();
801 this->a20 = out[0].getZ() *
in[0].getX() + out[1].getZ() *
in[1].getX() + out[2].getZ() *
in[2].getX();
802 this->a21 = out[0].getZ() *
in[0].getY() + out[1].getZ() *
in[1].getY() + out[2].getZ() *
in[2].getY();
803 this->a22 = out[0].getZ() *
in[0].getZ() + out[1].getZ() *
in[1].getZ() + out[2].getZ() *
in[2].getZ();
807 THROW(
JException,
"Module " << first.
getID() <<
" size " << N <<
" < " << NUMBER_OF_DIMENSIONS);
812 THROW(
JException,
"Module " << first.
getID() <<
" size " << first.size() <<
" != " << second.size());
833 for (
size_t i = index + 1; i !=
in.size(); ++i) {
834 if (
in[i].getLengthSquared() >
in[pos].getLengthSquared()) {
839 const double u =
in[pos].getLength();
847 swap(
in [pos],
in [index]);
848 swap(out[pos], out[index]);
851 for (
size_t i = index + 1; i !=
in.size(); ++i) {
853 const double dot =
in[index].getDot(
in[i]);
855 in [i] -= dot *
in [index];
856 out[i] -= dot * out[index];
Exception for opening of file.
double GetYrotationG4(const JVersor3D dir)
Get rotation over Y axis in Geant4 coordinate system.
Data structure for a composite optical module.
void setLocation(const JLocation &location)
Set location.
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.
JTOOLS::JRange< double > JTimeRange
Type definition for time range (unit [ns]).
#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.
then JShowerPostfit f $INPUT_FILE o $OUTPUT_FILE N
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.
static const char *const BINARY_DETECTOR_FILE_FORMAT[]
JIO binary file format.
Data structure for time calibration.
Auxiliary class for a type holder.
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.
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. so-called .det file).
Lookup table for PMT addresses in optical module.
double putTime(const T &t1, const JCalibration &cal)
Get de-calibrated time.
JVector3D & sub(const JVector3D &vector)
Subtract vector.
double getMaximalDistance(const JDetector &detector)
Get maximal distance between modules in detector.
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.
Logical location of module.
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.
Data structure for PMT geometry, calibration and status.
void compile()
Compile module data.
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.
then usage $script[distance] fi case set_variable R
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
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.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Auxiliary class to define a range between two values.
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.
void setID(const int id)
Set identifier.
double getX() const
Get x position.
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.
do set_variable DETECTOR_TXT $WORKDIR detector
Data structure for normalised vector in three dimensions.
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] 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
Binary buffered file output.
double GetXrotationG4(const JVersor3D dir)
Get rotation over X axis in Geant4 coordinate system.
static const char *const GDML_DETECTOR_FILE_FORMAT
KM3Sim input format.
double getZ() const
Get z position.
Binary buffered file input.
bool has(const int index) const
Test whether index is available.
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.
const JModule & getModule(const JDetector &detector, const JModuleLocation &location)
find module with a given string and floor number
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.