File indexing completed on 2023-03-17 11:14:32
0001
0002
0003
0004
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
0033 template <class NumType>
0034 inline constexpr NumType convertUnits(NumType millimeters)
0035 {
0036 return (geant_units::operators::convertMmToCm(millimeters));
0037 }
0038 }
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
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());
0113 double rIn = convertUnits(tubs.rIn());
0114 double rOut = convertUnits(tubs.rOut());
0115 double startPhi = tubs.startPhi();
0116 double deltaPhi = tubs.deltaPhi();
0117 double cutAtStart = convertUnits(tubs.cutAtStart());
0118 double cutAtDelta = convertUnits(tubs.cutAtDelta());
0119 bool cutInside = tubs.cutInside();
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
0126
0127
0128
0129
0130
0131
0132
0133
0134
0135
0136
0137
0138
0139
0140
0141
0142
0143
0144
0145
0146
0147
0148
0149
0150
0151
0152
0153
0154
0155
0156
0157
0158
0159
0160
0161
0162
0163
0164
0165
0166
0167
0168
0169
0170
0171
0172
0173
0174
0175
0176
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
0203
0204
0205
0206
0207
0208
0209
0210
0211
0212
0213
0214
0215
0216
0217
0218
0219
0220
0221
0222
0223
0224
0225
0226
0227
0228
0229
0230
0231 Surface::PositionType &posResult = center_;
0232
0233
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
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
0256 if (debug) {
0257 cout << "Refplane pos " << refPlane->position() << endl;
0258
0259
0260 LocalVector globalZdir(0., 0., 1.);
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
0278
0279 if (expand && (i == phiplus || i == phiminus))
0280 continue;
0281
0282
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"