GeomBuilder/src/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 "volumeHandle.h"
0008 #include "DetectorDescription/Core/interface/DDExpandedView.h"
0009 #include "DetectorDescription/Core/interface/DDTranslation.h"
0010 #include "DetectorDescription/Core/interface/DDSolid.h"
0011 #include "DetectorDescription/Core/interface/DDValue.h"
0012 #include "DetectorDescription/Core/interface/DDMaterial.h"
0013
0014 #include "DataFormats/GeometrySurface/interface/Plane.h"
0015 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0016 #include "DataFormats/GeometrySurface/interface/Cone.h"
0017 #include "DataFormats/GeometryVector/interface/CoordinateSets.h"
0018
0019 #include <CLHEP/Units/SystemOfUnits.h>
0020
0021 #include <string>
0022 #include <iterator>
0023 #include <iomanip>
0024 #include <iostream>
0025
0026 using namespace SurfaceOrientation;
0027 using namespace std;
0028
0029 #include "DataFormats/Math/interface/GeantUnits.h"
0030
0031 namespace {
0032   // Old DD returns lengths in mm, but CMS code uses cm
0033   template <class NumType>
0034   inline constexpr NumType convertUnits(NumType millimeters)  // Millimeters -> centimeters
0035   {
0036     return (geant_units::operators::convertMmToCm(millimeters));
0037   }
0038 }  // namespace
0039
0040 MagGeoBuilderFromDDD::volumeHandle::volumeHandle(const DDExpandedView &fv, bool expand2Pi, bool debugVal)
0041     : magneticfield::BaseVolumeHandle(expand2Pi, debugVal) {
0042   name = fv.logicalPart().name().name();
0043   copyno = fv.copyno();
0044   solid = fv.logicalPart().solid();
0045   center_ =
0046       GlobalPoint(fv.translation().x() / CLHEP::cm, fv.translation().y() / CLHEP::cm, fv.translation().z() / CLHEP::cm);
0047
0048   // ASSUMPTION: volume names ends with "_NUM" where NUM is the volume number
0049   string volName = name;
0050   volName.erase(0, volName.rfind('_') + 1);
0051   volumeno = std::stoul(volName);
0052
0053   for (int i = 0; i < 6; ++i) {
0054     isAssigned[i] = false;
0055   }
0056
0057   if (debug) {
0058     cout.precision(7);
0059   }
0060
0061   referencePlane(fv);
0062
0063   if (solid.shape() == DDSolidShape::ddbox) {
0064     DDBox box(solid);
0065     double halfX = convertUnits(box.halfX());
0066     double halfY = convertUnits(box.halfY());
0067     double halfZ = convertUnits(box.halfZ());
0068     buildBox(halfX, halfY, halfZ);
0069   } else if (solid.shape() == DDSolidShape::ddtrap) {
0070     DDTrap trap(solid);
0071     double x1 = convertUnits(trap.x1());
0072     double x2 = convertUnits(trap.x2());
0073     double x3 = convertUnits(trap.x3());
0074     double x4 = convertUnits(trap.x4());
0075     double y1 = convertUnits(trap.y1());
0076     double y2 = convertUnits(trap.y2());
0077     double theta = trap.theta();
0078     double phi = trap.phi();
0079     double halfZ = convertUnits(trap.halfZ());
0080     double alpha1 = trap.alpha1();
0081     double alpha2 = trap.alpha2();
0082     buildTrap(x1, x2, x3, x4, y1, y2, theta, phi, halfZ, alpha1, alpha2);
0083   } else if (solid.shape() == DDSolidShape::ddcons) {
0084     DDCons cons(solid);
0085     double zhalf = convertUnits(cons.zhalf());
0086     double rInMinusZ = convertUnits(cons.rInMinusZ());
0087     double rOutMinusZ = convertUnits(cons.rOutMinusZ());
0088     double rInPlusZ = convertUnits(cons.rInPlusZ());
0089     double rOutPlusZ = convertUnits(cons.rOutPlusZ());
0090     double startPhi = cons.phiFrom();
0091     double deltaPhi = cons.deltaPhi();
0092     buildCons(zhalf, rInMinusZ, rOutMinusZ, rInPlusZ, rOutPlusZ, startPhi, deltaPhi);
0093   } else if (solid.shape() == DDSolidShape::ddtubs) {
0094     DDTubs tubs(solid);
0095     double zhalf = convertUnits(tubs.zhalf());
0096     double rIn = convertUnits(tubs.rIn());
0097     double rOut = convertUnits(tubs.rOut());
0098     double startPhi = tubs.startPhi();
0099     double deltaPhi = tubs.deltaPhi();
0100     buildTubs(zhalf, rIn, rOut, startPhi, deltaPhi);
0101   } else if (solid.shape() == DDSolidShape::ddpseudotrap) {
0102     DDPseudoTrap ptrap(solid);
0103     double x1 = convertUnits(ptrap.x1());
0104     double x2 = convertUnits(ptrap.x2());
0105     double y1 = convertUnits(ptrap.y1());
0106     double y2 = convertUnits(ptrap.y2());
0107     double halfZ = convertUnits(ptrap.halfZ());
0108     double radius = convertUnits(ptrap.radius());
0109     bool atMinusZ = ptrap.atMinusZ();
0110     buildPseudoTrap(x1, x2, y1, y2, halfZ, radius, atMinusZ);
0111   } else if (solid.shape() == DDSolidShape::ddtrunctubs) {
0112     DDTruncTubs tubs(solid);
0113     double zhalf = convertUnits(tubs.zHalf());            // half of the z-Axis
0114     double rIn = convertUnits(tubs.rIn());                // inner radius
0115     double rOut = convertUnits(tubs.rOut());              // outer radius
0116     double startPhi = tubs.startPhi();                    // angular start of the tube-section
0117     double deltaPhi = tubs.deltaPhi();                    // angular span of the tube-section
0118     double cutAtStart = convertUnits(tubs.cutAtStart());  // truncation at begin of the tube-section
0119     double cutAtDelta = convertUnits(tubs.cutAtDelta());  // truncation at end of the tube-section
0120     bool cutInside = tubs.cutInside();  // true, if truncation is on the inner side of the tube-section
0121     buildTruncTubs(zhalf, rIn, rOut, startPhi, deltaPhi, cutAtStart, cutAtDelta, cutInside);
0122   } else {
0123     cout << "volumeHandle ctor: Unexpected solid: " << DDSolidShapesName::name(solid.shape()) << endl;
0124   }
0125
0126   // NOTE: Table name and master sector are no longer taken from xml!
0127   //   DDsvalues_type sv(fv.mergedSpecifics());
0128
0129   //   { // Extract the name of associated field file.
0130   //     std::vector<std::string> temp;
0131   //     std::string pname = "table";
0132   //     DDValue val(pname);
0133   //     DDsvalues_type sv(fv.mergedSpecifics());
0134   //     if (DDfetch(&sv,val)) {
0135   //       temp = val.strings();
0136   //       if (temp.size() != 1) {
0137   //    cout << "*** WARNING: volume has > 1 SpecPar " << pname << endl;
0138   //       }
0139   //       magFile = temp[0];
0140
0141   //       string find="[copyNo]";
0142   //       std::size_t j;
0143   //       for ( ; (j = magFile.find(find)) != string::npos ; ) {
0144   //    stringstream conv;
0145   //    conv << setfill('0') << setw(2) << copyno;
0146   //    string repl;
0147   //    conv >> repl;
0148   //    magFile.replace(j, find.length(), repl);
0149   //       }
0150
0151   //     } else {
0152   //       cout << "*** WARNING: volume does not have a SpecPar " << pname << endl;
0153   //       cout << " DDsvalues_type:  " << fv.mergedSpecifics() << endl;
0154   //     }
0155   //   }
0156
0157   //   { // Extract the number of the master sector.
0158   //     std::vector<double> temp;
0159   //     const std::string pname = "masterSector";
0160   //     DDValue val(pname);
0161   //     if (DDfetch(&sv,val)) {
0162   //       temp = val.doubles();
0163   //       if (temp.size() != 1) {
0164   //    cout << "*** WARNING: volume has > 1 SpecPar " << pname << endl;
0165   //       }
0166   //       masterSector = int(temp[0]+.5);
0167   //     } else {
0168   //       if (MagGeoBuilderFromDDD::debug) {
0169   //    cout << "Volume does not have a SpecPar " << pname
0170   //         << " using: " << copyno << endl;
0171   //    cout << " DDsvalues_type:  " << fv.mergedSpecifics() << endl;
0172   //       }
0173   //       masterSector = copyno;
0174   //     }
0175   //   }
0176
0177   // Get material for this volume
0178   if (fv.logicalPart().material().name().name() == "Iron")
0179     isIronFlag = true;
0180
0181   if (debug) {
0182     cout << " RMin =  " << theRMin << endl;
0183     cout << " RMax =  " << theRMax << endl;
0184
0185     if (theRMin < 0 || theRN < theRMin || theRMax < theRN)
0186       cout << "*** WARNING: wrong RMin/RN/RMax , shape: " << DDSolidShapesName::name(shape()) << endl;
0187
0188     cout << "Summary: " << name << " " << copyno << " Shape= " << DDSolidShapesName::name(shape()) << " trasl "
0189          << center() << " R " << center().perp() << " phi " << center().phi() << " magFile " << magFile
0190          << " Material= " << fv.logicalPart().material().name() << " isIron= " << isIronFlag
0191          << " masterSector= " << masterSector << std::endl;
0192
0193     cout << " Orientation of surfaces:";
0194     std::string sideName[3] = {"positiveSide", "negativeSide", "onSurface"};
0195     for (int i = 0; i < 6; ++i) {
0196       cout << "  " << i << ":" << sideName[surfaces[i]->side(center_, 0.3)];
0197     }
0198     cout << endl;
0199   }
0200 }
0201
0202 void MagGeoBuilderFromDDD::volumeHandle::referencePlane(const DDExpandedView &fv) {
0203   // The refPlane is the "main plane" for the solid. It corresponds to the
0204   // x,y plane in the DDD local frame, and defines a frame where the local
0205   // coordinates are the same as in DDD.
0206   // In the geometry version 85l_030919, this plane is normal to the
0207   // beam line for all volumes but pseudotraps, so that global R is along Y,
0208   // global phi is along -X and global Z along Z:
0209   //
0210   //   Global(for vol at pi/2)    Local
0211   //   +R (+Y)                    +Y
0212   //   +phi(-X)                   -X
0213   //   +Z                         +Z
0214   //
0215   // For pseudotraps the refPlane is parallel to beam line and global R is
0216   // along Z, global phi is along +-X and and global Z along Y:
0217   //
0218   //   Global(for vol at pi/2)    Local
0219   //   +R (+Y)                    +Z
0220   //   +phi(-X)                   +X
0221   //   +Z                         +Y
0222   //
0223   // Note that the frame is centered in the DDD volume center, which is
0224   // inside the volume for DDD boxes and (pesudo)trapezoids, on the beam line
0225   // for tubs, cons and trunctubs.
0226
0227   // In geometry version 1103l, trapezoids have X and Z in the opposite direction
0228   // than the above.  Boxes are either oriented as described above or in some case
0229   // have opposite direction for Y and X.
0230
0231   // The global position
0232   Surface::PositionType &posResult = center_;
0233
0234   // The reference plane rotation
0235   DD3Vector x, y, z;
0236   fv.rotation().GetComponents(x, y, z);
0237   if (debug) {
0238     if (x.Cross(y).Dot(z) < 0.5) {
0239       cout << "*** WARNING: Rotation is not RH " << endl;
0240     }
0241   }
0242
0243   // The global rotation
0244   Surface::RotationType rotResult(float(x.X()),
0245                                   float(x.Y()),
0246                                   float(x.Z()),
0247                                   float(y.X()),
0248                                   float(y.Y()),
0249                                   float(y.Z()),
0250                                   float(z.X()),
0251                                   float(z.Y()),
0252                                   float(z.Z()));
0253
0254   refPlane = new GloballyPositioned<float>(posResult, rotResult);
0255
0256   // Check correct orientation
0257   if (debug) {
0258     cout << "Refplane pos  " << refPlane->position() << endl;
0259
0260     // See comments above for the conventions for orientation.
0261     LocalVector globalZdir(0., 0., 1.);  // Local direction of the axis along global Z
0262     if (solid.shape() == DDSolidShape::ddpseudotrap) {
0263       globalZdir = LocalVector(0., 1., 0.);
0264     }
0265     if (refPlane->toGlobal(globalZdir).z() < 0.) {
0266       globalZdir = -globalZdir;
0267     }
0268
0269     float chk = refPlane->toGlobal(globalZdir).dot(GlobalVector(0, 0, 1));
0270     if (chk < .999)
0271       cout << "*** WARNING RefPlane check failed!***" << chk << endl;
0272   }
0273 }
0274
0275 std::vector<VolumeSide> MagGeoBuilderFromDDD::volumeHandle::sides() const {
0276   std::vector<VolumeSide> result;
0277   for (int i = 0; i < 6; ++i) {
0278     // If this is just a master volume out of wich a 2pi volume
0279     // should be built (e.g. central cylinder), skip the phi boundaries.
0280     if (expand && (i == phiplus || i == phiminus))
0281       continue;
0282
0283     // FIXME: Skip null inner degenerate cylindrical surface
0284     if (solid.shape() == DDSolidShape::ddtubs && i == SurfaceOrientation::inner && theRMin < 0.001)
0285       continue;
0286
0287     ReferenceCountingPointer<Surface> s = const_cast<Surface *>(surfaces[i].get());
0288     result.push_back(VolumeSide(s, GlobalFace(i), surfaces[i]->side(center_, 0.3)));
0289   }
0290   return result;
0291 }
0292
0293 using volumeHandle = MagGeoBuilderFromDDD::volumeHandle;
0294 using namespace magneticfield;
0295
0296 #include "buildBox.icc"
0297 #include "buildTrap.icc"
0298 #include "buildTubs.icc"
0299 #include "buildCons.icc"
0300 #include "buildPseudoTrap.icc"
0301 #include "buildTruncTubs.icc"