GeomBuilder/src/DD4hep_volumeHandle.cc

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author N. Amapane - INFN Torino (original developer)
0005  */
0006
0007 #include "DD4hep_volumeHandle.h"
0008
0009 #include "DataFormats/GeometrySurface/interface/Plane.h"
0010 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0011 #include "DataFormats/GeometrySurface/interface/Cone.h"
0012 #include "DataFormats/GeometryVector/interface/CoordinateSets.h"
0013 #include "DataFormats/Math/interface/GeantUnits.h"
0014 #include "DataFormats/Math/interface/Vector3D.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016
0017 #include <DD4hep/DD4hepUnits.h>
0018
0019 #include <string>
0020 #include <iterator>
0021
0022 using namespace SurfaceOrientation;
0023 using namespace std;
0024 using namespace magneticfield;
0025 using namespace edm;
0026
0027 using DDBox = dd4hep::Box;
0028 using DDTrap = dd4hep::Trap;
0029 using DDTubs = dd4hep::Tube;
0030 using DDCons = dd4hep::ConeSegment;
0031 using DDTruncTubs = dd4hep::TruncatedTube;
0032
0033 volumeHandle::volumeHandle(const cms::DDFilteredView &fv, bool expand2Pi, bool debugVal)
0034     : BaseVolumeHandle(expand2Pi, debugVal), theShape(fv.legacyShape(fv.shape())), solid(fv) {
0035   name = fv.name();
0036   copyno = fv.copyNum();
0037   const auto *const transArray = fv.trans();
0038
0039   // Convert from DD4hep units to cm
0040   center_ = GlobalPoint(transArray[0] / dd4hep::cm, transArray[1] / dd4hep::cm, transArray[2] / dd4hep::cm);
0041
0042   // ASSUMPTION: volume names ends with "_NUM" where NUM is the volume number
0043   string volName = name;
0044   volName.erase(0, volName.rfind('_') + 1);
0045   volumeno = static_cast<unsigned short>(std::atoi(volName.c_str()));
0046
0047   for (int i = 0; i < 6; ++i) {
0048     isAssigned[i] = false;
0049   }
0050   referencePlane(fv);
0051   switch (theShape) {
0052     case DDSolidShape::ddbox: {
0053       DDBox box(solid.solid());
0054       // Convert from DD4hep units to cm
0055       double halfX = box.x() / dd4hep::cm;
0056       double halfY = box.y() / dd4hep::cm;
0057       double halfZ = box.z() / dd4hep::cm;
0058       buildBox(halfX, halfY, halfZ);
0059     } break;
0060     case DDSolidShape::ddtrap: {
0061       DDTrap trap(solid.solid());
0062       double x1 = trap.bottomLow1() / dd4hep::cm;
0063       double x2 = trap.topLow1() / dd4hep::cm;
0064       double x3 = trap.bottomLow2() / dd4hep::cm;
0065       double x4 = trap.topLow2() / dd4hep::cm;
0066       double y1 = trap.high1() / dd4hep::cm;
0067       double y2 = trap.high2() / dd4hep::cm;
0068       double theta = trap.theta();
0069       double phi = trap.phi();
0070       double halfZ = trap.dZ() / dd4hep::cm;
0071       double alpha1 = trap.alpha1();
0072       double alpha2 = trap.alpha2();
0073       buildTrap(x1, x2, x3, x4, y1, y2, theta, phi, halfZ, alpha1, alpha2);
0074     } break;
0075
0076     case DDSolidShape::ddcons: {
0077       DDCons cons(solid.solid());
0078       double zhalf = cons.dZ() / dd4hep::cm;
0079       double rInMinusZ = cons.rMin1() / dd4hep::cm;
0080       double rOutMinusZ = cons.rMax1() / dd4hep::cm;
0081       double rInPlusZ = cons.rMin2() / dd4hep::cm;
0082       double rOutPlusZ = cons.rMax2() / dd4hep::cm;
0083       double startPhi = cons.startPhi();
0084       double deltaPhi = reco::deltaPhi(cons.endPhi(), startPhi);
0085       buildCons(zhalf, rInMinusZ, rOutMinusZ, rInPlusZ, rOutPlusZ, startPhi, deltaPhi);
0086     } break;
0087     case DDSolidShape::ddtubs: {
0088       DDTubs tubs(solid.solid());
0089       double zhalf = tubs.dZ() / dd4hep::cm;
0090       double rIn = tubs.rMin() / dd4hep::cm;
0091       double rOut = tubs.rMax() / dd4hep::cm;
0092       double startPhi = tubs.startPhi();
0093       double deltaPhi = tubs.endPhi() - startPhi;
0094       buildTubs(zhalf, rIn, rOut, startPhi, deltaPhi);
0095     } break;
0096     case DDSolidShape::ddpseudotrap: {
0097       vector<double> d = solid.parameters();
0098
0099       // The pseudo-trapezoid parameters are:
0100       // d[0] -- x1
0101       // d[1] -- x2
0102       // d[2] -- y1
0103       // d[3] -- y2
0104       // d[4] -- halfZ
0105       // d[5] -- radius
0106       // d[6] -- atMinusZ (0 or 1)
0107       // Note all are lengths except for the last one. The lengths come from
0108       // DD4hep in DD4hep units and must be converted to cm for use.
0109
0110       if (d.size() >= 7)
0111         LogTrace("MagGeoBuilder") << " Pseudo trap params raw =  " << d[0] << ", " << d[1] << ", " << d[2] << ", "
0112                                   << d[3] << ", " << d[4] << ", " << d[5] << ", " << d[6];
0113
0114       // Convert all but last parameter to cm (last one is a boolean).
0115       transform(d.begin(), --(d.end()), d.begin(), [](double val) { return val / dd4hep::cm; });
0116
0117       if (d.size() >= 7)
0118         LogTrace("MagGeoBuilder") << " Pseudo trap params converted =  " << d[0] << ", " << d[1] << ", " << d[2] << ", "
0119                                   << d[3] << ", " << d[4] << ", " << d[5] << ", " << d[6];
0120
0121       buildPseudoTrap(d[0], d[1], d[2], d[3], d[4], d[5], static_cast<bool>(d[6]));
0122     } break;
0123     case DDSolidShape::ddtrunctubs: {
0124       DDTruncTubs tubs(solid.solid());
0125       double zhalf = tubs.dZ() / dd4hep::cm;               // half of the z-Axis
0126       double rIn = tubs.rMin() / dd4hep::cm;               // inner radius
0127       double rOut = tubs.rMax() / dd4hep::cm;              // outer radius
0128       double startPhi = tubs.startPhi();                   // angular start of the tube-section
0129       double deltaPhi = tubs.deltaPhi();                   // angular span of the tube-section
0130       double cutAtStart = tubs.cutAtStart() / dd4hep::cm;  // truncation at begin of the tube-section
0131       double cutAtDelta = tubs.cutAtDelta() / dd4hep::cm;  // truncation at end of the tube-section
0132       bool cutInside = tubs.cutInside();  // true, if truncation is on the inner side of the tube-section
0133       buildTruncTubs(zhalf, rIn, rOut, startPhi, deltaPhi, cutAtStart, cutAtDelta, cutInside);
0134     } break;
0135     default:
0136       LogError("magneticfield::volumeHandle")
0137           << "ctor: Unexpected shape # " << static_cast<int>(theShape) << " for vol " << name;
0138   }
0139
0140   // The only materials used in the geometry are: materials:Air, d=0.001214; materials:Iron, d=7.87
0141   if (fv.volume().material().density() > 3.)
0142     isIronFlag = true;
0143
0144   if (debug) {
0145     LogTrace("MagGeoBuilder") << " RMin =  " << theRMin << newln << " RMax =  " << theRMax;
0146
0147     if (theRMin < 0 || theRN < theRMin || theRMax < theRN)
0148       LogTrace("MagGeoBuilder") << "*** WARNING: wrong RMin/RN/RMax";
0149
0150     LogTrace("MagGeoBuilder") << "Summary: " << name << " " << copyno << " shape = " << theShape << " trasl "
0151                               << center() << " R " << center().perp() << " phi " << center().phi() << " magFile "
0152                               << magFile << " Material= " << fv.materialName() << " isIron= " << isIronFlag
0153                               << " masterSector= " << masterSector;
0154
0155     LogTrace("MagGeoBuilder") << " Orientation of surfaces:";
0156     std::string sideName[3] = {"positiveSide", "negativeSide", "onSurface"};
0157     for (int i = 0; i < 6; ++i) {
0158       if (surfaces[i] != nullptr)
0159         LogTrace("MagGeoBuilder") << "  " << i << ":" << sideName[surfaces[i]->side(center_, 0.3)];
0160     }
0161   }
0162 }
0163
0164 void volumeHandle::referencePlane(const cms::DDFilteredView &fv) {
0165   // The refPlane is the "main plane" for the solid. It corresponds to the
0166   // x,y plane in the DDD local frame, and defines a frame where the local
0167   // coordinates are the same as in DDD.
0168   // In the geometry version 85l_030919, this plane is normal to the
0169   // beam line for all volumes but pseudotraps, so that global R is along Y,
0170   // global phi is along -X and global Z along Z:
0171   //
0172   //   Global(for vol at pi/2)    Local
0173   //   +R (+Y)                    +Y
0174   //   +phi(-X)                   -X
0175   //   +Z                         +Z
0176   //
0177   // For pseudotraps the refPlane is parallel to beam line and global R is
0178   // along Z, global phi is along +-X and and global Z along Y:
0179   //
0180   //   Global(for vol at pi/2)    Local
0181   //   +R (+Y)                    +Z
0182   //   +phi(-X)                   +X
0183   //   +Z                         +Y
0184   //
0185   // Note that the frame is centered in the DDD volume center, which is
0186   // inside the volume for DDD boxes and (pesudo)trapezoids, on the beam line
0187   // for tubs, cons and trunctubs.
0188
0189   // In geometry version 1103l, trapezoids have X and Z in the opposite direction
0190   // than the above.  Boxes are either oriented as described above or in some case
0191   // have opposite direction for Y and X.
0192
0193   // The global position
0194   Surface::PositionType &posResult = center_;
0195
0196   // The reference plane rotation
0197   math::XYZVector x, y, z;
0198   dd4hep::Rotation3D refRot;
0199   fv.rot(refRot);
0200   refRot.GetComponents(x, y, z);
0201   if (debug) {
0202     if (x.Cross(y).Dot(z) < 0.5) {
0203       LogTrace("MagGeoBuilder") << "*** WARNING: Rotation is not RH ";
0204     }
0205   }
0206
0207   // The global rotation
0208   Surface::RotationType rotResult(float(x.X()),
0209                                   float(x.Y()),
0210                                   float(x.Z()),
0211                                   float(y.X()),
0212                                   float(y.Y()),
0213                                   float(y.Z()),
0214                                   float(z.X()),
0215                                   float(z.Y()),
0216                                   float(z.Z()));
0217
0218   refPlane = new GloballyPositioned<float>(posResult, rotResult);
0219
0220   // Check correct orientation
0221   if (debug) {
0222     LogTrace("MagGeoBuilder") << "Refplane pos  " << refPlane->position();
0223
0224     // See comments above for the conventions for orientation.
0225     LocalVector globalZdir(0., 0., 1.);  // Local direction of the axis along global Z
0226
0227     if (theShape == DDSolidShape::ddpseudotrap) {
0228       globalZdir = LocalVector(0., 1., 0.);
0229     }
0230
0231     if (refPlane->toGlobal(globalZdir).z() < 0.) {
0232       globalZdir = -globalZdir;
0233     }
0234     float chk = refPlane->toGlobal(globalZdir).dot(GlobalVector(0, 0, 1));
0235     if (chk < .999)
0236       LogTrace("MagGeoBuilder") << "*** WARNING RefPlane check failed!***" << chk;
0237   }
0238 }
0239
0240 std::vector<VolumeSide> volumeHandle::sides() const {
0241   std::vector<VolumeSide> result;
0242   for (int i = 0; i < 6; ++i) {
0243     // If this is just a master volume out of wich a 2pi volume
0244     // should be built (e.g. central cylinder), skip the phi boundaries.
0245     if (expand && (i == phiplus || i == phiminus))
0246       continue;
0247
0248     // FIXME: Skip null inner degenerate cylindrical surface
0249     if (theShape == DDSolidShape::ddtubs && i == SurfaceOrientation::inner && theRMin < 0.001)
0250       continue;
0251
0252     ReferenceCountingPointer<Surface> s = const_cast<Surface *>(surfaces[i].get());
0253     result.push_back(VolumeSide(s, GlobalFace(i), surfaces[i]->side(center_, 0.3)));
0254   }
0255   return result;
0256 }
0257
0258 #include "buildBox.icc"
0259 #include "buildTrap.icc"
0260 #include "buildTubs.icc"
0261 #include "buildCons.icc"
0262 #include "buildPseudoTrap.icc"
0263 #include "buildTruncTubs.icc"