File indexing completed on 2021-02-14 13:08:05
0001 #ifndef GEOMETRY_ECALGEOMETRYLOADER_ICC
0002 #define GEOMETRY_ECALGEOMETRYLOADER_ICC 1
0003
0004 #include "Geometry/EcalTestBeam/test/ee/CaloGeometryLoaderTest.h"
0005 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0006 #include "DetectorDescription/Core/interface/DDCompactView.h"
0007 #include "DetectorDescription/Core/interface/DDSolid.h"
0008 #include "DetectorDescription/Core/interface/DDSpecifics.h"
0009
0010 #include "CLHEP/Units/GlobalSystemOfUnits.h"
0011 #include "CLHEP/Geometry/Plane3D.h"
0012
0013 #include <vector>
0014 #include <fstream>
0015
0016 typedef CaloCellGeometry::CCGFloat CCGFloat;
0017
0018 template <class T>
0019 const double CaloGeometryLoaderTest<T>::k_ScaleFromDDDtoGeant(0.1);
0020
0021 template <class T>
0022 CaloGeometryLoaderTest<T>::CaloGeometryLoaderTest() {
0023 m_filter.setCriteria(DDValue("Volume",
0024 "EndcapSC",
0025
0026 0),
0027 DDCompOp::equals);
0028 }
0029
0030 template <class T>
0031 typename CaloGeometryLoaderTest<T>::PtrType CaloGeometryLoaderTest<T>::load(const DDCompactView* cpv,
0032 const Alignments* alignments,
0033 const Alignments* globals) {
0034 PtrType geom = std::make_unique<T>();
0035
0036 makeGeometry(cpv, dynamic_cast<T*>(geom.get()), alignments, globals);
0037
0038 return geom;
0039 }
0040
0041 template <class T>
0042 void CaloGeometryLoaderTest<T>::makeGeometry(const DDCompactView* cpv,
0043 T* geom,
0044 const Alignments* alignments,
0045 const Alignments* globals) {
0046 geom->allocateCorners(T::k_NumberOfCellsForCorners);
0047 geom->allocatePar(T::k_NumberOfParametersPerShape * T::k_NumberOfShapes, T::k_NumberOfParametersPerShape);
0048
0049 DDFilteredView fv(*cpv, m_filter);
0050
0051 unsigned int counter(0);
0052 for (bool doSubDets = fv.firstChild(); doSubDets; doSubDets = fv.nextSibling()) {
0053 ++counter;
0054 const DDSolid& solid(fv.logicalPart().solid());
0055
0056 const ParmVec& parameters(solid.parameters());
0057
0058 DD3Vector x, y, z;
0059 fv.rotation().GetComponents(x, y, z);
0060 const CLHEP::HepRep3x3 temp(x.X(), y.X(), z.X(), x.Y(), y.Y(), z.Y(), x.Z(), y.Z(), z.Z());
0061 const CLHEP::HepRotation hr(temp);
0062 const CLHEP::Hep3Vector h3v(fv.translation().X(), fv.translation().Y(), fv.translation().Z());
0063 const HepGeom::Transform3D ht3d(hr,
0064 k_ScaleFromDDDtoGeant * h3v);
0065
0066 const unsigned int id(getDetIdForDDDNode(fv));
0067
0068 myFillGeom(geom, parameters, ht3d, id);
0069 }
0070 }
0071
0072 template <class T>
0073 void CaloGeometryLoaderTest<T>::myFillGeom(T* geom,
0074 const ParmVec& vv,
0075 const HepGeom::Transform3D& tr,
0076 const unsigned int id) {
0077 static const double kRadToDeg(180 / M_PI);
0078 static const HepGeom::Point3D<double> PIVOT(13533.35, -000.02, 0);
0079 static const HepGeom::Point3D<double> FLOOR(13533.35, -000.02, -152.23);
0080
0081 static const HepGeom::Point3D<double> IP(0, 0, 0);
0082 HepGeom::Point3D<double> AA0;
0083 HepGeom::Point3D<double> BB0;
0084 HepGeom::Point3D<double> CC0;
0085 HepGeom::Point3D<double> DD0;
0086 HepGeom::Point3D<double> EE0;
0087 HepGeom::Point3D<double> FF0;
0088 HepGeom::Point3D<double> GG0;
0089 HepGeom::Point3D<double> HH0;
0090
0091 std::vector<CCGFloat> pv;
0092 pv.reserve(vv.size());
0093 for (unsigned int i(0); i != vv.size(); ++i) {
0094 const double factor(1 == i || 2 == i || 6 == i || 10 == i ? 1 : k_ScaleFromDDDtoGeant);
0095 pv.push_back(factor * vv[i]);
0096 }
0097
0098 std::vector<GlobalPoint> corners(8);
0099
0100 TruncatedPyramid::createCorners(pv, tr, corners);
0101
0102 if ((160 == id || 164 == id || 134 == id || 133 == id || 153 == id || 150 == id || 130 == id) &&
0103 (0 < corners[0].x() && 0 < corners[0].y() && 0 < corners[0].z())) {
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115 if (160 == id)
0116 AA0 = HepGeom::Point3D<double>(corners[3].x(), corners[3].y(), corners[3].z());
0117 if (164 == id)
0118 BB0 = HepGeom::Point3D<double>(corners[2].x(), corners[2].y(), corners[2].z());
0119 if (134 == id)
0120 CC0 = HepGeom::Point3D<double>(corners[1].x(), corners[1].y(), corners[1].z());
0121 if (130 == id)
0122 DD0 = HepGeom::Point3D<double>(corners[0].x(), corners[0].y(), corners[0].z());
0123 if (150 == id)
0124 EE0 = HepGeom::Point3D<double>(corners[2].x(), corners[2].y(), corners[2].z());
0125 if (153 == id)
0126 FF0 = HepGeom::Point3D<double>(corners[2].x(), corners[2].y(), corners[2].z());
0127 if (133 == id)
0128 GG0 = HepGeom::Point3D<double>(corners[2].x(), corners[2].y(), corners[2].z());
0129 if (130 == id)
0130 HH0 = HepGeom::Point3D<double>(corners[2].x(), corners[2].y(), corners[2].z());
0131 }
0132 static bool done(false);
0133 if (0.001 < AA0.mag() && 0.001 < BB0.mag() && 0.001 < CC0.mag() && 0.001 < DD0.mag() && !done) {
0134 done = true;
0135
0136 const HepGeom::Transform3D H4ToBeamLine(PIVOT + HepGeom::Point3D<double>(0, 1, 0),
0137 PIVOT + HepGeom::Point3D<double>(0, 0, 1),
0138 PIVOT + HepGeom::Point3D<double>(1, 0, 0),
0139 HepGeom::Point3D<double>(-1, 0, 0),
0140 HepGeom::Point3D<double>(0, -1, 0),
0141 HepGeom::Point3D<double>(0, 0, 1));
0142
0143
0144
0145
0146 std::cout << "\nH4ToBeamLine"
0147 << " Theta_x = " << H4ToBeamLine.getRotation().thetaX() * kRadToDeg
0148 << ", Theta_y = " << H4ToBeamLine.getRotation().thetaY() * kRadToDeg
0149 << ", Theta_z = " << H4ToBeamLine.getRotation().thetaZ() * kRadToDeg
0150 << ", eta = " << -log(tan(H4ToBeamLine.getRotation().thetaZ() / 2.)) << std::endl;
0151
0152 std::cout << "\nH4ToBeamLine"
0153 << " Phi_x = " << H4ToBeamLine.getRotation().phiX() * kRadToDeg
0154 << ", Phi_y = " << H4ToBeamLine.getRotation().phiY() * kRadToDeg
0155 << ", Phi_z = " << H4ToBeamLine.getRotation().phiZ() * kRadToDeg << std::endl;
0156
0157 std::cout << "pivot to beamline=" << H4ToBeamLine * PIVOT << std::endl;
0158 std::cout << "endcap to beamline=" << H4ToBeamLine * HepGeom::Point3D<double>(PIVOT.x() + 316, 0, 50) << std::endl;
0159
0160 std::cout << "\nA0 = " << AA0 << "\nB0 = " << BB0 << "\nC0 = " << CC0 << "\nD0 = " << DD0 << std::endl;
0161
0162 const HepGeom::Point3D<double>& A0(AA0);
0163 const HepGeom::Point3D<double>& B0(BB0);
0164 const HepGeom::Point3D<double>& C0(CC0);
0165 const HepGeom::Point3D<double>& D0(DD0);
0166
0167 const double angle((C0 - B0).angle(A0 - B0));
0168
0169 static const std::vector<HepGeom::Point3D<double> > extra = {HepGeom::Point3D<double>(13859.05, 005.79, -002.66),
0170 HepGeom::Point3D<double>(13860.98, -038.39, -002.11),
0171 HepGeom::Point3D<double>(13853.19, -039.25, -030.77),
0172 HepGeom::Point3D<double>(13851.23, 004.99, -031.39)};
0173
0174 static const std::vector<HepGeom::Point3D<double> > sur = {
0175 HepGeom::Point3D<double>(13862.30, 20.68, 11.60), HepGeom::Point3D<double>(13865.69, -52.85, 12.64),
0176 HepGeom::Point3D<double>(13851.53, -54.27, -44.75), HepGeom::Point3D<double>(13848.21, 19.18, -45.89),
0177
0178 HepGeom::Point3D<double>(13861.72, 28.62, 11.61), HepGeom::Point3D<double>(13866.88, -44.81, 12.64),
0179 HepGeom::Point3D<double>(13852.76, -46.56, -44.73), HepGeom::Point3D<double>(13847.67, 26.79, -45.88),
0180 HepGeom::Point3D<double>(13856.82, 63.41, 11.60), HepGeom::Point3D<double>(13869.76, -9.06, 12.64),
0181 HepGeom::Point3D<double>(13855.90, -12.31, -44.73), HepGeom::Point3D<double>(13843.04, 60.08, -45.89),
0182 HepGeom::Point3D<double>(13860.86, -36.96, 11.60), HepGeom::Point3D<double>(13851.38, -109.96, 12.63),
0183 HepGeom::Point3D<double>(13837.21, -108.89, -44.73), HepGeom::Point3D<double>(13846.73, -35.98, -45.88),
0184 HepGeom::Point3D<double>(13854.03, -76.21, 11.60), HepGeom::Point3D<double>(13835.82, -147.53, 12.64),
0185 HepGeom::Point3D<double>(13821.88, -144.76, -44.73), HepGeom::Point3D<double>(13840.13, -73.54, -45.88),
0186 HepGeom::Point3D<double>(13861.10, 35.48, 09.88), HepGeom::Point3D<double>(13867.78, -37.82, 11.07),
0187 HepGeom::Point3D<double>(13853.42, -40.01, -46.22), HepGeom::Point3D<double>(13846.80, 33.21, -47.52),
0188 HepGeom::Point3D<double>(13860.68, 35.23, 19.89), HepGeom::Point3D<double>(13867.35, -38.08, 20.23),
0189 HepGeom::Point3D<double>(13854.67, -39.43, -37.48), HepGeom::Point3D<double>(13848.07, 33.79, -37.92),
0190 HepGeom::Point3D<double>(13859.95, 34.97, 29.87), HepGeom::Point3D<double>(13866.62, -38.34, 29.33),
0191 HepGeom::Point3D<double>(13855.61, -38.88, -28.76), HepGeom::Point3D<double>(13849.00, 34.34, -28.32),
0192 HepGeom::Point3D<double>(13858.93, 34.41, 39.81), HepGeom::Point3D<double>(13865.65, -38.90, 38.40),
0193 HepGeom::Point3D<double>(13856.34, -38.58, -19.97), HepGeom::Point3D<double>(13849.70, 34.65, -18.68),
0194 HepGeom::Point3D<double>(13857.67, 33.69, 49.76), HepGeom::Point3D<double>(13864.43, -39.60, 47.48),
0195 HepGeom::Point3D<double>(13856.81, -38.43, -11.13), HepGeom::Point3D<double>(13850.14, 34.78, -8.95)};
0196
0197 double fsumx(0);
0198 double fsumy(0);
0199 double fsumz(0);
0200 double fsumx2(0);
0201 double fsumy2(0);
0202 double fsumz2(0);
0203
0204 double psumx(0);
0205 double psumy(0);
0206 double psumz(0);
0207 double psumx2(0);
0208 double psumy2(0);
0209 double psumz2(0);
0210
0211 double sumx(0);
0212 double sumy(0);
0213 double sumz(0);
0214 double sumx2(0);
0215 double sumy2(0);
0216 double sumz2(0);
0217 unsigned int nsur(sur.size() / 4);
0218
0219 const std::string names[10] = {"H0", "H1", "H2", "H3", "H4", "V1", "V2", "V3", "V4", "V5"};
0220
0221 HepGeom::Transform3D V3Transform;
0222
0223 for (unsigned int i(0); i != nsur; ++i) {
0224 const HepGeom::Point3D<double> SA(sur[4 * i]);
0225 const HepGeom::Point3D<double> SB(sur[4 * i + 1]);
0226 const HepGeom::Point3D<double> SC(sur[4 * i + 2]);
0227 const HepGeom::Point3D<double> SD(sur[4 * i + 3]);
0228
0229
0230
0231
0232
0233 const HepGeom::Point3D<double>& aa(SA);
0234 const HepGeom::Point3D<double> bb(SA + (SB - SA).unit() * (B0 - A0).mag());
0235
0236 const HepGeom::Plane3D<double> plane(aa, bb, SC);
0237
0238 const HepGeom::Point3D<double> cc(
0239 bb +
0240 (HepGeom::Rotate3D(-angle, plane.normal()) * HepGeom::Vector3D<double>(aa - bb)).unit() * (C0 - B0).mag());
0241
0242 const HepGeom::Transform3D trform(A0, B0, C0, aa, bb, cc);
0243
0244 const HepGeom::Transform3D trforminv(aa, bb, cc, A0, B0, C0);
0245
0246 const HepGeom::Point3D<double> dd(trform * D0);
0247
0248 const HepGeom::Transform3D BeamLineToCMS(H4ToBeamLine * aa, H4ToBeamLine * bb, H4ToBeamLine * cc, A0, B0, C0);
0249
0250 if (7 == i)
0251 V3Transform = BeamLineToCMS;
0252
0253 const CLHEP::HepRotation blcmsrot(BeamLineToCMS.getRotation());
0254 const HepGeom::Vector3D<double> blcmstra(BeamLineToCMS.getTranslation());
0255
0256 std::cout << "\n********* Name is " << names[i];
0257
0258 std::cout << "\n\n B discrepancy is " << (bb - SB) << ", " << (bb - SB).mag() << "\n\n C discrepancy is "
0259 << (cc - SC) << ", " << (cc - SC).mag() << "\n\n D discrepancy is " << (dd - SD) << ", "
0260 << (dd - SD).mag() << std::endl;
0261
0262
0263
0264
0265
0266
0267
0268
0269
0270 const double dxdd(dd.x() - SD.x());
0271 const double dydd(dd.y() - SD.y());
0272 const double dzdd(dd.z() - SD.z());
0273
0274 sumx += dxdd;
0275 sumy += dydd;
0276 sumz += dzdd;
0277 sumx2 += dxdd * dxdd;
0278 sumy2 += dydd * dydd;
0279 sumz2 += dzdd * dzdd;
0280
0281 if (0 == i) {
0282 const HepGeom::Point3D<double> ee(trform * EE0);
0283 const HepGeom::Point3D<double> ff(trform * FF0);
0284 const HepGeom::Point3D<double> gg(trform * GG0);
0285 const HepGeom::Point3D<double> hh(trform * HH0);
0286
0287 std::cout << "\n\n E discrepancy is " << (ee - extra[0]) << ", " << (ee - extra[0]).mag()
0288 << "\n\n F discrepancy is " << (ff - extra[1]) << ", " << (ff - extra[1]).mag()
0289 << "\n\n G discrepancy is " << (gg - extra[2]) << ", " << (gg - extra[2]).mag()
0290 << "\n\n H discrepancy is " << (hh - extra[3]) << ", " << (hh - extra[3]).mag() << std::endl;
0291
0292 std::cout << "\n Point SA in CMS coord = " << trforminv * SA << std::endl;
0293 }
0294 const HepGeom::Point3D<double> cmsPiv(trforminv * PIVOT);
0295 psumx += cmsPiv.x();
0296 psumy += cmsPiv.y();
0297 psumz += cmsPiv.z();
0298 psumx2 += cmsPiv.x() * cmsPiv.x();
0299 psumy2 += cmsPiv.y() * cmsPiv.y();
0300 psumz2 += cmsPiv.z() * cmsPiv.z();
0301 const HepGeom::Point3D<double> cmsFlr(trforminv * FLOOR);
0302 fsumx += cmsFlr.x();
0303 fsumy += cmsFlr.y();
0304 fsumz += cmsFlr.z();
0305 fsumx2 += cmsFlr.x() * cmsFlr.x();
0306 fsumy2 += cmsFlr.y() * cmsFlr.y();
0307 fsumz2 += cmsFlr.z() * cmsFlr.z();
0308
0309 const CLHEP::HepRotation rotInv(trforminv.getRotation());
0310
0311 std::cout << "\n CMS IP in H4 coord = " << trform * IP << std::endl;
0312 std::cout << "\n Pivot point in CMS coord = " << cmsPiv << std::endl;
0313 std::cout << "\n Floor point in CMS coord = " << cmsFlr << std::endl;
0314
0315
0316
0317
0318
0319
0320
0321
0322
0323
0324
0325
0326
0327
0328 std::cout << "\n"
0329 << names[i] << " Theta_x = " << blcmsrot.thetaX() * kRadToDeg
0330 << ", Theta_y = " << blcmsrot.thetaY() * kRadToDeg << ", Theta_z = " << blcmsrot.thetaZ() * kRadToDeg
0331 << ", eta = " << -log(tan(blcmsrot.thetaZ() / 2.)) << std::endl;
0332
0333 std::cout << "\n"
0334 << names[i] << " Phi_x = " << blcmsrot.phiX() * kRadToDeg
0335 << ", Phi_y = " << blcmsrot.phiY() * kRadToDeg << ", Phi_z = " << blcmsrot.phiZ() * kRadToDeg
0336 << std::endl;
0337
0338 std::cout << "\n\n BeamLine to CMS translation" << blcmstra << std::endl;
0339 }
0340 const double ntot(1.0 * nsur);
0341 sumx /= ntot;
0342 sumy /= ntot;
0343 sumz /= ntot;
0344 sumx2 = sqrt(sumx2 / ntot - sumx * sumx);
0345 sumy2 = sqrt(sumy2 / ntot - sumy * sumy);
0346 sumz2 = sqrt(sumz2 / ntot - sumz * sumz);
0347
0348 std::cout << "\n Number of survey points = " << ntot;
0349 std::cout << "\n====== Mean x deviation = " << sumx << " +- " << sumx2;
0350 std::cout << "\n====== Mean y deviation = " << sumy << " +- " << sumy2;
0351 std::cout << "\n====== Mean z deviation = " << sumz << " +- " << sumz2 << std::endl;
0352
0353 psumx /= ntot;
0354 psumy /= ntot;
0355 psumz /= ntot;
0356 psumx2 = sqrt(psumx2 / ntot - psumx * psumx);
0357 psumy2 = sqrt(psumy2 / ntot - psumy * psumy);
0358 psumz2 = sqrt(psumz2 / ntot - psumz * psumz);
0359
0360 std::cout << "\n====== Mean x of pivot = " << psumx << " +- " << psumx2;
0361 std::cout << "\n====== Mean y of pivot = " << psumy << " +- " << psumy2;
0362 std::cout << "\n====== Mean z of pivot = " << psumz << " +- " << psumz2 << std::endl;
0363
0364 fsumx /= ntot;
0365 fsumy /= ntot;
0366 fsumz /= ntot;
0367 fsumx2 = sqrt(fsumx2 / ntot - fsumx * fsumx);
0368 fsumy2 = sqrt(fsumy2 / ntot - fsumy * fsumy);
0369 fsumz2 = sqrt(fsumz2 / ntot - fsumz * fsumz);
0370
0371 std::cout << "\n====== Mean x of floor = " << fsumx << " +- " << fsumx2;
0372 std::cout << "\n====== Mean y of floor = " << fsumy << " +- " << fsumy2;
0373 std::cout << "\n====== Mean z of floor = " << fsumz << " +- " << fsumz2 << std::endl;
0374
0375 const HepGeom::Vector3D<double> hpl[] = {HepGeom::Plane3D<double>(sur[16], sur[0], sur[8]).normal(),
0376 HepGeom::Plane3D<double>(sur[17], sur[1], sur[9]).normal(),
0377 HepGeom::Plane3D<double>(sur[18], sur[2], sur[10]).normal(),
0378 HepGeom::Plane3D<double>(sur[19], sur[3], sur[11]).normal(),
0379 HepGeom::Plane3D<double>(sur[16], sur[12], sur[8]).normal(),
0380 HepGeom::Plane3D<double>(sur[17], sur[13], sur[9]).normal(),
0381 HepGeom::Plane3D<double>(sur[18], sur[14], sur[10]).normal(),
0382 HepGeom::Plane3D<double>(sur[19], sur[15], sur[11]).normal()};
0383
0384 const HepGeom::Vector3D<double> hplavr(0.125 *
0385 (hpl[0] + hpl[1] + hpl[2] + hpl[3] + hpl[4] + hpl[5] + hpl[6] + hpl[7]));
0386
0387 std::cout << "\n\n+++ HNormal for A =" << hpl[0].theta() * kRadToDeg << ", " << hpl[0].phi() * kRadToDeg
0388 << "\n+++ HNormal for B =" << hpl[1].theta() * kRadToDeg << ", " << hpl[1].phi() * kRadToDeg
0389 << "\n+++ HNormal for C =" << hpl[2].theta() * kRadToDeg << ", " << hpl[2].phi() * kRadToDeg
0390 << "\n+++ HNormal for D =" << hpl[3].theta() * kRadToDeg << ", " << hpl[3].phi() * kRadToDeg
0391 << "\n+++ HNormal for A =" << hpl[4].theta() * kRadToDeg << ", " << hpl[4].phi() * kRadToDeg
0392 << "\n+++ HNormal for B =" << hpl[5].theta() * kRadToDeg << ", " << hpl[5].phi() * kRadToDeg
0393 << "\n+++ HNormal for C =" << hpl[6].theta() * kRadToDeg << ", " << hpl[6].phi() * kRadToDeg
0394 << "\n+++ HNormal for D =" << hpl[7].theta() * kRadToDeg << ", " << hpl[7].phi() * kRadToDeg
0395 << "\n+++ HAverage =" << hplavr.theta() * kRadToDeg << ", " << hplavr.phi() * kRadToDeg << std::endl;
0396
0397 const HepGeom::Vector3D<double> vpl[] = {HepGeom::Plane3D<double>(sur[20], sur[28], sur[36]).normal(),
0398 HepGeom::Plane3D<double>(sur[21], sur[29], sur[37]).normal(),
0399 HepGeom::Plane3D<double>(sur[22], sur[30], sur[38]).normal(),
0400 HepGeom::Plane3D<double>(sur[23], sur[31], sur[39]).normal(),
0401
0402 HepGeom::Plane3D<double>(sur[20], sur[24], sur[36]).normal(),
0403 HepGeom::Plane3D<double>(sur[21], sur[25], sur[37]).normal(),
0404 HepGeom::Plane3D<double>(sur[22], sur[26], sur[38]).normal(),
0405 HepGeom::Plane3D<double>(sur[23], sur[27], sur[39]).normal(),
0406
0407 HepGeom::Plane3D<double>(sur[20], sur[32], sur[36]).normal(),
0408 HepGeom::Plane3D<double>(sur[21], sur[33], sur[37]).normal(),
0409 HepGeom::Plane3D<double>(sur[22], sur[34], sur[38]).normal(),
0410 HepGeom::Plane3D<double>(sur[23], sur[35], sur[39]).normal(),
0411
0412 HepGeom::Plane3D<double>(sur[20], sur[32], sur[36]).normal(),
0413 HepGeom::Plane3D<double>(sur[21], sur[33], sur[37]).normal(),
0414 HepGeom::Plane3D<double>(sur[22], sur[34], sur[38]).normal(),
0415 HepGeom::Plane3D<double>(sur[23], sur[35], sur[39]).normal()};
0416
0417 const HepGeom::Vector3D<double> vplavr(0.0625 * (vpl[0] + vpl[1] + vpl[2] + vpl[3] + vpl[4] + vpl[5] + vpl[6] +
0418 vpl[7] + vpl[8] + vpl[9] + vpl[10] + vpl[11] + vpl[12] + vpl[13] +
0419 vpl[14] + vpl[15]));
0420
0421 std::cout << "\n\n+++ VNormal for A =" << vpl[0].theta() * kRadToDeg << ", " << vpl[0].phi() * kRadToDeg
0422 << "\n+++ VNormal for B =" << vpl[1].theta() * kRadToDeg << ", " << vpl[1].phi() * kRadToDeg
0423 << "\n+++ VNormal for C =" << vpl[2].theta() * kRadToDeg << ", " << vpl[2].phi() * kRadToDeg
0424 << "\n+++ VNormal for D =" << vpl[3].theta() * kRadToDeg << ", " << vpl[3].phi() * kRadToDeg
0425 << "\n+++ VNormal for A =" << vpl[4].theta() * kRadToDeg << ", " << vpl[4].phi() * kRadToDeg
0426 << "\n+++ VNormal for B =" << vpl[5].theta() * kRadToDeg << ", " << vpl[5].phi() * kRadToDeg
0427 << "\n+++ VNormal for C =" << vpl[6].theta() * kRadToDeg << ", " << vpl[6].phi() * kRadToDeg
0428 << "\n+++ VNormal for D =" << vpl[7].theta() * kRadToDeg << ", " << vpl[7].phi() * kRadToDeg
0429 << "\n+++ VNormal for A =" << vpl[8].theta() * kRadToDeg << ", " << vpl[8].phi() * kRadToDeg
0430 << "\n+++ VNormal for B =" << vpl[9].theta() * kRadToDeg << ", " << vpl[9].phi() * kRadToDeg
0431 << "\n+++ VNormal for C =" << vpl[10].theta() * kRadToDeg << ", " << vpl[10].phi() * kRadToDeg
0432 << "\n+++ VNormal for D =" << vpl[11].theta() * kRadToDeg << ", " << vpl[11].phi() * kRadToDeg
0433 << "\n+++ VNormal for A =" << vpl[12].theta() * kRadToDeg << ", " << vpl[12].phi() * kRadToDeg
0434 << "\n+++ VNormal for B =" << vpl[13].theta() * kRadToDeg << ", " << vpl[13].phi() * kRadToDeg
0435 << "\n+++ VNormal for C =" << vpl[14].theta() * kRadToDeg << ", " << vpl[14].phi() * kRadToDeg
0436 << "\n+++ VNormal for D =" << vpl[15].theta() * kRadToDeg << ", " << vpl[15].phi() * kRadToDeg
0437 << "\n+++ VAverage =" << vplavr.theta() * kRadToDeg << ", " << vplavr.phi() * kRadToDeg << std::endl;
0438
0439 unsigned int Hinput(0);
0440 unsigned int Vinput(0);
0441
0442 std::cout << "\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n"
0443 << "\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n"
0444 << "\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n"
0445 << "\n+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
0446 << "\n+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
0447 << "\n+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
0448 << "\n+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
0449 << "\n+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
0450 << "\n+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
0451 << "\n+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
0452
0453 bool keepGoing(true);
0454 while (keepGoing) {
0455 std::cout << "\n\n+++ Preparing to make two files for the table position you specify";
0456
0457 std::cout << "\n\n**Enter desired H position (integer, 25000-300000) [0 to quit]: " << std::ends;
0458 std::cin >> Hinput;
0459
0460 keepGoing = (0 != Hinput);
0461
0462 if (keepGoing) {
0463 std::cout << "\n\n**Enter desired V position (integer, 5000-30000): " << std::ends;
0464 std::cin >> Vinput;
0465
0466 std::cout << "\n\nH input value =" << Hinput << ", V input value =" << Vinput << std::endl;
0467
0468 if (25000 > Hinput || 300000 < Hinput || 5000 > Vinput || 30000 < Vinput) {
0469 std::cout << "!!!!!!!! Invalid input(s), try again !!!!!!!" << std::endl;
0470 } else {
0471
0472
0473 const HepGeom::Plane3D<double> plaV(sur[20], sur[28], sur[36]);
0474
0475 const HepGeom::Point3D<double> projPiv(plaV.point(PIVOT));
0476
0477 const double angleVPerUnit((plaV.point(sur[36]) - projPiv).angle(plaV.point(sur[20]) - projPiv) / 16000.);
0478
0479 const HepGeom::Plane3D<double> plaH(sur[16], sur[0], sur[8]);
0480
0481 const HepGeom::Point3D<double> projPivH(plaH.point(PIVOT));
0482
0483 const double angleHPerUnit((plaH.point(sur[8]) - projPivH).angle(plaH.point(sur[16]) - projPivH) / 176458);
0484
0485
0486
0487
0488
0489
0490
0491
0492
0493
0494
0495
0496
0497
0498
0499
0500
0501
0502
0503
0504
0505
0506
0507
0508
0509
0510
0511
0512
0513
0514 const double hAngle((1. * Hinput - 191000.) * angleHPerUnit);
0515
0516 const double vAngle((1. * Vinput - 16000.) * angleVPerUnit);
0517
0518 const HepGeom::Rotate3D vRotate(-vAngle, H4ToBeamLine * vpl[0]);
0519 const HepGeom::RotateY3D hRotate(hAngle);
0520
0521 const HepGeom::Transform3D myTrf(V3Transform * vRotate * hRotate);
0522
0523 const CLHEP::HepRotation myRot(myTrf.getRotation());
0524 const HepGeom::Vector3D<double> myTran(myTrf.getTranslation());
0525
0526 const HepGeom::RotateZ3D R1(myRot.phiZ() - M_PI);
0527 const HepGeom::Point3D<double> xUnit(0, 1, 0);
0528 const HepGeom::Point3D<double> yUnit(-1, 0, 0);
0529 const HepGeom::Point3D<double> zUnit(0, 0, 1);
0530
0531 const HepGeom::Transform3D RXRZ(HepGeom::Rotate3D(-myRot.thetaZ(), R1 * xUnit) * R1);
0532
0533 const HepGeom::Point3D<double> zzUnit(RXRZ * zUnit);
0534 const HepGeom::Point3D<double> yyUnit(RXRZ * yUnit);
0535 const HepGeom::Point3D<double> xxUnit(RXRZ * xUnit);
0536
0537
0538
0539
0540
0541
0542
0543
0544
0545
0546
0547
0548
0549
0550
0551 const HepGeom::Transform3D RPSI(HepGeom::Transform3D(myRot, HepGeom::Point3D<double>(0, 0, 0)) *
0552 (RXRZ.inverse()));
0553
0554 CLHEP::Hep3Vector axis;
0555 double psi;
0556 RPSI.getRotation().getAngleAxis(psi, axis);
0557
0558 const double zdiff(axis.angle(zzUnit - HepGeom::Point3D<double>(0, 0, 0)));
0559
0560 if (0.95 * M_PI < fabs(zdiff))
0561 psi = -psi;
0562
0563
0564
0565
0566
0567
0568
0569
0570
0571
0572
0573
0574
0575
0576
0577
0578
0579
0580
0581
0582
0583
0584
0585 const double TthetaX(myRot.thetaX());
0586 const double TthetaY(myRot.thetaY());
0587 const double TthetaZ(myRot.thetaZ());
0588 const double TphiX(myRot.phiX());
0589 const double TphiY(myRot.phiY());
0590 const double TphiZ(myRot.phiZ());
0591 const double Txtra(myTran.x());
0592 const double Tytra(myTran.y());
0593 const double Tztra(myTran.z());
0594 const double TetaZ(-log(tan(TthetaZ / 2.)));
0595
0596 const double Tpsi(psi);
0597
0598
0599
0600 std::stringstream ssTpsir;
0601 ssTpsir << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << Tpsi;
0602 const std::string strTpsir(ssTpsir.str());
0603
0604 std::stringstream ssTxtra;
0605 ssTxtra << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << Txtra << "*cm";
0606 const std::string strTxtra(ssTxtra.str());
0607 std::stringstream ssTytra;
0608 ssTytra << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << Tytra << "*cm";
0609 const std::string strTytra(ssTytra.str());
0610 std::stringstream ssTztra;
0611 ssTztra << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << Tztra << "*cm";
0612 const std::string strTztra(ssTztra.str());
0613
0614 std::stringstream ssTetaZ;
0615 ssTetaZ << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << TetaZ;
0616 const std::string strTetaZ(ssTetaZ.str());
0617
0618 std::stringstream ssTthetaX;
0619 ssTthetaX << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << TthetaX * kRadToDeg
0620 << "*deg";
0621 const std::string strTthetaX(ssTthetaX.str());
0622 std::stringstream ssTphiX;
0623 ssTphiX << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << TphiX * kRadToDeg
0624 << "*deg";
0625 const std::string strTphiX(ssTphiX.str());
0626
0627
0628
0629
0630
0631
0632
0633
0634
0635 std::stringstream ssTthetaY;
0636 ssTthetaY << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << TthetaY * kRadToDeg
0637 << "*deg";
0638 const std::string strTthetaY(ssTthetaY.str());
0639 std::stringstream ssTphiY;
0640 ssTphiY << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << TphiY * kRadToDeg
0641 << "*deg";
0642 const std::string strTphiY(ssTphiY.str());
0643
0644
0645
0646
0647
0648
0649
0650
0651
0652 std::stringstream ssTthetaZ;
0653 ssTthetaZ << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << TthetaZ * kRadToDeg
0654 << "*deg";
0655 const std::string strTthetaZ(ssTthetaZ.str());
0656 std::stringstream ssTphiZ;
0657 ssTphiZ << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << TphiZ * kRadToDeg
0658 << "*deg";
0659 const std::string strTphiZ(ssTphiZ.str());
0660 std::stringstream ssTphiZr;
0661 ssTphiZr << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << TphiZ;
0662 const std::string strTphiZr(ssTphiZr.str());
0663
0664 std::stringstream ssTxpos;
0665 ssTxpos << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << psumy;
0666 const std::string strTxpos(ssTxpos.str());
0667
0668 std::stringstream ssTypos;
0669 ssTypos << std::setw(7) << setiosflags(std::ios::fixed) << std::setprecision(5) << -psumx;
0670 const std::string strTypos(ssTypos.str());
0671
0672 std::stringstream ssTzpos;
0673 ssTzpos << std::setw(11) << setiosflags(std::ios::fixed) << std::setprecision(2) << -26733.5 + psumz;
0674 const std::string strTzpos(ssTzpos.str());
0675
0676 std::stringstream ss;
0677 ss << "tbrot_H" << std::setw(6) << std::setfill('0') << Hinput << "_V" << std::setw(5) << std::setfill('0')
0678 << Vinput << ".xml" << std::ends;
0679 const std::string myname(ss.str());
0680
0681 std::ofstream tbrotfile(myname.c_str());
0682
0683 tbrotfile << "<?xml version=\"1.0\"?>"
0684 << "<DDDefinition xmlns=\"http://www.cern.ch/cms/DDL\""
0685 << " xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\""
0686 << " xsi:schemaLocation=\"http://www.cern.ch/cms/DDL ../../"
0687 << "DetectorDescription/Schema/DDLSchema.xsd\">"
0688 << "\n\n<ConstantsSection label=\"tbrot.xml\" eval=\"true\" >"
0689 << "\n\n<Constant name=\"eta\" value=\"" << strTetaZ
0690 << "\"/>"
0691
0692 << "\n\n<Constant name=\"theta\" value=\"" << strTthetaZ << "\"/>"
0693 << "\n\n<Constant name=\"phi\" value=\"" << strTphiZ << "\"/>"
0694 << "\n\n<Constant name=\"theta_X\" value=\"" << strTthetaX << "\"/>"
0695 << "\n\n<Constant name=\"phi_X\" value=\"" << strTphiX << "\"/>"
0696 << "\n\n<Constant name=\"theta_Y\" value=\"" << strTthetaY << "\"/>"
0697 << "\n\n<Constant name=\"phi_Y\" value=\"" << strTphiY << "\"/>"
0698 << "\n\n<Constant name=\"xtran\" value=\"" << strTxtra << "\"/>"
0699 << "\n<Constant name=\"ytran\" value=\"" << strTytra << "\"/>"
0700 << "\n<Constant name=\"ztran\" value=\"" << strTztra << "\"/>"
0701 << "\n\n</ConstantsSection>"
0702 << "\n\n</DDDefinition>" << std::endl;
0703
0704 tbrotfile.close();
0705
0706
0707
0708 std::stringstream ss2;
0709 ss2 << "ee_beam_H" << std::setw(6) << std::setfill('0') << Hinput << "_V" << std::setw(5) << std::setfill('0')
0710 << Vinput << "_cff.py" << std::ends;
0711 const std::string myname2(ss2.str());
0712
0713 std::ofstream beamfile(myname2.c_str());
0714
0715 beamfile << "import FWCore.ParameterSet.Config as cms"
0716 << "\n\ncommon_beam_direction_parameters = cms.PSet("
0717 << "\n MinEta = cms.double(" << strTetaZ << "),"
0718 << "\n MaxEta = cms.double(" << strTetaZ << "),"
0719 << "\n MinPhi = cms.double(" << strTphiZr << "),"
0720 << "\n MaxPhi = cms.double(" << strTphiZr << "),"
0721 << "\n Psi = cms.double(" << strTpsir << "),"
0722 << "\n BeamMeanX = cms.double(" << strTxpos << "),"
0723 << "\n BeamMeanY = cms.double(" << strTypos << "),"
0724 << "\n BeamPosition = cms.double(" << strTzpos << ")\n )" << std::endl;
0725
0726 beamfile.close();
0727 }
0728 }
0729 }
0730 }
0731 }
0732
0733 template <class T>
0734 unsigned int CaloGeometryLoaderTest<T>::getDetIdForDDDNode(const DDFilteredView& fv) {
0735
0736
0737
0738 const DDGeoHistory& parents(fv.geoHistory());
0739 const DDGeoHistory::size_type psize(parents.size());
0740
0741 EcalBaseNumber baseNumber;
0742 baseNumber.setSize(psize);
0743
0744 for (unsigned int i = 1; i <= psize; ++i) {
0745 baseNumber.addLevel(parents[psize - i].logicalPart().name().name(), parents[psize - i].copyno());
0746 }
0747
0748 return baseNumber.getCopyNumber(0);
0749
0750
0751 }
0752
0753 #endif