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