Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 13:03:44 

0001 #include "Geometry/HGCalTBCommonData/interface/HGCalTBDDDConstants.h"
0002 
0003 #include "DataFormats/Math/interface/GeantUnits.h"
0004 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "FWCore/Utilities/interface/Exception.h"
0007 #include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h"
0008 #include "Geometry/HGCalCommonData/interface/HGCalTypes.h"
0009 #include "Geometry/HGCalTBCommonData/interface/HGCalTBGeomParameters.h"
0010 
0011 #include <algorithm>
0012 #include <bitset>
0013 #include <iterator>
0014 #include <functional>
0015 #include <numeric>
0016 
0017 //#define EDM_ML_DEBUG
0018 using namespace geant_units::operators;
0019 
0020 HGCalTBDDDConstants::HGCalTBDDDConstants(const HGCalTBParameters* hp, const std::string& name)
0021     : hgpar_(hp), sqrt3_(std::sqrt(3.0)), mode_(hgpar_->mode_) {
0022 #ifdef EDM_ML_DEBUG
0023   edm::LogVerbatim("HGCalGeom") << "Mode " << mode_;
0024 #endif
0025   if (waferHexagon6()) {
0026     rmax_ = (HGCalTBParameters::k_ScaleFromDDD * (hgpar_->waferR_) * std::cos(30._deg));
0027     rmaxT_ = rmax_ + 0.5 * hgpar_->sensorSeparation_;
0028     hexside_ = 2.0 * rmax_ * tan30deg_;
0029     hexsideT_ = 2.0 * rmaxT_ * tan30deg_;
0030 #ifdef EDM_ML_DEBUG
0031     edm::LogVerbatim("HGCalGeom") << "rmax_ " << rmax_ << ":" << rmaxT_ << ":" << hexside_ << ":" << hexsideT_
0032                                   << " CellSize " << 0.5 * HGCalTBParameters::k_ScaleFromDDD * hgpar_->cellSize_[0]
0033                                   << ":" << 0.5 * HGCalTBParameters::k_ScaleFromDDD * hgpar_->cellSize_[1];
0034 #endif
0035   }
0036   // init maps and constants
0037   modHalf_ = 0;
0038   maxWafersPerLayer_ = 0;
0039   for (int simreco = 0; simreco < 2; ++simreco) {
0040     tot_layers_[simreco] = layersInit((bool)simreco);
0041     max_modules_layer_[simreco].resize(tot_layers_[simreco] + 1);
0042     for (unsigned int layer = 1; layer <= tot_layers_[simreco]; ++layer) {
0043       max_modules_layer_[simreco][layer] = modulesInit(layer, (bool)simreco);
0044       if (simreco == 1) {
0045         modHalf_ += max_modules_layer_[simreco][layer];
0046         maxWafersPerLayer_ = std::max(maxWafersPerLayer_, max_modules_layer_[simreco][layer]);
0047 #ifdef EDM_ML_DEBUG
0048         edm::LogVerbatim("HGCalGeom") << "Layer " << layer << " with " << max_modules_layer_[simreco][layer] << ":"
0049                                       << modHalf_ << " modules in RECO";
0050       } else {
0051         edm::LogVerbatim("HGCalGeom") << "Layer " << layer << " with " << max_modules_layer_[simreco][layer]
0052                                       << " modules in SIM";
0053 #endif
0054       }
0055     }
0056 #ifdef EDM_ML_DEBUG
0057     edm::LogVerbatim("HGCalGeom") << "SimReco " << simreco << " with " << tot_layers_[simreco] << " Layers";
0058 #endif
0059   }
0060   tot_wafers_ = wafers();
0061 
0062 #ifdef EDM_ML_DEBUG
0063   edm::LogVerbatim("HGCalGeom") << "HGCalTBDDDConstants initialized for " << name << " with " << layers(false) << ":"
0064                                 << layers(true) << " layers, " << wafers() << ":" << 2 * modHalf_
0065                                 << " wafers with maximum " << maxWafersPerLayer_ << " per layer and "
0066                                 << "maximum of " << maxCells(false) << ":" << maxCells(true) << " cells";
0067 #endif
0068   if (waferHexagon6()) {
0069     int wminT(9999999), wmaxT(-9999999), kount1(0), kount2(0);
0070     for (unsigned int i = 0; i < getTrFormN(); ++i) {
0071       int lay0 = getTrForm(i).lay;
0072       int wmin(9999999), wmax(-9999999), kount(0);
0073       for (int wafer = 0; wafer < sectors(); ++wafer) {
0074         bool waferIn = waferInLayer(wafer, lay0, true);
0075         if (waferIn) {
0076           if (wafer < wmin)
0077             wmin = wafer;
0078           if (wafer > wmax)
0079             wmax = wafer;
0080           ++kount;
0081         }
0082       }
0083       if (wminT > wmin)
0084         wminT = wmin;
0085       if (wmaxT < wmax)
0086         wmaxT = wmax;
0087       if (kount1 < kount)
0088         kount1 = kount;
0089       kount2 += kount;
0090 #ifdef EDM_ML_DEBUG
0091       int lay1 = getIndex(lay0, true).first;
0092       edm::LogVerbatim("HGCalGeom") << "Index " << i << " Layer " << lay0 << ":" << lay1 << " Wafer " << wmin << ":"
0093                                     << wmax << ":" << kount;
0094 #endif
0095       HGCWaferParam a1{{wmin, wmax, kount}};
0096       waferLayer_[lay0] = a1;
0097     }
0098     waferMax_ = std::array<int, 4>{{wminT, wmaxT, kount1, kount2}};
0099 #ifdef EDM_ML_DEBUG
0100     edm::LogVerbatim("HGCalGeom") << "Overall wafer statistics: " << wminT << ":" << wmaxT << ":" << kount1 << ":"
0101                                   << kount2;
0102 #endif
0103   }
0104 }
0105 
0106 std::pair<int, int> HGCalTBDDDConstants::assignCell(float x, float y, int lay, int subSec, bool reco) const {
0107   const auto& index = getIndex(lay, reco);
0108   if (index.first < 0)
0109     return std::make_pair(-1, -1);
0110   if (waferHexagon6()) {
0111     float xx = (reco) ? x : HGCalTBParameters::k_ScaleFromDDD * x;
0112     float yy = (reco) ? y : HGCalTBParameters::k_ScaleFromDDD * y;
0113 
0114     // First the wafer
0115     int wafer = cellHex(xx, yy, rmax_, hgpar_->waferPosX_, hgpar_->waferPosY_);
0116     if (wafer < 0 || wafer >= static_cast<int>(hgpar_->waferTypeT_.size())) {
0117       edm::LogWarning("HGCalGeom") << "Wafer no. out of bound for " << wafer << ":" << (hgpar_->waferTypeT_).size()
0118                                    << ":" << (hgpar_->waferPosX_).size() << ":" << (hgpar_->waferPosY_).size()
0119                                    << " ***** ERROR *****";
0120       return std::make_pair(-1, -1);
0121     } else {
0122       // Now the cell
0123       xx -= hgpar_->waferPosX_[wafer];
0124       yy -= hgpar_->waferPosY_[wafer];
0125       if (hgpar_->waferTypeT_[wafer] == 1)
0126         return std::make_pair(wafer,
0127                               cellHex(xx,
0128                                       yy,
0129                                       0.5 * HGCalTBParameters::k_ScaleFromDDD * hgpar_->cellSize_[0],
0130                                       hgpar_->cellFineX_,
0131                                       hgpar_->cellFineY_));
0132       else
0133         return std::make_pair(wafer,
0134                               cellHex(xx,
0135                                       yy,
0136                                       0.5 * HGCalTBParameters::k_ScaleFromDDD * hgpar_->cellSize_[1],
0137                                       hgpar_->cellCoarseX_,
0138                                       hgpar_->cellCoarseY_));
0139     }
0140   } else {
0141     return std::make_pair(-1, -1);
0142   }
0143 }
0144 
0145 double HGCalTBDDDConstants::cellSizeHex(int type) const {
0146   int indx = (type == 1) ? 1 : 0;
0147   double cell = 0.5 * HGCalTBParameters::k_ScaleFromDDD * hgpar_->cellSize_[indx];
0148   return cell;
0149 }
0150 
0151 double HGCalTBDDDConstants::cellThickness(int layer, int wafer) const {
0152   double thick(-1);
0153   int type = waferType(layer, wafer);
0154   if (type >= 0)
0155     thick = 100.0 * (type + 1);  // type = 1,2,3 for 100,200,300 micron
0156   return thick;
0157 }
0158 
0159 double HGCalTBDDDConstants::distFromEdgeHex(double x, double y, double z) const {
0160   // Assming the point is within a hexagonal plane of the wafer, calculate
0161   // the shortest distance from the edge
0162   if (z < 0)
0163     x = -x;
0164   double dist(0);
0165   // Input x, y in Geant4 unit and transformed to CMSSW standard
0166   double xx = HGCalTBParameters::k_ScaleFromDDD * x;
0167   double yy = HGCalTBParameters::k_ScaleFromDDD * y;
0168   int sizew = static_cast<int>(hgpar_->waferPosX_.size());
0169   int wafer = sizew;
0170   // Transform to the local coordinate frame of the wafer first
0171   for (int k = 0; k < sizew; ++k) {
0172     double dx = std::abs(xx - hgpar_->waferPosX_[k]);
0173     double dy = std::abs(yy - hgpar_->waferPosY_[k]);
0174     if ((dx <= rmax_) && (dy <= hexside_) && ((dy <= 0.5 * hexside_) || (dx * tan30deg_ <= (hexside_ - dy)))) {
0175       wafer = k;
0176       xx -= hgpar_->waferPosX_[k];
0177       yy -= hgpar_->waferPosY_[k];
0178       break;
0179     }
0180   }
0181   // Look at only one quarter (both x,y are positive)
0182   if (wafer < sizew) {
0183     if (std::abs(yy) < 0.5 * hexside_) {
0184       dist = rmax_ - std::abs(xx);
0185     } else {
0186       dist = 0.5 * ((rmax_ - std::abs(xx)) - sqrt3_ * (std::abs(yy) - 0.5 * hexside_));
0187     }
0188   } else {
0189     dist = 0;
0190   }
0191   dist *= HGCalTBParameters::k_ScaleToDDD;
0192 #ifdef EDM_ML_DEBUG
0193   edm::LogVerbatim("HGCalGeom") << "DistFromEdgeHex: Local " << xx << ":" << yy << " wafer " << wafer << " flag "
0194                                 << (wafer < sizew) << " Distance " << rmax_ << ":" << (rmax_ - std::abs(xx)) << ":"
0195                                 << (std::abs(yy) - 0.5 * hexside_) << ":" << 0.5 * hexside_ << ":" << dist;
0196 #endif
0197   return dist;
0198 }
0199 
0200 int HGCalTBDDDConstants::getLayer(double z, bool reco) const {
0201   // Get the layer # from the gloabl z coordinate
0202   unsigned int k = 0;
0203   double zz = (reco ? std::abs(z) : HGCalTBParameters::k_ScaleFromDDD * std::abs(z));
0204   const auto& zLayerHex = hgpar_->zLayerHex_;
0205   auto itr = std::find_if(zLayerHex.begin() + 1, zLayerHex.end(), [&k, &zz, &zLayerHex](double zLayer) {
0206     ++k;
0207     return zz < 0.5 * (zLayerHex[k - 1] + zLayerHex[k]);
0208   });
0209   int lay = (itr == zLayerHex.end()) ? static_cast<int>(zLayerHex.size()) : k;
0210   if (waferHexagon6() && reco) {
0211     int indx = layerIndex(lay, false);
0212     if (indx >= 0)
0213       lay = hgpar_->layerGroupO_[indx];
0214   } else {
0215     lay += (hgpar_->firstLayer_ - 1);
0216   }
0217   return lay;
0218 }
0219 
0220 HGCalTBParameters::hgtrap HGCalTBDDDConstants::getModule(unsigned int indx, bool hexType, bool reco) const {
0221   HGCalTBParameters::hgtrap mytr;
0222   if (hexType) {
0223     if (indx >= hgpar_->waferTypeL_.size())
0224       edm::LogWarning("HGCalGeom") << "Wafer no. out bound for index " << indx << ":" << (hgpar_->waferTypeL_).size()
0225                                    << ":" << (hgpar_->waferPosX_).size() << ":" << (hgpar_->waferPosY_).size()
0226                                    << " ***** ERROR *****";
0227     unsigned int type =
0228         ((indx < hgpar_->waferTypeL_.size()) ? hgpar_->waferTypeL_[indx] - 1 : HGCalTBParameters::HGCalCoarseThick);
0229     mytr = hgpar_->getModule(type, reco);
0230   } else {
0231     mytr = hgpar_->getModule(indx, reco);
0232   }
0233   return mytr;
0234 }
0235 
0236 std::vector<HGCalTBParameters::hgtrap> HGCalTBDDDConstants::getModules() const {
0237   std::vector<HGCalTBParameters::hgtrap> mytrs;
0238   for (unsigned int k = 0; k < hgpar_->moduleLayR_.size(); ++k)
0239     mytrs.emplace_back(hgpar_->getModule(k, true));
0240   return mytrs;
0241 }
0242 
0243 std::vector<HGCalTBParameters::hgtrform> HGCalTBDDDConstants::getTrForms() const {
0244   std::vector<HGCalTBParameters::hgtrform> mytrs;
0245   for (unsigned int k = 0; k < hgpar_->trformIndex_.size(); ++k)
0246     mytrs.emplace_back(hgpar_->getTrForm(k));
0247   return mytrs;
0248 }
0249 
0250 bool HGCalTBDDDConstants::isHalfCell(int waferType, int cell) const {
0251   if (waferType < 1 || cell < 0)
0252     return false;
0253   return waferType == 2 ? hgpar_->cellCoarseHalf_[cell] : hgpar_->cellFineHalf_[cell];
0254 }
0255 
0256 bool HGCalTBDDDConstants::isValidHex(int lay, int mod, int cell, bool reco) const {
0257   // Check validity for a layer|wafer|cell of pre-TDR version
0258   bool result(false), resultMod(false);
0259   int cellmax(0);
0260   if (waferHexagon6()) {
0261     int32_t copyNumber = hgpar_->waferCopy_[mod];
0262     result = ((lay > 0 && lay <= static_cast<int>(layers(reco))));
0263     if (result) {
0264       const int32_t lay_idx = reco ? (lay - 1) * 3 + 1 : lay;
0265       const auto& the_modules = hgpar_->copiesInLayers_[lay_idx];
0266       auto moditr = the_modules.find(copyNumber);
0267       result = resultMod = (moditr != the_modules.end());
0268 #ifdef EDM_ML_DEBUG
0269       if (!result)
0270         edm::LogVerbatim("HGCalGeom") << "HGCalTBDDDConstants: Layer " << lay << ":" << lay_idx << " Copy "
0271                                       << copyNumber << ":" << mod << " Flag " << result;
0272 #endif
0273       if (result) {
0274         if (moditr->second >= 0) {
0275           if (mod >= static_cast<int>(hgpar_->waferTypeT_.size()))
0276             edm::LogWarning("HGCalGeom") << "Module no. out of bound for " << mod << " to be compared with "
0277                                          << (hgpar_->waferTypeT_).size() << " ***** ERROR *****";
0278           cellmax = ((hgpar_->waferTypeT_[mod] - 1 == HGCalTBParameters::HGCalFine)
0279                          ? static_cast<int>(hgpar_->cellFineX_.size())
0280                          : static_cast<int>(hgpar_->cellCoarseX_.size()));
0281           result = (cell >= 0 && cell <= cellmax);
0282         } else {
0283           result = isValidCell(lay_idx, mod, cell);
0284         }
0285       }
0286     }
0287   }
0288 
0289 #ifdef EDM_ML_DEBUG
0290   if (!result)
0291     edm::LogVerbatim("HGCalGeom") << "HGCalTBDDDConstants: Layer " << lay << ":"
0292                                   << (lay > 0 && (lay <= static_cast<int>(layers(reco)))) << " Module " << mod << ":"
0293                                   << resultMod << " Cell " << cell << ":" << cellmax << ":"
0294                                   << (cell >= 0 && cell <= cellmax) << ":" << maxCells(reco);
0295 #endif
0296   return result;
0297 }
0298 
0299 int HGCalTBDDDConstants::lastLayer(bool reco) const {
0300   return (hgpar_->firstLayer_ + tot_layers_[static_cast<int>(reco)] - 1);
0301 }
0302 
0303 int HGCalTBDDDConstants::layerIndex(int lay, bool reco) const {
0304   int ll = lay - hgpar_->firstLayer_;
0305   if (ll < 0 || ll >= static_cast<int>(hgpar_->layerIndex_.size()))
0306     return -1;
0307   if (waferHexagon6()) {
0308     if (reco && ll >= static_cast<int>(hgpar_->depthIndex_.size()))
0309       return -1;
0310     return (reco ? hgpar_->depthLayerF_[ll] : hgpar_->layerIndex_[ll]);
0311   }
0312   return -1;
0313 }
0314 
0315 unsigned int HGCalTBDDDConstants::layers(bool reco) const { return tot_layers_[static_cast<int>(reco)]; }
0316 
0317 unsigned int HGCalTBDDDConstants::layersInit(bool reco) const {
0318   return (reco ? hgpar_->depthIndex_.size() : hgpar_->layerIndex_.size());
0319 }
0320 
0321 std::pair<float, float> HGCalTBDDDConstants::locateCell(int cell, int lay, int type, bool reco) const {
0322   // type refers to wafer # for hexagon cell
0323   float x(999999.), y(999999.);
0324   const auto& index = getIndex(lay, reco);
0325   int i = index.first;
0326   if (i < 0)
0327     return std::make_pair(x, y);
0328   if (waferHexagon6()) {
0329     x = hgpar_->waferPosX_[type];
0330     y = hgpar_->waferPosY_[type];
0331 #ifdef EDM_ML_DEBUG
0332     float x0(x), y0(y);
0333 #endif
0334     if (hgpar_->waferTypeT_[type] - 1 == HGCalTBParameters::HGCalFine) {
0335       x += hgpar_->cellFineX_[cell];
0336       y += hgpar_->cellFineY_[cell];
0337     } else {
0338       x += hgpar_->cellCoarseX_[cell];
0339       y += hgpar_->cellCoarseY_[cell];
0340     }
0341 #ifdef EDM_ML_DEBUG
0342     edm::LogVerbatim("HGCalGeom") << "LocateCell (Wafer) " << x0 << ":" << y0 << " Final " << x << ":" << y;
0343 #endif
0344     if (!reco) {
0345       x *= HGCalTBParameters::k_ScaleToDDD;
0346       y *= HGCalTBParameters::k_ScaleToDDD;
0347     }
0348   }
0349   return std::make_pair(x, y);
0350 }
0351 
0352 std::pair<float, float> HGCalTBDDDConstants::locateCellHex(int cell, int wafer, bool reco) const {
0353   float x(0), y(0);
0354   if (hgpar_->waferTypeT_[wafer] - 1 == HGCalTBParameters::HGCalFine) {
0355     x = hgpar_->cellFineX_[cell];
0356     y = hgpar_->cellFineY_[cell];
0357   } else {
0358     x = hgpar_->cellCoarseX_[cell];
0359     y = hgpar_->cellCoarseY_[cell];
0360   }
0361   if (!reco) {
0362     x *= HGCalTBParameters::k_ScaleToDDD;
0363     y *= HGCalTBParameters::k_ScaleToDDD;
0364   }
0365   return std::make_pair(x, y);
0366 }
0367 
0368 int HGCalTBDDDConstants::maxCells(bool reco) const {
0369   int cells(0);
0370   for (unsigned int i = 0; i < layers(reco); ++i) {
0371     int lay = reco ? hgpar_->depth_[i] : hgpar_->layer_[i];
0372     if (cells < maxCells(lay, reco))
0373       cells = maxCells(lay, reco);
0374   }
0375   return cells;
0376 }
0377 
0378 int HGCalTBDDDConstants::maxCells(int lay, bool reco) const {
0379   const auto& index = getIndex(lay, reco);
0380   if ((index.first < 0) || (!waferHexagon6()))
0381     return 0;
0382   unsigned int cells(0);
0383   for (unsigned int k = 0; k < hgpar_->waferTypeT_.size(); ++k) {
0384     if (waferInLayerTest(k, index.first)) {
0385       unsigned int cell = (hgpar_->waferTypeT_[k] - 1 == HGCalTBParameters::HGCalFine) ? (hgpar_->cellFineX_.size())
0386                                                                                        : (hgpar_->cellCoarseX_.size());
0387       if (cell > cells)
0388         cells = cell;
0389     }
0390   }
0391   return static_cast<int>(cells);
0392 }
0393 
0394 int HGCalTBDDDConstants::maxRows(int lay, bool reco) const {
0395   int kymax(0);
0396   const auto& index = getIndex(lay, reco);
0397   int i = index.first;
0398   if ((i >= 0) && waferHexagon6()) {
0399     for (unsigned int k = 0; k < hgpar_->waferCopy_.size(); ++k) {
0400       if (waferInLayerTest(k, i)) {
0401         int ky = ((hgpar_->waferCopy_[k]) / 100) % 100;
0402         if (ky > kymax)
0403           kymax = ky;
0404       }
0405     }
0406   }
0407   return kymax;
0408 }
0409 
0410 int HGCalTBDDDConstants::modifyUV(int uv, int type1, int type2) const {
0411   // Modify u/v for transition of type1 to type2
0412   return (((type1 == type2) || (type1 * type2 != 0)) ? uv : ((type1 == 0) ? (2 * uv + 1) / 3 : (3 * uv) / 2));
0413 }
0414 
0415 int HGCalTBDDDConstants::modules(int lay, bool reco) const {
0416   if (getIndex(lay, reco).first < 0)
0417     return 0;
0418   else
0419     return max_modules_layer_[static_cast<int>(reco)][lay];
0420 }
0421 
0422 int HGCalTBDDDConstants::modulesInit(int lay, bool reco) const {
0423   int nmod(0);
0424   const auto& index = getIndex(lay, reco);
0425   if (index.first < 0)
0426     return nmod;
0427   for (unsigned int k = 0; k < hgpar_->waferPosX_.size(); ++k) {
0428     if (waferInLayerTest(k, index.first))
0429       ++nmod;
0430   }
0431   return nmod;
0432 }
0433 
0434 double HGCalTBDDDConstants::mouseBite(bool reco) const {
0435   return (reco ? hgpar_->mouseBite_ : HGCalTBParameters::k_ScaleToDDD * hgpar_->mouseBite_);
0436 }
0437 
0438 int HGCalTBDDDConstants::numberCells(bool reco) const {
0439   int cells(0);
0440   unsigned int nlayer = (reco) ? hgpar_->depth_.size() : hgpar_->layer_.size();
0441   for (unsigned k = 0; k < nlayer; ++k) {
0442     std::vector<int> ncells = numberCells(((reco) ? hgpar_->depth_[k] : hgpar_->layer_[k]), reco);
0443     cells = std::accumulate(ncells.begin(), ncells.end(), cells);
0444   }
0445   return cells;
0446 }
0447 
0448 std::vector<int> HGCalTBDDDConstants::numberCells(int lay, bool reco) const {
0449   const auto& index = getIndex(lay, reco);
0450   int i = index.first;
0451   std::vector<int> ncell;
0452   if ((i >= 0) && (waferHexagon6())) {
0453     for (unsigned int k = 0; k < hgpar_->waferTypeT_.size(); ++k) {
0454       if (waferInLayerTest(k, i)) {
0455         unsigned int cell = (hgpar_->waferTypeT_[k] - 1 == HGCalTBParameters::HGCalFine)
0456                                 ? (hgpar_->cellFineX_.size())
0457                                 : (hgpar_->cellCoarseX_.size());
0458         ncell.emplace_back(static_cast<int>(cell));
0459       }
0460     }
0461   }
0462   return ncell;
0463 }
0464 
0465 int HGCalTBDDDConstants::numberCellsHexagon(int wafer) const {
0466   if (wafer >= 0 && wafer < static_cast<int>(hgpar_->waferTypeT_.size())) {
0467     if (hgpar_->waferTypeT_[wafer] - 1 == 0)
0468       return static_cast<int>(hgpar_->cellFineX_.size());
0469     else
0470       return static_cast<int>(hgpar_->cellCoarseX_.size());
0471   } else {
0472     return 0;
0473   }
0474 }
0475 
0476 std::pair<double, double> HGCalTBDDDConstants::rangeR(double z, bool reco) const {
0477   double rmin(0), rmax(0);
0478   if (!reco) {
0479     rmin *= HGCalTBParameters::k_ScaleToDDD;
0480     rmax *= HGCalTBParameters::k_ScaleToDDD;
0481   }
0482 #ifdef EDM_ML_DEBUG
0483   edm::LogVerbatim("HGCalGeom") << "HGCalTBDDDConstants:rangeR: " << z << ":0"
0484                                 << " R " << rmin << ":" << rmax;
0485 #endif
0486   return std::make_pair(rmin, rmax);
0487 }
0488 
0489 std::pair<double, double> HGCalTBDDDConstants::rangeRLayer(int lay, bool reco) const {
0490   double rmin(0), rmax(0);
0491   const auto& index = getIndex(lay, reco);
0492   if (index.first >= 0 && index.first < static_cast<int>(hgpar_->rMinLayHex_.size())) {
0493     rmin = hgpar_->rMinLayHex_[index.first];
0494     rmax = hgpar_->rMaxLayHex_[index.first];
0495   }
0496   if (!reco) {
0497     rmin *= HGCalTBParameters::k_ScaleToDDD;
0498     rmax *= HGCalTBParameters::k_ScaleToDDD;
0499   }
0500 #ifdef EDM_ML_DEBUG
0501   edm::LogVerbatim("HGCalGeom") << "HGCalTBDDDConstants:rangeR: " << lay << ":" << index.first << " R " << rmin << ":"
0502                                 << rmax;
0503 #endif
0504   return std::make_pair(rmin, rmax);
0505 }
0506 
0507 std::pair<double, double> HGCalTBDDDConstants::rangeZ(bool reco) const {
0508   double zmin = (hgpar_->zLayerHex_[0] - hgpar_->waferThick_);
0509   double zmax = (hgpar_->zLayerHex_[hgpar_->zLayerHex_.size() - 1] + hgpar_->waferThick_);
0510 #ifdef EDM_ML_DEBUG
0511   edm::LogVerbatim("HGCalGeom") << "HGCalTBDDDConstants:rangeZ: " << zmin << ":" << zmax << ":" << hgpar_->waferThick_;
0512 #endif
0513   if (!reco) {
0514     zmin *= HGCalTBParameters::k_ScaleToDDD;
0515     zmax *= HGCalTBParameters::k_ScaleToDDD;
0516   }
0517   return std::make_pair(zmin, zmax);
0518 }
0519 
0520 std::pair<int, int> HGCalTBDDDConstants::rowColumnWafer(int wafer) const {
0521   int row(0), col(0);
0522   if (wafer < static_cast<int>(hgpar_->waferCopy_.size())) {
0523     int copy = hgpar_->waferCopy_[wafer];
0524     col = HGCalTypes::getUnpackedU(copy);
0525     row = HGCalTypes::getUnpackedV(copy);
0526   }
0527   return std::make_pair(row, col);
0528 }
0529 
0530 std::pair<int, int> HGCalTBDDDConstants::simToReco(int cell, int lay, int mod, bool half) const {
0531   return std::make_pair(cell, lay);
0532 }
0533 
0534 int HGCalTBDDDConstants::waferFromCopy(int copy) const {
0535   const int ncopies = hgpar_->waferCopy_.size();
0536   int wafer(ncopies);
0537   bool result(false);
0538   for (int k = 0; k < ncopies; ++k) {
0539     if (copy == hgpar_->waferCopy_[k]) {
0540       wafer = k;
0541       result = true;
0542       break;
0543     }
0544   }
0545   if (!result) {
0546     wafer = -1;
0547 #ifdef EDM_ML_DEBUG
0548     edm::LogVerbatim("HGCalGeom") << "Cannot find " << copy << " in a list of " << ncopies << " members";
0549     for (int k = 0; k < ncopies; ++k)
0550       edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << hgpar_->waferCopy_[k];
0551 #endif
0552   }
0553 #ifdef EDM_ML_DEBUG
0554   edm::LogVerbatim("HGCalGeom") << "WaferFromCopy " << copy << ":" << wafer << ":" << result;
0555 #endif
0556   return wafer;
0557 }
0558 
0559 void HGCalTBDDDConstants::waferFromPosition(const double x, const double y, int& wafer, int& icell, int& celltyp) const {
0560   // Input x, y in Geant4 unit and transformed to CMSSW standard
0561   double xx = HGCalTBParameters::k_ScaleFromDDD * x;
0562   double yy = HGCalTBParameters::k_ScaleFromDDD * y;
0563   int size_ = static_cast<int>(hgpar_->waferCopy_.size());
0564   wafer = size_;
0565   for (int k = 0; k < size_; ++k) {
0566     double dx = std::abs(xx - hgpar_->waferPosX_[k]);
0567     double dy = std::abs(yy - hgpar_->waferPosY_[k]);
0568     if (dx <= rmax_ && dy <= hexside_) {
0569       if ((dy <= 0.5 * hexside_) || (dx * tan30deg_ <= (hexside_ - dy))) {
0570         wafer = k;
0571         celltyp = hgpar_->waferTypeT_[k];
0572         xx -= hgpar_->waferPosX_[k];
0573         yy -= hgpar_->waferPosY_[k];
0574         break;
0575       }
0576     }
0577   }
0578   if (wafer < size_) {
0579     if (celltyp - 1 == HGCalTBParameters::HGCalFine)
0580       icell = cellHex(xx,
0581                       yy,
0582                       0.5 * HGCalTBParameters::k_ScaleFromDDD * hgpar_->cellSize_[0],
0583                       hgpar_->cellFineX_,
0584                       hgpar_->cellFineY_);
0585     else
0586       icell = cellHex(xx,
0587                       yy,
0588                       0.5 * HGCalTBParameters::k_ScaleFromDDD * hgpar_->cellSize_[1],
0589                       hgpar_->cellCoarseX_,
0590                       hgpar_->cellCoarseY_);
0591   } else {
0592     wafer = -1;
0593 #ifdef EDM_ML_DEBUG
0594     edm::LogWarning("HGCalGeom") << "Cannot get wafer type corresponding to " << x << ":" << y << "    " << xx << ":"
0595                                  << yy;
0596 #endif
0597   }
0598 #ifdef EDM_ML_DEBUG
0599   edm::LogVerbatim("HGCalGeom") << "Position " << x << ":" << y << " Wafer " << wafer << ":" << size_ << " XX " << xx
0600                                 << ":" << yy << " Cell " << icell << " Type " << celltyp;
0601 #endif
0602 }
0603 
0604 bool HGCalTBDDDConstants::waferInLayer(int wafer, int lay, bool reco) const {
0605   const auto& indx = getIndex(lay, reco);
0606   if (indx.first < 0)
0607     return false;
0608   return waferInLayerTest(wafer, indx.first);
0609 }
0610 
0611 bool HGCalTBDDDConstants::waferFullInLayer(int wafer, int lay, bool reco) const {
0612   const auto& indx = getIndex(lay, reco);
0613   if (indx.first < 0)
0614     return false;
0615   return waferInLayerTest(wafer, indx.first);
0616 }
0617 
0618 std::pair<double, double> HGCalTBDDDConstants::waferParameters(bool reco) const {
0619   if (reco)
0620     return std::make_pair(rmax_, hexside_);
0621   else
0622     return std::make_pair(HGCalTBParameters::k_ScaleToDDD * rmax_, HGCalTBParameters::k_ScaleToDDD * hexside_);
0623 }
0624 
0625 std::pair<double, double> HGCalTBDDDConstants::waferPosition(int wafer, bool reco) const {
0626   double xx(0), yy(0);
0627   if (wafer >= 0 && wafer < static_cast<int>(hgpar_->waferPosX_.size())) {
0628     xx = hgpar_->waferPosX_[wafer];
0629     yy = hgpar_->waferPosY_[wafer];
0630   }
0631   if (!reco) {
0632     xx *= HGCalTBParameters::k_ScaleToDDD;
0633     yy *= HGCalTBParameters::k_ScaleToDDD;
0634   }
0635   return std::make_pair(xx, yy);
0636 }
0637 
0638 int HGCalTBDDDConstants::wafers() const { return static_cast<int>(hgpar_->moduleLayR_.size()); }
0639 
0640 int HGCalTBDDDConstants::wafers(int layer, int type) const {
0641   int wafer(0);
0642   auto itr = waferLayer_.find(layer);
0643   if (itr != waferLayer_.end()) {
0644     unsigned ity = (type > 0 && type <= 2) ? type : 0;
0645     wafer = (itr->second)[ity];
0646   }
0647   return wafer;
0648 }
0649 
0650 int HGCalTBDDDConstants::waferType(DetId const& id) const {
0651   return waferType(HGCalDetId(id).layer(), HGCalDetId(id).wafer());
0652 }
0653 
0654 int HGCalTBDDDConstants::waferType(int layer, int wafer) const {
0655   int type(HGCalTBParameters::HGCalCoarseThick);
0656   if ((wafer >= 0) && (wafer < static_cast<int>(hgpar_->waferTypeL_.size())))
0657     type = (hgpar_->waferTypeL_[wafer] - 1);
0658   return type;
0659 }
0660 
0661 bool HGCalTBDDDConstants::waferVirtual(int layer, int wafer) const {
0662   int wl = HGCalWaferIndex::waferIndex(layer, wafer, 0, true);
0663   return (hgpar_->waferTypes_.find(wl) != hgpar_->waferTypes_.end());
0664 }
0665 
0666 double HGCalTBDDDConstants::waferZ(int lay, bool reco) const {
0667   const auto& index = getIndex(lay, reco);
0668   if (index.first < 0)
0669     return 0;
0670   else
0671     return (reco ? hgpar_->zLayerHex_[index.first] : HGCalTBParameters::k_ScaleToDDD * hgpar_->zLayerHex_[index.first]);
0672 }
0673 
0674 int HGCalTBDDDConstants::cellHex(
0675     double xx, double yy, const double& cellR, const std::vector<double>& posX, const std::vector<double>& posY) const {
0676   int num(0);
0677   const double tol(0.00001);
0678   double cellY = 2.0 * cellR * tan30deg_;
0679   for (unsigned int k = 0; k < posX.size(); ++k) {
0680     double dx = std::abs(xx - posX[k]);
0681     double dy = std::abs(yy - posY[k]);
0682     if (dx <= (cellR + tol) && dy <= (cellY + tol)) {
0683       double xmax = (dy <= 0.5 * cellY) ? cellR : (cellR - (dy - 0.5 * cellY) / tan30deg_);
0684       if (dx <= (xmax + tol)) {
0685         num = k;
0686         break;
0687       }
0688     }
0689   }
0690   return num;
0691 }
0692 
0693 std::pair<int, float> HGCalTBDDDConstants::getIndex(int lay, bool reco) const {
0694   int indx = layerIndex(lay, reco);
0695   if (indx < 0)
0696     return std::make_pair(-1, 0);
0697   float cell(0);
0698   if (waferHexagon6()) {
0699     cell = (reco ? hgpar_->moduleCellR_[0] : hgpar_->moduleCellS_[0]);
0700   }
0701   return std::make_pair(indx, cell);
0702 }
0703 
0704 int HGCalTBDDDConstants::layerFromIndex(int index, bool reco) const {
0705   int ll(-1);
0706   if (waferHexagon6() && reco) {
0707     ll = static_cast<int>(std::find(hgpar_->depthLayerF_.begin(), hgpar_->depthLayerF_.end(), index) -
0708                           hgpar_->depthLayerF_.begin());
0709     if (ll == static_cast<int>(hgpar_->depthLayerF_.size()))
0710       ll = -1;
0711   }
0712 #ifdef EDM_ML_DEBUG
0713   edm::LogVerbatim("HGCalGeom") << "LayerFromIndex for " << index << ":" << reco << ":" << waferHexagon6() << " is"
0714                                 << ll << ":" << (ll + hgpar_->firstLayer_);
0715 #endif
0716   return ((ll < 0) ? ll : (ll + hgpar_->firstLayer_));
0717 }
0718 
0719 bool HGCalTBDDDConstants::isValidCell(int lay, int wafer, int cell) const {
0720   // Calculate the position of the cell
0721   // Works for options HGCalHexagon/HGCalHexagonFull
0722   double x = hgpar_->waferPosX_[wafer];
0723   double y = hgpar_->waferPosY_[wafer];
0724   if (hgpar_->waferTypeT_[wafer] - 1 == HGCalTBParameters::HGCalFine) {
0725     x += hgpar_->cellFineX_[cell];
0726     y += hgpar_->cellFineY_[cell];
0727   } else {
0728     x += hgpar_->cellCoarseX_[cell];
0729     y += hgpar_->cellCoarseY_[cell];
0730   }
0731   double rr = sqrt(x * x + y * y);
0732   bool result = ((rr >= hgpar_->rMinLayHex_[lay - 1]) && (rr <= hgpar_->rMaxLayHex_[lay - 1]) &&
0733                  (wafer < static_cast<int>(hgpar_->waferPosX_.size())));
0734 #ifdef EDM_ML_DEBUG
0735   if (!result)
0736     edm::LogVerbatim("HGCalGeom") << "Input " << lay << ":" << wafer << ":" << cell << " Position " << x << ":" << y
0737                                   << ":" << rr << " Compare Limits " << hgpar_->rMinLayHex_[lay - 1] << ":"
0738                                   << hgpar_->rMaxLayHex_[lay - 1] << " Flag " << result;
0739 #endif
0740   return result;
0741 }
0742 
0743 int32_t HGCalTBDDDConstants::waferIndex(int wafer, int index) const {
0744   int layer = layerFromIndex(index, true);
0745   int indx = HGCalWaferIndex::waferIndex(layer, wafer, 0);
0746 #ifdef EDM_ML_DEBUG
0747   edm::LogVerbatim("HGCalGeom") << "WaferIndex for " << wafer << ":" << index << " (" << layer << ") " << indx;
0748 #endif
0749   return indx;
0750 }
0751 
0752 #include "FWCore/Utilities/interface/typelookup.h"
0753 
0754 TYPELOOKUP_DATA_REG(HGCalTBDDDConstants);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!