File indexing completed on 2025-01-18 09:14:42
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015 #include "DD4hep/Printout.h"
0016 #include "DD4hep/MatrixHelpers.h"
0017 #include "DD4hep/DetFactoryHelper.h"
0018 #include "XML/VolumeBuilder.h"
0019 #include "XML/Utilities.h"
0020 #include <iostream>
0021 #include <cmath>
0022
0023 using namespace std;
0024 using namespace dd4hep;
0025 using namespace dd4hep::detail;
0026
0027 namespace {
0028 Rotation3D makeRotReflect(double thetaX, double phiX, double thetaY, double phiY, double thetaZ, double phiZ) {
0029
0030 Position x(cos(phiX) * sin(thetaX), sin(phiX) * sin(thetaX), cos(thetaX));
0031 Position y(cos(phiY) * sin(thetaY), sin(phiY) * sin(thetaY), cos(thetaY));
0032 Position z(cos(phiZ) * sin(thetaZ), sin(phiZ) * sin(thetaZ), cos(thetaZ));
0033 printout(INFO,"ReflectionMatrices",
0034 "+++ Adding rotation: (theta/phi)[rad] X: %6.3f %6.3f Y: %6.3f %6.3f Z: %6.3f %6.3f",
0035 thetaX,phiX,thetaY,phiY,thetaZ,phiZ);
0036 Rotation3D rotation(x.x(), y.x(), z.x(), x.y(), y.y(), z.y(), x.z(), y.z(), z.z());
0037 return rotation;
0038 }
0039 Transform3D transform_3d(xml_h element) {
0040 xml_dim_t xrot(element.child(_U(rotation)));
0041 xml_dim_t xpos(element.child(_U(position)));
0042 double thetaX = xrot.attr<double>(Unicode("thetaX"));
0043 double phiX = xrot.attr<double>(Unicode("phiX"));
0044 double thetaY = xrot.attr<double>(Unicode("thetaY"));
0045 double phiY = xrot.attr<double>(Unicode("phiY"));
0046 double thetaZ = xrot.attr<double>(Unicode("thetaZ"));
0047 double phiZ = xrot.attr<double>(Unicode("phiZ"));
0048 printout(INFO, "ReflectionMatrices",
0049 "+++ Adding reflection rotation \"%s\": "
0050 "(theta/phi)[rad] X: %6.3f %6.3f Y: %6.3f %6.3f Z: %6.3f %6.3f",
0051 element.attr<string>(_U(name)).c_str(), thetaX, phiX, thetaY, phiY, thetaZ, phiZ);
0052 Rotation3D rot = makeRotReflect(thetaX, phiX, thetaY, phiY, thetaZ, phiZ);
0053 Position pos = Position(xpos.x(),xpos.y(),xpos.z());
0054 return Transform3D(rot, pos);
0055 }
0056
0057 struct ReflectionMatrices : public xml::tools::VolumeBuilder {
0058 int num_left = 0, num_right = 0;
0059
0060 using VolumeBuilder::VolumeBuilder;
0061
0062 void dump(const string& nam, const Transform3D& tr) {
0063 stringstream str;
0064 Position x, y, z, p;
0065 Rotation3D r;
0066 tr.GetDecomposition(r, p);
0067 r.GetComponents(x,y,z);
0068 double det1 = (x.Cross(y)).Dot(z);
0069 double det2 = detail::matrix::determinant(r);
0070 double det3 = detail::matrix::determinant(tr);
0071 str << "+++ Using rotation: " << nam << " "
0072 << (const char*)((det1 >= 0) ? "RIGHT" : "LEFT") << "-handed "
0073 << " Determinant: " << det1 << " " << det2 << " " << det3 << endl
0074 << "Pos: " << p << " Rotation:" << x << " " << y << " " << z
0075 << r;
0076 PrintLevel lvl = det1 != det2 || det1 != det3 ? ERROR : INFO;
0077 if ( det1 < 0e0 ) ++num_left;
0078 if ( det1 > 0e0 ) ++num_right;
0079 if ( det1 < 0e0 && nam.find("LEFT-handed") == string::npos ) {
0080 lvl = ERROR;
0081 }
0082 printout(lvl, "ReflectionMatrices",str.str().c_str());
0083 }
0084
0085 void dump(xml_elt_t elt, char reflection) {
0086 Position pos;
0087 RotationZYX rot;
0088 xml_dim_t e(elt);
0089 xml_dim_t xpos = e.child(_U(position), false);
0090 xml_dim_t xrot = e.child(_U(rotation), false);
0091 if ( !xpos.ptr() ) xpos = e;
0092 if ( xpos.ptr() ) pos = Position(xpos.x(0),xpos.y(0),xpos.z(0));
0093 if ( xrot.ptr() ) rot = RotationZYX(xrot.x(0),xrot.y(0),xrot.z(0));
0094
0095 TGeoHMatrix mat;
0096 detail::matrix::_transform(mat, pos, rot);
0097 switch(reflection) {
0098 case 'X':
0099 mat.ReflectX(kTRUE, kTRUE);
0100 break;
0101 case 'Y':
0102 mat.ReflectY(kTRUE, kTRUE);
0103 break;
0104 case 'Z':
0105 default:
0106 mat.ReflectZ(kTRUE, kTRUE);
0107 break;
0108 }
0109 dump(e.nameStr(), detail::matrix::_transform(mat));
0110 }
0111 DetElement create() {
0112 xml_dim_t x_box(x_det.dimensions());
0113 double bx = x_box.x();
0114 double by = x_box.y();
0115 double bz = x_box.z();
0116
0117 for(xml_coll_t c(x_det,_U(reflect_x)); c; ++c)
0118 dump(c, 'X');
0119 for(xml_coll_t c(x_det,_U(reflect_y)); c; ++c)
0120 dump(c, 'Y');
0121 for(xml_coll_t c(x_det,_U(reflect_z)); c; ++c)
0122 dump(c, 'Z');
0123 for(xml_coll_t c(x_det,_U(reflect)); c; ++c)
0124 dump(xml_dim_t(c).nameStr(), transform_3d(c));
0125
0126 auto pv = placeDetector(Volume("envelope",Box(bx, by, bz),description.air()), x_det);
0127 pv.addPhysVolID("system",x_det.id());
0128 return detector;
0129 }
0130 };
0131
0132 Ref_t create_test(Detector& description, xml_dim_t x_det, SensitiveDetector sens) {
0133 ReflectionMatrices builder(description, x_det, sens);
0134 DetElement det = builder.create();
0135 printout(ALWAYS,"ReflectionMatrices",
0136 "+++ Analysed %d right handed and %d left handed matrices.",
0137 builder.num_right, builder.num_left);
0138 return det;
0139 }
0140 }
0141
0142 DECLARE_DETELEMENT(DD4hep_Test_ReflectionMatrices,create_test)
