Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:15:04 

0001 #include "Geometry/HGCalCommonData/interface/HGCalGeomTools.h"
0002 #include "Geometry/HGCalCommonData/interface/HGCalTypes.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "Geometry/HGCalCommonData/interface/HGCalParameters.h"
0005 
0006 #include <sstream>
0007 #include <string>
0008 
0009 //#define EDM_ML_DEBUG
0010 HGCalGeomTools::HGCalGeomTools() : factor_(1.0 / std::sqrt(3.0)) {}
0011 
0012 void HGCalGeomTools::radius(double zf,
0013                             double zb,
0014                             std::vector<double> const& zFront1,
0015                             std::vector<double> const& rFront1,
0016                             std::vector<double> const& slope1,
0017                             std::vector<double> const& zFront2,
0018                             std::vector<double> const& rFront2,
0019                             std::vector<double> const& slope2,
0020                             int flag,
0021                             std::vector<double>& zz,
0022                             std::vector<double>& rin,
0023                             std::vector<double>& rout) {
0024 #ifdef EDM_ML_DEBUG
0025   std::ostringstream st1;
0026   st1 << "Z = " << zf << ":" << zb << "; zFront1 =";
0027   for (const auto& v : zFront1)
0028     st1 << " " << v;
0029   st1 << "; rFront1 =";
0030   for (const auto& v : rFront1)
0031     st1 << " " << v;
0032   st1 << "; slope1 =";
0033   for (const auto& v : slope1)
0034     st1 << " " << v;
0035   st1 << "; zFront2 =";
0036   for (const auto& v : zFront2)
0037     st1 << " " << v;
0038   st1 << "; rFront2 =";
0039   for (const auto& v : rFront2)
0040     st1 << " " << v;
0041   st1 << "; slope2 =";
0042   for (const auto& v : slope2)
0043     st1 << " " << v;
0044   edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::radiusX: " << st1.str();
0045   double dz2f(0), dz2b(0);
0046 #endif
0047   double dz1f(0), dz1b(0);
0048   auto zf1 = std::lower_bound(zFront1.begin(), zFront1.end(), zf);
0049   if (zf1 != zFront1.begin())
0050     --zf1;
0051   if (((zf1 + 1) != zFront1.end()) && (std::abs(*(zf1 + 1) - zf) < tol_)) {
0052     ++zf1;
0053     dz1f = 2 * tol_;
0054   }
0055   auto zf2 = std::lower_bound(zFront2.begin(), zFront2.end(), zf);
0056   if (zf2 != zFront2.begin())
0057     --zf2;
0058   if (((zf2 + 1) != zFront2.end()) && (std::abs(*(zf2 + 1) - zf) < tol_)) {
0059     if (static_cast<int>(zf2 - zFront2.begin()) >= 2)
0060       ++zf2;
0061 #ifdef EDM_ML_DEBUG
0062     dz2f = 2 * tol_;
0063 #endif
0064   }
0065   auto zb1 = std::lower_bound(zFront1.begin(), zFront1.end(), zb);
0066   if (zb1 != zFront1.begin())
0067     --zb1;
0068   if ((zb1 != zFront1.begin()) && (std::abs(*zb1 - zb) < tol_)) {
0069     --zb1;
0070     dz1b = -2 * tol_;
0071   }
0072   if (((zb1 + 1) != zFront1.end()) && (std::abs(*(zb1 + 1) - zb) < tol_)) {
0073     dz1b = -2 * tol_;
0074   }
0075   auto zb2 = std::lower_bound(zFront2.begin(), zFront2.end(), zb);
0076   if (zb2 != zFront2.begin()) {
0077     --zb2;
0078   }
0079   if ((zb2 != zFront2.begin()) && (std::abs(*zb2 - zb) < tol_)) {
0080     if (static_cast<int>(zb2 - zFront2.begin()) > 2)
0081       --zb2;
0082 #ifdef EDM_ML_DEBUG
0083     dz2b = -2 * tol_;
0084 #endif
0085   }
0086 #ifdef EDM_ML_DEBUG
0087   edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::radiusX:zf " << zf << " : "
0088                                 << static_cast<int>(zf1 - zFront1.begin()) << ":" << *zf1 << " : "
0089                                 << static_cast<int>(zf2 - zFront2.begin()) << ":" << *zf2 << " zb " << zb << ":"
0090                                 << static_cast<int>(zb1 - zFront1.begin()) << " : " << *zb1 << " : "
0091                                 << static_cast<int>(zb2 - zFront2.begin()) << ":" << *zb2 << " Flag " << flag << ":"
0092                                 << (zf1 == zb1) << ":" << (zf2 == zb2) << " dz " << dz1f << ":" << dz2f << ":" << dz1b
0093                                 << ":" << dz2b;
0094 #endif
0095   if ((zf1 == zb1) && (zf2 == zb2)) {
0096 #ifdef EDM_ML_DEBUG
0097     edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::radiusX:Try 1: " << zf << ":" << zb;
0098 #endif
0099     zz.emplace_back(zf);
0100     rin.emplace_back(radius(zf + dz1f, zFront1, rFront1, slope1));
0101     rout.emplace_back(radius(zf, zFront2, rFront2, slope2));
0102     zz.emplace_back(zb);
0103     rin.emplace_back(radius(zb + dz1b, zFront1, rFront1, slope1));
0104     rout.emplace_back(radius(zb, zFront2, rFront2, slope2));
0105   } else if (zf1 == zb1) {
0106 #ifdef EDM_ML_DEBUG
0107     double z1 = std::max(*zf2, *zb2);
0108     edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::radiusX:Try 2:" << zf << ":" << *zb2 << " (" << z1 << ") : " << zb
0109                                   << " Test " << (slope(*zb2, zFront2, slope2) < tol_) << ":"
0110                                   << ((std::abs(*zb2 - zb) > tol_) && (std::abs(*zb2 - zf) > tol_));
0111 #endif
0112     zz.emplace_back(zf);
0113     rin.emplace_back(radius(zf, zFront1, rFront1, slope1));
0114     rout.emplace_back(radius(zf, zFront2, rFront2, slope2));
0115     if (slope(*zb2, zFront2, slope2) < tol_) {
0116       zz.emplace_back(*zb2);
0117       rin.emplace_back(radius(*zb2, zFront1, rFront1, slope1));
0118       rout.emplace_back(radius(*zb2 - tol_, zFront2, rFront2, slope2));
0119     }
0120     if ((std::abs(*zb2 - zb) > tol_) && (std::abs(*zb2 - zf) > tol_)) {
0121       zz.emplace_back(*zb2);
0122       rin.emplace_back(radius(*zb2, zFront1, rFront1, slope1));
0123       rout.emplace_back(radius(*zb2, zFront2, rFront2, slope2));
0124     }
0125     zz.emplace_back(zb);
0126     rin.emplace_back(radius(zb, zFront1, rFront1, slope1));
0127     rout.emplace_back(radius(zb, zFront2, rFront2, slope2));
0128   } else if (zf2 == zb2) {
0129 #ifdef EDM_ML_DEBUG
0130     edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::radiusX:Try 3: " << zf << ":" << *zb1 << ":" << zb << " Test "
0131                                   << (slope(*zb1, zFront1, slope1) < tol_) << ":"
0132                                   << ((std::abs(*zb1 - zb) > tol_) && (std::abs(*zb1 - zf) > tol_));
0133 #endif
0134     zz.emplace_back(zf);
0135     rin.emplace_back(radius(zf, zFront1, rFront1, slope1));
0136     rout.emplace_back(radius(zf, zFront2, rFront2, slope2));
0137     if (slope(*zb1, zFront1, slope1) < tol_) {
0138       zz.emplace_back(*zb1);
0139       rin.emplace_back(radius(*zb1 - tol_, zFront1, rFront1, slope1));
0140       rout.emplace_back(radius(*zb1, zFront2, rFront2, slope2));
0141     }
0142     if ((std::abs(*zb1 - zb) > tol_) && (std::abs(*zb1 - zf) > tol_)) {
0143       zz.emplace_back(*zb1);
0144       rin.emplace_back(radius(*zb1, zFront1, rFront1, slope1));
0145       rout.emplace_back(radius(*zb1, zFront2, rFront2, slope2));
0146     }
0147     zz.emplace_back(zb);
0148     rin.emplace_back(radius(zb, zFront1, rFront1, slope1));
0149     rout.emplace_back(radius(zb, zFront2, rFront2, slope2));
0150   } else {
0151     double z1 = std::min(*zf2, *zb1);
0152     double z2 = std::max(*zf2, *zb1);
0153 #ifdef EDM_ML_DEBUG
0154     edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::radiusX:Try 4: " << zf << ":" << z1 << " : " << z2 << ":" << zb;
0155 #endif
0156     zz.emplace_back(zf);
0157     rin.emplace_back(radius(zf, zFront1, rFront1, slope1));
0158     rout.emplace_back(radius(zf, zFront2, rFront2, slope2));
0159     zz.emplace_back(z1);
0160     rin.emplace_back(radius(z1, zFront1, rFront1, slope1));
0161     rout.emplace_back(radius(z1, zFront2, rFront2, slope2));
0162     zz.emplace_back(z2);
0163     rin.emplace_back(radius(z2, zFront1, rFront1, slope1));
0164     rout.emplace_back(radius(z2, zFront2, rFront2, slope2));
0165     zz.emplace_back(zb);
0166     rin.emplace_back(radius(zb, zFront1, rFront1, slope1));
0167     rout.emplace_back(radius(zb, zFront2, rFront2, slope2));
0168   }
0169   double rmin = *(std::min_element(rout.begin(), rout.end()));
0170   if (flag > 1) {
0171     for (auto& rr : rout)
0172       rr = rmin;
0173   }
0174 #ifdef EDM_ML_DEBUG
0175   edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::radiusX has " << zz.size() << " sections: " << rmin;
0176   for (unsigned int k = 0; k < zz.size(); ++k)
0177     edm::LogVerbatim("HGCalGeom") << "[" << k << "] Z = " << zz[k] << " R = " << rin[k] << ":" << rout[k];
0178 #endif
0179 }
0180 
0181 double HGCalGeomTools::radius(double z,
0182                               std::vector<double> const& zFront,
0183                               std::vector<double> const& rFront,
0184                               std::vector<double> const& slope) {
0185 #ifdef EDM_ML_DEBUG
0186   std::ostringstream st1;
0187   st1 << "Z = " << z << "; zFront =";
0188   for (const auto& v : zFront)
0189     st1 << " " << v;
0190   st1 << "; rFront =";
0191   for (const auto& v : rFront)
0192     st1 << " " << v;
0193   st1 << "; slope =";
0194   for (const auto& v : slope)
0195     st1 << " " << v;
0196   edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::radius: " << st1.str();
0197 #endif
0198   auto itrz = std::lower_bound(zFront.begin(), zFront.end(), z);
0199   if (itrz != zFront.begin())
0200     --itrz;
0201   unsigned int ik = static_cast<unsigned int>(itrz - zFront.begin());
0202   if ((ik + 1) < zFront.size() && std::abs(z - zFront[ik + 1]) < tol_)
0203     ++ik;
0204   double r = rFront[ik] + (z - zFront[ik]) * slope[ik];
0205 #ifdef EDM_ML_DEBUG
0206   edm::LogVerbatim("HGCalGeom") << "DDHGCalGeomTools: Z " << z << " k " << ik << " R " << r;
0207 #endif
0208   return r;
0209 }
0210 
0211 double HGCalGeomTools::radius(
0212     double z, int layer0, int layerf, std::vector<double> const& zFront, std::vector<double> const& rFront) {
0213   double r = rFront[0];
0214 #ifdef EDM_ML_DEBUG
0215   unsigned int ik(0);
0216 #endif
0217   for (unsigned int k = 0; k < rFront.size(); ++k) {
0218     int k1 = layerf - layer0 + (int)(k);
0219     if (k1 < (int)(zFront.size())) {
0220       r = rFront[k];
0221 #ifdef EDM_ML_DEBUG
0222       ik = k;
0223 #endif
0224       if (z < zFront[k1] + tol_)
0225         break;
0226     }
0227   }
0228 #ifdef EDM_ML_DEBUG
0229   edm::LogVerbatim("HGCalGeom") << "DDHGCalGeomTools: Z " << z << ":" << ik << " R " << r;
0230 #endif
0231   return r;
0232 }
0233 
0234 std::pair<double, double> HGCalGeomTools::shiftXY(int waferPosition, double waferSize) const {
0235   double dx(0), dy(0);
0236   switch (waferPosition) {
0237     case (HGCalTypes::CornerCenterYp): {
0238       dy = factor_ * waferSize;
0239       break;
0240     }
0241     case (HGCalTypes::CornerCenterYm): {
0242       dy = -factor_ * waferSize;
0243       break;
0244     }
0245     case (HGCalTypes::CornerCenterXp): {
0246       dx = factor_ * waferSize;
0247       break;
0248     }
0249     case (HGCalTypes::CornerCenterXm): {
0250       dx = -factor_ * waferSize;
0251       break;
0252     }
0253   }
0254 #ifdef EDM_ML_DEBUG
0255   edm::LogVerbatim("HGCalGeom") << "Shift for " << waferPosition << " is (" << dx << ":" << dy << ")";
0256 #endif
0257   return std::make_pair(dx, dy);
0258 }
0259 
0260 double HGCalGeomTools::slope(double z, std::vector<double> const& zFront, std::vector<double> const& slope) {
0261   auto itrz = std::lower_bound(zFront.begin(), zFront.end(), z);
0262   if (itrz != zFront.begin())
0263     --itrz;
0264   unsigned int ik = static_cast<unsigned int>(itrz - zFront.begin());
0265   //  if (ik < zFront.size() && std::abs(z-zFront[ik+1]) < tol_) ++ik;
0266 #ifdef EDM_ML_DEBUG
0267   edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::slope:z " << z << " k " << ik << " Slope " << slope[ik];
0268 #endif
0269   return slope[ik];
0270 }
0271 
0272 std::pair<double, double> HGCalGeomTools::zradius(double z1,
0273                                                   double z2,
0274                                                   std::vector<double> const& zF,
0275                                                   std::vector<double> const& rF) {
0276   double z(-1), r(-1);
0277   for (unsigned int k = 0; k < rF.size(); ++k) {
0278     if ((z1 > zF[k] - tol_) && (z2 < zF[k] + tol_)) {
0279       z = zF[k];
0280       r = rF[k];
0281       break;
0282     }
0283   }
0284   return std::make_pair(z, r);
0285 }
0286 
0287 std::pair<int32_t, int32_t> HGCalGeomTools::waferCorner(
0288     double xpos, double ypos, double r, double R, double rMin, double rMax, bool oldBug) {
0289   double xc[HGCalParameters::k_CornerSize], yc[HGCalParameters::k_CornerSize];
0290   xc[0] = xpos;
0291   yc[0] = ypos + R;
0292   xc[1] = xpos - r;
0293   yc[1] = ypos + 0.5 * R;
0294   if (oldBug) {
0295     xc[2] = xpos + r;
0296     yc[2] = ypos - 0.5 * R;
0297   } else {
0298     xc[2] = xpos - r;
0299     yc[2] = ypos - 0.5 * R;
0300   }
0301   xc[3] = xpos;
0302   yc[3] = ypos - R;
0303   xc[4] = xpos + r;
0304   yc[4] = ypos - 0.5 * R;
0305   xc[5] = xpos + r;
0306   yc[5] = ypos + 0.5 * R;
0307   int32_t nCorner(0), firstCorner(-1), firstMiss(-1);
0308 #ifdef EDM_ML_DEBUG
0309   std::vector<uint32_t> corners;
0310 #endif
0311   for (uint32_t k = 0; k < HGCalParameters::k_CornerSize; ++k) {
0312     double rpos = sqrt(xc[k] * xc[k] + yc[k] * yc[k]);
0313     if ((rpos <= rMax) && (rpos >= rMin)) {
0314 #ifdef EDM_ML_DEBUG
0315       corners.emplace_back(k);
0316 #endif
0317       if (firstCorner < 0)
0318         firstCorner = k;
0319       ++nCorner;
0320     } else {
0321       if (firstMiss < 0)
0322         firstMiss = k;
0323     }
0324   }
0325   if ((nCorner > 1) && (firstCorner == 0) && (firstMiss < nCorner)) {
0326     firstCorner = firstMiss + HGCalParameters::k_CornerSize - nCorner;
0327   }
0328 #ifdef EDM_ML_DEBUG
0329   edm::LogVerbatim("HGCalGeom") << "waferCorner:: R " << rMin << ":" << rMax << nCorner << " corners; first corner "
0330                                 << firstCorner;
0331   for (uint32_t k = 0; k < HGCalParameters::k_CornerSize; ++k) {
0332     double rpos = std::sqrt(xc[k] * xc[k] + yc[k] * yc[k]);
0333     std::string ok = (std::find(corners.begin(), corners.end(), k) != corners.end()) ? " In" : " Out";
0334     edm::LogVerbatim("HGCalGeom") << "Corner[" << k << "] x " << xc[k] << " y " << yc[k] << " R " << rpos << ok;
0335   }
0336 #endif
0337   return std::make_pair(nCorner, firstCorner);
0338 }
[ 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!