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());
827 bool orthonormalise(
const size_t index,
const double precision = std::numeric_limits<double>::epsilon())
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.
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 source JAcoustics sh $DETECTOR_ID CHECK_EXIT_CODE typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS 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.