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");
86 for (JDetector::const_iterator i1 = detector.begin(); i1 != detector.end(); ++i1) {
87 for (JDetector::const_iterator i2 = detector.begin(); i2 != i1; ++i2) {
89 if (option || (i1->getFloor() != 0 && i2->getFloor() != 0)) {
91 const double ds =
getDistance(i1->getPosition(), i2->getPosition());
114 const double phi = atan2(dir.
getDY(), dir.
getDZ())*(180.0/
PI);
135 return asin(-dir.
getDX())*(180.0/
PI);
156 const double DEFAULT_CAN_MARGIN_M = 350.0;
157 const double DEFAULT_WORLD_MARGIN_M = 50.0;
159 const JCylinder3D cylinder(detector.begin(), detector.end());
166 cerr <<
"CAN_MARGIN_M not defined! Setting standard CAN_MARGIN_M " << DEFAULT_CAN_MARGIN_M <<
" m." << endl;
167 can_margin = DEFAULT_CAN_MARGIN_M;
170 const double WorldCylinderHeight = 2*(cylinder.getZmax() - cylinder.getZmin() + can_margin + DEFAULT_WORLD_MARGIN_M);
171 const double WorldRadius = cylinder.getLength() + cylinder.getRadius() + can_margin + DEFAULT_WORLD_MARGIN_M;
173 const double Crust_Z_size = WorldCylinderHeight/2;
174 const double Crust_Z_position = -WorldCylinderHeight/4;
176 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=\"";
178 cerr <<
"GDML_SCHEMA_DIR NOT DEFINED! Setting default path." << endl;
185 out <<
"<rotation name=\"identity\"/>\n<position name=\"zero\"/>\n";
189 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
191 if(module->empty())
continue;
193 const JVector3D center = module->getCenter();
195 out <<
"<position name=\"PosString" << module->getString() <<
"_Dom" << module->getID() <<
"\" unit=\"m\" x=\"" << module->
getX() <<
"\" y=\"" << module->getY() <<
"\" z=\"" << module->getZ() <<
"\"/>\n";
197 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
200 out <<
"<position name=\"CathodPosition" << pmt->getID() <<
"_" << module->getID() <<
"\" unit=\"m\" x=\"" << vec.
getX() <<
"\" y=\"" << vec.
getY() <<
"\" z=\"" << vec.
getZ() <<
"\"/>\n";
201 out <<
"<rotation name=\"CathodRotation" << pmt->getID() <<
"_" << module->getID() <<
"\" unit=\"deg\" x=\"" <<
GetXrotationG4(*pmt) <<
"\" y=\"" <<
GetYrotationG4(*pmt) <<
"\" z=\"0.000000\"/>\n";
202 out <<
"<constant name=\"CathodID_" << PMTs_NO <<
"\" value=\"" << pmt->getID() <<
"\"/>\n";
208 out <<
"<position name=\"OMShift\" unit=\"m\" x=\"0\" y=\"0\" z=\"0.0392\"/>\n";
210 out <<
"<!-- end of DU position definitions -->\n<position name=\"CrustPosition\" unit=\"m\" x=\"0\" y=\"0\" z=\"" << Crust_Z_position <<
"\"/>\n";
212 out <<
"</define>\n";
213 out <<
"<materials>\n";
214 out <<
"</materials>\n";
216 out <<
" <box name=\"CrustBox\" lunit=\"m\" x=\"2200\" y=\"2200\" z=\"" << Crust_Z_size <<
"\"/>\n";
217 out <<
" <box name=\"StoreyBox\" lunit=\"m\" x=\"0.6\" y=\"0.6\" z=\"0.6\"/>\n";
218 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";
219 out <<
" <tube name=\"CathodTube\" lunit=\"cm\" aunit=\"deg\" rmin=\"0.0\" rmax=\"4.7462\" z=\"0.5\" startphi=\"0.0\" deltaphi=\"360.0\"/>\n";
220 out <<
" <tube name=\"WorldTube\" lunit=\"m\" aunit=\"deg\" rmin=\"0.0\" rmax=\"" << WorldRadius <<
"\" z=\"" << WorldCylinderHeight <<
"\" startphi=\"0.0\" deltaphi=\"360.0\"/>\n";
221 out <<
"</solids>\n\n\n";
223 out <<
"<structure>\n";
224 out <<
" <volume name=\"CathodVolume\">\n";
225 out <<
" <materialref ref=\"Cathod\"/>\n";
226 out <<
" <solidref ref=\"CathodTube\"/>\n";
227 out <<
" </volume>\n";
229 out <<
"<!-- OMVolume(s) construction -->\n";
231 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
233 if(module->empty())
continue;
234 out <<
" <volume name=\"OMVolume" << module->getID() <<
"\">\n";
235 out <<
" <materialref ref=\"Water\"/>\n";
236 out <<
" <solidref ref=\"OMSphere\"/>\n";
238 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
239 out <<
" <physvol>\n";
240 out <<
" <volumeref ref=\"CathodVolume\"/>\n";
241 out <<
" <positionref ref=\"CathodPosition" << pmt->getID() <<
"_" << module->getID() <<
"\"/>\n";
242 out <<
" <rotationref ref=\"CathodRotation" << pmt->getID() <<
"_" << module->getID() <<
"\"/>\n";
243 out <<
" </physvol>\n";
246 out <<
" </volume>\n";
249 out <<
"<!-- StoreyVolume(s) construction -->\n";
251 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
252 if(module->empty())
continue;
253 out <<
" <volume name=\"StoreyVolume" << module->getID() <<
"\">\n";
254 out <<
" <materialref ref=\"Water\"/>\n";
255 out <<
" <solidref ref=\"StoreyBox\"/>\n";
256 out <<
" <physvol>\n";
257 out <<
" <volumeref ref=\"OMVolume" << module->getID() <<
"\"/>\n";
258 out <<
" <positionref ref=\"OMShift\"/>\n";
259 out <<
" <rotationref ref=\"identity\"/>\n";
260 out <<
" </physvol>\n";
261 out <<
" </volume>\n";
264 out <<
"<!-- Crust Volume construction-->\n";
265 out <<
"<volume name=\"CrustVolume\">\n";
266 out <<
" <materialref ref=\"Crust\"/>\n";
267 out <<
" <solidref ref=\"CrustBox\"/>\n";
268 out <<
"</volume>\n";
270 out <<
"<!-- World Volume construction-->\n";
271 out <<
"<volume name=\"WorldVolume\">\n";
272 out <<
" <materialref ref=\"Water\"/>\n";
273 out <<
" <solidref ref=\"WorldTube\"/>\n";
275 out <<
" <physvol>\n";
276 out <<
" <volumeref ref=\"CrustVolume\"/>\n";
277 out <<
" <positionref ref=\"CrustPosition\"/>\n";
278 out <<
" <rotationref ref=\"identity\"/>\n";
279 out <<
" </physvol>\n";
281 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
282 if(module->empty())
continue;
283 out <<
" <physvol>\n";
284 out <<
" <volumeref ref=\"StoreyVolume" << module->getID() <<
"\"/>\n";
285 out <<
" <positionref ref=\"PosString" << module->getString() <<
"_Dom" << module->getID() <<
"\"/>\n";
286 out <<
" <rotationref ref=\"identity\"/>\n";
287 out <<
" </physvol>\n";
290 out <<
"</volume>\n";
292 out <<
"</structure>\n";
293 out <<
"<setup name=\"Default\" version=\"1.0\">\n";
294 out <<
"<world ref=\"WorldVolume\"/>\n";
345 if (!module.empty()) {
349 for (JModule::const_iterator pmt = module.begin(); pmt != module.end(); ++pmt) {
374 return module.size();
386 int number_of_pmts = 0;
388 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
392 return number_of_pmts;
406 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
407 buffer.insert(module->getString());
424 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
425 buffer.insert(module->getFloor());
428 return buffer.size();
448 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
449 buffer.
include(module->getFloor());
466 for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
467 buffer.insert(module->getID());
470 return buffer.size();
495 ifstream
in(file_name.c_str());
505 detector.swap(buffer);
520 }
catch(
const exception& error) {
538 ifstream
in(file_name.c_str());
598 std::ofstream out(file_name.c_str());
622 std::ofstream out(file_name.c_str());
658 module.resize(memo.size());
711 template<
class JDetector_t>
716 return getModule(getModuleAddressMap<JDetector_t>(
id),
id, location);
738 if (module.empty()) {
742 const double R = 0.5;
744 const double st =
sin(0.75*
PI);
745 const double ct = cos(0.75*
PI);
747 for (
int i = 0; i != 3; ++i) {
749 const double phi = (2.0*
PI*i) / 3.0;
750 const double cp = cos(phi);
751 const double sp =
sin(phi);
768 template<
class JDetector_t>
783 static const size_t NUMBER_OF_DIMENSIONS = 3;
797 if (first.size() == second.size()) {
799 const size_t N = first.size();
801 if (N >= NUMBER_OF_DIMENSIONS) {
806 for (
size_t i = 0; i !=
N; ++i) {
811 for (
size_t i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
812 if (!orthonormalise(i)) {
817 this->a00 = out[0].getX() *
in[0].getX() + out[1].getX() *
in[1].getX() + out[2].getX() *
in[2].getX();
818 this->a01 = out[0].getX() *
in[0].getY() + out[1].getX() *
in[1].getY() + out[2].getX() *
in[2].getY();
819 this->a02 = out[0].getX() *
in[0].getZ() + out[1].getX() *
in[1].getZ() + out[2].getX() *
in[2].getZ();
821 this->a10 = out[0].getY() *
in[0].getX() + out[1].getY() *
in[1].getX() + out[2].getY() *
in[2].getX();
822 this->a11 = out[0].getY() *
in[0].getY() + out[1].getY() *
in[1].getY() + out[2].getY() *
in[2].getY();
823 this->a12 = out[0].getY() *
in[0].getZ() + out[1].getY() *
in[1].getZ() + out[2].getY() *
in[2].getZ();
825 this->a20 = out[0].getZ() *
in[0].getX() + out[1].getZ() *
in[1].getX() + out[2].getZ() *
in[2].getX();
826 this->a21 = out[0].getZ() *
in[0].getY() + out[1].getZ() *
in[1].getY() + out[2].getZ() *
in[2].getY();
827 this->a22 = out[0].getZ() *
in[0].getZ() + out[1].getZ() *
in[1].getZ() + out[2].getZ() *
in[2].getZ();
831 THROW(
JException,
"Module " << first.
getID() <<
" size " << N <<
" < " << NUMBER_OF_DIMENSIONS);
836 THROW(
JException,
"Module " << first.
getID() <<
" size " << first.size() <<
" != " << second.size());
857 for (
size_t i = index + 1; i !=
in.size(); ++i) {
858 if (
in[i].getLengthSquared() >
in[pos].getLengthSquared()) {
863 const double u =
in[pos].getLength();
871 swap(
in [pos],
in [index]);
872 swap(out[pos], out[index]);
875 for (
size_t i = index + 1; i !=
in.size(); ++i) {
877 const double dot =
in[index].getDot(
in[i]);
879 in [i] -= dot *
in [index];
880 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.
#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.
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.
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
then JCookie sh JDataQuality D $DETECTOR_ID R
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.
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.
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.