Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-05-06 02:07:23 

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     edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::radiusX:Try 2X:" << (slope(*zb2, zFront2, slope2) < tol_) << ":"
0112                                   << (std::abs(*zb2 - zf) < tol_);
0113 #endif
0114     zz.emplace_back(zf);
0115     rin.emplace_back(radius(zf, zFront1, rFront1, slope1));
0116     rout.emplace_back(radius(zf, zFront2, rFront2, slope2));
0117     if (slope(*zb2, zFront2, slope2) < tol_) {
0118       if (std::abs(*zb2 - zf) > tol_) {
0119         zz.emplace_back(*zb2);
0120         rin.emplace_back(radius(*zb2, zFront1, rFront1, slope1));
0121         rout.emplace_back(radius(*zb2 - tol_, zFront2, rFront2, slope2));
0122       }
0123     }
0124     if ((std::abs(*zb2 - zb) > tol_) && (std::abs(*zb2 - zf) > tol_)) {
0125       zz.emplace_back(*zb2);
0126       rin.emplace_back(radius(*zb2, zFront1, rFront1, slope1));
0127       rout.emplace_back(radius(*zb2, zFront2, rFront2, slope2));
0128     }
0129     zz.emplace_back(zb);
0130     rin.emplace_back(radius(zb, zFront1, rFront1, slope1));
0131     rout.emplace_back(radius(zb, zFront2, rFront2, slope2));
0132   } else if (zf2 == zb2) {
0133 #ifdef EDM_ML_DEBUG
0134     edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::radiusX:Try 3: " << zf << ":" << *zb1 << ":" << zb << " Test "
0135                                   << (slope(*zb1, zFront1, slope1) < tol_) << ":"
0136                                   << ((std::abs(*zb1 - zb) > tol_) && (std::abs(*zb1 - zf) > tol_));
0137 #endif
0138     zz.emplace_back(zf);
0139     rin.emplace_back(radius(zf, zFront1, rFront1, slope1));
0140     rout.emplace_back(radius(zf, zFront2, rFront2, slope2));
0141     if (slope(*zb1, zFront1, slope1) < tol_) {
0142       zz.emplace_back(*zb1);
0143       rin.emplace_back(radius(*zb1 - tol_, zFront1, rFront1, slope1));
0144       rout.emplace_back(radius(*zb1, zFront2, rFront2, slope2));
0145     }
0146     if ((std::abs(*zb1 - zb) > tol_) && (std::abs(*zb1 - zf) > tol_)) {
0147       zz.emplace_back(*zb1);
0148       rin.emplace_back(radius(*zb1, zFront1, rFront1, slope1));
0149       rout.emplace_back(radius(*zb1, zFront2, rFront2, slope2));
0150     }
0151     zz.emplace_back(zb);
0152     rin.emplace_back(radius(zb, zFront1, rFront1, slope1));
0153     rout.emplace_back(radius(zb, zFront2, rFront2, slope2));
0154   } else {
0155     double z1 = std::min(*zf2, *zb1);
0156     double z2 = std::max(*zf2, *zb1);
0157 #ifdef EDM_ML_DEBUG
0158     edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::radiusX:Try 4: " << zf << ":" << z1 << " : " << z2 << ":" << zb;
0159 #endif
0160     zz.emplace_back(zf);
0161     rin.emplace_back(radius(zf, zFront1, rFront1, slope1));
0162     rout.emplace_back(radius(zf, zFront2, rFront2, slope2));
0163     zz.emplace_back(z1);
0164     rin.emplace_back(radius(z1, zFront1, rFront1, slope1));
0165     rout.emplace_back(radius(z1, zFront2, rFront2, slope2));
0166     zz.emplace_back(z2);
0167     rin.emplace_back(radius(z2, zFront1, rFront1, slope1));
0168     rout.emplace_back(radius(z2, zFront2, rFront2, slope2));
0169     zz.emplace_back(zb);
0170     rin.emplace_back(radius(zb, zFront1, rFront1, slope1));
0171     rout.emplace_back(radius(zb, zFront2, rFront2, slope2));
0172   }
0173   double rmin = *(std::min_element(rout.begin(), rout.end()));
0174   if (flag > 1) {
0175     for (auto& rr : rout)
0176       rr = rmin;
0177   }
0178 #ifdef EDM_ML_DEBUG
0179   edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::radiusX has " << zz.size() << " sections: " << rmin;
0180   for (unsigned int k = 0; k < zz.size(); ++k)
0181     edm::LogVerbatim("HGCalGeom") << "[" << k << "] Z = " << zz[k] << " R = " << rin[k] << ":" << rout[k];
0182 #endif
0183 }
0184 
0185 double HGCalGeomTools::radius(double z,
0186                               std::vector<double> const& zFront,
0187                               std::vector<double> const& rFront,
0188                               std::vector<double> const& slope) {
0189 #ifdef EDM_ML_DEBUG
0190   std::ostringstream st1;
0191   st1 << "Z = " << z << "; zFront =";
0192   for (const auto& v : zFront)
0193     st1 << " " << v;
0194   st1 << "; rFront =";
0195   for (const auto& v : rFront)
0196     st1 << " " << v;
0197   st1 << "; slope =";
0198   for (const auto& v : slope)
0199     st1 << " " << v;
0200   edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::radius: " << st1.str();
0201 #endif
0202   auto itrz = std::lower_bound(zFront.begin(), zFront.end(), z);
0203   if (itrz != zFront.begin())
0204     --itrz;
0205   unsigned int ik = static_cast<unsigned int>(itrz - zFront.begin());
0206   if ((ik + 1) < zFront.size() && std::abs(z - zFront[ik + 1]) < tol_)
0207     ++ik;
0208   double r = rFront[ik] + (z - zFront[ik]) * slope[ik];
0209 #ifdef EDM_ML_DEBUG
0210   edm::LogVerbatim("HGCalGeom") << "DDHGCalGeomTools: Z " << z << " k " << ik << " R " << r;
0211 #endif
0212   return r;
0213 }
0214 
0215 double HGCalGeomTools::radius(
0216     double z, int layer0, int layerf, std::vector<double> const& zFront, std::vector<double> const& rFront) {
0217   double r = rFront[0];
0218 #ifdef EDM_ML_DEBUG
0219   unsigned int ik(0);
0220 #endif
0221   for (unsigned int k = 0; k < rFront.size(); ++k) {
0222     int k1 = layerf - layer0 + (int)(k);
0223     if (k1 < (int)(zFront.size())) {
0224       r = rFront[k];
0225 #ifdef EDM_ML_DEBUG
0226       ik = k;
0227 #endif
0228       if (z < zFront[k1] + tol_)
0229         break;
0230     }
0231   }
0232 #ifdef EDM_ML_DEBUG
0233   edm::LogVerbatim("HGCalGeom") << "DDHGCalGeomTools: Z " << z << ":" << ik << " R " << r;
0234 #endif
0235   return r;
0236 }
0237 
0238 std::pair<double, double> HGCalGeomTools::shiftXY(int waferPosition, double waferSize) const {
0239   double dx(0), dy(0);
0240   switch (waferPosition) {
0241     case (HGCalTypes::CornerCenterYp): {
0242       dy = factor_ * waferSize;
0243       break;
0244     }
0245     case (HGCalTypes::CornerCenterYm): {
0246       dy = -factor_ * waferSize;
0247       break;
0248     }
0249     case (HGCalTypes::CornerCenterXp): {
0250       dx = factor_ * waferSize;
0251       break;
0252     }
0253     case (HGCalTypes::CornerCenterXm): {
0254       dx = -factor_ * waferSize;
0255       break;
0256     }
0257   }
0258 #ifdef EDM_ML_DEBUG
0259   edm::LogVerbatim("HGCalGeom") << "Shift for " << waferPosition << " is (" << dx << ":" << dy << ")";
0260 #endif
0261   return std::make_pair(dx, dy);
0262 }
0263 
0264 double HGCalGeomTools::slope(double z, std::vector<double> const& zFront, std::vector<double> const& slope) {
0265   auto itrz = std::lower_bound(zFront.begin(), zFront.end(), z);
0266   if (itrz != zFront.begin())
0267     --itrz;
0268   unsigned int ik = static_cast<unsigned int>(itrz - zFront.begin());
0269   //  if (ik < zFront.size() && std::abs(z-zFront[ik+1]) < tol_) ++ik;
0270 #ifdef EDM_ML_DEBUG
0271   edm::LogVerbatim("HGCalGeom") << "HGCalGeomTools::slope:z " << z << " k " << ik << " Slope " << slope[ik];
0272 #endif
0273   return slope[ik];
0274 }
0275 
0276 std::pair<double, double> HGCalGeomTools::zradius(double z1,
0277                                                   double z2,
0278                                                   std::vector<double> const& zF,
0279                                                   std::vector<double> const& rF) {
0280   double z(-1), r(-1);
0281   for (unsigned int k = 0; k < rF.size(); ++k) {
0282     if ((z1 > zF[k] - tol_) && (z2 < zF[k] + tol_)) {
0283       z = zF[k];
0284       r = rF[k];
0285       break;
0286     }
0287   }
0288   return std::make_pair(z, r);
0289 }
0290 
0291 std::pair<int32_t, int32_t> HGCalGeomTools::waferCorner(
0292     double xpos, double ypos, double r, double R, double rMin, double rMax, bool oldBug) {
0293   double xc[HGCalParameters::k_CornerSize], yc[HGCalParameters::k_CornerSize];
0294   xc[0] = xpos;
0295   yc[0] = ypos + R;
0296   xc[1] = xpos - r;
0297   yc[1] = ypos + 0.5 * R;
0298   if (oldBug) {
0299     xc[2] = xpos + r;
0300     yc[2] = ypos - 0.5 * R;
0301   } else {
0302     xc[2] = xpos - r;
0303     yc[2] = ypos - 0.5 * R;
0304   }
0305   xc[3] = xpos;
0306   yc[3] = ypos - R;
0307   xc[4] = xpos + r;
0308   yc[4] = ypos - 0.5 * R;
0309   xc[5] = xpos + r;
0310   yc[5] = ypos + 0.5 * R;
0311   int32_t nCorner(0), firstCorner(-1), firstMiss(-1);
0312 #ifdef EDM_ML_DEBUG
0313   std::vector<uint32_t> corners;
0314 #endif
0315   for (uint32_t k = 0; k < HGCalParameters::k_CornerSize; ++k) {
0316     double rpos = sqrt(xc[k] * xc[k] + yc[k] * yc[k]);
0317     if ((rpos <= rMax) && (rpos >= rMin)) {
0318 #ifdef EDM_ML_DEBUG
0319       corners.emplace_back(k);
0320 #endif
0321       if (firstCorner < 0)
0322         firstCorner = k;
0323       ++nCorner;
0324     } else {
0325       if (firstMiss < 0)
0326         firstMiss = k;
0327     }
0328   }
0329   if ((nCorner > 1) && (firstCorner == 0) && (firstMiss < nCorner)) {
0330     firstCorner = firstMiss + HGCalParameters::k_CornerSize - nCorner;
0331   }
0332 #ifdef EDM_ML_DEBUG
0333   edm::LogVerbatim("HGCalGeom") << "waferCorner:: R " << rMin << ":" << rMax << nCorner << " corners; first corner "
0334                                 << firstCorner;
0335   for (uint32_t k = 0; k < HGCalParameters::k_CornerSize; ++k) {
0336     double rpos = std::sqrt(xc[k] * xc[k] + yc[k] * yc[k]);
0337     std::string ok = (std::find(corners.begin(), corners.end(), k) != corners.end()) ? " In" : " Out";
0338     edm::LogVerbatim("HGCalGeom") << "Corner[" << k << "] x " << xc[k] << " y " << yc[k] << " R " << rpos << ok;
0339   }
0340 #endif
0341   return std::make_pair(nCorner, firstCorner);
0342 }
[ 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!