File indexing completed on 2024-09-07 04:35:04
0001 #include "Calibration/Tools/interface/EcalRingCalibrationTools.h"
0002 #include "DataFormats/DetId/interface/DetId.h"
0003
0004
0005 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0006 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0007 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0008
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010
0011 #include <iostream>
0012
0013
0014 std::atomic<bool> EcalRingCalibrationTools::isInitializedFromGeometry_(false);
0015 short EcalRingCalibrationTools::endcapRingIndex_[EEDetId::IX_MAX][EEDetId::IY_MAX];
0016 std::once_flag EcalRingCalibrationTools::once_;
0017
0018 constexpr short EcalRingCalibrationTools::N_RING_TOTAL;
0019 constexpr short EcalRingCalibrationTools::N_RING_BARREL;
0020 constexpr short EcalRingCalibrationTools::N_RING_ENDCAP;
0021
0022 constexpr short EcalRingCalibrationTools::N_MODULES_BARREL;
0023
0024 short EcalRingCalibrationTools::getRingIndex(DetId id) {
0025 if (id.det() != DetId::Ecal)
0026 return -1;
0027
0028 if (id.subdetId() == EcalBarrel) {
0029 if (EBDetId(id).ieta() < 0)
0030 return EBDetId(id).ieta() + 85;
0031 else
0032 return EBDetId(id).ieta() + 84;
0033 }
0034 if (id.subdetId() == EcalEndcap) {
0035
0036 if (not isInitializedFromGeometry_)
0037 throw std::logic_error(
0038 "EcalRingCalibrationTools::initializeFromGeometry Ecal Endcap geometry is not initialized");
0039
0040 EEDetId eid(id);
0041 short endcapRingIndex = endcapRingIndex_[eid.ix() - 1][eid.iy() - 1] + N_RING_BARREL;
0042 if (eid.zside() == 1)
0043 endcapRingIndex += N_RING_ENDCAP / 2;
0044 return endcapRingIndex;
0045 }
0046 return -1;
0047 }
0048
0049 short EcalRingCalibrationTools::getModuleIndex(DetId id) {
0050 if (id.det() != DetId::Ecal)
0051 return -1;
0052
0053 if (id.subdetId() == EcalBarrel) {
0054 short module = 4 * (EBDetId(id).ism() - 1) + EBDetId(id).im() - 1;
0055
0056
0057 return module;
0058 }
0059 if (id.subdetId() == EcalEndcap) {
0060 return -1;
0061 }
0062
0063 return -1;
0064 }
0065
0066 std::vector<DetId> EcalRingCalibrationTools::getDetIdsInRing(short etaIndex) {
0067 std::vector<DetId> ringIds;
0068 if (etaIndex < 0)
0069 return ringIds;
0070
0071 if (etaIndex < N_RING_BARREL) {
0072 int k = 0;
0073 if (etaIndex < 85)
0074 k = -85 + etaIndex;
0075 else
0076 k = etaIndex - 84;
0077
0078 for (int iphi = EBDetId::MIN_IPHI; iphi <= EBDetId::MAX_IPHI; ++iphi)
0079
0080 if (EBDetId::validDetId(k, iphi))
0081 ringIds.push_back(EBDetId(k, iphi));
0082 }
0083
0084 else if (etaIndex < N_RING_TOTAL) {
0085
0086 if (not isInitializedFromGeometry_)
0087 throw std::logic_error(
0088 "EcalRingCalibrationTools::initializeFromGeometry Ecal Endcap geometry is not initialized");
0089
0090 int zside = (etaIndex < N_RING_BARREL + (N_RING_ENDCAP / 2)) ? -1 : 1;
0091 short eeEtaIndex = (etaIndex - N_RING_BARREL) % (N_RING_ENDCAP / 2);
0092
0093 for (int ix = 0; ix < EEDetId::IX_MAX; ++ix)
0094 for (int iy = 0; iy < EEDetId::IY_MAX; ++iy)
0095 if (endcapRingIndex_[ix][iy] == eeEtaIndex)
0096 ringIds.push_back(EEDetId(ix + 1, iy + 1, zside));
0097 }
0098
0099 return ringIds;
0100 }
0101
0102 std::vector<DetId> EcalRingCalibrationTools::getDetIdsInECAL() {
0103 std::vector<DetId> ringIds;
0104
0105 for (int ieta = -EBDetId::MAX_IETA; ieta <= EBDetId::MAX_IETA; ++ieta)
0106 for (int iphi = EBDetId::MIN_IPHI; iphi <= EBDetId::MAX_IPHI; ++iphi)
0107 if (EBDetId::validDetId(ieta, iphi))
0108 ringIds.push_back(EBDetId(ieta, iphi));
0109
0110
0111 if (not isInitializedFromGeometry_)
0112 throw std::logic_error("EcalRingCalibrationTools::initializeFromGeometry Ecal Endcap geometry is not initialized");
0113
0114 for (int ix = 0; ix < EEDetId::IX_MAX; ++ix)
0115 for (int iy = 0; iy < EEDetId::IY_MAX; ++iy)
0116 for (int zside = -1; zside < 2; zside += 2)
0117 if (EEDetId::validDetId(ix + 1, iy + 1, zside))
0118 ringIds.push_back(EEDetId(ix + 1, iy + 1, zside));
0119
0120
0121
0122 return ringIds;
0123 }
0124
0125 std::vector<DetId> EcalRingCalibrationTools::getDetIdsInModule(short moduleIndex) {
0126 std::vector<DetId> ringIds;
0127 if (moduleIndex < 0)
0128 return ringIds;
0129
0130 short moduleBound[5] = {1, 26, 46, 66, 86};
0131 if (moduleIndex < EcalRingCalibrationTools::N_MODULES_BARREL) {
0132
0133 short sm, moduleInSm, zsm;
0134
0135 short minModuleiphi, maxModuleiphi, minModuleieta = 360, maxModuleieta = 0;
0136
0137
0138 sm = moduleIndex / 4 + 1;
0139
0140
0141
0142
0143 moduleInSm = moduleIndex % 4;
0144
0145
0146
0147 if (moduleIndex > 71)
0148 zsm = -1;
0149 else
0150 zsm = 1;
0151
0152 minModuleiphi = ((sm - 1) % 18 + 1) * 20 - 19;
0153 maxModuleiphi = ((sm - 1) % 18 + 1) * 20;
0154
0155 if (zsm == 1) {
0156 minModuleieta = moduleBound[moduleInSm];
0157 maxModuleieta = moduleBound[moduleInSm + 1] - 1;
0158 } else if (zsm == -1) {
0159 minModuleieta = -moduleBound[moduleInSm + 1] + 1;
0160 maxModuleieta = -moduleBound[moduleInSm];
0161 }
0162
0163
0164 std::cout << "Called moduleIndex " << moduleIndex << std::endl;
0165 std::cout << "minModuleieta " << minModuleieta << " maxModuleieta " << maxModuleieta << " minModuleiphi "
0166 << minModuleiphi << " maxModuleiphi " << maxModuleiphi << std::endl;
0167
0168 for (int ieta = minModuleieta; ieta <= maxModuleieta; ++ieta) {
0169 for (int iphi = minModuleiphi; iphi <= maxModuleiphi; ++iphi) {
0170 ringIds.push_back(EBDetId(ieta, iphi));
0171
0172
0173
0174 }
0175 }
0176 }
0177
0178 return ringIds;
0179 }
0180
0181 void EcalRingCalibrationTools::setCaloGeometry(const CaloGeometry* geometry) {
0182 std::call_once(once_, EcalRingCalibrationTools::initializeFromGeometry, geometry);
0183 }
0184
0185 void EcalRingCalibrationTools::initializeFromGeometry(CaloGeometry const* geometry) {
0186 if (not geometry)
0187 throw std::invalid_argument("EcalRingCalibrationTools::initializeFromGeometry called with a nullptr argument");
0188
0189 float cellPosEta[EEDetId::IX_MAX][EEDetId::IY_MAX];
0190 for (int ix = 0; ix < EEDetId::IX_MAX; ++ix)
0191 for (int iy = 0; iy < EEDetId::IY_MAX; ++iy) {
0192 cellPosEta[ix][iy] = -1.;
0193 endcapRingIndex_[ix][iy] = -9;
0194 }
0195
0196 CaloSubdetectorGeometry const* endcapGeometry = geometry->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0197 if (not endcapGeometry)
0198 throw std::logic_error("EcalRingCalibrationTools::initializeFromGeometry Ecal Endcap geometry not found");
0199
0200 std::vector<DetId> const& endcapCells = geometry->getValidDetIds(DetId::Ecal, EcalEndcap);
0201
0202 for (std::vector<DetId>::const_iterator endcapIt = endcapCells.begin(); endcapIt != endcapCells.end(); ++endcapIt) {
0203 EEDetId ee(*endcapIt);
0204 if (ee.zside() == -1)
0205 continue;
0206 auto cellGeometry = endcapGeometry->getGeometry(*endcapIt);
0207 int ics = ee.ix() - 1;
0208 int ips = ee.iy() - 1;
0209 cellPosEta[ics][ips] = fabs(cellGeometry->getPosition().eta());
0210
0211 }
0212
0213 float eta_ring[N_RING_ENDCAP / 2];
0214 for (int ring = 0; ring < N_RING_ENDCAP / 2; ++ring)
0215 eta_ring[ring] = cellPosEta[ring][50];
0216
0217 double etaBoundary[N_RING_ENDCAP / 2 + 1];
0218 etaBoundary[0] = 1.47;
0219 etaBoundary[N_RING_ENDCAP / 2] = 4.0;
0220
0221 for (int ring = 1; ring < N_RING_ENDCAP / 2; ++ring)
0222 etaBoundary[ring] = (eta_ring[ring] + eta_ring[ring - 1]) / 2.;
0223
0224 for (int ring = 0; ring < N_RING_ENDCAP / 2; ++ring) {
0225
0226 for (int ix = 0; ix < EEDetId::IX_MAX; ix++)
0227 for (int iy = 0; iy < EEDetId::IY_MAX; iy++)
0228 if (cellPosEta[ix][iy] > etaBoundary[ring] && cellPosEta[ix][iy] < etaBoundary[ring + 1]) {
0229 endcapRingIndex_[ix][iy] = ring;
0230
0231 }
0232 }
0233
0234 std::vector<DetId> const& barrelCells = geometry->getValidDetIds(DetId::Ecal, EcalBarrel);
0235
0236 for (std::vector<DetId>::const_iterator barrelIt = barrelCells.begin(); barrelIt != barrelCells.end(); ++barrelIt) {
0237 EBDetId eb(*barrelIt);
0238 }
0239
0240
0241
0242 isInitializedFromGeometry_ = true;
0243 }