HGCalCommonData/src/HGCalDDDConstants.cc

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