Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:14:32 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!