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