File indexing completed on 2024-04-06 12:22:30
0001
0002
0003
0004
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
0040 center_ = GlobalPoint(transArray[0] / dd4hep::cm, transArray[1] / dd4hep::cm, transArray[2] / dd4hep::cm);
0041
0042
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
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
0100
0101
0102
0103
0104
0105
0106
0107
0108
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
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;
0126 double rIn = tubs.rMin() / dd4hep::cm;
0127 double rOut = tubs.rMax() / dd4hep::cm;
0128 double startPhi = tubs.startPhi();
0129 double deltaPhi = tubs.deltaPhi();
0130 double cutAtStart = tubs.cutAtStart() / dd4hep::cm;
0131 double cutAtDelta = tubs.cutAtDelta() / dd4hep::cm;
0132 bool cutInside = tubs.cutInside();
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
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
0166
0167
0168
0169
0170
0171
0172
0173
0174
0175
0176
0177
0178
0179
0180
0181
0182
0183
0184
0185
0186
0187
0188
0189
0190
0191
0192
0193
0194 Surface::PositionType &posResult = center_;
0195
0196
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
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
0221 if (debug) {
0222 LogTrace("MagGeoBuilder") << "Refplane pos " << refPlane->position();
0223
0224
0225 LocalVector globalZdir(0., 0., 1.);
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
0244
0245 if (expand && (i == phiplus || i == phiminus))
0246 continue;
0247
0248
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"