1 #ifndef __JDETECTOR__JDETECTORTOOLKIT__
2 #define __JDETECTOR__JDETECTORTOOLKIT__
67 static const char*
const GDML_SCHEMA = getenv(
"GDML_SCHEMA_DIR");
68 static const char*
const CAN_MARGIN = getenv(
"CAN_MARGIN");
82 for (JDetector::const_iterator i =
detector.begin(); i !=
detector.end(); ++i) {
83 for (JDetector::const_iterator
j =
detector.begin();
j != i; ++
j) {
84 if (i->getDistance(*
j) > dmax) {
85 dmax = i->getDistance(*
j);
102 const double phi = atan2(dir.
getDY(), dir.
getDZ())*(180.0/
PI);
121 return asin(-dir.
getDX())*(180.0/
PI);
143 std::cout<<
"CAN_MARGIN not defined! Setting standard CAN_MARGIN = 280 m"<<std::endl;
158 const double WorldCylinderHeight = 2*(cylinder.getZmax() - cylinder.getZmin() + can_margin);
159 const double WorldRadius = cylinder.getRadius() + can_margin;
161 std::cout<<WorldCylinderHeight/2<<std::endl;
171 const double Crust_Z_size = WorldCylinderHeight/2;
172 const double Crust_Z_position= -WorldCylinderHeight/4;
174 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=\"";
176 std::cout<<
"GDML_SCHEMA_DIR NOT DEFINED! Setting default path"<<std::endl;
184 out<<
"<rotation name=\"identity\"/>\n<position name=\"zero\"/>\n";
187 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
192 out<<
"<position name=\"PosString"<<module->getString()<<
"_Dom"<<module->getID()<<
"\" unit=\"m\" x=\""<<module->
getX()<<
"\" y=\""<<module->getY()<<
"\" z=\""<<module->getZ()<<
"\"/>\n";
194 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
198 const JVector3D vec = static_cast<JVector3D>(*pmt).
sub(center);
199 out<<
"<position name=\"CathodPosition"<<pmt->getID()<<
"_"<<module->getID()<<
"\" unit=\"m\" x=\""<<vec.
getX()<<
"\" y=\""<<vec.
getY()<<
"\" z=\""<<vec.
getZ()<<
"\"/>\n";
200 out<<
"<rotation name=\"CathodRotation"<<pmt->getID()<<
"_"<<module->getID()<<
"\" unit=\"deg\" x=\""<<
GetXrotationG4(*pmt)<<
"\" y=\""<<
GetYrotationG4(*pmt)<<
"\" z=\"0.000000\"/>\n";
201 out<<
"<constant name=\"CathodID_"<<PMTs_NO<<
"\" value=\""<<pmt->getID()<<
"\"/>\n";
207 out<<
"<position name=\"OMShift\" unit=\"m\" x=\"0\" y=\"0\" z=\"0.0392\"/>\n";
209 out<<
"<!-- end of DU position definitions -->\n<position name=\"CrustPosition\" unit=\"m\" x=\"0\" y=\"0\" z=\""<<Crust_Z_position<<
"\"/>\n";
212 out<<
"<materials>\n";
213 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";
229 out<<
"<!-- OMVolume(s) construction -->\n";
231 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
233 out<<
" <volume name=\"OMVolume"<<module->getID()<<
"\">\n";
234 out<<
" <materialref ref=\"Water\"/>\n";
235 out<<
" <solidref ref=\"OMSphere\"/>\n";
237 for (JModule::const_iterator pmt = module->begin(); pmt != module->end(); ++pmt) {
239 out<<
" <volumeref ref=\"CathodVolume\"/>\n";
240 out<<
" <positionref ref=\"CathodPosition"<<pmt->getID()<<
"_"<<module->getID()<<
"\"/>\n";
241 out<<
" <rotationref ref=\"CathodRotation"<<pmt->getID()<<
"_"<<module->getID()<<
"\"/>\n";
242 out<<
" </physvol>\n";
248 out<<
"<!-- StoreyVolume(s) construction -->\n";
250 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
251 out<<
" <volume name=\"StoreyVolume"<<module->getID()<<
"\">\n";
252 out<<
" <materialref ref=\"Water\"/>\n";
253 out<<
" <solidref ref=\"StoreyBox\"/>\n";
255 out<<
" <volumeref ref=\"OMVolume"<<module->getID()<<
"\"/>\n";
256 out<<
" <positionref ref=\"OMShift\"/>\n";
258 out<<
" <rotationref ref=\"identity\"/>\n";
259 out<<
" </physvol>\n";
263 out<<
"<!-- Crust Volume construction-->\n";
264 out<<
"<volume name=\"CrustVolume\">\n";
265 out<<
" <materialref ref=\"Crust\"/>\n";
266 out<<
" <solidref ref=\"CrustBox\"/>\n";
269 out<<
"<!-- World Volume construction-->\n";
270 out<<
"<volume name=\"WorldVolume\">\n";
271 out<<
" <materialref ref=\"Water\"/>\n";
273 out<<
" <solidref ref=\"WorldTube\"/>\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) {
283 out<<
" <volumeref ref=\"StoreyVolume"<<module->getID()<<
"\"/>\n";
284 out<<
" <positionref ref=\"PosString"<< module->getString() <<
"_Dom"<< module->getID() <<
"\"/>\n";
286 out<<
" <rotationref ref=\"identity\"/>\n";
287 out<<
" </physvol>\n";
292 out<<
"</structure>\n";
293 out<<
"<setup name=\"Default\" version=\"1.0\">\n";
294 out<<
"<world ref=\"WorldVolume\"/>\n";
341 if (!module.empty()) {
345 for (JModule::const_iterator pmt = module.begin(); pmt != module.end(); ++pmt) {
370 return module.size();
382 int number_of_pmts = 0;
384 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
388 return number_of_pmts;
402 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
403 buffer.insert(module->getString());
406 return buffer.size();
420 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
421 buffer.insert(module->getString());
438 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
439 buffer.insert(module->getFloor());
442 return buffer.size();
456 for (JDetector::const_iterator module =
detector.begin(); module !=
detector.end(); ++module) {
457 buffer.insert(module->getID());
460 return buffer.size();
485 ifstream in(file_name.c_str());
511 ifstream in(file_name.c_str());
567 std::ofstream out(file_name.c_str());
583 std::ofstream out(file_name.c_str());
615 module.resize(memo.size());
678 if (module.empty()) {
682 const double R = 0.5;
684 const double st = sin(0.75*
PI);
685 const double ct = cos(0.75*
PI);
687 for (
int i = 0; i != 3; ++i) {
689 const double phi = (2.0*
PI*i) / 3.0;
690 const double cp = cos(phi);
691 const double sp = sin(phi);
714 os <<
"(" << setfill(
'0') << setw(3) << location.
getString() <<
"," << setfill(
'0') << setw(2) << location.
getFloor() <<
")";
727 static const size_t NUMBER_OF_DIMENSIONS = 3;
741 if (first.size() == second.size()) {
743 const size_t N = first.size();
745 if (N >= NUMBER_OF_DIMENSIONS) {
750 for (
size_t i = 0; i != N; ++i) {
755 for (
size_t i = 0; i != NUMBER_OF_DIMENSIONS; ++i) {
756 if (!orthonormalise(i)) {
761 this->a00 = out[0].getX() * in[0].getX() + out[1].getX() * in[1].getX() + out[2].getX() * in[2].getX();
762 this->a01 = out[0].getX() * in[0].getY() + out[1].getX() * in[1].getY() + out[2].getX() * in[2].getY();
763 this->a02 = out[0].getX() * in[0].getZ() + out[1].getX() * in[1].getZ() + out[2].getX() * in[2].getZ();
765 this->a10 = out[0].getY() * in[0].getX() + out[1].getY() * in[1].getX() + out[2].getY() * in[2].getX();
766 this->a11 = out[0].getY() * in[0].getY() + out[1].getY() * in[1].getY() + out[2].getY() * in[2].getY();
767 this->a12 = out[0].getY() * in[0].getZ() + out[1].getY() * in[1].getZ() + out[2].getY() * in[2].getZ();
769 this->a20 = out[0].getZ() * in[0].getX() + out[1].getZ() * in[1].getX() + out[2].getZ() * in[2].getX();
770 this->a21 = out[0].getZ() * in[0].getY() + out[1].getZ() * in[1].getY() + out[2].getZ() * in[2].getY();
771 this->a22 = out[0].getZ() * in[0].getZ() + out[1].getZ() * in[1].getZ() + out[2].getZ() * in[2].getZ();
775 THROW(
JException,
"Module size" << N <<
" < " << NUMBER_OF_DIMENSIONS);
780 THROW(
JException,
"Module size" << first.size() <<
" != " << second.size());
800 for (
size_t i = index + 1; i != in.size(); ++i) {
801 if (in[i].getLengthSquared() > in[pos].getLengthSquared()) {
806 const double u = in[pos].getLength();
814 swap(in [pos], in [index]);
815 swap(out[pos], out[index]);
818 for (
size_t i = index + 1; i != in.size(); ++i) {
820 const double dot = in[index].getDot(in[i]);
822 in [i] -= dot * in [index];
823 out[i] -= dot * out[index];