File indexing completed on 2024-04-06 12:14:51
0001 #include "Geometry/HcalTowerAlgo/interface/HcalFlexiHardcodeGeometryLoader.h"
0002 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0003 #include "Geometry/CaloGeometry/interface/IdealObliquePrism.h"
0004 #include "Geometry/CaloGeometry/interface/IdealZPrism.h"
0005 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0006 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0007 #include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009
0010 #include <algorithm>
0011 #include <sstream>
0012 #include <vector>
0013
0014 using CCGFloat = CaloCellGeometry::CCGFloat;
0015
0016
0017
0018
0019 HcalFlexiHardcodeGeometryLoader::HcalFlexiHardcodeGeometryLoader() {
0020 MAX_HCAL_PHI = 72;
0021 DEGREE2RAD = M_PI / 180.;
0022 }
0023
0024 CaloSubdetectorGeometry* HcalFlexiHardcodeGeometryLoader::load(const HcalTopology& fTopology,
0025 const HcalDDDRecConstants& hcons) {
0026 HcalGeometry* hcalGeometry = new HcalGeometry(fTopology);
0027 if (nullptr == hcalGeometry->cornersMgr())
0028 hcalGeometry->allocateCorners(fTopology.ncells() + fTopology.getHFSize());
0029 if (nullptr == hcalGeometry->parMgr())
0030 hcalGeometry->allocatePar(hcalGeometry->numberOfShapes(), HcalGeometry::k_NumberOfParametersPerShape);
0031 isBH_ = hcons.isBH();
0032 #ifdef EDM_ML_DEBUG
0033 edm::LogVerbatim("HCalGeom") << "FlexiGeometryLoader initialize with ncells " << fTopology.ncells() << " and shapes "
0034 << hcalGeometry->numberOfShapes() << ":" << HcalGeometry::k_NumberOfParametersPerShape
0035 << " with BH Flag " << isBH_;
0036 #endif
0037 if (fTopology.mode() == HcalTopologyMode::H2) {
0038 fillHBHO(hcalGeometry, makeHBCells(hcons), true);
0039 fillHBHO(hcalGeometry, makeHOCells(), false);
0040 fillHE(hcalGeometry, makeHECells_H2());
0041 } else {
0042 fillHBHO(hcalGeometry, makeHBCells(hcons), true);
0043 fillHBHO(hcalGeometry, makeHOCells(), false);
0044 fillHF(hcalGeometry, makeHFCells(hcons));
0045 fillHE(hcalGeometry, makeHECells(hcons));
0046 }
0047
0048 hcalGeometry->sortValidIds();
0049 return hcalGeometry;
0050 }
0051
0052
0053 std::vector<HcalFlexiHardcodeGeometryLoader::HBHOCellParameters> HcalFlexiHardcodeGeometryLoader::makeHBCells(
0054 const HcalDDDRecConstants& hcons) {
0055 std::vector<HcalFlexiHardcodeGeometryLoader::HBHOCellParameters> result;
0056 std::vector<std::pair<double, double> > gconsHB = hcons.getConstHBHE(0);
0057 std::vector<HcalDDDRecConstants::HcalEtaBin> etabins = hcons.getEtaBins(0);
0058
0059 #ifdef EDM_ML_DEBUG
0060 edm::LogVerbatim("HCalGeom") << "FlexiGeometryLoader called for " << etabins.size() << " Eta Bins";
0061 for (unsigned int k = 0; k < gconsHB.size(); ++k) {
0062 edm::LogVerbatim("HCalGeom") << "gconsHB[" << k << "] = " << gconsHB[k].first << " +- " << gconsHB[k].second;
0063 }
0064 #endif
0065 for (auto& etabin : etabins) {
0066 int iring = (etabin.zside >= 0) ? etabin.ieta : -etabin.ieta;
0067 int depth = etabin.depthStart;
0068 double dphi =
0069 (etabin.phis.size() > 1) ? (etabin.phis[1].second - etabin.phis[0].second) : ((2.0 * M_PI) / MAX_HCAL_PHI);
0070 for (unsigned int k = 0; k < etabin.layer.size(); ++k) {
0071 int layf = etabin.layer[k].first - 1;
0072 int layl = etabin.layer[k].second - 1;
0073 double rmin = gconsHB[layf].first - gconsHB[layf].second;
0074 double rmax = gconsHB[layl].first + gconsHB[layl].second;
0075 for (unsigned int j = 0; j < etabin.phis.size(); ++j) {
0076 #ifdef EDM_ML_DEBUG
0077 edm::LogVerbatim("HCalGeom") << "HBRing " << iring << " eta " << etabin.etaMin << ":" << etabin.etaMax
0078 << " depth " << depth << " R " << rmin << ":" << rmax << " Phi "
0079 << etabin.phis[j].first << ":" << etabin.phis[j].second << ":" << dphi << " layer["
0080 << k << "]: " << etabin.layer[k].first - 1 << ":" << etabin.layer[k].second;
0081 #endif
0082 result.emplace_back(HcalFlexiHardcodeGeometryLoader::HBHOCellParameters(
0083 iring, depth, etabin.phis[j].first, etabin.phis[j].second, dphi, rmin, rmax, etabin.etaMin, etabin.etaMax));
0084 }
0085 depth++;
0086 }
0087 }
0088 return result;
0089 }
0090
0091
0092 std::vector<HcalFlexiHardcodeGeometryLoader::HBHOCellParameters> HcalFlexiHardcodeGeometryLoader::makeHOCells() {
0093 const double HORMIN0 = 390.0;
0094 const double HORMIN1 = 412.6;
0095 const double HORMAX = 413.6;
0096 const int nCells = 15;
0097 const int nPhi = 72;
0098 const double etamin[nCells] = {
0099 0.000, 0.087, 0.174, 0.261, 0.3395, 0.435, 0.522, 0.609, 0.696, 0.783, 0.873, 0.957, 1.044, 1.131, 1.218};
0100 const double etamax[nCells] = {
0101 0.087, 0.174, 0.261, 0.3075, 0.435, 0.522, 0.609, 0.696, 0.783, 0.8494, 0.957, 1.044, 1.131, 1.218, 1.305};
0102 std::vector<HcalFlexiHardcodeGeometryLoader::HBHOCellParameters> result;
0103 result.reserve(nCells * nPhi);
0104 double dphi = ((2.0 * M_PI) / nPhi);
0105 for (int i = 0; i < nCells; ++i) {
0106 double rmin = ((i < 4) ? HORMIN0 : HORMIN1);
0107 for (int iside = -1; iside <= 1; iside += 2) {
0108 for (int j = 0; j < nPhi; ++j) {
0109 double phi = (j + 0.5) * dphi;
0110
0111 result.emplace_back(HcalFlexiHardcodeGeometryLoader::HBHOCellParameters(
0112 iside * (i + 1), 4, j + 1, phi, dphi, rmin, HORMAX, etamin[i], etamax[i]));
0113 }
0114 }
0115 }
0116 return result;
0117 }
0118
0119
0120
0121
0122 void HcalFlexiHardcodeGeometryLoader::fillHBHO(
0123 HcalGeometry* fGeometry, const std::vector<HcalFlexiHardcodeGeometryLoader::HBHOCellParameters>& fCells, bool fHB) {
0124 fGeometry->increaseReserve(fCells.size());
0125 for (const auto& param : fCells) {
0126 HcalDetId hid(fHB ? HcalBarrel : HcalOuter, param.ieta, param.iphi, param.depth);
0127 float phiCenter = param.phi;
0128 float etaCenter = 0.5 * (param.etaMin + param.etaMax);
0129 float x = param.rMin * cos(phiCenter);
0130 float y = param.rMin * sin(phiCenter);
0131 float z = (param.ieta < 0) ? -(param.rMin * sinh(etaCenter)) : (param.rMin * sinh(etaCenter));
0132
0133 GlobalPoint refPoint(x, y, z);
0134 std::vector<CCGFloat> cellParams;
0135 cellParams.reserve(5);
0136 cellParams.emplace_back(0.5 * (param.etaMax - param.etaMin));
0137 cellParams.emplace_back(0.5 * param.dphi);
0138 cellParams.emplace_back(0.5 * (param.rMax - param.rMin) * cosh(etaCenter));
0139 cellParams.emplace_back(fabs(refPoint.eta()));
0140 cellParams.emplace_back(fabs(refPoint.z()));
0141 #ifdef EDM_ML_DEBUG
0142 edm::LogVerbatim("HCalGeom") << "HcalFlexiHardcodeGeometryLoader::fillHBHO-> " << hid << " " << hid.rawId() << " "
0143 << std::hex << hid.rawId() << std::dec << " " << hid << " " << refPoint << '/'
0144 << cellParams[0] << '/' << cellParams[1] << '/' << cellParams[2];
0145 #endif
0146 fGeometry->newCellFast(refPoint,
0147 refPoint,
0148 refPoint,
0149 CaloCellGeometry::getParmPtr(cellParams, fGeometry->parMgr(), fGeometry->parVecVec()),
0150 hid);
0151 }
0152 }
0153
0154
0155 std::vector<HcalFlexiHardcodeGeometryLoader::HECellParameters> HcalFlexiHardcodeGeometryLoader::makeHECells(
0156 const HcalDDDRecConstants& hcons) {
0157 std::vector<HcalFlexiHardcodeGeometryLoader::HECellParameters> result;
0158 std::vector<std::pair<double, double> > gconsHE = hcons.getConstHBHE(1);
0159 #ifdef EDM_ML_DEBUG
0160 edm::LogVerbatim("HCalGeom") << "HcalFlexiHardcodeGeometryLoader:HE with " << gconsHE.size() << " cells";
0161 #endif
0162 if (!gconsHE.empty()) {
0163 std::vector<HcalDDDRecConstants::HcalEtaBin> etabins = hcons.getEtaBins(1);
0164
0165 #ifdef EDM_ML_DEBUG
0166 edm::LogVerbatim("HCalGeom") << "FlexiGeometryLoader called for HE with " << etabins.size() << " Eta Bins and "
0167 << gconsHE.size() << " depths";
0168 std::ostringstream st1;
0169 for (unsigned int i = 0; i < gconsHE.size(); ++i)
0170 st1 << " Depth[" << i << "] = " << gconsHE[i].first << " +- " << gconsHE[i].second;
0171 edm::LogVerbatim("HCalGeom") << st1.str();
0172 #endif
0173 for (auto& etabin : etabins) {
0174 int iring = (etabin.zside >= 0) ? etabin.ieta : -etabin.ieta;
0175 int depth = etabin.depthStart;
0176 double dphi =
0177 (etabin.phis.size() > 1) ? (etabin.phis[1].second - etabin.phis[0].second) : ((2.0 * M_PI) / MAX_HCAL_PHI);
0178 #ifdef EDM_ML_DEBUG
0179 edm::LogVerbatim("HCalGeom") << "FlexiGeometryLoader::Ring " << iring << " nphi " << etabin.phis.size()
0180 << " dstart " << depth << " dphi " << dphi << " layers " << etabin.layer.size();
0181 std::ostringstream st2;
0182 for (unsigned int j = 0; j < etabin.phis.size(); ++j)
0183 st2 << " [" << j << "] " << etabin.phis[j].first << ":" << etabin.phis[j].second;
0184 edm::LogVerbatim("HCalGeom") << st2.str();
0185 #endif
0186 for (unsigned int k = 0; k < etabin.layer.size(); ++k) {
0187 int layf = etabin.layer[k].first - 1;
0188 int layl = etabin.layer[k].second - 1;
0189 double zmin = gconsHE[layf].first - gconsHE[layf].second;
0190 double zmax = gconsHE[layl].first + gconsHE[layl].second;
0191 if (zmin < 1.0) {
0192 for (int k2 = layf; k2 <= layl; ++k2) {
0193 if (gconsHE[k2].first > 10) {
0194 zmin = gconsHE[k2].first - gconsHE[k2].second;
0195 break;
0196 }
0197 }
0198 }
0199 if (zmin >= zmax)
0200 zmax = zmin + 10.;
0201 for (unsigned int j = 0; j < etabin.phis.size(); ++j) {
0202 #ifdef EDM_ML_DEBUG
0203 edm::LogVerbatim("HCalGeom") << "HERing " << iring << " eta " << etabin.etaMin << ":" << etabin.etaMax
0204 << " depth " << depth << " Z " << zmin << ":" << zmax
0205 << " Phi :" << etabin.phis[j].first << ":" << etabin.phis[j].second << ":"
0206 << dphi << " layer[" << k << "]: " << etabin.layer[k].first - 1 << ":"
0207 << etabin.layer[k].second - 1;
0208 #endif
0209 result.emplace_back(HcalFlexiHardcodeGeometryLoader::HECellParameters(
0210 iring, depth, etabin.phis[j].first, etabin.phis[j].second, dphi, zmin, zmax, etabin.etaMin, etabin.etaMax));
0211 }
0212 depth++;
0213 }
0214 }
0215 }
0216 return result;
0217 }
0218
0219
0220 std::vector<HcalFlexiHardcodeGeometryLoader::HECellParameters> HcalFlexiHardcodeGeometryLoader::makeHECells_H2() {
0221 const double HEZMIN_H2 = 400.715;
0222 const double HEZMID_H2 = 436.285;
0223 const double HEZMAX_H2 = 541.885;
0224 const int nEtas = 10;
0225 const int nDepth[nEtas] = {1, 2, 2, 2, 2, 2, 2, 2, 3, 3};
0226 const int dStart[nEtas] = {3, 1, 1, 1, 1, 1, 1, 1, 1, 1};
0227 const int nPhis[nEtas] = {8, 8, 8, 8, 8, 8, 4, 4, 4, 4};
0228 const double etas[nEtas + 1] = {1.305, 1.373, 1.444, 1.521, 1.603, 1.693, 1.790, 1.880, 1.980, 2.090, 2.210};
0229 const double zval[4 * nEtas] = {
0230 409.885, 462.685, 0., 0., HEZMIN_H2, 427.485, 506.685, 0.0, HEZMIN_H2, HEZMID_H2,
0231 524.285, 0., HEZMIN_H2, HEZMID_H2, HEZMAX_H2, 0., HEZMIN_H2, HEZMID_H2, HEZMAX_H2, 0.,
0232 HEZMIN_H2, HEZMID_H2, HEZMAX_H2, 0., HEZMIN_H2, HEZMID_H2, HEZMAX_H2, 0., HEZMIN_H2, HEZMID_H2,
0233 HEZMAX_H2, 0., HEZMIN_H2, 418.685, HEZMID_H2, HEZMAX_H2, HEZMIN_H2, 418.685, HEZMID_H2, HEZMAX_H2};
0234 std::vector<HcalFlexiHardcodeGeometryLoader::HECellParameters> result;
0235
0236 for (int i = 0; i < nEtas; ++i) {
0237 int ieta = 16 + i;
0238 for (int k = 0; k < nDepth[i]; ++k) {
0239 int depth = dStart[i] + k;
0240 for (int j = 0; j < nPhis[i]; ++j) {
0241 int iphi = (nPhis[i] == 8) ? (j + 1) : (2 * j + 1);
0242 double dphi = (40.0 * DEGREE2RAD) / nPhis[i];
0243 double phi0 = (j + 0.5) * dphi;
0244
0245 result.emplace_back(HcalFlexiHardcodeGeometryLoader::HECellParameters(
0246 ieta, depth, iphi, phi0, dphi, zval[4 * i + k + 1], zval[4 * i + k + 2], etas[i], etas[i + 1]));
0247 }
0248 }
0249 }
0250 return result;
0251 }
0252
0253
0254 std::vector<HcalFlexiHardcodeGeometryLoader::HFCellParameters> HcalFlexiHardcodeGeometryLoader::makeHFCells(
0255 const HcalDDDRecConstants& hcons) {
0256 const float HFZMIN1 = 1115.;
0257 const float HFZMIN2 = 1137.;
0258 const float HFZMAX = 1280.1;
0259 std::vector<HcalDDDRecConstants::HFCellParameters> cells = hcons.getHFCellParameters();
0260 unsigned int nCells = cells.size();
0261 std::vector<HcalFlexiHardcodeGeometryLoader::HFCellParameters> result;
0262 result.reserve(nCells);
0263 for (unsigned int i = 0; i < nCells; ++i) {
0264 HcalFlexiHardcodeGeometryLoader::HFCellParameters cell1(cells[i].ieta,
0265 cells[i].depth,
0266 cells[i].firstPhi,
0267 cells[i].stepPhi,
0268 cells[i].nPhi,
0269 5 * cells[i].stepPhi,
0270 HFZMIN1,
0271 HFZMAX,
0272 cells[i].rMin,
0273 cells[i].rMax);
0274 result.emplace_back(cell1);
0275 HcalFlexiHardcodeGeometryLoader::HFCellParameters cell2(cells[i].ieta,
0276 1 + cells[i].depth,
0277 cells[i].firstPhi,
0278 cells[i].stepPhi,
0279 cells[i].nPhi,
0280 5 * cells[i].stepPhi,
0281 HFZMIN2,
0282 HFZMAX,
0283 cells[i].rMin,
0284 cells[i].rMax);
0285 result.emplace_back(cell2);
0286 }
0287 return result;
0288 }
0289
0290 void HcalFlexiHardcodeGeometryLoader::fillHE(
0291 HcalGeometry* fGeometry, const std::vector<HcalFlexiHardcodeGeometryLoader::HECellParameters>& fCells) {
0292 fGeometry->increaseReserve(fCells.size());
0293 for (const auto& param : fCells) {
0294 HcalDetId hid(HcalEndcap, param.ieta, param.iphi, param.depth);
0295 float phiCenter = param.phi;
0296 float etaCenter = 0.5 * (param.etaMin + param.etaMax);
0297 int iside = (param.ieta >= 0) ? 1 : -1;
0298 float z = (isBH_) ? (iside * 0.5 * (param.zMin + param.zMax)) : (iside * param.zMin);
0299 float perp = std::abs(z) / sinh(etaCenter);
0300 float x = perp * cos(phiCenter);
0301 float y = perp * sin(phiCenter);
0302
0303 GlobalPoint refPoint(x, y, z);
0304 std::vector<CCGFloat> cellParams;
0305 cellParams.reserve(5);
0306 cellParams.emplace_back(0.5 * (param.etaMax - param.etaMin));
0307 cellParams.emplace_back(0.5 * param.dphi);
0308 cellParams.emplace_back(-0.5 * (param.zMax - param.zMin) / tanh(etaCenter));
0309 cellParams.emplace_back(fabs(refPoint.eta()));
0310 cellParams.emplace_back(fabs(refPoint.z()));
0311 #ifdef EDM_ML_DEBUG
0312 edm::LogVerbatim("HCalGeom") << "HcalFlexiHardcodeGeometryLoader::fillHE-> " << hid << " " << hid.rawId() << " "
0313 << std::hex << hid.rawId() << std::dec << " " << hid << refPoint << '/'
0314 << cellParams[0] << '/' << cellParams[1] << '/' << cellParams[2];
0315 #endif
0316 fGeometry->newCellFast(refPoint,
0317 refPoint,
0318 refPoint,
0319 CaloCellGeometry::getParmPtr(cellParams, fGeometry->parMgr(), fGeometry->parVecVec()),
0320 hid);
0321 }
0322 }
0323
0324 void HcalFlexiHardcodeGeometryLoader::fillHF(
0325 HcalGeometry* fGeometry, const std::vector<HcalFlexiHardcodeGeometryLoader::HFCellParameters>& fCells) {
0326 fGeometry->increaseReserve(fCells.size());
0327 for (const auto& param : fCells) {
0328 for (int kPhi = 0; kPhi < param.nPhi; ++kPhi) {
0329 int iPhi = param.phiFirst + kPhi * param.phiStep;
0330 HcalDetId hid(HcalForward, param.eta, iPhi, param.depth);
0331
0332 float phiCenter = ((iPhi - 1) * 360. / MAX_HCAL_PHI + 0.5 * param.dphi) * DEGREE2RAD;
0333 GlobalPoint inner(param.rMin, 0, param.zMin);
0334 GlobalPoint outer(param.rMax, 0, param.zMin);
0335 float iEta = inner.eta();
0336 float oEta = outer.eta();
0337 float etaCenter = 0.5 * (iEta + oEta);
0338
0339 float perp = param.zMin / sinh(etaCenter);
0340 float x = perp * cos(phiCenter);
0341 float y = perp * sin(phiCenter);
0342 float z = (param.eta > 0) ? param.zMin : -param.zMin;
0343
0344 GlobalPoint refPoint(x, y, z);
0345 std::vector<CCGFloat> cellParams;
0346 cellParams.reserve(5);
0347 cellParams.emplace_back(0.5 * (iEta - oEta));
0348 cellParams.emplace_back(0.5 * param.dphi * DEGREE2RAD);
0349 cellParams.emplace_back(0.5 * (param.zMax - param.zMin));
0350 cellParams.emplace_back(fabs(refPoint.eta()));
0351 cellParams.emplace_back(fabs(refPoint.z()));
0352 #ifdef EDM_ML_DEBUG
0353 edm::LogVerbatim("HCalGeom") << "HcalFlexiHardcodeGeometryLoader::fillHF-> " << hid << " " << hid.rawId() << " "
0354 << std::hex << hid.rawId() << std::dec << " " << hid << " " << refPoint << '/'
0355 << cellParams[0] << '/' << cellParams[1] << '/' << cellParams[2];
0356 #endif
0357 fGeometry->newCellFast(refPoint,
0358 refPoint,
0359 refPoint,
0360 CaloCellGeometry::getParmPtr(cellParams, fGeometry->parMgr(), fGeometry->parVecVec()),
0361 hid);
0362 }
0363 }
0364 }