Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 #include "Geometry/HGCalCommonData/interface/HGCalCellUV.h"
0002 #include "Geometry/HGCalCommonData/interface/HGCalTypes.h"
0003 #include "Geometry/HGCalCommonData/interface/HGCalWaferMask.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include <iostream>
0006 #include <array>
0007 #include <algorithm>
0008 #include <cassert>
0009 
0010 HGCalCellUV::HGCalCellUV(double waferSize, double separation, int32_t nFine, int32_t nCoarse) : waferSize_(waferSize) {
0011   hgcalcell_ = std::make_unique<HGCalCell>(waferSize, nFine, nCoarse);
0012   assert(nFine > 0 && nCoarse > 0);
0013   ncell_[0] = nFine;
0014   ncell_[1] = nCoarse;
0015   for (int k = 0; k < 2; ++k) {
0016     cellX_[k] = waferSize / (3 * ncell_[k]);
0017     cellY_[k] = 0.5 * sqrt3_ * cellX_[k];
0018     cellXTotal_[k] = (waferSize + separation) / (3 * ncell_[k]);
0019     cellYTotal_[k] = 0.5 * sqrt3_ * cellXTotal_[k];
0020   }
0021   // Fill up the placement vectors
0022   for (int placement = 0; placement < HGCalCell::cellPlacementTotal; ++placement) {
0023     // Fine cells
0024     for (int iu = 0; iu < 2 * ncell_[0]; ++iu) {
0025       for (int iv = 0; iv < 2 * ncell_[0]; ++iv) {
0026         int u = (placement < HGCalCell::cellPlacementExtra) ? iv : iu;
0027         int v = (placement < HGCalCell::cellPlacementExtra) ? iu : iv;
0028         if (((v - u) < ncell_[0]) && (u - v) <= ncell_[0]) {
0029           cellPosFine_[placement][std::pair<int, int>(u, v)] = hgcalcell_->cellUV2XY1(u, v, placement, 0);
0030         }
0031       }
0032     }
0033     // Coarse cells
0034     for (int iu = 0; iu < 2 * ncell_[1]; ++iu) {
0035       for (int iv = 0; iv < 2 * ncell_[1]; ++iv) {
0036         int u = (placement < HGCalCell::cellPlacementExtra) ? iv : iu;
0037         int v = (placement < HGCalCell::cellPlacementExtra) ? iu : iv;
0038         if (((v - u) < ncell_[1]) && (u - v) <= ncell_[1]) {
0039           cellPosCoarse_[placement][std::pair<int, int>(u, v)] = hgcalcell_->cellUV2XY1(u, v, placement, 1);
0040         }
0041       }
0042     }
0043   }
0044 }
0045 
0046 std::pair<int32_t, int32_t> HGCalCellUV::cellUVFromXY1(
0047     double xloc, double yloc, int32_t placement, int32_t type, bool extend, bool debug) const {
0048   if (type != 0)
0049     type = 1;
0050   //--- Reverse transform to placement=0, if placement index ≠ 6
0051   double xloc1 = (placement >= HGCalCell::cellPlacementExtra) ? xloc : -xloc;
0052   double cellY = (extend) ? cellYTotal_[type] : cellY_[type];
0053   int rot = placement % HGCalCell::cellPlacementExtra;
0054   static constexpr std::array<double, 6> fcos = {{1.0, cos60_, -cos60_, -1.0, -cos60_, cos60_}};
0055   static constexpr std::array<double, 6> fsin = {{0.0, sin60_, sin60_, 0.0, -sin60_, -sin60_}};
0056   double x = xloc1 * fcos[rot] - yloc * fsin[rot];
0057   double y = xloc1 * fsin[rot] + yloc * fcos[rot];
0058 
0059   //--- Calculate coordinates in u,v,w system
0060   double u = x * sin60_ + y * cos60_;
0061   double v = -x * sin60_ + y * cos60_;
0062   double w = y;
0063 
0064   //--- Rounding in u, v, w coordinates
0065   int iu = std::floor(u / cellY) + 3 * (ncell_[type]) - 1;
0066   int iv = std::floor(v / cellY) + 3 * (ncell_[type]);
0067   int iw = std::floor(w / cellY) + 1;
0068 
0069   int isv = (iu + iw) / 3;
0070   int isu = (iv + iw) / 3;
0071 
0072   //--- Taking care of extending cells
0073   if ((iu + iw) < 0) {
0074     isu = (iv + iw + 1) / 3;
0075     isv = 0;
0076   } else if (isv - isu > ncell_[type] - 1) {
0077     isu = (iv + iw + 1) / 3;
0078     isv = (iu + iw - 1) / 3;
0079   } else if (isu > 2 * ncell_[type] - 1) {
0080     isu = 2 * ncell_[type] - 1;
0081     isv = (iu + iw - 1) / 3;
0082   }
0083   if (debug)
0084     edm::LogVerbatim("HGCalGeom") << "cellUVFromXY1: Input " << xloc << ":" << yloc << ":" << extend << " Output "
0085                                   << isu << ":" << isv;
0086   return std::make_pair(isu, isv);
0087 }
0088 
0089 std::pair<int32_t, int32_t> HGCalCellUV::cellUVFromXY2(
0090     double xloc, double yloc, int32_t placement, int32_t type, bool extend, bool debug) const {
0091   //--- Using multiple inequalities to find (u, v)
0092   //--- Reverse transform to placement=0, if placement index ≠ 7
0093   if (type != 0)
0094     type = 1;
0095   double cellY = (extend) ? cellYTotal_[type] : cellY_[type];
0096   double xloc1 = (placement >= HGCalCell::cellPlacementExtra) ? xloc : -1 * xloc;
0097   int rot = placement % HGCalCell::cellPlacementExtra;
0098   static constexpr std::array<double, 6> fcos = {{cos60_, 1.0, cos60_, -cos60_, -1.0, -cos60_}};
0099   static constexpr std::array<double, 6> fsin = {{-sin60_, 0.0, sin60_, sin60_, 0.0, -sin60_}};
0100   double x = xloc1 * fcos[rot] - yloc * fsin[rot];
0101   double y = xloc1 * fsin[rot] + yloc * fcos[rot];
0102 
0103   int32_t u(-100), v(-100);
0104   int ncell = ncell_[type];
0105   double r = cellY;
0106   double l1 = (y / r) + ncell - 1.0;
0107   int l2 = std::floor((0.5 * y + 0.5 * x / sqrt3_) / r + ncell - 4.0 / 3.0);
0108   int l3 = std::floor((x / sqrt3_) / r + ncell - 4.0 / 3.0);
0109   double l4 = (y + sqrt3_ * x) / (2 * r) + 2 * ncell - 2;
0110   double l5 = (y - sqrt3_ * x) / (2 * r) - ncell;
0111   double u1 = (y / r) + ncell + 1.0;
0112   int u2 = std::ceil((0.5 * y + 0.5 * x / sqrt3_) / r + ncell + 2.0 / 3.0);
0113   int u3 = std::ceil((x / sqrt3_) / r + ncell);
0114   double u4 = l4 + 2;
0115   double u5 = l5 + 2;
0116 
0117   for (int ui = l2 + 1; ui < u2; ui++) {
0118     for (int vi = l3 + 1; vi < u3; vi++) {
0119       int c1 = 2 * ui - vi;
0120       int c4 = ui + vi;
0121       int c5 = ui - 2 * vi;
0122       if ((c1 < u1) && (c1 > l1) && (c4 < u4) && (c4 > l4) && (c5 < u5) && (c5 > l5)) {
0123         u = ui;
0124         v = vi;
0125       }
0126     }
0127   }
0128 
0129   //--- Taking care of extending cells
0130   if (v == -1) {
0131     if (y < (2 * u - v - ncell) * r) {
0132       ++v;
0133     } else {
0134       ++u;
0135       ++v;
0136     }
0137   }
0138   if (v - u == ncell) {
0139     if ((y + sqrt3_ * x) < ((u + v - 2 * ncell + 1) * 2 * r)) {
0140       --v;
0141     } else {
0142       ++u;
0143     }
0144   }
0145   if (u == 2 * ncell) {
0146     if ((y - sqrt3_ * x) < ((u - 2 * v + ncell - 1) * 2 * r)) {
0147       --u;
0148     } else {
0149       --u;
0150       --v;
0151     }
0152   }
0153   if (debug)
0154     edm::LogVerbatim("HGCalGeom") << "cellUVFromXY2: Input " << xloc << ":" << yloc << ":" << extend << " Output " << u
0155                                   << ":" << v;
0156   return std::make_pair(u, v);
0157 }
0158 
0159 std::pair<int32_t, int32_t> HGCalCellUV::cellUVFromXY3(
0160     double xloc, double yloc, int32_t placement, int32_t type, bool extend, bool debug) const {
0161   //--- Using Cube coordinates to find the (u, v)
0162   //--- Reverse transform to placement=0, if placement index ≠ 6
0163   if (type != 0)
0164     type = 1;
0165   double cellX = (extend) ? cellXTotal_[type] : cellX_[type];
0166   double cellY = (extend) ? cellYTotal_[type] : cellY_[type];
0167   double xloc1 = (placement >= HGCalCell::cellPlacementExtra) ? xloc : -xloc;
0168   int rot = placement % HGCalCell::cellPlacementExtra;
0169   static constexpr std::array<double, 6> fcos = {{1.0, cos60_, -cos60_, -1.0, -cos60_, cos60_}};
0170   static constexpr std::array<double, 6> fsin = {{0.0, sin60_, sin60_, 0.0, -sin60_, -sin60_}};
0171   double xprime = xloc1 * fcos[rot] - yloc * fsin[rot];
0172   double yprime = xloc1 * fsin[rot] + yloc * fcos[rot];
0173   double x = xprime + cellX;
0174   double y = yprime;
0175 
0176   x = x / cellX;
0177   y = y / cellY;
0178 
0179   double cu = 2 * x / 3;
0180   double cv = -x / 3 + y / 2;
0181   double cw = -x / 3 - y / 2;
0182 
0183   int iu = std::round(cu);
0184   int iv = std::round(cv);
0185   int iw = std::round(cw);
0186 
0187   if (iu + iv + iw != 0) {
0188     double arr[] = {std::abs(cu - iu), std::abs(cv - iv), std::abs(cw - iw)};
0189     int i = std::distance(arr, std::max_element(arr, arr + 3));
0190 
0191     if (i == 1)
0192       iv = (std::round(cv) == std::floor(cv)) ? std::ceil(cv) : std::floor(cv);
0193     else if (i == 2)
0194       iw = (std::round(cw) == std::floor(cw)) ? std::ceil(cw) : std::floor(cw);
0195   }
0196 
0197   //--- Taking care of extending cells
0198   int u(ncell_[type] + iv), v(ncell_[type] - 1 - iw);
0199   double xcell = (1.5 * (v - u) + 0.5) * cellX_[type];
0200   double ycell = (v + u - 2 * ncell_[type] + 1) * cellY;
0201   if (v == -1) {
0202     if ((yprime - sqrt3_ * xprime) < (ycell - sqrt3_ * xcell)) {
0203       ++v;
0204     } else {
0205       ++u;
0206       ++v;
0207     }
0208   }
0209   if (v - u == ncell_[type]) {
0210     if (yprime < ycell) {
0211       --v;
0212     } else {
0213       ++u;
0214     }
0215   }
0216   if (u == 2 * ncell_[type]) {
0217     if ((yprime + sqrt3_ * xprime) > (ycell + sqrt3_ * xcell)) {
0218       --u;
0219     } else {
0220       --u;
0221       --v;
0222     }
0223   }
0224 
0225   if (debug)
0226     edm::LogVerbatim("HGCalGeom") << "cellUVFromXY3: Input " << xloc << ":" << yloc << ":" << extend << " Output " << u
0227                                   << ":" << v;
0228   return std::make_pair(u, v);
0229 }
0230 
0231 std::pair<int, int> HGCalCellUV::cellUVFromXY4(
0232     double xloc, double yloc, int32_t placement, int32_t type, bool extend, bool debug) {
0233   if (type != 0)
0234     return cellUVFromXY4(
0235         xloc, yloc, ncell_[1], cellX_[1], cellY_[1], cellXTotal_[1], cellY_[1], cellPosCoarse_[placement], extend, debug);
0236   else
0237     return cellUVFromXY4(
0238         xloc, yloc, ncell_[0], cellX_[0], cellY_[0], cellXTotal_[0], cellY_[0], cellPosFine_[placement], extend, debug);
0239 }
0240 
0241 std::pair<int, int> HGCalCellUV::cellUVFromXY4(double xloc,
0242                                                double yloc,
0243                                                int n,
0244                                                double cellX,
0245                                                double cellY,
0246                                                double cellXTotal,
0247                                                double cellYTotal,
0248                                                std::map<std::pair<int, int>, std::pair<double, double> >& cellPos,
0249                                                bool extend,
0250                                                bool debug) {
0251   std::pair<int, int> uv = std::make_pair(-1, -1);
0252   std::map<std::pair<int, int>, std::pair<double, double> >::const_iterator itr;
0253   for (itr = cellPos.begin(); itr != cellPos.end(); ++itr) {
0254     double delX = std::abs(xloc - (itr->second).first);
0255     double delY = std::abs(yloc - (itr->second).second);
0256     if ((delX < cellX) && (delY < cellY)) {
0257       if ((delX < (0.5 * cellX)) || (delY < (2.0 * cellY - sqrt3_ * delX))) {
0258         uv = itr->first;
0259         break;
0260       }
0261     }
0262   }
0263   if ((uv.first < 0) && extend) {
0264     for (itr = cellPos.begin(); itr != cellPos.end(); ++itr) {
0265       double delX = std::abs(xloc - (itr->second).first);
0266       double delY = std::abs(yloc - (itr->second).second);
0267       if ((delX < cellXTotal) && (delY < cellYTotal)) {
0268         if ((delX < (0.5 * cellXTotal)) || (delY < (2.0 * cellYTotal - sqrt3_ * delX))) {
0269           uv = itr->first;
0270           break;
0271         }
0272       }
0273     }
0274   }
0275   if (debug)
0276     edm::LogVerbatim("HGCalGeom") << "cellUVFromXY4: Input " << xloc << ":" << yloc << ":" << extend << " Output "
0277                                   << uv.first << ":" << uv.second;
0278   return uv;
0279 }
0280 
0281 std::pair<int32_t, int32_t> HGCalCellUV::cellUVFromXY1(  // for v17
0282     double xloc,
0283     double yloc,
0284     int32_t placement,
0285     int32_t type,
0286     int32_t partial,
0287     bool extend,
0288     bool debug) const {
0289   if (type != 0)
0290     type = 1;
0291   double cellX = (extend) ? cellXTotal_[type] : cellX_[type];
0292   double cellY = (extend) ? cellYTotal_[type] : cellY_[type];
0293   std::pair<int, int> uv = HGCalCellUV::cellUVFromXY1(xloc, yloc, placement, type, extend, debug);
0294   int u = uv.first;
0295   int v = uv.second;
0296   if (partial == HGCalTypes::WaferLDTop) {
0297     if (u * HGCalTypes::edgeWaferLDTop[0] + v * HGCalTypes::edgeWaferLDTop[1] == HGCalTypes::edgeWaferLDTop[2] + 1) {
0298       double xloc1 = (placement >= HGCalCell::cellPlacementExtra) ? xloc : -xloc;
0299       int rot = placement % HGCalCell::cellPlacementExtra;
0300       static constexpr std::array<double, 6> fcos = {{1.0, cos60_, -cos60_, -1.0, -cos60_, cos60_}};
0301       static constexpr std::array<double, 6> fsin = {{0.0, sin60_, sin60_, 0.0, -sin60_, -sin60_}};
0302       double xprime = -1 * (xloc1 * fcos[rot] - yloc * fsin[rot]);
0303       double yprime = xloc1 * fsin[rot] + yloc * fcos[rot];
0304       double xcell = -1 * (1.5 * (v - u) + 0.5) * cellX;
0305       double ycell = (v + u - 2 * ncell_[type] + 1) * cellY;
0306       if ((yprime - sqrt3_ * xprime) > (ycell - sqrt3_ * xcell)) {
0307         --u;
0308         if ((v - u) >= ncell_[1])
0309           --v;
0310       } else {
0311         --u;
0312         --v;
0313         v = std::max(v, 0);
0314       }
0315     }
0316   } else if (partial == HGCalTypes::WaferHDBottom) {
0317     if (u * HGCalTypes::edgeWaferHDBottom[0] + v * HGCalTypes::edgeWaferHDBottom[1] ==
0318         HGCalTypes::edgeWaferHDBottom[2] + 1) {
0319       double xloc1 = (placement >= HGCalCell::cellPlacementExtra) ? xloc : -xloc;
0320       int rot = placement % HGCalCell::cellPlacementExtra;
0321       static constexpr std::array<double, 6> fcos = {{1.0, cos60_, -cos60_, -1.0, -cos60_, cos60_}};
0322       static constexpr std::array<double, 6> fsin = {{0.0, sin60_, sin60_, 0.0, -sin60_, -sin60_}};
0323       double xprime = -1 * (xloc1 * fcos[rot] - yloc * fsin[rot]);
0324       double yprime = xloc1 * fsin[rot] + yloc * fcos[rot];
0325       double xcell = -1 * (1.5 * (v - u) + 0.5) * cellX;
0326       double ycell = (v + u - 2 * ncell_[type] + 1) * cellY;
0327       if ((yprime - sqrt3_ * xprime) > (ycell - sqrt3_ * xcell)) {
0328         ++u;
0329         ++v;
0330       } else {
0331         ++u;
0332       }
0333     }
0334   }
0335   if (debug)
0336     edm::LogVerbatim("HGCalGeom") << "cellUVFromXY5: Input " << xloc << ":" << yloc << ":" << extend << " Output " << u
0337                                   << ":" << v;
0338   return std::make_pair(u, v);
0339 }
0340 
0341 std::pair<int32_t, int32_t> HGCalCellUV::cellUVFromXY2(  // for v18
0342     double xloc,
0343     double yloc,
0344     int32_t placement,
0345     int32_t type,
0346     int32_t partial,
0347     bool extend,
0348     bool debug) const {
0349   if (type != 0)
0350     type = 1;
0351   std::pair<int, int> uv = HGCalCellUV::cellUVFromXY1(xloc, yloc, placement, type, extend, debug);
0352   int u = uv.first;
0353   int v = uv.second;
0354   if ((partial != HGCalTypes::WaferFull) && (type == 1)) {
0355     if (u == 1 && v == 8) {
0356       std::array<double, 4> criterion =
0357           HGCalWaferMask::maskCut(HGCalTypes::WaferLDThree, placement, waferSize_, 0.0, false);
0358       if ((criterion[0] * yloc) + (criterion[1] * xloc) < -criterion[2]) {
0359         ++u;
0360         ++v;
0361       }
0362     }
0363     if (u == 15 && v == 15) {
0364       std::array<double, 4> criterion =
0365           HGCalWaferMask::maskCut(HGCalTypes::WaferLDThree, placement, waferSize_, 0.0, false);
0366       if ((criterion[0] * yloc) + (criterion[1] * xloc) < -criterion[2]) {
0367         --u;
0368       }
0369     }
0370     if (u * HGCalTypes::edgeWaferLDTop[0] + v * HGCalTypes::edgeWaferLDTop[1] == HGCalTypes::edgeWaferLDTop[2] + 1) {
0371       std::array<double, 4> criterion =
0372           HGCalWaferMask::maskCut(HGCalTypes::WaferLDTop, placement, waferSize_, 0.0, false);
0373       if ((criterion[0] * yloc) + (criterion[1] * xloc) < -criterion[2]) {
0374         std::pair<double, double> xy1 = hgcalcell_->cellUV2XY1(u, v, placement, 1);
0375         std::pair<double, double> xy2 = hgcalcell_->cellUV2XY1(u - 2, v - 1, placement, 1);
0376         if (((placement >= HGCalCell::cellPlacementExtra) &&
0377              ((((xloc - xy1.first) / (xy2.first - xy1.first)) - ((yloc - xy1.second) / (xy2.second - xy1.second))) >
0378               0.0)) ||
0379             ((placement < HGCalCell::cellPlacementExtra) &&
0380              ((((xloc - xy1.first) / (xy2.first - xy1.first)) - ((yloc - xy1.second) / (xy2.second - xy1.second))) <
0381               0.0))) {
0382           --u;
0383           if ((v - u) >= ncell_[1])
0384             --v;
0385         } else {
0386           --u;
0387           --v;
0388           v = std::max(v, 0);
0389         }
0390       }
0391     }
0392   } else if ((partial != HGCalTypes::WaferFull) && (type == 0)) {
0393     if (u == 10 && v == 0) {
0394       std::array<double, 4> criterion =
0395           HGCalWaferMask::maskCut(HGCalTypes::WaferHDTop, placement, waferSize_, 0.0, false);
0396       if ((criterion[0] * yloc) + (criterion[1] * xloc) < -criterion[2]) {
0397         --u;
0398       }
0399     }
0400     if (u == 10 && v == 21) {
0401       std::array<double, 4> criterion =
0402           HGCalWaferMask::maskCut(HGCalTypes::WaferHDTop, placement, waferSize_, 0.0, false);
0403       if ((criterion[0] * yloc) + (criterion[1] * xloc) < -criterion[2]) {
0404         --u;
0405         --v;
0406       }
0407     }
0408     if (u * HGCalTypes::edgeWaferHDBottom[0] + v * HGCalTypes::edgeWaferHDBottom[1] ==
0409         HGCalTypes::edgeWaferHDBottom[2] + 1) {
0410       std::array<double, 4> criterion =
0411           HGCalWaferMask::maskCut(HGCalTypes::WaferHDBottom, placement, waferSize_, 0.0, false);
0412       if ((criterion[0] * yloc) + (criterion[1] * xloc) < -criterion[2]) {
0413         std::pair<double, double> xy1 = hgcalcell_->cellUV2XY1(u, v, placement, 0);
0414         std::pair<double, double> xy2 = hgcalcell_->cellUV2XY1(u - 2, v - 1, placement, 0);
0415         if (((placement >= HGCalCell::cellPlacementExtra) &&
0416              ((((xloc - xy1.first) / (xy2.first - xy1.first)) - ((yloc - xy1.second) / (xy2.second - xy1.second))) >
0417               0.0)) ||
0418             ((placement < HGCalCell::cellPlacementExtra) &&
0419              ((((xloc - xy1.first) / (xy2.first - xy1.first)) - ((yloc - xy1.second) / (xy2.second - xy1.second))) <
0420               0.0))) {
0421           ++u;
0422           ++v;
0423         } else {
0424           ++u;
0425         }
0426       }
0427     }
0428   }
0429   if (debug)
0430     edm::LogVerbatim("HGCalGeom") << "cellUVFromXY5: Input " << xloc << ":" << yloc << ":" << extend << " Output " << u
0431                                   << ":" << v;
0432   return std::make_pair(u, v);
0433 }