DDG4/src/Geant4Helpers.cpp

0001 //==========================================================================
0002 //  AIDA Detector description implementation
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 // Author     : M.Frank
0011 //
0012 //==========================================================================
0013
0014 // Framework include files
0015 #include <DD4hep/Printout.h>
0016 #include <DD4hep/DD4hepUnits.h>
0017 #include <DDG4/Geant4Helpers.h>
0018
0019 #include <CLHEP/Units/SystemOfUnits.h>
0020
0021 // ROOT include files
0022 #include <TGeoMatrix.h>
0023
0024 namespace  {
0025   static constexpr const double CM_2_MM = (CLHEP::centimeter/dd4hep::centimeter);
0026
0027   /// Overload to access protected constructor
0028   struct MyTransform3D : public G4Transform3D {
0029 #if 0
0030     MyTransform3D(double XX, double XY, double XZ, double DX,
0031           double YX, double YY, double YZ, double DY,
0032           double ZX, double ZY, double ZZ, double DZ)
0033       : G4Transform3D(XX, XY, XZ, DX, YX, YY, YZ, DY, ZX, ZY, ZZ, DZ) {
0034     }
0035 #endif
0036     MyTransform3D(const double* t, const double* r)
0037       : G4Transform3D(r[0],r[1],r[2],t[0]*CM_2_MM,r[3],r[4],r[5],t[1]*CM_2_MM,r[6],r[7],r[8],t[2]*CM_2_MM)  {
0038     }
0039     MyTransform3D(const double* t)
0040       : G4Transform3D(1.0, 0.0, 0.0, t[0]*CM_2_MM, 0.0, 1.0, 0.0, t[1]*CM_2_MM, 0.0, 0.0, 1.0, t[2]*CM_2_MM)  {
0041     }
0042     MyTransform3D(Transform3D&& copy) : Transform3D(std::move(copy)) {}
0043   };
0044   /// Overload to access protected constructor
0045   class MyG4RotationMatrix : public G4RotationMatrix   {
0046   public:
0047     MyG4RotationMatrix() : G4RotationMatrix() {}
0048     MyG4RotationMatrix(const double* r) : G4RotationMatrix(r[0],r[1],r[2],r[3],r[4],r[5],r[6],r[7],r[8])  {
0049     }
0050   };
0051 }
0052
0053 G4RotationMatrix dd4hep::sim::g4Rotation(const TGeoMatrix* matrix)   {
0054   return matrix->IsRotation() ? MyG4RotationMatrix(matrix->GetRotationMatrix()) : MyG4RotationMatrix();
0055 }
0056
0057 G4RotationMatrix dd4hep::sim::g4Rotation(const TGeoMatrix& matrix)   {
0058   return g4Rotation(&matrix);
0059 }
0060
0061 G4RotationMatrix dd4hep::sim::g4Rotation(const TGeoRotation* matrix)   {
0062   return matrix->IsRotation() ? MyG4RotationMatrix(matrix->GetRotationMatrix()) : MyG4RotationMatrix();
0063 }
0064
0065 G4RotationMatrix dd4hep::sim::g4Rotation(const TGeoRotation& matrix)   {
0066   return g4Rotation(&matrix);
0067 }
0068
0069 G4RotationMatrix dd4hep::sim::g4Rotation(const Rotation3D& rot)   {
0070   double r[9];
0071   rot.GetComponents(r, r+9);
0072   return MyG4RotationMatrix(r);
0073 }
0074
0075 G4RotationMatrix dd4hep::sim::g4Rotation(const RotationZYX& rot)   {
0076   return g4Rotation(Rotation3D(rot));
0077 }
0078
0079 G4Transform3D dd4hep::sim::g4Transform(const double translation[])    {
0080   return MyTransform3D(translation);
0081 }
0082
0083 G4Transform3D dd4hep::sim::g4Transform(const double* translation, const double* rotation)    {
0084   return MyTransform3D(translation, rotation);
0085 }
0086
0087 void dd4hep::sim::g4Transform(const double translation[], G4Transform3D& transform)    {
0088   transform = MyTransform3D(translation);
0089 }
0090
0091 void dd4hep::sim::g4Transform(const double* translation, const double* rotation, G4Transform3D& transform)    {
0092   transform = MyTransform3D(translation, rotation);
0093 }
0094
0095 G4Transform3D dd4hep::sim::g4Transform(const TGeoMatrix& matrix)   {
0096   return g4Transform(&matrix);
0097 }
0098
0099 G4Transform3D dd4hep::sim::g4Transform(const TGeoMatrix* matrix)   {
0100   return matrix->IsRotation()
0101     ? g4Transform(matrix->GetTranslation(), matrix->GetRotationMatrix())
0102     : g4Transform(matrix->GetTranslation());
0103 }
0104
0105 void dd4hep::sim::g4Transform(const TGeoMatrix* matrix, G4Transform3D& transform)   {
0106   matrix->IsRotation()
0107     ? g4Transform(matrix->GetTranslation(), matrix->GetRotationMatrix(), transform)
0108     : g4Transform(matrix->GetTranslation(), transform);
0109 }
0110
0111 void dd4hep::sim::g4Transform(const TGeoMatrix& matrix, G4Transform3D& transform)   {
0112   g4Transform(&matrix, transform);
0113 }
0114
0115 G4Transform3D dd4hep::sim::g4Transform(const TGeoTranslation& translation, const TGeoRotation& rotation)   {
0116   return g4Transform(&translation, &rotation);
0117 }
0118
0119 G4Transform3D dd4hep::sim::g4Transform(const TGeoTranslation* translation, const TGeoRotation* rotation)   {
0120   return rotation->IsRotation()
0121     ? g4Transform(translation->GetTranslation(), rotation->GetRotationMatrix())
0122     : g4Transform(translation->GetTranslation());
0123 }
0124
0125 G4Transform3D dd4hep::sim::g4Transform(const Position& pos, const Rotation3D& rot)   {
0126   double r[9], t[3] = {pos.X(), pos.Y(), pos.Z()};
0127   rot.GetComponents(r, r+9);
0128   return g4Transform(t, r);
0129 }
0130
0131 void dd4hep::sim::g4Transform(const Position& pos, const Rotation3D& rot, G4Transform3D& transform)   {
0132   double r[9], t[3] = {pos.X(), pos.Y(), pos.Z()};
0133   rot.GetComponents(r, r+9);
0134   g4Transform(t, r, transform);
0135 }
0136
0137 G4Transform3D dd4hep::sim::g4Transform(const Position& pos, const RotationZYX& rot)   {
0138   return g4Transform(pos, Rotation3D(rot));
0139 }
0140
0141 G4Transform3D dd4hep::sim::g4Transform(const Transform3D& matrix)   {
0142   Position pos;
0143   Rotation3D rot;
0144   matrix.GetDecomposition(rot, pos);
0145   return g4Transform(pos, rot);
0146 }
0147
0148 void dd4hep::sim::g4Transform(const Transform3D& matrix, G4Transform3D& transform)   {
0149   Position pos;
0150   Rotation3D rot;
0151   matrix.GetDecomposition(rot, pos);
0152   g4Transform(pos, rot, transform);
0153 }
0154
0155 /// Generate parameterised placements in 2 dimension according to transformation delta
0156 G4Transform3D
0157 dd4hep::sim::generate_placements(const G4Transform3D& start,
0158                  const G4Transform3D& delta,
0159                  std::size_t count,
0160                  const std::function<void(const G4Transform3D& delta)>& callback)
0161 {
0162   G4Transform3D transform(start);
0163   for( std::size_t i = 0; i < count; ++i )    {
0164     callback(transform);
0165     transform = transform * delta;
0166   }
0167   return transform;
0168 }
0169
0170 /// Generate parameterised placements in 2 dimensions according to transformation delta
0171 G4Transform3D
0172 dd4hep::sim::generate_placements(const G4Transform3D& start,
0173                  const G4Transform3D& delta1,
0174                  std::size_t count1,
0175                  const G4Transform3D& delta2,
0176                  std::size_t count2,
0177                  const std::function<void(const G4Transform3D& delta)>& callback)
0178 {
0179   G4Transform3D transform2 = start;
0180   for(std::size_t j = 0; j < count2; ++j)    {
0181     G4Transform3D transform1 = transform2;
0182     generate_placements(transform1, delta1, count1, callback);
0183     transform2 = transform2 * delta2;
0184   }
0185   return transform2;
0186 }
0187
0188 /// Generate parameterised placements in 3 dimensions according to transformation delta
0189 G4Transform3D
0190 dd4hep::sim::generate_placements(const G4Transform3D& start,
0191                  const G4Transform3D& delta1,
0192                  std::size_t count1,
0193                  const G4Transform3D& delta2,
0194                  std::size_t count2,
0195                  const G4Transform3D& delta3,
0196                  std::size_t count3,
0197                  const std::function<void(const G4Transform3D& delta)>& callback)
0198 {
0199   G4Transform3D transform3 = start;
0200   for(std::size_t k = 0; k < count3; ++k)    {
0201     G4Transform3D transform2 = transform3;
0202     generate_placements(transform2, delta1, count1, delta2, count2, callback);
0203     transform3 = transform3 * delta3;
0204   }
0205   return transform3;
0206 }
0207
0208 std::pair<double, EAxis> dd4hep::sim::extract_axis(const Transform3D& trafo)    {
0209   constexpr double eps = 2e-14;
0210   G4Transform3D tr = g4Transform(trafo);
0211   if ( fabs(tr.xx()) > (1e0-eps) &&
0212        fabs(tr.yy()) > (1e0-eps) &&
0213        fabs(tr.zz()) > (1e0-eps) )   {
0214     if ( fabs(tr.dx()) > eps )
0215       return { tr.dx(), kYAxis };
0216     else if ( fabs(tr.dy()) > eps )
0217       return { tr.dy(), kYAxis };
0218     else if ( fabs(tr.dz()) > eps )
0219       return { tr.dz(), kZAxis };
0220   }
0221   else if ( tr.getTranslation().mag() > eps && fabs(tr.dz()) < eps )  {
0222     return { tr.getTranslation().rho(), kRho };
0223   }
0224   else if ( fabs(tr.dz()) < eps )  {
0225     return { tr.getRotation().phi(), kPhi };   // Is this correct ?
0226   }
0227   except("Geant4Converter","Invalid volume parametrization matrix. Unknown Axis!");
0228   return { 0e0, kUndefined };
0229 }