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