Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-01-09 03:27:01

0001 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0002 
0003 #include "DataFormats/Math/interface/GeantUnits.h"
0004 #include "DataFormats/ForwardDetId/interface/HFNoseDetId.h"
0005 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0006 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0007 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "FWCore/Utilities/interface/Exception.h"
0010 #include "Geometry/HGCalCommonData/interface/HGCalGeomParameters.h"
0011 #include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h"
0012 #include "Geometry/HGCalCommonData/interface/HGCalTypes.h"
0013 #include "Geometry/HGCalCommonData/interface/HGCalWaferMask.h"
0014 #include "Geometry/HGCalCommonData/interface/HGCalWaferType.h"
0015 
0016 #include <algorithm>
0017 #include <bitset>
0018 #include <iterator>
0019 #include <functional>
0020 #include <numeric>
0021 
0022 //#define EDM_ML_DEBUG
0023 using namespace geant_units::operators;
0024 
0025 HGCalDDDConstants::HGCalDDDConstants(const HGCalParameters* hp, const std::string& name)
0026     : hgpar_(hp), sqrt3_(std::sqrt(3.0)), mode_(hgpar_->mode_), fullAndPart_(waferHexagon8File()) {
0027 #ifdef EDM_ML_DEBUG
0028   edm::LogVerbatim("HGCalGeom") << "Mode " << mode_ << " FullAndPart " << fullAndPart_;
0029 #endif
0030   if (waferHexagon6() || waferHexagon8()) {
0031     rmax_ = (HGCalParameters::k_ScaleFromDDD * (hgpar_->waferR_) * std::cos(30._deg));
0032     rmaxT_ = rmax_ + 0.5 * hgpar_->sensorSeparation_;
0033     hexside_ = 2.0 * rmax_ * tan30deg_;
0034     hexsideT_ = 2.0 * rmaxT_ * tan30deg_;
0035     hgcell_ = std::make_unique<HGCalCell>(2.0 * rmaxT_, hgpar_->nCellsFine_, hgpar_->nCellsCoarse_);
0036     hgcellUV_ = std::make_unique<HGCalCellUV>(
0037         2.0 * rmax_, hgpar_->sensorSeparation_, hgpar_->nCellsFine_, hgpar_->nCellsCoarse_);
0038 #ifdef EDM_ML_DEBUG
0039     edm::LogVerbatim("HGCalGeom") << "rmax_ " << rmax_ << ":" << rmaxT_ << ":" << hexside_ << ":" << hexsideT_
0040                                   << " CellSize " << 0.5 * HGCalParameters::k_ScaleFromDDD * hgpar_->cellSize_[0] << ":"
0041                                   << 0.5 * HGCalParameters::k_ScaleFromDDD * hgpar_->cellSize_[1];
0042 #endif
0043   }
0044   if (cassetteMode()) {
0045     hgcassette_.setParameter(hgpar_->cassettes_, hgpar_->cassetteShift_);
0046 #ifdef EDM_ML_DEBUG
0047     edm::LogVerbatim("HGCalGeom") << "Setup HGCalCassette for " << hgpar_->cassettes_ << " cassettes";
0048 #endif
0049   }
0050   // init maps and constants
0051   modHalf_ = 0;
0052   maxWafersPerLayer_ = 0;
0053   for (int simreco = 0; simreco < 2; ++simreco) {
0054     tot_layers_[simreco] = layersInit((bool)simreco);
0055     max_modules_layer_[simreco].resize(tot_layers_[simreco] + 1);
0056     for (unsigned int layer = 1; layer <= tot_layers_[simreco]; ++layer) {
0057       max_modules_layer_[simreco][layer] = modulesInit(layer, (bool)simreco);
0058       if (simreco == 1) {
0059         modHalf_ += max_modules_layer_[simreco][layer];
0060         maxWafersPerLayer_ = std::max(maxWafersPerLayer_, max_modules_layer_[simreco][layer]);
0061 #ifdef EDM_ML_DEBUG
0062         edm::LogVerbatim("HGCalGeom") << "Layer " << layer << " with " << max_modules_layer_[simreco][layer] << ":"
0063                                       << modHalf_ << " modules in RECO";
0064       } else {
0065         edm::LogVerbatim("HGCalGeom") << "Layer " << layer << " with " << max_modules_layer_[simreco][layer]
0066                                       << " modules in SIM";
0067 #endif
0068       }
0069     }
0070 #ifdef EDM_ML_DEBUG
0071     edm::LogVerbatim("HGCalGeom") << "SimReco " << simreco << " with " << tot_layers_[simreco] << " Layers";
0072 #endif
0073   }
0074   tot_wafers_ = wafers();
0075 
0076 #ifdef EDM_ML_DEBUG
0077   edm::LogVerbatim("HGCalGeom") << "HGCalDDDConstants initialized for " << name << " with " << layers(false) << ":"
0078                                 << layers(true) << " layers, " << wafers() << ":" << 2 * modHalf_
0079                                 << " wafers with maximum " << maxWafersPerLayer_ << " per layer and "
0080                                 << "maximum of " << maxCells(false) << ":" << maxCells(true) << " cells";
0081 #endif
0082   if (waferHexagon6() || waferHexagon8()) {
0083     int wminT(9999999), wmaxT(-9999999), kount1(0), kount2(0);
0084     for (unsigned int i = 0; i < getTrFormN(); ++i) {
0085       int lay0 = getTrForm(i).lay;
0086       int wmin(9999999), wmax(-9999999), kount(0);
0087       for (int wafer = 0; wafer < sectors(); ++wafer) {
0088         bool waferIn = waferInLayer(wafer, lay0, true);
0089         if (waferHexagon8()) {
0090           int kndx = HGCalWaferIndex::waferIndex(lay0,
0091                                                  HGCalWaferIndex::waferU(hgpar_->waferCopy_[wafer]),
0092                                                  HGCalWaferIndex::waferV(hgpar_->waferCopy_[wafer]));
0093           waferIn_[kndx] = waferIn;
0094         }
0095         if (waferIn) {
0096           int waferU = ((waferHexagon6()) ? wafer : HGCalWaferIndex::waferU(hgpar_->waferCopy_[wafer]));
0097           if (waferU < wmin)
0098             wmin = waferU;
0099           if (waferU > wmax)
0100             wmax = waferU;
0101           ++kount;
0102         }
0103       }
0104       if (wminT > wmin)
0105         wminT = wmin;
0106       if (wmaxT < wmax)
0107         wmaxT = wmax;
0108       if (kount1 < kount)
0109         kount1 = kount;
0110       kount2 += kount;
0111 #ifdef EDM_ML_DEBUG
0112       int lay1 = getIndex(lay0, true).first;
0113       edm::LogVerbatim("HGCalGeom") << "Index " << i << " Layer " << lay0 << ":" << lay1 << " Wafer " << wmin << ":"
0114                                     << wmax << ":" << kount;
0115 #endif
0116       HGCWaferParam a1{{wmin, wmax, kount}};
0117       waferLayer_[lay0] = a1;
0118     }
0119     waferMax_ = std::array<int, 4>{{wminT, wmaxT, kount1, kount2}};
0120 #ifdef EDM_ML_DEBUG
0121     edm::LogVerbatim("HGCalGeom") << "Overall wafer statistics: " << wminT << ":" << wmaxT << ":" << kount1 << ":"
0122                                   << kount2;
0123 #endif
0124   }
0125 }
0126 
0127 std::pair<int, int> HGCalDDDConstants::assignCell(float x, float y, int lay, int subSec, bool reco) const {
0128   const auto& index = getIndex(lay, reco);
0129   if (index.first < 0)
0130     return std::make_pair(-1, -1);
0131   if (waferHexagon6()) {
0132     float xx = (reco) ? x : HGCalParameters::k_ScaleFromDDD * x;
0133     float yy = (reco) ? y : HGCalParameters::k_ScaleFromDDD * y;
0134 
0135     // First the wafer
0136     int wafer = cellHex(xx, yy, rmax_, hgpar_->waferPosX_, hgpar_->waferPosY_);
0137     if (wafer < 0 || wafer >= static_cast<int>(hgpar_->waferTypeT_.size())) {
0138       edm::LogWarning("HGCalGeom") << "Wafer no. out of bound for " << wafer << ":" << (hgpar_->waferTypeT_).size()
0139                                    << ":" << (hgpar_->waferPosX_).size() << ":" << (hgpar_->waferPosY_).size()
0140                                    << " ***** ERROR *****";
0141       return std::make_pair(-1, -1);
0142     } else {
0143       // Now the cell
0144       xx -= hgpar_->waferPosX_[wafer];
0145       yy -= hgpar_->waferPosY_[wafer];
0146       if (hgpar_->waferTypeT_[wafer] == 1)
0147         return std::make_pair(wafer,
0148                               cellHex(xx,
0149                                       yy,
0150                                       0.5 * HGCalParameters::k_ScaleFromDDD * hgpar_->cellSize_[0],
0151                                       hgpar_->cellFineX_,
0152                                       hgpar_->cellFineY_));
0153       else
0154         return std::make_pair(wafer,
0155                               cellHex(xx,
0156                                       yy,
0157                                       0.5 * HGCalParameters::k_ScaleFromDDD * hgpar_->cellSize_[1],
0158                                       hgpar_->cellCoarseX_,
0159                                       hgpar_->cellCoarseY_));
0160     }
0161   } else {
0162     return std::make_pair(-1, -1);
0163   }
0164 }
0165 
0166 std::array<int, 5> HGCalDDDConstants::assignCellHex(
0167     float x, float y, int zside, int lay, bool reco, bool extend, bool debug) const {
0168   int waferU(0), waferV(0), waferType(-1), cellU(0), cellV(0);
0169   if (waferHexagon8()) {
0170     double xx = (reco) ? HGCalParameters::k_ScaleToDDD * x : x;
0171     double yy = (reco) ? HGCalParameters::k_ScaleToDDD * y : y;
0172     double wt(1.0);
0173 #ifdef EDM_ML_DEBUG
0174     edm::LogVerbatim("HGCalGeom") << "assignCellHex x " << x << ":" << xx << " y " << y << ":" << yy << " Lay " << lay;
0175 #endif
0176     waferFromPosition(xx, yy, zside, lay, waferU, waferV, cellU, cellV, waferType, wt, extend, debug);
0177   }
0178   return std::array<int, 5>{{waferU, waferV, waferType, cellU, cellV}};
0179 }
0180 
0181 std::array<int, 3> HGCalDDDConstants::assignCellTrap(float x, float y, float z, int layer, bool reco) const {
0182   int irad(-1), iphi(-1), type(-1);
0183   const auto& indx = getIndex(layer, reco);
0184   if (indx.first < 0)
0185     return std::array<int, 3>{{irad, iphi, type}};
0186   int zside = (z > 0) ? 1 : -1;
0187   double xx = (reco) ? (zside * x) : (zside * HGCalParameters::k_ScaleFromDDD * x);
0188   double yy = (reco) ? y : HGCalParameters::k_ScaleFromDDD * y;
0189   int ll = layer - hgpar_->firstLayer_;
0190   xx -= hgpar_->xLayerHex_[ll];
0191   yy -= hgpar_->yLayerHex_[ll];
0192   double phi = (((yy == 0.0) && (xx == 0.0)) ? 0. : std::atan2(yy, xx));
0193   if (phi < 0)
0194     phi += (2.0 * M_PI);
0195   if (indx.second != 0)
0196     iphi = (1 + static_cast<int>(phi / indx.second)) % hgpar_->scintCells(layer);
0197   if (iphi == 0)
0198     iphi = hgpar_->scintCells(layer);
0199   if (cassetteMode()) {
0200     int cassette = HGCalTileIndex::tileCassette(iphi, hgpar_->phiOffset_, hgpar_->nphiCassette_, hgpar_->cassettes_);
0201     auto cshift = hgcassette_.getShift(layer, -1, cassette);
0202 #ifdef EDM_ML_DEBUG
0203     std::ostringstream st1;
0204     st1 << "Cassette " << cassette << " Shift " << cshift.first << ":" << cshift.second << " Original " << xx << ":"
0205         << yy;
0206 #endif
0207     xx += (zside * cshift.first);
0208     yy -= cshift.second;
0209 #ifdef EDM_ML_DEBUG
0210     st1 << " Shifted " << xx << ":" << yy;
0211     edm::LogVerbatim("HGCalGeomT") << st1.str();
0212 #endif
0213   }
0214   type = hgpar_->scintType(layer);
0215   double r = std::sqrt(xx * xx + yy * yy);
0216   auto ir = std::lower_bound(hgpar_->radiusLayer_[type].begin(), hgpar_->radiusLayer_[type].end(), r);
0217   irad = static_cast<int>(ir - hgpar_->radiusLayer_[type].begin());
0218   irad = std::clamp(irad, hgpar_->iradMinBH_[indx.first], hgpar_->iradMaxBH_[indx.first]);
0219 #ifdef EDM_ML_DEBUG
0220   edm::LogVerbatim("HGCalGeomT") << "assignCellTrap Input " << x << ":" << y << ":" << z << ":" << layer << ":" << reco
0221                                  << " x|y|r " << xx << ":" << yy << ":" << r << " phi " << phi << ":"
0222                                  << convertRadToDeg(phi) << " o/p " << irad << ":" << iphi << ":" << type;
0223 #endif
0224   if (!tileExist(zside, layer, irad, iphi)) {
0225     if (tileRingEdge(r, layer, irad)) {
0226       if (std::abs(r - hgpar_->radiusLayer_[type][irad - 1]) < tol_) {
0227         --irad;
0228         if (irad <= hgpar_->iradMinBH_[indx.first])
0229           irad = hgpar_->iradMinBH_[indx.first];
0230       } else {
0231         ++irad;
0232         if (irad > hgpar_->iradMaxBH_[indx.first])
0233           irad = hgpar_->iradMaxBH_[indx.first];
0234       }
0235 #ifdef EDM_ML_DEBUG
0236       edm::LogVerbatim("HGCalGeomT") << "assignCellTrap: ring # modified to " << irad << ":"
0237                                      << hgpar_->iradMinBH_[indx.first] << ":" << hgpar_->iradMaxBH_[indx.first];
0238       ;
0239 #endif
0240     } else if (tilePhiEdge(phi, layer, iphi)) {
0241       if (std::abs(phi - hgpar_->scintCellSize(layer) * (iphi - 1)) < tol_) {
0242         --iphi;
0243         if (iphi <= 0)
0244           iphi = 1;
0245       } else {
0246         ++iphi;
0247         if (iphi > hgpar_->scintCells(layer))
0248           iphi = 1;
0249       }
0250 #ifdef EDM_ML_DEBUG
0251       edm::LogVerbatim("HGCalGeomT") << "assignCellTrap: iphi # modified to " << iphi << ":"
0252                                      << hgpar_->scintCells(layer);
0253 #endif
0254     }
0255   }
0256   return std::array<int, 3>{{irad, iphi, type}};
0257 }
0258 
0259 bool HGCalDDDConstants::cassetteShiftSilicon(int zside, int layer, int waferU, int waferV) const {
0260   bool shift(false);
0261   if (cassetteMode()) {
0262     int indx = HGCalWaferIndex::waferIndex(layer, waferU, waferV);
0263     auto ktr = hgpar_->waferInfoMap_.find(indx);
0264     if (ktr != hgpar_->waferInfoMap_.end()) {
0265       auto cshift = hgcassette_.getShift(layer, zside, (ktr->second).cassette);
0266       if ((cshift.first != 0) || (cshift.second != 0))
0267         shift = true;
0268     }
0269   }
0270   return shift;
0271 }
0272 
0273 bool HGCalDDDConstants::cassetteShiftScintillator(int zside, int layer, int iphi) const {
0274   bool shift(false);
0275   if (cassetteMode()) {
0276     auto cshift = hgcassette_.getShift(layer, zside, cassetteTile(iphi));
0277     if ((cshift.first != 0) || (cshift.second != 0))
0278       shift = true;
0279   }
0280   return shift;
0281 }
0282 
0283 std::pair<double, double> HGCalDDDConstants::cellEtaPhiTrap(int type, int irad) const {
0284   double dr(0), df(0);
0285   if (tileTrapezoid()) {
0286     double r = 0.5 * ((hgpar_->radiusLayer_[type][irad - 1] + hgpar_->radiusLayer_[type][irad]));
0287     dr = (hgpar_->radiusLayer_[type][irad] - hgpar_->radiusLayer_[type][irad - 1]);
0288     df = r * hgpar_->cellSize_[type];
0289   }
0290   return std::make_pair(dr, df);
0291 }
0292 
0293 bool HGCalDDDConstants::cellInLayer(int waferU, int waferV, int cellU, int cellV, int lay, int zside, bool reco) const {
0294   const auto& indx = getIndex(lay, true);
0295   if (indx.first >= 0) {
0296     if (cassetteMode()) {
0297       int indx = HGCalWaferIndex::waferIndex(lay, waferU, waferV);
0298       auto ktr = hgpar_->waferInfoMap_.find(indx);
0299       int part = (ktr != hgpar_->waferInfoMap_.end()) ? (ktr->second).part : HGCalTypes::WaferFull;
0300       return HGCalWaferMask::goodCell(cellU, cellV, part);
0301     } else if (mode_ == HGCalGeometryMode::Hexagon8Module) {
0302       int indx = HGCalWaferIndex::waferIndex(lay, waferU, waferV);
0303       auto ktr = hgpar_->waferInfoMap_.find(indx);
0304       int thck(HGCalTypes::WaferFineThin), part(HGCalTypes::WaferFull), rotn(HGCalTypes::WaferOrient0);
0305       if (ktr != hgpar_->waferInfoMap_.end()) {
0306         thck = (ktr->second).type;
0307         part = (ktr->second).part;
0308         rotn = (ktr->second).orient;
0309       }
0310       int ncell = (thck == HGCalTypes::WaferFineThin) ? hgpar_->nCellsFine_ : hgpar_->nCellsCoarse_;
0311       return HGCalWaferMask::goodCell(cellU, cellV, ncell, part, rotn);
0312     } else if (waferHexagon8() || waferHexagon6()) {
0313       const auto& xy =
0314           ((waferHexagon8()) ? locateCell(zside, lay, waferU, waferV, cellU, cellV, reco, true, false, false)
0315                              : locateCell(cellU, lay, waferU, reco));
0316       double rpos = sqrt(xy.first * xy.first + xy.second * xy.second);
0317       return ((rpos >= hgpar_->rMinLayHex_[indx.first]) && (rpos <= hgpar_->rMaxLayHex_[indx.first]));
0318     } else {
0319       return true;
0320     }
0321   } else {
0322     return false;
0323   }
0324 }
0325 
0326 double HGCalDDDConstants::cellThickness(int layer, int waferU, int waferV) const {
0327   double thick(-1);
0328   int type = waferType(layer, waferU, waferV, false);
0329   if (type >= 0) {
0330     if (waferHexagon8()) {
0331       thick = 10000.0 * hgpar_->cellThickness_[type];  // cm to micron
0332     } else if (waferHexagon6()) {
0333       thick = 100.0 * (type + 1);  // type = 1,2,3 for 100,200,300 micron
0334     }
0335   }
0336   return thick;
0337 }
0338 
0339 double HGCalDDDConstants::cellSizeHex(int type) const {
0340   int indx = ((waferHexagon8()) ? ((type >= 1) ? 1 : 0) : ((type == 1) ? 1 : 0));
0341   double cell = (tileTrapezoid() ? 0.5 * hgpar_->cellSize_[indx]
0342                                  : 0.5 * HGCalParameters::k_ScaleFromDDD * hgpar_->cellSize_[indx]);
0343   return cell;
0344 }
0345 
0346 int32_t HGCalDDDConstants::cellType(int type, int cellU, int cellV, int iz, int fwdBack, int orient) const {
0347   int placement = (orient < 0) ? HGCalCell::cellPlacementOld : HGCalCell::cellPlacementIndex(iz, fwdBack, orient);
0348   int ncell = (type == 0) ? hgpar_->nCellsFine_ : hgpar_->nCellsCoarse_;
0349   auto cellType = HGCalCell::cellType(cellU, cellV, ncell, placement);
0350   return cellType.first;
0351 }
0352 
0353 double HGCalDDDConstants::distFromEdgeHex(double x, double y, double z) const {
0354   // Assming the point is within a hexagonal plane of the wafer, calculate
0355   // the shortest distance from the edge
0356   if (z < 0)
0357     x = -x;
0358   double dist(0);
0359   // Input x, y in Geant4 unit and transformed to CMSSW standard
0360   double xx = HGCalParameters::k_ScaleFromDDD * x;
0361   double yy = HGCalParameters::k_ScaleFromDDD * y;
0362   if (waferHexagon8()) {
0363     int ll = layerIndex(getLayer(z, false), false);
0364     xx -= hgpar_->xLayerHex_[ll];
0365     yy -= hgpar_->yLayerHex_[ll];
0366   }
0367   int sizew = static_cast<int>(hgpar_->waferPosX_.size());
0368   int wafer = sizew;
0369   // Transform to the local coordinate frame of the wafer first
0370   for (int k = 0; k < sizew; ++k) {
0371     double dx = std::abs(xx - hgpar_->waferPosX_[k]);
0372     double dy = std::abs(yy - hgpar_->waferPosY_[k]);
0373     if ((dx <= rmax_) && (dy <= hexside_) && ((dy <= 0.5 * hexside_) || (dx * tan30deg_ <= (hexside_ - dy)))) {
0374       wafer = k;
0375       xx -= hgpar_->waferPosX_[k];
0376       yy -= hgpar_->waferPosY_[k];
0377       break;
0378     }
0379   }
0380   // Look at only one quarter (both x,y are positive)
0381   if (wafer < sizew) {
0382     if (std::abs(yy) < 0.5 * hexside_) {
0383       dist = rmax_ - std::abs(xx);
0384     } else {
0385       dist = 0.5 * ((rmax_ - std::abs(xx)) - sqrt3_ * (std::abs(yy) - 0.5 * hexside_));
0386     }
0387   } else {
0388     dist = 0;
0389   }
0390   dist *= HGCalParameters::k_ScaleToDDD;
0391 #ifdef EDM_ML_DEBUG
0392   edm::LogVerbatim("HGCalGeom") << "DistFromEdgeHex: Local " << xx << ":" << yy << " wafer " << wafer << " flag "
0393                                 << (wafer < sizew) << " Distance " << rmax_ << ":" << (rmax_ - std::abs(xx)) << ":"
0394                                 << (std::abs(yy) - 0.5 * hexside_) << ":" << 0.5 * hexside_ << ":" << dist;
0395 #endif
0396   return dist;
0397 }
0398 
0399 double HGCalDDDConstants::distFromEdgeTrap(double x, double y, double z) const {
0400   // Assming the point is within the eta-phi plane of the scintillator tile,
0401   // calculate the shortest distance from the edge
0402   int lay = getLayer(z, false);
0403   double xx = (z < 0) ? -x : x;
0404   int indx = layerIndex(lay, false);
0405   double r = HGCalParameters::k_ScaleFromDDD * std::sqrt(x * x + y * y);
0406   double phi = (r == 0. ? 0. : std::atan2(y, xx));
0407   if (phi < 0)
0408     phi += (2.0 * M_PI);
0409   int type = hgpar_->scintType(lay);
0410   double cell = hgpar_->scintCellSize(lay);
0411   // Compare with the center of the tile find distances along R and also phi
0412   // Take the smaller value
0413   auto ir = std::lower_bound(hgpar_->radiusLayer_[type].begin(), hgpar_->radiusLayer_[type].end(), r);
0414   int irad = static_cast<int>(ir - hgpar_->radiusLayer_[type].begin());
0415   irad = std::clamp(irad, hgpar_->iradMinBH_[indx], hgpar_->iradMaxBH_[indx]);
0416   int iphi = 1 + static_cast<int>(phi / cell);
0417   double dphi = std::max(0.0, (0.5 * cell - std::abs(phi - (iphi - 0.5) * cell)));
0418   double dist = std::min((r - hgpar_->radiusLayer_[type][irad - 1]), (hgpar_->radiusLayer_[type][irad] - r));
0419 #ifdef EDM_ML_DEBUG
0420   edm::LogVerbatim("HGCalGeom") << "DistFromEdgeTrap: Global " << x << ":" << y << ":" << z << " Layer " << lay
0421                                 << " Index " << indx << ":" << type << " xx " << xx << " R " << r << ":" << irad << ":"
0422                                 << hgpar_->radiusLayer_[type][irad - 1] << ":" << hgpar_->radiusLayer_[type][irad]
0423                                 << " Phi " << phi << ":" << iphi << ":" << (iphi - 0.5) * cell << " cell " << cell
0424                                 << " Dphi " << dphi << " Dist " << dist << ":" << r * dphi;
0425 #endif
0426   return HGCalParameters::k_ScaleToDDD * std::min(r * dphi, dist);
0427 }
0428 
0429 int HGCalDDDConstants::getLayer(double z, bool reco) const {
0430   // Get the layer # from the gloabl z coordinate
0431   unsigned int k = 0;
0432   double zz = (reco ? std::abs(z) : HGCalParameters::k_ScaleFromDDD * std::abs(z));
0433   const auto& zLayerHex = hgpar_->zLayerHex_;
0434   auto itr = std::find_if(zLayerHex.begin() + 1, zLayerHex.end(), [&k, &zz, &zLayerHex](double zLayer) {
0435     ++k;
0436     return zz < 0.5 * (zLayerHex[k - 1] + zLayerHex[k]);
0437   });
0438   int lay = (itr == zLayerHex.end()) ? static_cast<int>(zLayerHex.size()) : k;
0439   if (waferHexagon6() && reco) {
0440     int indx = layerIndex(lay, false);
0441     if (indx >= 0)
0442       lay = hgpar_->layerGroupO_[indx];
0443   } else {
0444     lay += (hgpar_->firstLayer_ - 1);
0445   }
0446   return lay;
0447 }
0448 
0449 HGCalParameters::hgtrap HGCalDDDConstants::getModule(unsigned int indx, bool hexType, bool reco) const {
0450   HGCalParameters::hgtrap mytr;
0451   if (hexType) {
0452     if (indx >= hgpar_->waferTypeL_.size())
0453       edm::LogWarning("HGCalGeom") << "Wafer no. out bound for index " << indx << ":" << (hgpar_->waferTypeL_).size()
0454                                    << ":" << (hgpar_->waferPosX_).size() << ":" << (hgpar_->waferPosY_).size()
0455                                    << " ***** ERROR *****";
0456     unsigned int type =
0457         ((indx < hgpar_->waferTypeL_.size()) ? hgpar_->waferTypeL_[indx] - 1 : HGCSiliconDetId::HGCalCoarseThick);
0458     mytr = hgpar_->getModule(type, reco);
0459   } else {
0460     mytr = hgpar_->getModule(indx, reco);
0461   }
0462   return mytr;
0463 }
0464 
0465 std::vector<HGCalParameters::hgtrap> HGCalDDDConstants::getModules() const {
0466   std::vector<HGCalParameters::hgtrap> mytrs;
0467   for (unsigned int k = 0; k < hgpar_->moduleLayR_.size(); ++k)
0468     mytrs.emplace_back(hgpar_->getModule(k, true));
0469   return mytrs;
0470 }
0471 
0472 int HGCalDDDConstants::getPhiBins(int lay) const { return (tileTrapezoid() ? hgpar_->scintCells(lay) : 0); }
0473 
0474 std::pair<double, double> HGCalDDDConstants::getRangeR(int lay, bool reco) const {
0475   int indx = layerIndex(lay, false);
0476   if ((indx >= 0) && (indx < static_cast<int>(hgpar_->rMinLayHex_.size())))
0477     return std::make_pair(hgpar_->rMinLayHex_[indx], hgpar_->rMaxLayHex_[indx]);
0478   else
0479     return std::make_pair(0, -1.);
0480 }
0481 
0482 std::pair<int, int> HGCalDDDConstants::getREtaRange(int lay) const {
0483   int irmin(0), irmax(0);
0484   if (tileTrapezoid()) {
0485     int indx = layerIndex(lay, false);
0486     if ((indx >= 0) && (indx < static_cast<int>(hgpar_->iradMinBH_.size()))) {
0487       irmin = hgpar_->iradMinBH_[indx];
0488       irmax = hgpar_->iradMaxBH_[indx];
0489     }
0490   }
0491   return std::make_pair(irmin, irmax);
0492 }
0493 
0494 std::vector<HGCalParameters::hgtrform> HGCalDDDConstants::getTrForms() const {
0495   std::vector<HGCalParameters::hgtrform> mytrs;
0496   for (unsigned int k = 0; k < hgpar_->trformIndex_.size(); ++k)
0497     mytrs.emplace_back(hgpar_->getTrForm(k));
0498   return mytrs;
0499 }
0500 
0501 int HGCalDDDConstants::getTypeTrap(int layer) const {
0502   // Get the module type for scinitllator
0503   if (tileTrapezoid()) {
0504     return hgpar_->scintType(layer);
0505   } else {
0506     return -1;
0507   }
0508 }
0509 
0510 int HGCalDDDConstants::getTypeHex(int layer, int waferU, int waferV) const {
0511   // Get the module type for a silicon wafer
0512   if (waferHexagon8()) {
0513     auto itr = hgpar_->typesInLayers_.find(HGCalWaferIndex::waferIndex(layer, waferU, waferV));
0514     return ((itr == hgpar_->typesInLayers_.end() ? 2 : hgpar_->waferTypeL_[itr->second]));
0515   } else {
0516     return -1;
0517   }
0518 }
0519 
0520 std::pair<double, double> HGCalDDDConstants::getXY(int layer, double x, double y, bool forwd) const {
0521   int ll = layer - hgpar_->firstLayer_;
0522   double x0(x), y0(y);
0523   if ((!hgpar_->layerType_.empty()) && (ll < static_cast<int>(hgpar_->layerRotV_.size()))) {
0524     if (forwd) {
0525       x0 = x * hgpar_->layerRotV_[ll].first + y * hgpar_->layerRotV_[ll].second;
0526       y0 = y * hgpar_->layerRotV_[ll].first - x * hgpar_->layerRotV_[ll].second;
0527     } else {
0528       x0 = x * hgpar_->layerRotV_[ll].first - y * hgpar_->layerRotV_[ll].second;
0529       y0 = y * hgpar_->layerRotV_[ll].first + x * hgpar_->layerRotV_[ll].second;
0530     }
0531   }
0532 #ifdef EDM_ML_DEBUG
0533   double x1(x0), y1(y0);
0534   if (ll < static_cast<int>(hgpar_->layerRotV_.size())) {
0535     if (forwd) {
0536       x1 = x0 * hgpar_->layerRotV_[ll].first - y0 * hgpar_->layerRotV_[ll].second;
0537       y1 = y0 * hgpar_->layerRotV_[ll].first + x0 * hgpar_->layerRotV_[ll].second;
0538     } else {
0539       x1 = x0 * hgpar_->layerRotV_[ll].first + y0 * hgpar_->layerRotV_[ll].second;
0540       y1 = y0 * hgpar_->layerRotV_[ll].first - x0 * hgpar_->layerRotV_[ll].second;
0541     }
0542   }
0543   edm::LogVerbatim("HGCalGeom") << "HGCalDDDConstants: Layer " << layer << ":" << ll << " mode " << forwd << " x " << x
0544                                 << ":" << x0 << ":" << x1 << " y " << y << ":" << y0 << ":" << y1;
0545 #endif
0546   return std::make_pair(x0, y0);
0547 }
0548 
0549 double HGCalDDDConstants::guardRingOffset(bool reco) const {
0550   return (reco ? hgpar_->guardRingOffset_ : HGCalParameters::k_ScaleToDDD * hgpar_->guardRingOffset_);
0551 }
0552 
0553 bool HGCalDDDConstants::isHalfCell(int waferType, int cell) const {
0554   if (waferType < 1 || cell < 0)
0555     return false;
0556   return waferType == 2 ? hgpar_->cellCoarseHalf_[cell] : hgpar_->cellFineHalf_[cell];
0557 }
0558 
0559 bool HGCalDDDConstants::isValidHex(int lay, int mod, int cell, bool reco) const {
0560   // Check validity for a layer|wafer|cell of pre-TDR version
0561   bool result(false), resultMod(false);
0562   int cellmax(0);
0563   if (waferHexagon6()) {
0564     int32_t copyNumber = hgpar_->waferCopy_[mod];
0565     result = ((lay > 0 && lay <= static_cast<int>(layers(reco))));
0566     if (result) {
0567       const int32_t lay_idx = reco ? (lay - 1) * 3 + 1 : lay;
0568       const auto& the_modules = hgpar_->copiesInLayers_[lay_idx];
0569       auto moditr = the_modules.find(copyNumber);
0570       result = resultMod = (moditr != the_modules.end());
0571 #ifdef EDM_ML_DEBUG
0572       if (!result)
0573         edm::LogVerbatim("HGCalGeom") << "HGCalDDDConstants: Layer " << lay << ":" << lay_idx << " Copy " << copyNumber
0574                                       << ":" << mod << " Flag " << result;
0575 #endif
0576       if (result) {
0577         if (moditr->second >= 0) {
0578           if (mod >= static_cast<int>(hgpar_->waferTypeT_.size()))
0579             edm::LogWarning("HGCalGeom") << "Module no. out of bound for " << mod << " to be compared with "
0580                                          << (hgpar_->waferTypeT_).size() << " ***** ERROR *****";
0581           cellmax = ((hgpar_->waferTypeT_[mod] - 1 == HGCSiliconDetId::HGCalFine)
0582                          ? static_cast<int>(hgpar_->cellFineX_.size())
0583                          : static_cast<int>(hgpar_->cellCoarseX_.size()));
0584           result = (cell >= 0 && cell <= cellmax);
0585         } else {
0586           result = isValidCell(lay_idx, mod, cell);
0587         }
0588       }
0589     }
0590   }
0591 
0592 #ifdef EDM_ML_DEBUG
0593   if (!result)
0594     edm::LogVerbatim("HGCalGeom") << "HGCalDDDConstants: Layer " << lay << ":"
0595                                   << (lay > 0 && (lay <= static_cast<int>(layers(reco)))) << " Module " << mod << ":"
0596                                   << resultMod << " Cell " << cell << ":" << cellmax << ":"
0597                                   << (cell >= 0 && cell <= cellmax) << ":" << maxCells(reco);
0598 #endif
0599   return result;
0600 }
0601 
0602 bool HGCalDDDConstants::isValidHex8(int layer, int modU, int modV, bool fullAndPart) const {
0603   // Check validity for a layer|wafer|cell of post-TDR version
0604   int indx = HGCalWaferIndex::waferIndex(layer, modU, modV);
0605   auto itr = hgpar_->typesInLayers_.find(indx);
0606 #ifdef EDM_ML_DEBUG
0607   edm::LogVerbatim("HGCalGeom") << "HGCalDDDConstants::isValidHex8:WaferType " << layer << ":" << modU << ":" << modV
0608                                 << ":" << indx << " Test " << (itr != hgpar_->typesInLayers_.end());
0609 #endif
0610   if (itr == hgpar_->typesInLayers_.end()) {
0611 #ifdef EDM_ML_DEBUG
0612     edm::LogVerbatim("HGCalGeom") << "HGCalDDDConstants:: Cannot find " << layer << ":" << modU << ":" << modV
0613                                   << " in wadferIndex";
0614 #endif
0615     return false;
0616   }
0617 
0618   if (fullAndPart_) {
0619     auto ktr = hgpar_->waferInfoMap_.find(indx);
0620 #ifdef EDM_ML_DEBUG
0621     edm::LogVerbatim("HGCalGeom") << "HGCalDDDConstants::isValidHex8:WaferInfoMap " << layer << ":" << modU << ":"
0622                                   << modV << ":" << indx << " Test " << (ktr != hgpar_->waferInfoMap_.end());
0623 #endif
0624     if (ktr == hgpar_->waferInfoMap_.end()) {
0625 #ifdef EDM_ML_DEBUG
0626       edm::LogVerbatim("HGCalGeom") << "HGCalDDDConstants:: Cannot find " << layer << ":" << modU << ":" << modV
0627                                     << " in wadferInfoMap";
0628 #endif
0629       return false;
0630     }
0631   } else {
0632     auto jtr = waferIn_.find(indx);
0633 #ifdef EDM_ML_DEBUG
0634     edm::LogVerbatim("HGCalGeom") << "HGCalDDDConstants::isValidHex8:WaferIn " << jtr->first << ":" << jtr->second;
0635 #endif
0636     if (!(jtr->second)) {
0637 #ifdef EDM_ML_DEBUG
0638       edm::LogVerbatim("HGCalGeom") << "HGCalDDDConstants:: Cannot find " << layer << ":" << modU << ":" << modV
0639                                     << " in wadferIn";
0640 #endif
0641       return false;
0642     }
0643   }
0644 
0645   if (fullAndPart || fullAndPart_) {
0646     auto ktr = hgpar_->waferTypes_.find(indx);
0647     if (ktr != hgpar_->waferTypes_.end()) {
0648       if (hgpar_->waferMaskMode_ > 0) {
0649         if (ktr->second.first == HGCalTypes::WaferOut) {
0650 #ifdef EDM_ML_DEBUG
0651           edm::LogVerbatim("HGCalGeom") << "HGCalDDDConstants:: Cannot find " << layer << ":" << modU << ":" << modV
0652                                         << " due to WaferOut";
0653 #endif
0654           return false;
0655         }
0656       } else {
0657         if (ktr->second.first < HGCalTypes::WaferCornerMin) {
0658 #ifdef EDM_ML_DEBUG
0659           edm::LogVerbatim("HGCalGeom") << "HGCalDDDConstants:: Cannot find " << layer << ":" << modU << ":" << modV
0660                                         << " due to WaferCornerMin";
0661 #endif
0662           return false;
0663         }
0664       }
0665     }
0666   }
0667   return true;
0668 }
0669 
0670 bool HGCalDDDConstants::isValidHex8(int layer, int modU, int modV, int cellU, int cellV, bool fullAndPart) const {
0671   // First check validity for a layer|wafer| of post TDR version
0672 #ifdef EDM_ML_DEBUG
0673   edm::LogVerbatim("HGCalGeom") << "HGCalDDDConstants:: layer|wafer " << layer << ":" << modU << ":" << modV << ":"
0674                                 << fullAndPart << " Valid " << isValidHex8(layer, modU, modV, fullAndPart);
0675 #endif
0676   if (!isValidHex8(layer, modU, modV, fullAndPart))
0677     return false;
0678   int indx = HGCalWaferIndex::waferIndex(layer, modU, modV);
0679   auto itr = hgpar_->typesInLayers_.find(indx);
0680   int type = hgpar_->waferTypeL_[itr->second];
0681   int N = ((hgpar_->waferTypeL_[itr->second] == 0) ? hgpar_->nCellsFine_ : hgpar_->nCellsCoarse_);
0682 #ifdef EDM_ML_DEBUG
0683   edm::LogVerbatim("HGCalGeom") << "HGCalDDDConstants::isValidHex8:Cell " << cellU << ":" << cellV << ":" << N
0684                                 << " Tests " << (cellU >= 0) << ":" << (cellU < 2 * N) << ":" << (cellV >= 0) << ":"
0685                                 << (cellV < 2 * N) << ":" << ((cellV - cellU) < N) << ":" << ((cellU - cellV) <= N);
0686 #endif
0687   if ((cellU < 0) || (cellU >= 2 * N) || (cellV < 0) || (cellV >= 2 * N)) {
0688 #ifdef EDM_ML_DEBUG
0689     edm::LogVerbatim("HGCalGeom") << "HGCalDDDConstants:: Cannot statisfy Cell 1 condition " << cellU << ":" << cellV
0690                                   << ":" << N;
0691 #endif
0692     return false;
0693   }
0694   if (((cellV - cellU) >= N) || ((cellU - cellV) > N)) {
0695 #ifdef EDM_ML_DEBUG
0696     edm::LogVerbatim("HGCalGeom") << "HGCalDDDConstants:: Cannot statisfy Cell 2 condition " << cellU << ":" << cellV
0697                                   << ":" << N;
0698 #endif
0699     return false;
0700   }
0701   return isValidCell8(layer, modU, modV, cellU, cellV, type);
0702 }
0703 
0704 bool HGCalDDDConstants::isValidTrap(int zside, int layer, int irad, int iphi) const {
0705   // Check validity for a layer|eta|phi of scintillator
0706   const auto& indx = getIndex(layer, true);
0707   if (indx.first < 0)
0708     return false;
0709   bool ok = ((irad >= hgpar_->iradMinBH_[indx.first]) && (irad <= (hgpar_->iradMaxBH_[indx.first] + 1)) && (iphi > 0) &&
0710              (iphi <= hgpar_->scintCells(layer)));
0711   bool valid = ((ok && trapezoidFile()) ? tileExist(zside, layer, irad, iphi) : ok);
0712 #ifdef EDM_ML_DEBUG
0713   bool tileEx = trapezoidFile() ? tileExist(zside, layer, irad, iphi) : true;
0714   edm::LogVerbatim("HGCalGeomT") << "HGCalDDDConstants::isValidityTrap: Input " << zside << ":" << layer << ":" << irad
0715                                  << ":" << iphi << " Range on Ring " << hgpar_->iradMinBH_[indx.first] << ":"
0716                                  << (hgpar_->iradMaxBH_[indx.first] + 1)
0717                                  << " Range on phi 0:" << hgpar_->scintCells(layer) << " tileExist " << tileEx
0718                                  << " Valid " << ok << ":" << valid;
0719 #endif
0720   return valid;
0721 }
0722 
0723 int HGCalDDDConstants::lastLayer(bool reco) const {
0724   return (hgpar_->firstLayer_ + tot_layers_[static_cast<int>(reco)] - 1);
0725 }
0726 
0727 unsigned int HGCalDDDConstants::layers(bool reco) const { return tot_layers_[static_cast<int>(reco)]; }
0728 
0729 int HGCalDDDConstants::layerIndex(int lay, bool reco) const {
0730   int ll = lay - hgpar_->firstLayer_;
0731   if (ll < 0 || ll >= static_cast<int>(hgpar_->layerIndex_.size()))
0732     return -1;
0733   if (waferHexagon6()) {
0734     if (reco && ll >= static_cast<int>(hgpar_->depthIndex_.size()))
0735       return -1;
0736     return (reco ? hgpar_->depthLayerF_[ll] : hgpar_->layerIndex_[ll]);
0737   } else {
0738     return (hgpar_->layerIndex_[ll]);
0739   }
0740 }
0741 
0742 unsigned int HGCalDDDConstants::layersInit(bool reco) const {
0743   return (reco ? hgpar_->depthIndex_.size() : hgpar_->layerIndex_.size());
0744 }
0745 
0746 std::pair<float, float> HGCalDDDConstants::localToGlobal8(
0747     int zside, int lay, int waferU, int waferV, double localX, double localY, bool reco, bool debug) const {
0748   double x(localX), y(localY);
0749   bool rotx =
0750       ((!hgpar_->layerType_.empty()) && (hgpar_->layerType_[lay - hgpar_->firstLayer_] == HGCalTypes::WaferCenterR));
0751   if (debug)
0752     edm::LogVerbatim("HGCalGeom") << "LocalToGlobal8 " << lay << ":" << (lay - hgpar_->firstLayer_) << ":" << rotx
0753                                   << " Local (" << x << ":" << y << ") Reco " << reco;
0754   if (!reco) {
0755     x *= HGCalParameters::k_ScaleToDDD;
0756     y *= HGCalParameters::k_ScaleToDDD;
0757   }
0758   const auto& xy = waferPositionNoRot(lay, waferU, waferV, reco, debug);
0759   x += xy.first;
0760   y += xy.second;
0761   int indx = HGCalWaferIndex::waferIndex(lay, waferU, waferV);
0762   auto ktr = hgpar_->waferInfoMap_.find(indx);
0763   if (cassetteMode() && (ktr != hgpar_->waferInfoMap_.end())) {
0764     auto cshift = hgcassette_.getShift(lay, -1, (ktr->second).cassette);
0765     std::ostringstream st1;
0766     if (debug)
0767       st1 << "Cassette " << (ktr->second).cassette << " Shift " << cshift.first << ":" << cshift.second << " Original "
0768           << x << ":" << y;
0769     if (!reco) {
0770       x -= ((HGCalParameters::k_ScaleToDDD)*zside * cshift.first);
0771       y += ((HGCalParameters::k_ScaleToDDD)*cshift.second);
0772     } else {
0773       x -= (zside * cshift.first);
0774       y += cshift.second;
0775     }
0776     if (debug) {
0777       st1 << " Final " << x << ":" << y;
0778       edm::LogVerbatim("HGCalGeom") << st1.str();
0779     }
0780   }
0781   if (debug)
0782     edm::LogVerbatim("HGCalGeom") << "With wafer " << x << ":" << y << " by adding " << xy.first << ":" << xy.second;
0783   return (rotx ? getXY(lay, x, y, false) : std::make_pair(x, y));
0784 }
0785 
0786 std::pair<float, float> HGCalDDDConstants::locateCell(int cell, int lay, int type, bool reco) const {
0787   // type refers to wafer # for hexagon cell
0788   float x(999999.), y(999999.);
0789   const auto& index = getIndex(lay, reco);
0790   int i = index.first;
0791   if (i < 0)
0792     return std::make_pair(x, y);
0793   if (waferHexagon6()) {
0794     x = hgpar_->waferPosX_[type];
0795     y = hgpar_->waferPosY_[type];
0796 #ifdef EDM_ML_DEBUG
0797     float x0(x), y0(y);
0798 #endif
0799     if (hgpar_->waferTypeT_[type] - 1 == HGCSiliconDetId::HGCalFine) {
0800       x += hgpar_->cellFineX_[cell];
0801       y += hgpar_->cellFineY_[cell];
0802     } else {
0803       x += hgpar_->cellCoarseX_[cell];
0804       y += hgpar_->cellCoarseY_[cell];
0805     }
0806 #ifdef EDM_ML_DEBUG
0807     edm::LogVerbatim("HGCalGeom") << "LocateCell (Wafer) " << x0 << ":" << y0 << " Final " << x << ":" << y;
0808 #endif
0809     if (!reco) {
0810       x *= HGCalParameters::k_ScaleToDDD;
0811       y *= HGCalParameters::k_ScaleToDDD;
0812     }
0813   }
0814   return std::make_pair(x, y);
0815 }
0816 
0817 std::pair<float, float> HGCalDDDConstants::locateCell(
0818     int zside, int lay, int waferU, int waferV, int cellU, int cellV, bool reco, bool all, bool norot, bool debug)
0819     const {
0820   double x(0), y(0);
0821   int indx = HGCalWaferIndex::waferIndex(lay, waferU, waferV);
0822   auto itr = hgpar_->typesInLayers_.find(indx);
0823   int type = ((itr == hgpar_->typesInLayers_.end()) ? 2 : hgpar_->waferTypeL_[itr->second]);
0824   int layertype = layerType(lay);
0825   bool rotx = (norot) ? false : (layertype == HGCalTypes::WaferCenterR);
0826   if (debug) {
0827     edm::LogVerbatim("HGCalGeom") << "LocateCell " << lay << ":" << (lay - hgpar_->firstLayer_) << ":" << layertype
0828                                   << ":" << rotx << ":" << waferU << ":" << waferV << ":" << indx << ":"
0829                                   << (itr == hgpar_->typesInLayers_.end()) << ":" << type << " Flags " << reco << ":"
0830                                   << all;
0831   }
0832   auto ktr = hgpar_->waferInfoMap_.end();
0833   int place(HGCalCell::cellPlacementOld);
0834   if (waferHexagon8File()) {
0835     if (cassetteMode()) {
0836       ktr = hgpar_->waferInfoMap_.find(indx);
0837       if (ktr != hgpar_->waferInfoMap_.end())
0838         place = HGCalCell::cellPlacementIndex(1, HGCalTypes::layerFrontBack(layertype), (ktr->second).orient);
0839     }
0840     auto xy = hgcell_->cellUV2XY2(cellU, cellV, place, type);
0841     x = xy.first;
0842     y = xy.second;
0843     if (debug)
0844       edm::LogVerbatim("HGCalGeom") << "Type " << type << " Place " << place << " Cell " << cellU << ":" << cellV
0845                                     << " Position " << x << ":" << y;
0846   } else {
0847     int kndx = cellV * 100 + cellU;
0848     if (type == 0) {
0849       auto jtr = hgpar_->cellFineIndex_.find(kndx);
0850       if (jtr != hgpar_->cellFineIndex_.end()) {
0851         x = hgpar_->cellFineX_[jtr->second];
0852         y = hgpar_->cellFineY_[jtr->second];
0853       }
0854       if (debug)
0855         edm::LogVerbatim("HGCalGeom") << "Fine " << cellU << ":" << cellV << ":" << kndx << ":" << x << ":" << y << ":"
0856                                       << (jtr != hgpar_->cellFineIndex_.end());
0857     } else {
0858       auto jtr = hgpar_->cellCoarseIndex_.find(kndx);
0859       if (jtr != hgpar_->cellCoarseIndex_.end()) {
0860         x = hgpar_->cellCoarseX_[jtr->second];
0861         y = hgpar_->cellCoarseY_[jtr->second];
0862       }
0863       if (debug)
0864         edm::LogVerbatim("HGCalGeom") << "Coarse " << cellU << ":" << cellV << ":" << kndx << ":" << x << ":" << y
0865                                       << ":" << (jtr != hgpar_->cellCoarseIndex_.end());
0866     }
0867   }
0868   if (!reco) {
0869     x *= HGCalParameters::k_ScaleToDDD;
0870     y *= HGCalParameters::k_ScaleToDDD;
0871   }
0872   if (all) {
0873     const auto& xy = waferPositionNoRot(lay, waferU, waferV, reco, debug);
0874     x += xy.first;
0875     y += xy.second;
0876     if (cassetteMode() && (ktr != hgpar_->waferInfoMap_.end())) {
0877       auto cshift = hgcassette_.getShift(lay, -1, (ktr->second).cassette);
0878       std::ostringstream st1;
0879       if (debug)
0880         st1 << "Cassette " << (ktr->second).cassette << " Shift " << cshift.first << ":" << cshift.second
0881             << " Original " << x << ":" << y;
0882       if (!reco) {
0883         x -= ((HGCalParameters::k_ScaleToDDD)*cshift.first);
0884         y += ((HGCalParameters::k_ScaleToDDD)*cshift.second);
0885       } else {
0886         x -= cshift.first;
0887         y += cshift.second;
0888       }
0889       if (debug) {
0890         st1 << " Final " << x << ":" << y;
0891         edm::LogVerbatim("HGCalGeom") << st1.str();
0892       }
0893     }
0894     if (debug)
0895       edm::LogVerbatim("HGCalGeom") << "With wafer " << x << ":" << y << " by adding " << xy.first << ":" << xy.second;
0896   }
0897   return (rotx ? getXY(lay, x, y, false) : std::make_pair(x, y));
0898 }
0899 
0900 std::pair<float, float> HGCalDDDConstants::locateCell(const HGCSiliconDetId& id, bool debug) const {
0901   return locateCell(id.zside(), id.layer(), id.waferU(), id.waferV(), id.cellU(), id.cellV(), true, true, false, debug);
0902 }
0903 
0904 std::pair<float, float> HGCalDDDConstants::locateCell(const HGCScintillatorDetId& id, bool debug) const {
0905   return locateCellTrap(id.zside(), id.layer(), id.iradius(), id.iphi(), true, debug);
0906 }
0907 
0908 std::pair<float, float> HGCalDDDConstants::locateCellHex(int cell, int wafer, bool reco) const {
0909   float x(0), y(0);
0910   if (hgpar_->waferTypeT_[wafer] - 1 == HGCSiliconDetId::HGCalFine) {
0911     x = hgpar_->cellFineX_[cell];
0912     y = hgpar_->cellFineY_[cell];
0913   } else {
0914     x = hgpar_->cellCoarseX_[cell];
0915     y = hgpar_->cellCoarseY_[cell];
0916   }
0917   if (!reco) {
0918     x *= HGCalParameters::k_ScaleToDDD;
0919     y *= HGCalParameters::k_ScaleToDDD;
0920   }
0921   return std::make_pair(x, y);
0922 }
0923 
0924 std::pair<float, float> HGCalDDDConstants::locateCellTrap(
0925     int zside, int lay, int irad, int iphi, bool reco, bool debug) const {
0926   float x(0), y(0);
0927   const auto& indx = getIndex(lay, reco);
0928   if (indx.first >= 0) {
0929     int ir = std::abs(irad);
0930     int type = hgpar_->scintType(lay);
0931     double phi = (iphi - 0.5) * indx.second;
0932     double z = hgpar_->zLayerHex_[indx.first];
0933     double r = 0.5 * (hgpar_->radiusLayer_[type][ir - 1] + hgpar_->radiusLayer_[type][ir]);
0934     std::pair<double, double> range = rangeR(z, true);
0935     if (debug)
0936       edm::LogVerbatim("HGCalGeom") << "locateCellTrap:: Input " << lay << ":" << irad << ":" << iphi << ":" << reco
0937                                     << " IR " << ir << ":" << hgpar_->iradMinBH_[indx.first] << ":"
0938                                     << hgpar_->iradMaxBH_[indx.first] << " Type " << type << " Z " << indx.first << ":"
0939                                     << z << " phi " << phi << ":" << convertRadToDeg(phi) << " R " << r << ":"
0940                                     << range.first << ":" << range.second;
0941     if (!trapezoidFile())
0942       r = std::max(range.first, std::min(r, range.second));
0943     x = r * std::cos(phi);
0944     y = r * std::sin(phi);
0945     int ll = lay - hgpar_->firstLayer_;
0946     x += hgpar_->xLayerHex_[ll];
0947     y += hgpar_->yLayerHex_[ll];
0948     if (irad < 0)
0949       x = -x;
0950     if (cassetteMode()) {
0951       int cassette = HGCalTileIndex::tileCassette(iphi, hgpar_->phiOffset_, hgpar_->nphiCassette_, hgpar_->cassettes_);
0952       auto cshift = hgcassette_.getShift(lay, -1, cassette);
0953       std::ostringstream st1;
0954       if (debug)
0955         st1 << "Cassette " << cassette << " Shift " << cshift.first << ":" << cshift.second << " Original " << x << ":"
0956             << y;
0957       x -= cshift.first;
0958       y += cshift.second;
0959       if (debug) {
0960         st1 << " Final " << x << ":" << y;
0961         edm::LogVerbatim("HGCalGeom") << st1.str();
0962       }
0963     }
0964   }
0965   if (!reco) {
0966     x *= HGCalParameters::k_ScaleToDDD;
0967     y *= HGCalParameters::k_ScaleToDDD;
0968   }
0969   return std::make_pair(x, y);
0970 }
0971 
0972 bool HGCalDDDConstants::maskCell(const DetId& detId, int corners) const {
0973   bool mask(false);
0974   if (corners > 2 && corners <= static_cast<int>(HGCalParameters::k_CornerSize)) {
0975     if (waferHexagon8()) {
0976       int N(0), layer(0), waferU(0), waferV(0), u(0), v(0);
0977       if (detId.det() == DetId::Forward) {
0978         HFNoseDetId id(detId);
0979         N = getUVMax(id.type());
0980         layer = id.layer();
0981         waferU = id.waferU();
0982         waferV = id.waferV();
0983         u = id.cellU();
0984         v = id.cellV();
0985       } else {
0986         HGCSiliconDetId id(detId);
0987         N = getUVMax(id.type());
0988         layer = id.layer();
0989         waferU = id.waferU();
0990         waferV = id.waferV();
0991         u = id.cellU();
0992         v = id.cellV();
0993       }
0994       int wl = HGCalWaferIndex::waferIndex(layer, waferU, waferV);
0995       auto itr = hgpar_->waferTypes_.find(wl);
0996       auto ktr = hgpar_->waferInfoMap_.find(wl);
0997 #ifdef EDM_ML_DEBUG
0998       edm::LogVerbatim("HGCalGeom") << "MaskCell: Layer " << layer << " Wafer " << waferU << ":" << waferV << " Index "
0999                                     << wl << ":" << (itr != hgpar_->waferTypes_.end()) << ":"
1000                                     << (ktr != hgpar_->waferInfoMap_.end());
1001 #endif
1002       if (cassetteMode()) {
1003         int part = (ktr != hgpar_->waferInfoMap_.end()) ? (ktr->second).part : HGCalTypes::WaferFull;
1004         mask = !(HGCalWaferMask::goodCell(u, v, part));
1005       } else if (itr != hgpar_->waferTypes_.end()) {
1006         if ((itr->second).second <= HGCalTypes::k_OffsetRotation)
1007           mask = HGCalWaferMask::maskCell(u, v, N, (itr->second).first, (itr->second).second, corners);
1008         else
1009           mask = !(HGCalWaferMask::goodCell(
1010               u, v, N, (itr->second).first, ((itr->second).second - HGCalTypes::k_OffsetRotation)));
1011       }
1012     }
1013   }
1014   return mask;
1015 }
1016 
1017 int HGCalDDDConstants::maxCells(bool reco) const {
1018   int cells(0);
1019   for (unsigned int i = 0; i < layers(reco); ++i) {
1020     int lay = reco ? hgpar_->depth_[i] : hgpar_->layer_[i];
1021     if (cells < maxCells(lay, reco))
1022       cells = maxCells(lay, reco);
1023   }
1024   return cells;
1025 }
1026 
1027 int HGCalDDDConstants::maxCells(int lay, bool reco) const {
1028   const auto& index = getIndex(lay, reco);
1029   if (index.first < 0)
1030     return 0;
1031   if (waferHexagon6()) {
1032     unsigned int cells(0);
1033     for (unsigned int k = 0; k < hgpar_->waferTypeT_.size(); ++k) {
1034       if (waferInLayerTest(k, index.first, hgpar_->defineFull_)) {
1035         unsigned int cell = (hgpar_->waferTypeT_[k] - 1 == HGCSiliconDetId::HGCalFine) ? (hgpar_->cellFineX_.size())
1036                                                                                        : (hgpar_->cellCoarseX_.size());
1037         if (cell > cells)
1038           cells = cell;
1039       }
1040     }
1041     return static_cast<int>(cells);
1042   } else if (waferHexagon8()) {
1043     int cells(0);
1044     for (unsigned int k = 0; k < hgpar_->waferCopy_.size(); ++k) {
1045       if (waferInLayerTest(k, index.first, hgpar_->defineFull_)) {
1046         auto itr = hgpar_->typesInLayers_.find(HGCalWaferIndex::waferIndex(
1047             lay, HGCalWaferIndex::waferU(hgpar_->waferCopy_[k]), HGCalWaferIndex::waferV(hgpar_->waferCopy_[k])));
1048         int type = ((itr == hgpar_->typesInLayers_.end()) ? HGCSiliconDetId::HGCalCoarseThick
1049                                                           : hgpar_->waferTypeL_[itr->second]);
1050         int N = (type == HGCSiliconDetId::HGCalFine) ? hgpar_->nCellsFine_ : hgpar_->nCellsCoarse_;
1051         cells = std::max(cells, 3 * N * N);
1052       }
1053     }
1054     return cells;
1055   } else if (tileTrapezoid()) {
1056     return hgpar_->scintCells(index.first + hgpar_->firstLayer_);
1057   } else {
1058     return 0;
1059   }
1060 }
1061 
1062 int HGCalDDDConstants::maxRows(int lay, bool reco) const {
1063   int kymax(0);
1064   const auto& index = getIndex(lay, reco);
1065   int i = index.first;
1066   if (i < 0)
1067     return kymax;
1068   if (waferHexagon6()) {
1069     for (unsigned int k = 0; k < hgpar_->waferCopy_.size(); ++k) {
1070       if (waferInLayerTest(k, i, hgpar_->defineFull_)) {
1071         int ky = ((hgpar_->waferCopy_[k]) / 100) % 100;
1072         if (ky > kymax)
1073           kymax = ky;
1074       }
1075     }
1076   } else if (waferHexagon8()) {
1077     kymax = 1 + 2 * hgpar_->waferUVMaxLayer_[index.first];
1078   }
1079   return kymax;
1080 }
1081 
1082 int HGCalDDDConstants::modifyUV(int uv, int type1, int type2) const {
1083   // Modify u/v for transition of type1 to type2
1084   return (((type1 == type2) || (type1 * type2 != 0)) ? uv : ((type1 == 0) ? (2 * uv + 1) / 3 : (3 * uv) / 2));
1085 }
1086 
1087 int HGCalDDDConstants::modules(int lay, bool reco) const {
1088   if (getIndex(lay, reco).first < 0)
1089     return 0;
1090   else
1091     return max_modules_layer_[static_cast<int>(reco)][lay];
1092 }
1093 
1094 int HGCalDDDConstants::modulesInit(int lay, bool reco) const {
1095   int nmod(0);
1096   const auto& index = getIndex(lay, reco);
1097   if (index.first < 0)
1098     return nmod;
1099   if (!tileTrapezoid()) {
1100     for (unsigned int k = 0; k < hgpar_->waferPosX_.size(); ++k) {
1101       if (waferInLayerTest(k, index.first, hgpar_->defineFull_))
1102         ++nmod;
1103     }
1104   } else {
1105     nmod = 1 + hgpar_->lastModule_[index.first] - hgpar_->firstModule_[index.first];
1106   }
1107   return nmod;
1108 }
1109 
1110 double HGCalDDDConstants::mouseBite(bool reco) const {
1111   return (reco ? hgpar_->mouseBite_ : HGCalParameters::k_ScaleToDDD * hgpar_->mouseBite_);
1112 }
1113 
1114 int HGCalDDDConstants::numberCells(bool reco) const {
1115   int cells = (tileTrapezoid() && ((hgpar_->waferMaskMode_ == HGCalGeomParameters::scintillatorFile) ||
1116                                    (hgpar_->waferMaskMode_ == HGCalGeomParameters::scintillatorCassette)))
1117                   ? tileCount(0, -1)
1118                   : 0;
1119   if (cells == 0) {
1120     unsigned int nlayer = (reco) ? hgpar_->depth_.size() : hgpar_->layer_.size();
1121     for (unsigned k = 0; k < nlayer; ++k) {
1122       std::vector<int> ncells = numberCells(((reco) ? hgpar_->depth_[k] : hgpar_->layer_[k]), reco);
1123       cells = std::accumulate(ncells.begin(), ncells.end(), cells);
1124     }
1125   }
1126   return cells;
1127 }
1128 
1129 std::vector<int> HGCalDDDConstants::numberCells(int lay, bool reco) const {
1130   const auto& index = getIndex(lay, reco);
1131   int i = index.first;
1132   std::vector<int> ncell;
1133   if (i >= 0) {
1134     if (waferHexagon6()) {
1135       for (unsigned int k = 0; k < hgpar_->waferTypeT_.size(); ++k) {
1136         if (waferInLayerTest(k, i, hgpar_->defineFull_)) {
1137           unsigned int cell = (hgpar_->waferTypeT_[k] - 1 == HGCSiliconDetId::HGCalFine)
1138                                   ? (hgpar_->cellFineX_.size())
1139                                   : (hgpar_->cellCoarseX_.size());
1140           ncell.emplace_back(static_cast<int>(cell));
1141         }
1142       }
1143     } else if (tileTrapezoid()) {
1144       int nphi = hgpar_->scintCells(lay);
1145       for (int k = hgpar_->firstModule_[i]; k <= hgpar_->lastModule_[i]; ++k)
1146         ncell.emplace_back(nphi);
1147     } else {
1148       for (unsigned int k = 0; k < hgpar_->waferCopy_.size(); ++k) {
1149         if (waferInLayerTest(k, index.first, hgpar_->defineFull_)) {
1150           int cell = numberCellsHexagon(lay,
1151                                         HGCalWaferIndex::waferU(hgpar_->waferCopy_[k]),
1152                                         HGCalWaferIndex::waferV(hgpar_->waferCopy_[k]),
1153                                         true);
1154           ncell.emplace_back(cell);
1155         }
1156       }
1157     }
1158   }
1159   return ncell;
1160 }
1161 
1162 int HGCalDDDConstants::numberCellsHexagon(int wafer) const {
1163   if (wafer >= 0 && wafer < static_cast<int>(hgpar_->waferTypeT_.size())) {
1164     if (hgpar_->waferTypeT_[wafer] - 1 == HGCSiliconDetId::HGCalFine)
1165       return static_cast<int>(hgpar_->cellFineX_.size());
1166     else
1167       return static_cast<int>(hgpar_->cellCoarseX_.size());
1168   } else {
1169     return 0;
1170   }
1171 }
1172 
1173 int HGCalDDDConstants::numberCellsHexagon(int lay, int waferU, int waferV, bool flag) const {
1174   auto itr = hgpar_->typesInLayers_.find(HGCalWaferIndex::waferIndex(lay, waferU, waferV));
1175   int type =
1176       ((itr == hgpar_->typesInLayers_.end()) ? HGCSiliconDetId::HGCalCoarseThick : hgpar_->waferTypeL_[itr->second]);
1177   int N = (type == 0) ? hgpar_->nCellsFine_ : hgpar_->nCellsCoarse_;
1178   if (flag)
1179     return (3 * N * N);
1180   else
1181     return N;
1182 }
1183 
1184 std::pair<double, double> HGCalDDDConstants::rangeR(double z, bool reco) const {
1185   double rmin(0), rmax(0), zz(0);
1186   if (hgpar_->detectorType_ > 0) {
1187     zz = (reco ? std::abs(z) : HGCalParameters::k_ScaleFromDDD * std::abs(z));
1188     if (hgpar_->detectorType_ <= 2) {
1189       rmin = HGCalGeomTools::radius(zz, hgpar_->zFrontMin_, hgpar_->rMinFront_, hgpar_->slopeMin_);
1190     } else {
1191       rmin = HGCalGeomTools::radius(
1192           zz, hgpar_->firstLayer_, hgpar_->firstMixedLayer_, hgpar_->zLayerHex_, hgpar_->radiusMixBoundary_);
1193     }
1194     if ((hgpar_->detectorType_ == 2) && (zz >= hgpar_->zLayerHex_[hgpar_->firstMixedLayer_ - 1])) {
1195       rmax = HGCalGeomTools::radius(
1196           zz, hgpar_->firstLayer_, hgpar_->firstMixedLayer_, hgpar_->zLayerHex_, hgpar_->radiusMixBoundary_);
1197     } else {
1198       rmax = HGCalGeomTools::radius(zz, hgpar_->zFrontTop_, hgpar_->rMaxFront_, hgpar_->slopeTop_);
1199     }
1200   }
1201   if (!reco) {
1202     rmin *= HGCalParameters::k_ScaleToDDD;
1203     rmax *= HGCalParameters::k_ScaleToDDD;
1204   }
1205 #ifdef EDM_ML_DEBUG
1206   edm::LogVerbatim("HGCalGeom") << "HGCalDDDConstants:rangeR: " << z << ":" << zz << " R " << rmin << ":" << rmax;
1207 #endif
1208   return std::make_pair(rmin, rmax);
1209 }
1210 
1211 std::pair<double, double> HGCalDDDConstants::rangeRLayer(int lay, bool reco) const {
1212   double rmin(0), rmax(0);
1213   const auto& index = getIndex(lay, reco);
1214   if (index.first >= 0 && index.first < static_cast<int>(hgpar_->rMinLayHex_.size())) {
1215     rmin = hgpar_->rMinLayHex_[index.first];
1216     rmax = hgpar_->rMaxLayHex_[index.first];
1217   }
1218   if (!reco) {
1219     rmin *= HGCalParameters::k_ScaleToDDD;
1220     rmax *= HGCalParameters::k_ScaleToDDD;
1221   }
1222 #ifdef EDM_ML_DEBUG
1223   edm::LogVerbatim("HGCalGeom") << "HGCalDDDConstants:rangeR: " << lay << ":" << index.first << " R " << rmin << ":"
1224                                 << rmax;
1225 #endif
1226   return std::make_pair(rmin, rmax);
1227 }
1228 
1229 std::pair<double, double> HGCalDDDConstants::rangeZ(bool reco) const {
1230   double zmin = (hgpar_->zLayerHex_[0] - hgpar_->waferThick_);
1231   double zmax = (hgpar_->zLayerHex_[hgpar_->zLayerHex_.size() - 1] + hgpar_->waferThick_);
1232 #ifdef EDM_ML_DEBUG
1233   edm::LogVerbatim("HGCalGeom") << "HGCalDDDConstants:rangeZ: " << zmin << ":" << zmax << ":" << hgpar_->waferThick_;
1234 #endif
1235   if (!reco) {
1236     zmin *= HGCalParameters::k_ScaleToDDD;
1237     zmax *= HGCalParameters::k_ScaleToDDD;
1238   }
1239   return std::make_pair(zmin, zmax);
1240 }
1241 
1242 std::pair<int, int> HGCalDDDConstants::rowColumnWafer(int wafer) const {
1243   int row(0), col(0);
1244   if (wafer < static_cast<int>(hgpar_->waferCopy_.size())) {
1245     int copy = hgpar_->waferCopy_[wafer];
1246     col = HGCalTypes::getUnpackedU(copy);
1247     row = HGCalTypes::getUnpackedV(copy);
1248     ;
1249   }
1250   return std::make_pair(row, col);
1251 }
1252 
1253 double HGCalDDDConstants::sensorSizeOffset(bool reco) const {
1254   return (reco ? hgpar_->sensorSizeOffset_ : HGCalParameters::k_ScaleToDDD * hgpar_->sensorSizeOffset_);
1255 }
1256 
1257 std::pair<int, int> HGCalDDDConstants::simToReco(int cell, int lay, int mod, bool half) const {
1258   if (!waferHexagon6()) {
1259     return std::make_pair(cell, lay);
1260   } else {
1261     const auto& index = getIndex(lay, false);
1262     int i = index.first;
1263     if (i < 0) {
1264       edm::LogWarning("HGCalGeom") << "Wrong Layer # " << lay << " not in the list ***** ERROR *****";
1265       return std::make_pair(-1, -1);
1266     }
1267     if (mod >= static_cast<int>(hgpar_->waferTypeL_.size())) {
1268       edm::LogWarning("HGCalGeom") << "Invalid Wafer # " << mod << "should be < " << (hgpar_->waferTypeL_).size()
1269                                    << " ***** ERROR *****";
1270       return std::make_pair(-1, -1);
1271     }
1272     int depth(-1);
1273     int kx = cell;
1274     int type = hgpar_->waferTypeL_[mod];
1275     if (type == 1) {
1276       depth = hgpar_->layerGroup_[i];
1277     } else if (type == 2) {
1278       depth = hgpar_->layerGroupM_[i];
1279     } else {
1280       depth = hgpar_->layerGroupO_[i];
1281     }
1282     return std::make_pair(kx, depth);
1283   }
1284 }
1285 
1286 int HGCalDDDConstants::tileCount(int layer, int ring) const {
1287   int laymin(layer), laymax(layer), ringmin(ring), ringmax(ring), kount(0);
1288   if (layer == 0) {
1289     laymin = hgpar_->firstLayer_;
1290     laymax = lastLayer(true);
1291   }
1292 #ifdef EDM_ML_DEBUG
1293   edm::LogVerbatim("HGCalGeom") << "tileCount: layer " << layer << " ring " << ring << " layerMin/Max " << laymin << ":"
1294                                 << laymax;
1295 #endif
1296   for (int lay = laymin; lay <= laymax; ++lay) {
1297     if (ring < 0) {
1298       int ll = lay - hgpar_->firstLayer_;
1299       ringmin = hgpar_->tileRingRange_[ll].first;
1300       ringmax = hgpar_->tileRingRange_[ll].second;
1301     }
1302 #ifdef EDM_ML_DEBUG
1303     edm::LogVerbatim("HGCalGeom") << "tileCount: lay " << lay << ":" << (lay - hgpar_->firstLayer_) << " rings "
1304                                   << ringmin << ":" << ringmax;
1305 #endif
1306     for (int rin = ringmin; rin <= ringmax; ++rin) {
1307       int indx = HGCalTileIndex::tileIndex(lay, rin + 1, 0);
1308       auto itr = hgpar_->tileInfoMap_.find(indx);
1309 #ifdef EDM_ML_DEBUG
1310       edm::LogVerbatim("HGCalGeom") << "tileCount: rin  " << rin << " indx " << indx << " itr "
1311                                     << (itr != hgpar_->tileInfoMap_.end());
1312 #endif
1313       if (itr != hgpar_->tileInfoMap_.end()) {
1314         for (int k = 0; k < 4; ++k) {
1315           std::bitset<24> b(itr->second.hex[k]);
1316           kount += b.count();
1317         }
1318       }
1319 #ifdef EDM_ML_DEBUG
1320       edm::LogVerbatim("HGCalGeom") << "tileCount: lay|rin  " << lay << ":" << rin << " kount " << kount;
1321 #endif
1322     }
1323   }
1324   return (3 * kount);
1325 }
1326 
1327 bool HGCalDDDConstants::tileExist(int zside, int layer, int ring, int phi) const {
1328   int indx = HGCalTileIndex::tileIndex(layer, ring, 0);
1329   auto itr = hgpar_->tileInfoMap_.find(indx);
1330   bool ok = (itr == hgpar_->tileInfoMap_.end()) ? false : HGCalTileIndex::tileExist(itr->second.hex, zside, phi);
1331   return ok;
1332 }
1333 
1334 HGCalParameters::tileInfo HGCalDDDConstants::tileInfo(int zside, int layer, int ring) const {
1335   int indx = HGCalTileIndex::tileIndex(layer, ring, 0);
1336   auto itr = hgpar_->tileInfoMap_.find(indx);
1337   return ((itr == hgpar_->tileInfoMap_.end()) ? HGCalParameters::tileInfo() : itr->second);
1338 }
1339 
1340 bool HGCalDDDConstants::tilePhiEdge(double phi, int layer, int iphi) const {
1341   double dif1 = std::abs(phi - hgpar_->scintCellSize(layer) * (iphi - 1));
1342   double dif2 = std::abs(phi - hgpar_->scintCellSize(layer) * iphi);
1343 #ifdef EDM_ML_DEBUG
1344   edm::LogVerbatim("HGCalGeomT") << "HGCalDDDConstants::tilePhiEdge:: input: " << phi << ":" << layer << ":" << iphi
1345                                  << " Differences " << dif1 << ":" << dif2;
1346 #endif
1347   return ((dif1 < tol_) || (dif2 < tol_));
1348 }
1349 
1350 bool HGCalDDDConstants::tileRingEdge(double r, int layer, int ring) const {
1351   int type = hgpar_->scintType(layer);
1352   double dif1 = std::abs(r - hgpar_->radiusLayer_[type][ring - 1]);
1353   double dif2 = std::abs(r - hgpar_->radiusLayer_[type][ring]);
1354 #ifdef EDM_ML_DEBUG
1355   edm::LogVerbatim("HGCalGeomT") << "HGCalDDDConstants::tileRingEdge:: input: " << r << ":" << layer << ":" << ring
1356                                  << " Differences " << dif1 << ":" << dif2;
1357 #endif
1358   return ((dif1 < tol_) || (dif2 < tol_));
1359 }
1360 std::pair<int, int> HGCalDDDConstants::tileRings(int layer) const {
1361   if (trapezoidFile()) {
1362     int ll = layer - hgpar_->firstLayer_;
1363     if (ll >= 0 && ll < static_cast<int>(hgpar_->tileRingRange_.size()))
1364       return hgpar_->tileRingRange_[ll];
1365   }
1366   return std::make_pair(0, 0);
1367 }
1368 
1369 std::pair<int, int> HGCalDDDConstants::tileType(int layer, int ring, int phi) const {
1370   int indx = HGCalTileIndex::tileIndex(layer, ring, phi);
1371   int type(-1), sipm(-1);
1372   auto itr = hgpar_->tileInfoMap_.find(indx);
1373   if (itr != hgpar_->tileInfoMap_.end()) {
1374     type = 1 + (itr->second).type;
1375     sipm = ((itr->second).sipm == HGCalTypes::SiPMLarge) ? 0 : 1;
1376   }
1377   return std::make_pair(type, sipm);
1378 }
1379 
1380 int HGCalDDDConstants::waferFromCopy(int copy) const {
1381   const int ncopies = hgpar_->waferCopy_.size();
1382   int wafer(ncopies);
1383   bool result(false);
1384   for (int k = 0; k < ncopies; ++k) {
1385     if (copy == hgpar_->waferCopy_[k]) {
1386       wafer = k;
1387       result = true;
1388       break;
1389     }
1390   }
1391   if (!result) {
1392     wafer = -1;
1393 #ifdef EDM_ML_DEBUG
1394     edm::LogVerbatim("HGCalGeom") << "Cannot find " << copy << " in a list of " << ncopies << " members";
1395     for (int k = 0; k < ncopies; ++k)
1396       edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << hgpar_->waferCopy_[k];
1397 #endif
1398   }
1399 #ifdef EDM_ML_DEBUG
1400   edm::LogVerbatim("HGCalGeom") << "WaferFromCopy " << copy << ":" << wafer << ":" << result;
1401 #endif
1402   return wafer;
1403 }
1404 
1405 void HGCalDDDConstants::waferFromPosition(const double x, const double y, int& wafer, int& icell, int& celltyp) const {
1406   // Input x, y in Geant4 unit and transformed to CMSSW standard
1407   double xx = HGCalParameters::k_ScaleFromDDD * x;
1408   double yy = HGCalParameters::k_ScaleFromDDD * y;
1409   int size_ = static_cast<int>(hgpar_->waferCopy_.size());
1410   wafer = size_;
1411   for (int k = 0; k < size_; ++k) {
1412     double dx = std::abs(xx - hgpar_->waferPosX_[k]);
1413     double dy = std::abs(yy - hgpar_->waferPosY_[k]);
1414     if (dx <= rmax_ && dy <= hexside_) {
1415       if ((dy <= 0.5 * hexside_) || (dx * tan30deg_ <= (hexside_ - dy))) {
1416         wafer = k;
1417         celltyp = hgpar_->waferTypeT_[k];
1418         xx -= hgpar_->waferPosX_[k];
1419         yy -= hgpar_->waferPosY_[k];
1420         break;
1421       }
1422     }
1423   }
1424   if (wafer < size_) {
1425     if (celltyp - 1 == HGCSiliconDetId::HGCalFine)
1426       icell = cellHex(
1427           xx, yy, 0.5 * HGCalParameters::k_ScaleFromDDD * hgpar_->cellSize_[0], hgpar_->cellFineX_, hgpar_->cellFineY_);
1428     else
1429       icell = cellHex(xx,
1430                       yy,
1431                       0.5 * HGCalParameters::k_ScaleFromDDD * hgpar_->cellSize_[1],
1432                       hgpar_->cellCoarseX_,
1433                       hgpar_->cellCoarseY_);
1434   } else {
1435     wafer = -1;
1436 #ifdef EDM_ML_DEBUG
1437     edm::LogWarning("HGCalGeom") << "Cannot get wafer type corresponding to " << x << ":" << y << "    " << xx << ":"
1438                                  << yy;
1439 #endif
1440   }
1441 #ifdef EDM_ML_DEBUG
1442   edm::LogVerbatim("HGCalGeom") << "Position " << x << ":" << y << " Wafer " << wafer << ":" << size_ << " XX " << xx
1443                                 << ":" << yy << " Cell " << icell << " Type " << celltyp;
1444 #endif
1445 }
1446 
1447 void HGCalDDDConstants::waferFromPosition(const double x,
1448                                           const double y,
1449                                           const int zside,
1450                                           const int layer,
1451                                           int& waferU,
1452                                           int& waferV,
1453                                           int& cellU,
1454                                           int& cellV,
1455                                           int& celltype,
1456                                           double& wt,
1457                                           bool extend,
1458                                           bool debug) const {
1459   // Expect x, y as in SIM step
1460   bool waferin = ((waferU == 0) && (waferV == 0));
1461   if (waferin)
1462     waferU = waferV = 1 + hgpar_->waferUVMax_;
1463   cellU = cellV = celltype = 0;
1464   if ((hgpar_->xLayerHex_.empty()) || (hgpar_->yLayerHex_.empty()))
1465     return;
1466   int ll = layer - hgpar_->firstLayer_;
1467   int layertype = layerType(layer);
1468   bool rotx = ((!hgpar_->layerType_.empty()) && (layertype == HGCalTypes::WaferCenterR));
1469   double xx(0), yy(0);
1470   if (rotx) {
1471     std::pair<double, double> xy =
1472         getXY(layer, HGCalParameters::k_ScaleFromDDD * x, HGCalParameters::k_ScaleFromDDD * y, true);
1473     xx = xy.first - hgpar_->xLayerHex_[ll];
1474     yy = xy.second - hgpar_->yLayerHex_[ll];
1475   } else {
1476     xx = HGCalParameters::k_ScaleFromDDD * x - hgpar_->xLayerHex_[ll];
1477     yy = HGCalParameters::k_ScaleFromDDD * y - hgpar_->yLayerHex_[ll];
1478   }
1479   if (debug)
1480     edm::LogVerbatim("HGCalGeom") << "waferFromPosition:: Layer " << layer << ":" << ll << " Rot " << rotx << " X " << x
1481                                   << ":" << xx << " Y " << y << ":" << yy << " side " << zside << " extend " << extend
1482                                   << " initial wafer index " << waferU << ":" << waferV;
1483   ;
1484   double rmax = extend ? rmaxT_ : rmax_;
1485   double hexside = extend ? hexsideT_ : hexside_;
1486   if (waferin) {
1487     double tolmin(100.0);
1488     for (unsigned int k = 0; k < hgpar_->waferPosX_.size(); ++k) {
1489       double dx0(0), dy0(0);
1490       waferU = HGCalWaferIndex::waferU(hgpar_->waferCopy_[k]);
1491       waferV = HGCalWaferIndex::waferV(hgpar_->waferCopy_[k]);
1492       if (cassetteMode()) {
1493         int indx = HGCalWaferIndex::waferIndex(layer, waferU, waferV);
1494         auto ktr = hgpar_->waferInfoMap_.find(indx);
1495         if (ktr != hgpar_->waferInfoMap_.end()) {
1496           auto cshift = hgcassette_.getShift(layer, -1, (ktr->second).cassette);
1497           if (debug)
1498             edm::LogVerbatim("HGCalGeom")
1499                 << "Cassette " << (ktr->second).cassette << " Shift " << cshift.first << ":" << cshift.second;
1500           dx0 = -cshift.first;
1501           dy0 = cshift.second;
1502         }
1503       }
1504       double dx = std::abs(xx - dx0 - hgpar_->waferPosX_[k]);
1505       double dy = std::abs(yy - dy0 - hgpar_->waferPosY_[k]);
1506       constexpr double tolc = 0.01;
1507       if (debug) {
1508         edm::LogVerbatim("HGCalGeom") << "Wafer " << waferU << ":" << waferV << " position " << xx << ":" << yy
1509                                       << " Distance " << dx << ":" << dy << " diff0 " << (dx - rmax) << ":"
1510                                       << (dy - hexside) << " diff1 " << (dy - 0.5 * hexside) << ":"
1511                                       << (dx * tan30deg_ - (hexside - dy));
1512         if ((dx - rmax) <= tolc && (dy - hexside) <= tolc) {
1513           tolmin = std::min(tolmin, (dy - 0.5 * hexside));
1514           tolmin = std::min(tolmin, (dx * tan30deg_ - (hexside - dy)));
1515         }
1516       }
1517       if ((dx - rmax) <= tolc && (dy - hexside) <= tolc) {
1518         if (((dy - 0.5 * hexside) <= tolc) || ((dx * tan30deg_ - (hexside - dy)) <= tolc)) {
1519           if (waferHexagon8File()) {
1520             int index = HGCalWaferIndex::waferIndex(layer, waferU, waferV);
1521             celltype = HGCalWaferType::getType(index, hgpar_->waferInfoMap_);
1522             if (debug)
1523               edm::LogVerbatim("HGCalGeom")
1524                   << "Position (" << x << ", " << y << ") Wafer type:partial:orient:cassette " << celltype << ":"
1525                   << HGCalWaferType::getPartial(index, hgpar_->waferInfoMap_) << ":"
1526                   << HGCalWaferType::getOrient(index, hgpar_->waferInfoMap_) << ":"
1527                   << HGCalWaferType::getCassette(index, hgpar_->waferInfoMap_);
1528           } else {
1529             auto itr = hgpar_->typesInLayers_.find(HGCalWaferIndex::waferIndex(layer, waferU, waferV));
1530             celltype = ((itr == hgpar_->typesInLayers_.end()) ? HGCSiliconDetId::HGCalCoarseThick
1531                                                               : hgpar_->waferTypeL_[itr->second]);
1532           }
1533           if (debug)
1534             edm::LogVerbatim("HGCalGeom")
1535                 << "WaferFromPosition:: Input " << layer << ":" << ll << ":" << hgpar_->firstLayer_ << ":" << rotx
1536                 << ":" << x << ":" << y << ":" << hgpar_->xLayerHex_[ll] << ":" << hgpar_->yLayerHex_[ll] << ":" << xx
1537                 << ":" << yy << " compared with " << hgpar_->waferPosX_[k] << ":" << hgpar_->waferPosY_[k]
1538                 << " difference " << dx << ":" << dy << ":" << dx * tan30deg_ << ":" << (hexside_ - dy)
1539                 << " comparator " << rmax_ << ":" << rmaxT_ << ":" << hexside_ << ":" << hexsideT_ << " wafer "
1540                 << waferU << ":" << waferV << ":" << celltype;
1541           xx -= (dx0 + hgpar_->waferPosX_[k]);
1542           yy -= (dy0 + hgpar_->waferPosY_[k]);
1543           break;
1544         }
1545       }
1546     }
1547     if (debug)
1548       edm::LogVerbatim("HGCalGeom") << "Tolmin " << tolmin;
1549   } else {
1550     for (unsigned int k = 0; k < hgpar_->waferPosX_.size(); ++k) {
1551       double dx0(0), dy0(0);
1552       if ((waferU == HGCalWaferIndex::waferU(hgpar_->waferCopy_[k])) &&
1553           (waferV == HGCalWaferIndex::waferV(hgpar_->waferCopy_[k]))) {
1554         if (cassetteMode()) {
1555           int indx = HGCalWaferIndex::waferIndex(layer, waferU, waferV);
1556           auto ktr = hgpar_->waferInfoMap_.find(indx);
1557           if (ktr != hgpar_->waferInfoMap_.end()) {
1558             auto cshift = hgcassette_.getShift(layer, -1, (ktr->second).cassette);
1559             if (debug)
1560               edm::LogVerbatim("HGCalGeom")
1561                   << "Cassette " << (ktr->second).cassette << " Shift " << cshift.first << ":" << cshift.second;
1562             dx0 = -cshift.first;
1563             dy0 = cshift.second;
1564           }
1565         }
1566         if (waferHexagon8File()) {
1567           int index = HGCalWaferIndex::waferIndex(layer, waferU, waferV);
1568           celltype = HGCalWaferType::getType(index, hgpar_->waferInfoMap_);
1569           if (debug)
1570             edm::LogVerbatim("HGCalGeom") << "Position (" << x << ", " << y << ") Wafer type:partial:orient:cassette "
1571                                           << celltype << ":" << HGCalWaferType::getPartial(index, hgpar_->waferInfoMap_)
1572                                           << ":" << HGCalWaferType::getOrient(index, hgpar_->waferInfoMap_) << ":"
1573                                           << HGCalWaferType::getCassette(index, hgpar_->waferInfoMap_);
1574         } else {
1575           auto itr = hgpar_->typesInLayers_.find(HGCalWaferIndex::waferIndex(layer, waferU, waferV));
1576           celltype = ((itr == hgpar_->typesInLayers_.end()) ? HGCSiliconDetId::HGCalCoarseThick
1577                                                             : hgpar_->waferTypeL_[itr->second]);
1578         }
1579         xx -= (dx0 + hgpar_->waferPosX_[k]);
1580         yy -= (dy0 + hgpar_->waferPosY_[k]);
1581         break;
1582       }
1583     }
1584   }
1585   if ((std::abs(waferU) <= hgpar_->waferUVMax_) && (celltype >= 0)) {
1586     int place(HGCalCell::cellPlacementOld), part(HGCalTypes::WaferFull);
1587     if (cassetteMode()) {
1588       int indx = HGCalWaferIndex::waferIndex(layer, waferU, waferV);
1589       auto ktr = hgpar_->waferInfoMap_.find(indx);
1590       if (ktr != hgpar_->waferInfoMap_.end()) {
1591         place = HGCalCell::cellPlacementIndex(1, HGCalTypes::layerFrontBack(layertype), (ktr->second).orient);
1592         part = (ktr->second).part;
1593         if (debug)
1594           edm::LogVerbatim("HGCalGeom") << "waferFromPosition: frontback " << layertype << ":"
1595                                         << HGCalTypes::layerFrontBack(layertype) << " Orient " << (ktr->second).orient
1596                                         << " place " << place << " part " << part;
1597       }
1598     }
1599     cellHex(xx, yy, celltype, place, part, cellU, cellV, extend, debug);
1600     wt = (((celltype < 2) && (hgpar_->useSimWt_ > 0)) ? (hgpar_->cellThickness_[celltype] / hgpar_->waferThick_) : 1.0);
1601   } else {
1602     cellU = cellV = 2 * hgpar_->nCellsFine_;
1603     wt = 1.0;
1604     celltype = -1;
1605   }
1606   if ((celltype < 0) && debug) {
1607     double x1(xx);
1608     double y1(yy);
1609     edm::LogVerbatim("HGCalGeom") << "waferfFromPosition: Bad type for X " << x << ":" << x1 << ":" << xx << " Y " << y
1610                                   << ":" << y1 << ":" << yy << " Wafer " << waferU << ":" << waferV << " Cell " << cellU
1611                                   << ":" << cellV;
1612     for (unsigned int k = 0; k < hgpar_->waferPosX_.size(); ++k) {
1613       double dx = std::abs(x1 - hgpar_->waferPosX_[k]);
1614       double dy = std::abs(y1 - hgpar_->waferPosY_[k]);
1615       edm::LogVerbatim("HGCalGeom") << "Wafer [" << k << "] Position (" << hgpar_->waferPosX_[k] << ", "
1616                                     << hgpar_->waferPosY_[k] << ") difference " << dx << ":" << dy << ":"
1617                                     << dx * tan30deg_ << ":" << hexside - dy << " Paramerers " << rmax << ":"
1618                                     << hexside;
1619     }
1620   }
1621   edm::LogVerbatim("HGCalGeomX") << "Input x:y:layer " << x << ":" << y << ":" << layer << " Wafer " << waferU << ":"
1622                                  << waferV << " Cell " << cellU << ":" << cellV << ":" << celltype << " wt " << wt;
1623 }
1624 
1625 bool HGCalDDDConstants::waferInLayer(int wafer, int lay, bool reco) const {
1626   const auto& indx = getIndex(lay, reco);
1627   if (indx.first < 0)
1628     return false;
1629   return waferInLayerTest(wafer, indx.first, hgpar_->defineFull_);
1630 }
1631 
1632 bool HGCalDDDConstants::waferFullInLayer(int wafer, int lay, bool reco) const {
1633   const auto& indx = getIndex(lay, reco);
1634   if (indx.first < 0)
1635     return false;
1636   return waferInLayerTest(wafer, indx.first, false);
1637 }
1638 
1639 HGCalParameters::waferInfo HGCalDDDConstants::waferInfo(int lay, int waferU, int waferV) const {
1640   int indx = HGCalWaferIndex::waferIndex(lay, waferU, waferV);
1641   auto itr = hgpar_->waferInfoMap_.find(indx);
1642   return ((itr == hgpar_->waferInfoMap_.end()) ? HGCalParameters::waferInfo() : itr->second);
1643 }
1644 
1645 std::pair<double, double> HGCalDDDConstants::waferParameters(bool reco) const {
1646   if (reco)
1647     return std::make_pair(rmax_, hexside_);
1648   else
1649     return std::make_pair(HGCalParameters::k_ScaleToDDD * rmax_, HGCalParameters::k_ScaleToDDD * hexside_);
1650 }
1651 
1652 std::pair<double, double> HGCalDDDConstants::waferPosition(int wafer, bool reco) const {
1653   double xx(0), yy(0);
1654   if (wafer >= 0 && wafer < static_cast<int>(hgpar_->waferPosX_.size())) {
1655     xx = hgpar_->waferPosX_[wafer];
1656     yy = hgpar_->waferPosY_[wafer];
1657   }
1658   if (!reco) {
1659     xx *= HGCalParameters::k_ScaleToDDD;
1660     yy *= HGCalParameters::k_ScaleToDDD;
1661   }
1662   return std::make_pair(xx, yy);
1663 }
1664 
1665 std::pair<double, double> HGCalDDDConstants::waferPosition(
1666     int lay, int waferU, int waferV, bool reco, bool debug) const {
1667   int ll = lay - hgpar_->firstLayer_;
1668   bool rotx = ((!hgpar_->layerType_.empty()) && (hgpar_->layerType_[ll] == HGCalTypes::WaferCenterR));
1669 #ifdef EDM_ML_DEBUG
1670   if (debug)
1671     edm::LogVerbatim("HGCalGeom") << "Layer " << lay << ":" << ll << " Rotation " << rotx << " U:V " << waferU << ":"
1672                                   << waferV;
1673 #endif
1674   auto xy = waferPositionNoRot(lay, waferU, waferV, reco, debug);
1675   std::pair<double, double> xy0 = (rotx) ? getXY(lay, xy.first, xy.second, false) : xy;
1676 #ifdef EDM_ML_DEBUG
1677   if (debug)
1678     edm::LogVerbatim("HGCalGeom") << "Without and with rotation " << xy.first << ":" << xy.second << ":" << xy0.first
1679                                   << ":" << xy0.second;
1680 #endif
1681   return xy0;
1682 }
1683 
1684 int HGCalDDDConstants::waferFileIndex(unsigned int kk) const {
1685   if (kk < hgpar_->waferInfoMap_.size()) {
1686     auto itr = hgpar_->waferInfoMap_.begin();
1687     std::advance(itr, kk);
1688     return itr->first;
1689   } else
1690     return 0;
1691 }
1692 
1693 std::tuple<int, int, int, int> HGCalDDDConstants::waferFileInfo(unsigned int kk) const {
1694   if (kk < hgpar_->waferInfoMap_.size()) {
1695     auto itr = hgpar_->waferInfoMap_.begin();
1696     std::advance(itr, kk);
1697     return std::make_tuple(itr->second.type, itr->second.part, itr->second.orient, itr->second.cassette);
1698   } else
1699     return std::make_tuple(0, 0, 0, 0);
1700 }
1701 
1702 std::tuple<int, int, int, int> HGCalDDDConstants::waferFileInfoFromIndex(int kk) const {
1703   auto itr = hgpar_->waferInfoMap_.find(kk);
1704   if (itr != hgpar_->waferInfoMap_.end()) {
1705     return std::make_tuple(itr->second.type, itr->second.part, itr->second.orient, itr->second.cassette);
1706   } else
1707     return std::make_tuple(0, 0, 0, 0);
1708 }
1709 
1710 GlobalPoint HGCalDDDConstants::waferLocal2Global(
1711     HepGeom::Point3D<float>& loc, const DetId& id, bool useWafer, bool reco, bool debug) const {
1712   HGCSiliconDetId detid(id);
1713   double x(0), y(0);
1714   if (useWafer) {
1715     auto xyw = waferPositionNoRot(detid.layer(), detid.waferU(), detid.waferV(), reco, debug);
1716     x = xyw.first;
1717     y = xyw.second;
1718   }
1719   auto xy = getXY(detid.layer(), (x + loc.x()), (y + loc.y()), false);
1720   double zz = (detid.zside() < 0) ? -(loc.z() + waferZ(detid.layer(), reco)) : (loc.z() + waferZ(detid.layer(), reco));
1721   double xx = (detid.zside() < 0) ? -xy.first : xy.first;
1722   return GlobalPoint(xx, xy.second, zz);
1723 }
1724 
1725 int HGCalDDDConstants::wafers() const {
1726   int wafer(0);
1727   if (!tileTrapezoid()) {
1728     for (unsigned int i = 0; i < layers(true); ++i) {
1729       int lay = hgpar_->depth_[i];
1730       wafer += modules(lay, true);
1731     }
1732   } else {
1733     wafer = static_cast<int>(hgpar_->moduleLayR_.size());
1734   }
1735   return wafer;
1736 }
1737 
1738 int HGCalDDDConstants::wafers(int layer, int type) const {
1739   int wafer(0);
1740   if (!tileTrapezoid()) {
1741     auto itr = waferLayer_.find(layer);
1742     if (itr != waferLayer_.end()) {
1743       unsigned ity = (type > 0 && type <= 2) ? type : 0;
1744       wafer = (itr->second)[ity];
1745     }
1746   } else {
1747     const auto& index = getIndex(layer, true);
1748     wafer = 1 + hgpar_->lastModule_[index.first] - hgpar_->firstModule_[index.first];
1749   }
1750   return wafer;
1751 }
1752 
1753 int HGCalDDDConstants::waferType(DetId const& id, bool fromFile) const {
1754   int type(1);
1755   if (waferHexagon8()) {
1756     if (fromFile && (waferFileSize() > 0)) {
1757       int layer(0), waferU(0), waferV(0);
1758       if (id.det() != DetId::Forward) {
1759         HGCSiliconDetId hid(id);
1760         layer = hid.layer();
1761         waferU = hid.waferU();
1762         waferV = hid.waferV();
1763       } else {
1764         HFNoseDetId hid(id);
1765         layer = hid.layer();
1766         waferU = hid.waferU();
1767         waferV = hid.waferV();
1768       }
1769       auto itr = hgpar_->waferInfoMap_.find(HGCalWaferIndex::waferIndex(layer, waferU, waferV));
1770       if (itr != hgpar_->waferInfoMap_.end())
1771         type = (itr->second).type;
1772     } else {
1773       type = ((id.det() != DetId::Forward) ? HGCSiliconDetId(id).type() : HFNoseDetId(id).type());
1774     }
1775   } else if (waferHexagon6()) {
1776     type = waferTypeL(HGCalDetId(id).wafer()) - 1;
1777   }
1778   return type;
1779 }
1780 
1781 int HGCalDDDConstants::waferType(int layer, int waferU, int waferV, bool fromFile) const {
1782   int type(HGCSiliconDetId::HGCalCoarseThick);
1783   if (waferHexagon8()) {
1784     if (fromFile && (waferFileSize() > 0)) {
1785       auto itr = hgpar_->waferInfoMap_.find(HGCalWaferIndex::waferIndex(layer, waferU, waferV));
1786       if (itr != hgpar_->waferInfoMap_.end())
1787         type = (itr->second).type;
1788     } else {
1789       auto itr = hgpar_->typesInLayers_.find(HGCalWaferIndex::waferIndex(layer, waferU, waferV));
1790       if (itr != hgpar_->typesInLayers_.end())
1791         type = hgpar_->waferTypeL_[itr->second];
1792     }
1793   } else if (waferHexagon6()) {
1794     if ((waferU >= 0) && (waferU < static_cast<int>(hgpar_->waferTypeL_.size())))
1795       type = (hgpar_->waferTypeL_[waferU] - 1);
1796   }
1797   return type;
1798 }
1799 
1800 std::tuple<int, int, int> HGCalDDDConstants::waferType(HGCSiliconDetId const& id, bool fromFile) const {
1801   const auto& index = HGCalWaferIndex::waferIndex(id.layer(), id.waferU(), id.waferV());
1802   int type(-1), part(-1), orient(-1);
1803   if (fromFile && (waferFileSize() > 0)) {
1804     auto itr = hgpar_->waferInfoMap_.find(index);
1805     if (itr != hgpar_->waferInfoMap_.end()) {
1806       type = (itr->second).type;
1807       part = (itr->second).part;
1808       orient = (itr->second).orient;
1809     }
1810   } else {
1811     auto ktr = hgpar_->typesInLayers_.find(index);
1812     if (ktr != hgpar_->typesInLayers_.end())
1813       type = hgpar_->waferTypeL_[ktr->second];
1814     auto itr = hgpar_->waferTypes_.find(index);
1815     if (itr != hgpar_->waferTypes_.end()) {
1816       if ((itr->second).second < HGCalTypes::k_OffsetRotation) {
1817         orient = (itr->second).second;
1818         if ((itr->second).first == HGCalGeomTools::k_allCorners) {
1819           part = HGCalTypes::WaferFull;
1820         } else if ((itr->second).first == HGCalGeomTools::k_fiveCorners) {
1821           part = HGCalTypes::WaferFive;
1822         } else if ((itr->second).first == HGCalGeomTools::k_fourCorners) {
1823           part = HGCalTypes::WaferHalf;
1824         } else if ((itr->second).first == HGCalGeomTools::k_threeCorners) {
1825           part = HGCalTypes::WaferThree;
1826         }
1827       } else {
1828         part = (itr->second).first;
1829         orient = ((itr->second).second - HGCalTypes::k_OffsetRotation);
1830       }
1831     } else {
1832       part = HGCalTypes::WaferFull;
1833       orient = 0;
1834     }
1835   }
1836   return std::make_tuple(type, part, orient);
1837 }
1838 
1839 std::pair<int, int> HGCalDDDConstants::waferTypeRotation(
1840     int layer, int waferU, int waferV, bool fromFile, bool debug) const {
1841   int type(HGCalTypes::WaferOut), rotn(0);
1842   int wl = HGCalWaferIndex::waferIndex(layer, waferU, waferV);
1843   bool withinList(true);
1844   if (fromFile && (waferFileSize() > 0)) {
1845     auto itr = hgpar_->waferInfoMap_.find(wl);
1846     withinList = (itr != hgpar_->waferInfoMap_.end());
1847     if (withinList) {
1848       type = (itr->second).part;
1849       rotn = (itr->second).orient;
1850     }
1851   } else {
1852     auto itr = hgpar_->waferTypes_.find(wl);
1853     if (waferHexagon8()) {
1854       withinList = (itr != hgpar_->waferTypes_.end());
1855       if (withinList) {
1856         if ((itr->second).second < HGCalTypes::k_OffsetRotation) {
1857           rotn = (itr->second).second;
1858           if ((itr->second).first == HGCalGeomTools::k_allCorners) {
1859             type = HGCalTypes::WaferFull;
1860           } else if ((itr->second).first == HGCalGeomTools::k_fiveCorners) {
1861             type = HGCalTypes::WaferFive;
1862           } else if ((itr->second).first == HGCalGeomTools::k_fourCorners) {
1863             type = HGCalTypes::WaferHalf;
1864           } else if ((itr->second).first == HGCalGeomTools::k_threeCorners) {
1865             type = HGCalTypes::WaferThree;
1866           }
1867         } else {
1868           type = (itr->second).first;
1869           rotn = ((itr->second).second - HGCalTypes::k_OffsetRotation);
1870         }
1871       } else {
1872         type = HGCalTypes::WaferFull;
1873         rotn = HGCalTypes::WaferCorner0;
1874       }
1875     }
1876   }
1877 #ifdef EDM_ML_DEBUG
1878   if (debug)
1879     edm::LogVerbatim("HGCalGeom") << "waferTypeRotation: Layer " << layer << " Wafer " << waferU << ":" << waferV
1880                                   << " Index " << std::hex << wl << std::dec << ":" << withinList << " Type " << type
1881                                   << " Rotation " << rotn;
1882 #endif
1883   return std::make_pair(type, rotn);
1884 }
1885 
1886 bool HGCalDDDConstants::waferVirtual(int layer, int waferU, int waferV) const {
1887   bool type(false);
1888   if (waferHexagon8()) {
1889     int wl = HGCalWaferIndex::waferIndex(layer, waferU, waferV, false);
1890     type = (hgpar_->waferTypes_.find(wl) != hgpar_->waferTypes_.end());
1891   } else if (waferHexagon6()) {
1892     int wl = HGCalWaferIndex::waferIndex(layer, waferU, 0, true);
1893     type = (hgpar_->waferTypes_.find(wl) != hgpar_->waferTypes_.end());
1894   }
1895   return type;
1896 }
1897 
1898 double HGCalDDDConstants::waferZ(int lay, bool reco) const {
1899   const auto& index = getIndex(lay, reco);
1900   if (index.first < 0)
1901     return 0;
1902   else
1903     return (reco ? hgpar_->zLayerHex_[index.first] : HGCalParameters::k_ScaleToDDD * hgpar_->zLayerHex_[index.first]);
1904 }
1905 
1906 int HGCalDDDConstants::cellHex(
1907     double xx, double yy, const double& cellR, const std::vector<double>& posX, const std::vector<double>& posY) const {
1908   int num(0);
1909   const double tol(0.00001);
1910   double cellY = 2.0 * cellR * tan30deg_;
1911   for (unsigned int k = 0; k < posX.size(); ++k) {
1912     double dx = std::abs(xx - posX[k]);
1913     double dy = std::abs(yy - posY[k]);
1914     if (dx <= (cellR + tol) && dy <= (cellY + tol)) {
1915       double xmax = (dy <= 0.5 * cellY) ? cellR : (cellR - (dy - 0.5 * cellY) / tan30deg_);
1916       if (dx <= (xmax + tol)) {
1917         num = k;
1918         break;
1919       }
1920     }
1921   }
1922   return num;
1923 }
1924 
1925 void HGCalDDDConstants::cellHex(
1926     double xloc, double yloc, int cellType, int place, int part, int& cellU, int& cellV, bool extend, bool debug) const {
1927   if (cassetteMode()) {
1928     auto uv = (part == HGCalTypes::WaferFull)
1929                   ? hgcellUV_->cellUVFromXY3(xloc, yloc, place, cellType, true, debug)
1930                   : hgcellUV_->cellUVFromXY1(xloc, yloc, place, cellType, part, true, debug);
1931     cellU = uv.first;
1932     cellV = uv.second;
1933   } else if (waferHexagon8File()) {
1934     auto uv = hgcellUV_->cellUVFromXY3(xloc, yloc, place, cellType, extend, debug);
1935     cellU = uv.first;
1936     cellV = uv.second;
1937   } else {
1938     int ncell = (cellType == 0) ? hgpar_->nCellsFine_ : hgpar_->nCellsCoarse_;
1939     double delY = 2 * rmax_ / (3 * ncell);
1940     double delX = 0.5 * delY * sqrt3_;
1941     double delYT = (extend) ? (2 * rmaxT_ / (3 * ncell)) : delY;
1942     double delXT = 0.5 * delYT * sqrt3_;
1943     double v0 = ((xloc / delY - 1.0) / 1.5);
1944     int cv0 = (v0 > 0) ? (ncell + static_cast<int>(v0 + 0.5)) : (ncell - static_cast<int>(-v0 + 0.5));
1945     double u0 = (0.5 * yloc / delX + 0.5 * cv0);
1946     int cu0 = (u0 > 0) ? (ncell / 2 + static_cast<int>(u0 + 0.5)) : (ncell / 2 - static_cast<int>(-u0 + 0.5));
1947     cu0 = std::max(0, std::min(cu0, 2 * ncell - 1));
1948     cv0 = std::max(0, std::min(cv0, 2 * ncell - 1));
1949     if (cv0 - cu0 >= ncell)
1950       cv0 = cu0 + ncell - 1;
1951     if (debug)
1952       edm::LogVerbatim("HGCalGeom") << "cellHex: input " << xloc << ":" << yloc << ":" << cellType << " parameter "
1953                                     << delX << ":" << delY << " u0 " << u0 << ":" << cu0 << " v0 " << v0 << ":" << cv0;
1954     bool found(false);
1955     static constexpr int shift[3] = {0, 1, -1};
1956     for (int i1 = 0; i1 < 3; ++i1) {
1957       cellU = cu0 + shift[i1];
1958       for (int i2 = 0; i2 < 3; ++i2) {
1959         cellV = cv0 + shift[i2];
1960         if (((cellV - cellU) < ncell) && ((cellU - cellV) <= ncell) && (cellU >= 0) && (cellV >= 0) &&
1961             (cellU < 2 * ncell) && (cellV < 2 * ncell)) {
1962           double xc = (1.5 * (cellV - ncell) + 1.0) * delY;
1963           double yc = (2 * cellU - cellV - ncell) * delX;
1964           if ((std::abs(yloc - yc) <= delXT) && (std::abs(xloc - xc) <= delYT) &&
1965               ((std::abs(xloc - xc) <= 0.5 * delYT) ||
1966                (std::abs(yloc - yc) <= sqrt3_ * (delYT - std::abs(xloc - xc))))) {
1967             if (debug)
1968               edm::LogVerbatim("HGCalGeom")
1969                   << "cellHex: local " << xc << ":" << yc << " difference " << std::abs(xloc - xc) << ":"
1970                   << std::abs(yloc - yc) << ":" << sqrt3_ * (delY - std::abs(yloc - yc)) << " comparator " << delX
1971                   << ":" << delY << " (u,v) = (" << cellU << "," << cellV << ")";
1972             found = true;
1973             break;
1974           }
1975         }
1976       }
1977       if (found)
1978         break;
1979     }
1980     if (!found) {
1981       cellU = cu0;
1982       cellV = cv0;
1983     }
1984   }
1985 }
1986 
1987 std::pair<int, float> HGCalDDDConstants::getIndex(int lay, bool reco) const {
1988   int indx = layerIndex(lay, reco);
1989   if (indx < 0)
1990     return std::make_pair(-1, 0);
1991   float cell(0);
1992   if (waferHexagon6()) {
1993     cell = (reco ? hgpar_->moduleCellR_[0] : hgpar_->moduleCellS_[0]);
1994   } else {
1995     if (waferHexagon8()) {
1996       cell = (reco ? hgpar_->moduleCellR_[0] : hgpar_->moduleCellS_[0]);
1997     } else {
1998       cell = hgpar_->scintCellSize(lay);
1999     }
2000   }
2001   return std::make_pair(indx, cell);
2002 }
2003 
2004 int HGCalDDDConstants::layerFromIndex(int index, bool reco) const {
2005   int ll(-1);
2006   if (waferHexagon6() && reco) {
2007     ll = static_cast<int>(std::find(hgpar_->depthLayerF_.begin(), hgpar_->depthLayerF_.end(), index) -
2008                           hgpar_->depthLayerF_.begin());
2009     if (ll == static_cast<int>(hgpar_->depthLayerF_.size()))
2010       ll = -1;
2011   } else {
2012     ll = static_cast<int>(std::find(hgpar_->layerIndex_.begin(), hgpar_->layerIndex_.end(), index) -
2013                           hgpar_->layerIndex_.begin());
2014     if (ll == static_cast<int>(hgpar_->layerIndex_.size()))
2015       ll = -1;
2016   }
2017 #ifdef EDM_ML_DEBUG
2018   edm::LogVerbatim("HGCalGeom") << "LayerFromIndex for " << index << ":" << reco << ":" << waferHexagon6() << " is"
2019                                 << ll << ":" << (ll + hgpar_->firstLayer_);
2020 #endif
2021   return ((ll < 0) ? ll : (ll + hgpar_->firstLayer_));
2022 }
2023 
2024 bool HGCalDDDConstants::isValidCell(int lay, int wafer, int cell) const {
2025   // Calculate the position of the cell
2026   // Works for options HGCalHexagon/HGCalHexagonFull
2027   double x = hgpar_->waferPosX_[wafer];
2028   double y = hgpar_->waferPosY_[wafer];
2029   if (hgpar_->waferTypeT_[wafer] - 1 == HGCSiliconDetId::HGCalFine) {
2030     x += hgpar_->cellFineX_[cell];
2031     y += hgpar_->cellFineY_[cell];
2032   } else {
2033     x += hgpar_->cellCoarseX_[cell];
2034     y += hgpar_->cellCoarseY_[cell];
2035   }
2036   double rr = sqrt(x * x + y * y);
2037   bool result = ((rr >= hgpar_->rMinLayHex_[lay - 1]) && (rr <= hgpar_->rMaxLayHex_[lay - 1]) &&
2038                  (wafer < static_cast<int>(hgpar_->waferPosX_.size())));
2039 #ifdef EDM_ML_DEBUG
2040   if (!result)
2041     edm::LogVerbatim("HGCalGeom") << "Input " << lay << ":" << wafer << ":" << cell << " Position " << x << ":" << y
2042                                   << ":" << rr << " Compare Limits " << hgpar_->rMinLayHex_[lay - 1] << ":"
2043                                   << hgpar_->rMaxLayHex_[lay - 1] << " Flag " << result;
2044 #endif
2045   return result;
2046 }
2047 
2048 bool HGCalDDDConstants::isValidCell8(int lay, int waferU, int waferV, int cellU, int cellV, int type) const {
2049   bool result(false);
2050   auto partn = waferTypeRotation(lay, waferU, waferV, false, false);
2051 #ifdef EDM_ML_DEBUG
2052   edm::LogVerbatim("HGCalGeom") << "waferHexagon8 " << waferHexagon8File() << ":" << mode_ << ":" << cassetteMode()
2053                                 << " part " << partn.first << ":" << partn.second;
2054 #endif
2055   if (cassetteMode()) {
2056     result = HGCalWaferMask::goodCell(cellU, cellV, partn.first);
2057 #ifdef EDM_ML_DEBUG
2058     edm::LogVerbatim("HGCalGeom") << "Input " << lay << ":" << waferU << ":" << waferV << ":" << cellU << ":" << cellV
2059                                   << " Result " << result << " from goodCell";
2060 #endif
2061   } else {
2062     float x(0), y(0);
2063     int kndx = cellV * 100 + cellU;
2064     if (type == 0) {
2065       auto ktr = hgpar_->cellFineIndex_.find(kndx);
2066       if (ktr != hgpar_->cellFineIndex_.end()) {
2067         x = hgpar_->cellFineX_[ktr->second];
2068         y = hgpar_->cellFineY_[ktr->second];
2069       }
2070 #ifdef EDM_ML_DEBUG
2071       edm::LogVerbatim("HGCalGeom") << "Fine " << cellU << ":" << cellV << ":" << kndx << ":" << x << ":" << y << ":"
2072                                     << (ktr != hgpar_->cellFineIndex_.end());
2073 #endif
2074     } else {
2075       auto ktr = hgpar_->cellCoarseIndex_.find(kndx);
2076       if (ktr != hgpar_->cellCoarseIndex_.end()) {
2077         x = hgpar_->cellCoarseX_[ktr->second];
2078         y = hgpar_->cellCoarseY_[ktr->second];
2079       }
2080 #ifdef EDM_ML_DEBUG
2081       edm::LogVerbatim("HGCalGeom") << "Coarse " << cellU << ":" << cellV << ":" << kndx << ":" << x << ":" << y << ":"
2082                                     << (ktr != hgpar_->cellCoarseIndex_.end());
2083 #endif
2084     }
2085     const auto& xy = waferPositionNoRot(lay, waferU, waferV, true, false);
2086     x += xy.first;
2087     y += xy.second;
2088 #ifdef EDM_ML_DEBUG
2089     edm::LogVerbatim("HGCalGeom") << "With wafer (" << waferU << "," << waferV << ") " << x << ":" << y;
2090 #endif
2091     double rr = sqrt(x * x + y * y);
2092     int ll = lay - hgpar_->firstLayer_;
2093     double tol = waferHexagon8File() ? 0.5 : 0.0;
2094     result = (((rr + tol) >= hgpar_->rMinLayHex_[ll]) && (rr <= hgpar_->rMaxLayHex_[ll]));
2095 #ifdef EDM_ML_DEBUG
2096     edm::LogVerbatim("HGCalGeom") << "Input " << lay << ":" << ll << ":" << waferU << ":" << waferV << ":" << cellU
2097                                   << ":" << cellV << " Position " << x << ":" << y << ":" << rr << " Compare Limits "
2098                                   << hgpar_->rMinLayHex_[ll] << ":" << hgpar_->rMaxLayHex_[ll] << " Flag " << result
2099                                   << " from Radius Limits";
2100 #endif
2101     if (result && waferHexagon8File()) {
2102       int N = (type == 0) ? hgpar_->nCellsFine_ : hgpar_->nCellsCoarse_;
2103       result = HGCalWaferMask::goodCell(cellU, cellV, N, partn.first, partn.second);
2104 #ifdef EDM_ML_DEBUG
2105       edm::LogVerbatim("HGCalGeom") << "Input " << lay << ":" << waferU << ":" << waferV << ":" << cellU << ":" << cellV
2106                                     << " N " << N << " part " << partn.first << ":" << partn.second << " Result "
2107                                     << result << " from goodCell";
2108 #endif
2109     }
2110   }
2111   return result;
2112 }
2113 
2114 int32_t HGCalDDDConstants::waferIndex(int wafer, int index) const {
2115   int layer = layerFromIndex(index, true);
2116   int waferU = HGCalWaferIndex::waferU(hgpar_->waferCopy_[wafer]);
2117   int waferV = HGCalWaferIndex::waferV(hgpar_->waferCopy_[wafer]);
2118   int indx = HGCalWaferIndex::waferIndex(layer, waferU, waferV);
2119 #ifdef EDM_ML_DEBUG
2120   edm::LogVerbatim("HGCalGeom") << "WaferIndex for " << wafer << ":" << index << " (" << layer << ":" << waferU << ":"
2121                                 << waferV << ") " << indx;
2122 #endif
2123   return indx;
2124 }
2125 
2126 bool HGCalDDDConstants::waferInLayerTest(int wafer, int lay, bool full) const {
2127   bool in = (waferHexagon6()) ? true : false;
2128   if (!in) {
2129     double xpos = hgpar_->waferPosX_[wafer] + hgpar_->xLayerHex_[lay];
2130     double ypos = hgpar_->waferPosY_[wafer] + hgpar_->yLayerHex_[lay];
2131     std::pair<int, int> corner = HGCalGeomTools::waferCorner(
2132         xpos, ypos, rmax_, hexside_, hgpar_->rMinLayHex_[lay], hgpar_->rMaxLayHex_[lay], in);
2133     in = (full ? (corner.first > 0) : (corner.first == static_cast<int>(HGCalParameters::k_CornerSize)));
2134     if (in && fullAndPart_) {
2135       int indx = waferIndex(wafer, lay);
2136       in = (hgpar_->waferInfoMap_.find(indx) != hgpar_->waferInfoMap_.end());
2137 #ifdef EDM_ML_DEBUG
2138       if (!in)
2139         edm::LogVerbatim("HGCalGeom") << "WaferInLayerTest: Layer " << lay << " wafer " << wafer << " index " << indx
2140                                       << "( " << HGCalWaferIndex::waferLayer(indx) << ", "
2141                                       << HGCalWaferIndex::waferU(indx) << ", " << HGCalWaferIndex::waferV(indx)
2142                                       << ") in " << in;
2143 #endif
2144     }
2145 #ifdef EDM_ML_DEBUG
2146     edm::LogVerbatim("HGCalGeom") << "WaferInLayerTest: Layer " << lay << " wafer " << wafer << " R-limits "
2147                                   << hgpar_->rMinLayHex_[lay] << ":" << hgpar_->rMaxLayHex_[lay] << " Corners "
2148                                   << corner.first << ":" << corner.second << " In " << in;
2149 #endif
2150   }
2151   return in;
2152 }
2153 
2154 std::pair<double, double> HGCalDDDConstants::waferPositionNoRot(
2155     int lay, int waferU, int waferV, bool reco, bool debug) const {
2156   int ll = lay - hgpar_->firstLayer_;
2157   double x = hgpar_->xLayerHex_[ll];
2158   double y = hgpar_->yLayerHex_[ll];
2159 #ifdef EDM_ML_DEBUG
2160   if (debug)
2161     edm::LogVerbatim("HGCalGeom") << "Layer " << lay << ":" << ll << " Shift " << hgpar_->xLayerHex_[ll] << ":"
2162                                   << hgpar_->yLayerHex_[ll] << " U:V " << waferU << ":" << waferV;
2163 #endif
2164   if (!reco) {
2165     x *= HGCalParameters::k_ScaleToDDD;
2166     y *= HGCalParameters::k_ScaleToDDD;
2167   }
2168   const auto& xy = waferPosition(waferU, waferV, reco);
2169   x += xy.first;
2170   y += xy.second;
2171 #ifdef EDM_ML_DEBUG
2172   if (debug)
2173     edm::LogVerbatim("HGCalGeom") << "With wafer " << x << ":" << y << ":" << xy.first << ":" << xy.second;
2174 #endif
2175   return std::make_pair(x, y);
2176 }
2177 
2178 std::pair<double, double> HGCalDDDConstants::waferPosition(int waferU, int waferV, bool reco) const {
2179   double xx(0), yy(0);
2180   int indx = HGCalWaferIndex::waferIndex(0, waferU, waferV);
2181   auto itr = hgpar_->wafersInLayers_.find(indx);
2182   if (itr != hgpar_->wafersInLayers_.end()) {
2183     xx = hgpar_->waferPosX_[itr->second];
2184     yy = hgpar_->waferPosY_[itr->second];
2185   }
2186   if (!reco) {
2187     xx *= HGCalParameters::k_ScaleToDDD;
2188     yy *= HGCalParameters::k_ScaleToDDD;
2189   }
2190   return std::make_pair(xx, yy);
2191 }
2192 
2193 #include "FWCore/Utilities/interface/typelookup.h"
2194 
2195 TYPELOOKUP_DATA_REG(HGCalDDDConstants);