File indexing completed on 2023-10-25 09:45:58
0001 #include "FastSimulation/CaloGeometryTools/interface/CaloGeometryHelper.h"
0002
0003
0004
0005
0006 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0007 #include "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h"
0008 #include "Geometry/EcalAlgo/interface/EcalEndcapGeometry.h"
0009 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0010 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0011
0012 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0013 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0014 #include "DataFormats/EcalDetId/interface/ESDetId.h"
0015 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0016 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0017 #include "Geometry/EcalAlgo/interface/EcalPreshowerGeometry.h"
0018
0019 #include "FastSimulation/CaloGeometryTools/interface/DistanceToCell.h"
0020 #include "FastSimulation/CaloGeometryTools/interface/Crystal.h"
0021
0022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0023
0024 #include <algorithm>
0025
0026 CaloGeometryHelper::CaloGeometryHelper() : Calorimeter() {
0027 neighbourmapcalculated_ = false;
0028 psLayer1Z_ = 303;
0029 psLayer2Z_ = 307;
0030 }
0031
0032 CaloGeometryHelper::CaloGeometryHelper(const edm::ParameterSet& fastCalo) : Calorimeter(fastCalo) {
0033
0034 psLayer1Z_ = 303;
0035 psLayer2Z_ = 307;
0036 }
0037
0038 void CaloGeometryHelper::initialize(double bField) {
0039 buildCrystalArray();
0040 buildNeighbourArray();
0041 bfield_ = bField;
0042 preshowerPresent_ = (getEcalPreshowerGeometry() != nullptr);
0043
0044 if (preshowerPresent_) {
0045 ESDetId cps1(getEcalPreshowerGeometry()->getClosestCellInPlane(GlobalPoint(80., 80., 303.), 1));
0046 psLayer1Z_ = getEcalPreshowerGeometry()->getGeometry(cps1)->getPosition().z();
0047 ESDetId cps2(getEcalPreshowerGeometry()->getClosestCellInPlane(GlobalPoint(80., 80., 307.), 2));
0048 psLayer2Z_ = getEcalPreshowerGeometry()->getGeometry(cps2)->getPosition().z();
0049 LogDebug("CaloGeometryTools") << " Preshower layer positions " << psLayer1Z_ << " " << psLayer2Z_ << std::endl;
0050 } else
0051 LogDebug("CaloGeometryTools") << " No preshower present" << std::endl;
0052
0053
0054 }
0055
0056 CaloGeometryHelper::~CaloGeometryHelper() { ; }
0057
0058 DetId CaloGeometryHelper::getClosestCell(const XYZPoint& point, bool ecal, bool central) const {
0059 DetId result;
0060 if (ecal) {
0061 if (central) {
0062
0063 result = EcalBarrelGeometry_->getClosestCell(GlobalPoint(point.X(), point.Y(), point.Z()));
0064 #ifdef DEBUGGCC
0065 if (result.null())
0066 return result;
0067 GlobalPoint ip = GlobalPoint(point.X(), point.Y(), point.Z());
0068 GlobalPoint cc = EcalBarrelGeometry_->getGeometry(result)->getPosition();
0069 float deltaeta2 = ip.eta() - cc.eta();
0070 deltaeta2 *= deltaeta2;
0071 float deltaphi2 = acos(cos(ip.phi() - cc.phi()));
0072 deltaphi2 *= deltaphi2;
0073 Histos::instance()->fill("h100", point.eta(), sqrt(deltaeta2 + deltaphi2));
0074 #endif
0075 } else {
0076 result = EcalEndcapGeometry_->getClosestCell(GlobalPoint(point.X(), point.Y(), point.Z()));
0077 #ifdef DEBUGGCC
0078 if (result.null()) {
0079 return result;
0080 }
0081 GlobalPoint ip = GlobalPoint(point.X(), point.Y(), point.Z());
0082 GlobalPoint cc = EcalEndcapGeometry_->getGeometry(result)->getPosition();
0083 Histos::instance()->fill("h110", point.eta(), (ip - cc).perp());
0084 #endif
0085 }
0086 } else {
0087 result = static_cast<const HcalGeometry*>(HcalGeometry_)
0088 ->getClosestCell(GlobalPoint(point.X(), point.Y(), point.Z()), true);
0089 HcalDetId myDetId(result);
0090
0091
0092 if (myDetId.subdetId() == HcalForward) {
0093 int mylayer;
0094 if (fabs(point.Z()) > 1132.) {
0095 mylayer = 2;
0096 } else {
0097 mylayer = 1;
0098 }
0099 HcalDetId myDetId2((HcalSubdetector)myDetId.subdetId(), myDetId.ieta(), myDetId.iphi(), mylayer);
0100 result = myDetId2;
0101
0102 }
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131 #ifdef DEBUGGCC
0132 if (result.null()) {
0133 return result;
0134 }
0135 GlobalPoint ip = GlobalPoint(point.x(), point.y(), point.z());
0136 GlobalPoint cc = HcalGeometry_->getPosition(result);
0137 float deltaeta2 = ip.eta() - cc.eta();
0138 deltaeta2 *= deltaeta2;
0139 float deltaphi2 = acos(cos(ip.phi() - cc.phi()));
0140 deltaphi2 *= deltaphi2;
0141
0142 Histos::instance()->fill("h120", point.eta(), sqrt(deltaeta2 + deltaphi2));
0143 #endif
0144 }
0145 return result;
0146 }
0147
0148 void CaloGeometryHelper::getWindow(const DetId& pivot, int s1, int s2, std::vector<DetId>& vec) const {
0149
0150
0151
0152 vec = getEcalTopology(pivot.subdetId())->getWindow(pivot, s1, s2);
0153 DistanceToCell distance(getEcalGeometry(pivot.subdetId()), pivot);
0154 sort(vec.begin(), vec.end(), distance);
0155 }
0156
0157 void CaloGeometryHelper::buildCrystal(const DetId& cell, Crystal& xtal) const {
0158 if (cell.subdetId() == EcalBarrel) {
0159 xtal = Crystal(cell, &barrelCrystals_[EBDetId(cell).hashedIndex()]);
0160 return;
0161 }
0162 if (cell.subdetId() == EcalEndcap) {
0163 xtal = Crystal(cell, &endcapCrystals_[EEDetId(cell).hashedIndex()]);
0164 return;
0165 }
0166 }
0167
0168
0169 void CaloGeometryHelper::buildNeighbourArray() {
0170 static const CaloDirection orderedDir[8] = {SOUTHWEST, SOUTH, SOUTHEAST, WEST, EAST, NORTHWEST, NORTH, NORTHEAST};
0171
0172 const unsigned nbarrel = EBDetId::kSizeForDenseIndexing;
0173
0174 barrelNeighbours_.resize(nbarrel);
0175
0176
0177
0178 const std::vector<DetId>& vec(EcalBarrelGeometry_->getValidDetIds(DetId::Ecal, EcalBarrel));
0179 unsigned size = vec.size();
0180 for (unsigned ic = 0; ic < size; ++ic) {
0181
0182 std::vector<DetId> neighbours(EcalBarrelTopology_->getWindow(vec[ic], 3, 3));
0183
0184 unsigned nneighbours = neighbours.size();
0185
0186 unsigned hashedindex = EBDetId(vec[ic]).hashedIndex();
0187 if (hashedindex >= nbarrel) {
0188 LogDebug("CaloGeometryTools") << " Array overflow " << std::endl;
0189 }
0190
0191
0192
0193
0194
0195
0196 if (nneighbours == 9) {
0197
0198 unsigned int nn = 0;
0199 for (unsigned in = 0; in < nneighbours; ++in) {
0200
0201 if (neighbours[in] != vec[ic]) {
0202 barrelNeighbours_[hashedindex][nn] = (neighbours[in]);
0203 nn++;
0204
0205 }
0206 }
0207 } else {
0208 DetId central(vec[ic]);
0209
0210 for (unsigned idir = 0; idir < 8; ++idir) {
0211 DetId testid = central;
0212 bool status = move(testid, orderedDir[idir], false);
0213 if (status)
0214 barrelNeighbours_[hashedindex][idir] = testid;
0215 }
0216 }
0217 }
0218
0219
0220
0221
0222
0223
0224 const std::vector<DetId>& vece(EcalEndcapGeometry_->getValidDetIds(DetId::Ecal, EcalEndcap));
0225 size = vece.size();
0226
0227
0228 const unsigned nendcap = EEDetId::kSizeForDenseIndexing;
0229
0230 endcapNeighbours_.resize(nendcap);
0231 for (unsigned ic = 0; ic < size; ++ic) {
0232
0233 std::vector<DetId> neighbours(EcalEndcapTopology_->getWindow(vece[ic], 3, 3));
0234 unsigned nneighbours = neighbours.size();
0235
0236 unsigned hashedindex = EEDetId(vece[ic]).hashedIndex();
0237
0238 if (hashedindex >= nendcap) {
0239 LogDebug("CaloGeometryTools") << " Array overflow " << std::endl;
0240 }
0241
0242 if (nneighbours == 9) {
0243
0244 unsigned int nn = 0;
0245 for (unsigned in = 0; in < nneighbours; ++in) {
0246
0247 if (neighbours[in] != vece[ic]) {
0248 endcapNeighbours_[hashedindex][nn] = (neighbours[in]);
0249 nn++;
0250 }
0251 }
0252 } else {
0253 DetId central(vece[ic]);
0254
0255 for (unsigned idir = 0; idir < 8; ++idir) {
0256 DetId testid = central;
0257 bool status = move(testid, orderedDir[idir], false);
0258 if (status)
0259 endcapNeighbours_[hashedindex][idir] = testid;
0260 }
0261 }
0262 }
0263
0264 neighbourmapcalculated_ = true;
0265 }
0266
0267 const CaloGeometryHelper::NeiVect& CaloGeometryHelper::getNeighbours(const DetId& detid) const {
0268 return (detid.subdetId() == EcalBarrel) ? barrelNeighbours_[EBDetId(detid).hashedIndex()]
0269 : endcapNeighbours_[EEDetId(detid).hashedIndex()];
0270 }
0271
0272 bool CaloGeometryHelper::move(DetId& cell, const CaloDirection& dir, bool fast) const {
0273 DetId originalcell = cell;
0274 if (dir == NONE || cell == DetId(0))
0275 return false;
0276
0277
0278
0279
0280 static const int calodirections[9] = {-1, 1, 2, 0, 4, 3, 7, 5, 6};
0281
0282 if (fast && neighbourmapcalculated_) {
0283 DetId result = (originalcell.subdetId() == EcalBarrel)
0284 ? barrelNeighbours_[EBDetId(originalcell).hashedIndex()][calodirections[dir]]
0285 : endcapNeighbours_[EEDetId(originalcell).hashedIndex()][calodirections[dir]];
0286 bool status = !result.null();
0287 cell = result;
0288 return status;
0289 }
0290
0291 if (dir == NORTH || dir == SOUTH || dir == EAST || dir == WEST) {
0292 return simplemove(cell, dir);
0293 } else {
0294 if (dir == NORTHEAST || dir == NORTHWEST || dir == SOUTHEAST || dir == SOUTHWEST)
0295 return diagonalmove(cell, dir);
0296 }
0297
0298 cell = DetId(0);
0299 return false;
0300 }
0301
0302 bool CaloGeometryHelper::simplemove(DetId& cell, const CaloDirection& dir) const {
0303 std::vector<DetId> neighbours;
0304 if (cell.subdetId() == EcalBarrel)
0305 neighbours = EcalBarrelTopology_->getNeighbours(cell, dir);
0306 else if (cell.subdetId() == EcalEndcap)
0307 neighbours = EcalEndcapTopology_->getNeighbours(cell, dir);
0308
0309 if ((!neighbours.empty()) && (!neighbours[0].null())) {
0310 cell = neighbours[0];
0311 return true;
0312 } else {
0313 cell = DetId(0);
0314 return false;
0315 }
0316 }
0317
0318 bool CaloGeometryHelper::diagonalmove(DetId& cell, const CaloDirection& dir) const {
0319 bool result;
0320
0321 if (dir == NORTHEAST) {
0322 result = simplemove(cell, NORTH);
0323 if (result)
0324 return simplemove(cell, EAST);
0325 else {
0326 result = simplemove(cell, EAST);
0327 if (result)
0328 return simplemove(cell, NORTH);
0329 else
0330 return false;
0331 }
0332 } else if (dir == NORTHWEST) {
0333 result = simplemove(cell, NORTH);
0334 if (result)
0335 return simplemove(cell, WEST);
0336 else {
0337 result = simplemove(cell, WEST);
0338 if (result)
0339 return simplemove(cell, NORTH);
0340 else
0341 return false;
0342 }
0343 } else if (dir == SOUTHEAST) {
0344 result = simplemove(cell, SOUTH);
0345 if (result)
0346 return simplemove(cell, EAST);
0347 else {
0348 result = simplemove(cell, EAST);
0349 if (result)
0350 return simplemove(cell, SOUTH);
0351 else
0352 return false;
0353 }
0354 } else if (dir == SOUTHWEST) {
0355 result = simplemove(cell, SOUTH);
0356 if (result)
0357 return simplemove(cell, WEST);
0358 else {
0359 result = simplemove(cell, SOUTH);
0360 if (result)
0361 return simplemove(cell, WEST);
0362 else
0363 return false;
0364 }
0365 }
0366 cell = DetId(0);
0367 return false;
0368 }
0369
0370 bool CaloGeometryHelper::borderCrossing(const DetId& c1, const DetId& c2) const {
0371 if (c1.subdetId() != c2.subdetId())
0372 return false;
0373
0374 if (c1.subdetId() == EcalBarrel) {
0375
0376
0377 EBDetId cc1(c1);
0378 EBDetId cc2(c2);
0379 return (cc1.im() != cc2.im() || cc1.ism() != cc2.ism());
0380 }
0381
0382 if (c1.subdetId() == EcalEndcap) {
0383
0384
0385 return (EEDetId(c1).isc() != EEDetId(c2).isc());
0386 }
0387 return false;
0388 }
0389
0390 void CaloGeometryHelper::buildCrystalArray() {
0391 const unsigned nbarrel = EBDetId::kSizeForDenseIndexing;
0392
0393 barrelCrystals_.resize(nbarrel, BaseCrystal());
0394
0395
0396 const std::vector<DetId>& vec(EcalBarrelGeometry_->getValidDetIds(DetId::Ecal, EcalBarrel));
0397 unsigned size = vec.size();
0398 for (unsigned ic = 0; ic < size; ++ic) {
0399 unsigned hashedindex = EBDetId(vec[ic]).hashedIndex();
0400 auto geom = EcalBarrelGeometry_->getGeometry(vec[ic]);
0401 BaseCrystal xtal(vec[ic]);
0402 xtal.setCorners(geom->getCorners(), geom->getPosition());
0403 barrelCrystals_[hashedindex] = xtal;
0404 }
0405
0406
0407
0408
0409 const std::vector<DetId>& vece(EcalEndcapGeometry_->getValidDetIds(DetId::Ecal, EcalEndcap));
0410 size = vece.size();
0411
0412
0413 const unsigned nendcap = EEDetId::kSizeForDenseIndexing;
0414
0415 endcapCrystals_.resize(nendcap, BaseCrystal());
0416 for (unsigned ic = 0; ic < size; ++ic) {
0417 unsigned hashedindex = EEDetId(vece[ic]).hashedIndex();
0418 auto geom = EcalEndcapGeometry_->getGeometry(vece[ic]);
0419 BaseCrystal xtal(vece[ic]);
0420 xtal.setCorners(geom->getCorners(), geom->getPosition());
0421 endcapCrystals_[hashedindex] = xtal;
0422 }
0423
0424 }