1 #ifndef ECALBOUNDARYINFOCALCULATOR_H_
2 #define ECALBOUNDARYINFOCALCULATOR_H_
25 template <
class EcalDetId>
45 std::vector<int>& stati)
const {
47 EcalDetId hitdetid = EcalDetId(hit.
id());
52 int hitIeta = ebhitdetid.
ieta();
53 int hitIphi = ebhitdetid.
iphi();
55 for (
int ieta = -1; ieta <= 1; ieta++) {
56 for (
int iphi = -1; iphi <= 1; iphi++) {
57 if ((iphi == 0 && ieta == 0) || iphi * ieta != 0)
60 int neighbourIeta = hitIeta + ieta;
61 int neighbourIphi = hitIphi + iphi;
63 if (neighbourIphi < 1)
65 if (neighbourIphi > 360)
67 if (neighbourIeta == 0) {
68 neighbourIeta += ieta;
75 int status = (chit != ecalStatus.
end()) ? chit->getStatusCode() & 0x1F : -1;
79 for (std::vector<int>::const_iterator
s = stati.begin();
s != stati.end(); ++
s) {
86 stati.push_back(status);
94 int hitIx = eehitdetid.
ix();
95 int hitIy = eehitdetid.
iy();
96 int hitIz = eehitdetid.
zside();
98 for (
int ix = -1; ix <= 1; ix++) {
99 for (
int iy = -1; iy <= 1; iy++) {
100 if ((ix == 0 && iy == 0) || ix * iy != 0)
103 int neighbourIx = hitIx + ix;
104 int neighbourIy = hitIy + iy;
109 int status = (chit != ecalStatus.
end()) ? chit->getStatusCode() & 0x1F : -1;
112 bool present =
false;
113 for (std::vector<int>::const_iterator
s = stati.begin();
s != stati.end(); ++
s) {
120 stati.push_back(status);
127 edm::LogWarning(
"EcalBoundaryInfoCalculator") <<
"ERROR - RecHit belongs to wrong sub detector";
138 EcalDetId hitdetid = EcalDetId(hit.
id());
143 int hitIeta = ebhitdetid.
ieta();
144 int hitIphi = ebhitdetid.
iphi();
146 for (
int ieta = -1; ieta <= 1; ieta++) {
147 for (
int iphi = -1; iphi <= 1; iphi++) {
148 if ((iphi == 0 && ieta == 0) || iphi * ieta != 0)
151 int neighbourIeta = hitIeta + ieta;
152 int neighbourIphi = hitIphi + iphi;
154 if (neighbourIphi < 1)
155 neighbourIphi += 360;
156 if (neighbourIphi > 360)
157 neighbourIphi -= 360;
158 if (neighbourIeta == 0) {
159 neighbourIeta += ieta;
171 int hitIx = eehitdetid.
ix();
172 int hitIy = eehitdetid.
iy();
173 int hitIz = eehitdetid.
zside();
175 for (
int ix = -1; ix <= 1; ix++) {
176 for (
int iy = -1; iy <= 1; iy++) {
177 if ((ix == 0 && iy == 0) || ix * iy != 0)
180 int neighbourIx = hitIx + ix;
181 int neighbourIy = hitIy + iy;
190 edm::LogWarning(
"EcalBoundaryInfoCalculator") <<
"ERROR - RecHit belongs to wrong sub detector";
197 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"set Debug Mode!";
206 next = theNavi->
north();
210 next = theNavi->
east();
214 next = theNavi->
south();
218 next = theNavi->
west();
231 oppDirection = oIt->second;
234 EcalDetId currDetId = theEcalNav->
pos();
241 std::map<CdOrientation, CdOrientation> turnMap =
nextDirs;
242 if (reverseOrientation)
244 std::map<CdOrientation, CdOrientation>::iterator nIt = turnMap.find(currDirection);
246 if (nIt != turnMap.end())
247 nextDirection = (*nIt).second;
249 edm::LogWarning(
"EcalBoundaryInfoCalculator") <<
"ERROR - no Next Direction found!?!?";
250 return nextDirection;
255 std::map<CdOrientation, CdOrientation> turnMap =
prevDirs;
256 if (reverseOrientation)
258 std::map<CdOrientation, CdOrientation>::iterator nIt = turnMap.find(currDirection);
260 if (nIt != turnMap.end())
261 nextDirection = (*nIt).second;
263 edm::LogWarning(
"EcalBoundaryInfoCalculator") <<
"ERROR - no Next Direction found!?!?";
264 return nextDirection;
270 std::unique_ptr<CaloNavigator<EcalDetId>> theEcalNav;
272 theEcalNav = std::make_unique<CaloNavigator<EcalDetId>>(
275 theEcalNav = std::make_unique<CaloNavigator<EcalDetId>>(
279 <<
"initializeEcalNavigator not implemented for subDet: " << ecalSubDet;
290 template <
class EcalDetId>
304 oppositeDirs.clear();
313 template <
class EcalDetId>
316 template <
class EcalDetId>
323 std::vector<EcalRecHit> boundaryRecHits;
324 std::vector<DetId> boundaryDetIds;
325 std::vector<int> stati;
327 double boundaryEnergy = 0;
328 double boundaryET = 0;
329 int beCellCounter = 0;
330 bool nextToBorder =
false;
332 boundaryRecHits.push_back(*hit);
334 boundaryEnergy += hit->
energy();
335 EcalDetId hitdetid = (EcalDetId)hit->
id();
336 boundaryDetIds.push_back(hitdetid);
339 double eta = cellGeom->getPosition().eta();
340 boundaryET += hit->
energy() / cosh(eta);
343 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"Find Boundary RecHits...";
347 <<
"Starting at : (" << ((
EBDetId)hitdetid).ieta() <<
"," << ((
EBDetId)hitdetid).iphi() <<
")";
350 <<
"Starting at : (" << ((
EEDetId)hitdetid).ix() <<
"," << ((
EEDetId)hitdetid).iy() <<
")";
355 auto theEcalNav = initializeEcalNavigator(hitdetid, theCaloTopology, EcalDetId::subdet());
357 bool reverseOrientation =
false;
360 EcalDetId
start = hitdetid;
361 EcalDetId current =
start;
362 int current_status = 0;
365 bool startAlgo =
false;
368 next = makeStepInDirection(currDirection, theEcalNav.get());
369 theEcalNav->setHome(current);
372 int status = (chit != ecalStatus.
end()) ? chit->getStatusCode() & 0x1F : -1;
374 stati.push_back(status);
378 currDirection = turnLeft(currDirection, reverseOrientation);
381 throw cms::Exception(
"NoStartingDirection") <<
"EcalBoundaryInfoCalculator: No starting direction can be found: "
382 "This should never happen if RecHit has a dead neighbour!";
387 currDirection = turnRight(currDirection, reverseOrientation);
390 bool nextIsStart =
false;
391 bool atBorder =
false;
393 while (!nextIsStart) {
394 bool nextStepFound =
false;
397 while (!nextStepFound) {
398 next = makeStepInDirection(currDirection, theEcalNav.get());
399 theEcalNav->setHome(current);
402 status = (chit != ecalStatus.
end()) ? chit->getStatusCode() & 0x1F : -1;
405 bool present =
false;
406 for (std::vector<int>::const_iterator
s = stati.begin();
s != stati.end(); ++
s) {
413 stati.push_back(status);
416 nextStepFound =
true;
418 currDirection = turnRight(currDirection, reverseOrientation);
420 }
else if (next == EcalDetId(0)) {
422 currDirection = turnLeft(currDirection, reverseOrientation);
424 }
else if (status == 0) {
425 nextStepFound =
true;
430 <<
"EcalBoundaryInfoCalculator: No valid next direction can be found: This should never happen!";
435 next = makeStepInDirection(currDirection, theEcalNav.get());
440 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"Boundary path reached starting position!";
445 <<
"Next step: " << (EcalDetId)next <<
" Status: " << status <<
" Start: " << (EcalDetId)
start;
448 if ((!atBorder || status == 0) && !nextIsStart) {
449 boundaryDetIds.push_back(next);
450 if (RecHits.
find(next) != RecHits.
end() && status == 0) {
453 boundaryRecHits.push_back(nexthit);
454 boundaryEnergy += nexthit.
energy();
456 eta = cellGeom->getPosition().eta();
457 boundaryET += nexthit.
energy() / cosh(eta);
461 if (current_status == 0 && status == 0 && atBorder) {
463 currDirection = turnRight(currDirection, reverseOrientation);
466 if (status == 0 && atBorder) {
468 currDirection = turnLeft(currDirection, reverseOrientation);
472 currDirection = turnLeft(currDirection, reverseOrientation);
475 currDirection = turnRight(currDirection, reverseOrientation);
485 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"<<<<<<<<<<<<<<< Final Boundary object <<<<<<<<<<<<<<<";
486 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"no of neighbouring RecHits: " << boundaryRecHits.size();
487 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"no of neighbouring DetIds: " << boundaryDetIds.size();
488 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"boundary energy: " << boundaryEnergy;
489 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"boundary ET: " << boundaryET;
490 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"no of cells contributing to boundary energy: " << beCellCounter;
491 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"Channel stati: ";
492 for (std::vector<int>::iterator it = stati.begin(); it != stati.end(); ++it) {
493 edm::LogInfo(
"EcalBoundaryInfoCalculator") << *it <<
" ";
499 boundInfo.
subdet = hitdetid.subdet();
500 boundInfo.
detIds = boundaryDetIds;
501 boundInfo.
recHits = boundaryRecHits;
510 template <
class EcalDetId>
517 std::vector<EcalRecHit> gapRecHits;
518 std::vector<DetId> gapDetIds;
520 double gapEnergy = 0;
522 int gapCellCounter = 0;
523 bool nextToBorder =
false;
525 gapRecHits.push_back(*hit);
527 gapEnergy += hit->
energy();
528 EcalDetId hitdetid = (EcalDetId)hit->
id();
529 gapDetIds.push_back(hitdetid);
532 double eta = cellGeom->getPosition().eta();
533 gapET += hit->
energy() / cosh(eta);
536 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"Find Border RecHits...";
540 <<
"Starting at : (" << ((
EBDetId)hitdetid).ieta() <<
"," << ((
EBDetId)hitdetid).iphi() <<
")";
543 <<
"Starting at : (" << ((
EEDetId)hitdetid).ix() <<
"," << ((
EEDetId)hitdetid).iy() <<
")";
548 auto theEcalNav = initializeEcalNavigator(hitdetid, theCaloTopology, EcalDetId::subdet());
550 bool reverseOrientation =
false;
553 EcalDetId
start = hitdetid;
554 EcalDetId current =
start;
557 bool startAlgo =
false;
560 next = makeStepInDirection(currDirection, theEcalNav.get());
561 theEcalNav->setHome(start);
563 if (next == EcalDetId(0)) {
568 currDirection = turnLeft(currDirection, reverseOrientation);
571 throw cms::Exception(
"NoStartingDirection") <<
"EcalBoundaryInfoCalculator: No starting direction can be found: "
572 "This should never happen if RecHit is at border!";
578 currDirection = turnLeft(currDirection, reverseOrientation);
581 bool endIsFound =
false;
582 bool startIsEnd =
false;
584 while (!endIsFound) {
585 bool nextStepFound =
false;
588 while (!nextStepFound) {
589 next = makeStepInDirection(currDirection, theEcalNav.get());
590 theEcalNav->setHome(current);
593 status = (chit != ecalStatus.
end()) ? chit->getStatusCode() & 0x1F : -1;
598 }
else if (next == EcalDetId(0)) {
600 currDirection = turnLeft(currDirection, reverseOrientation);
601 }
else if (status == 0) {
602 if (RecHits.
find(next) != RecHits.
end()) {
603 nextStepFound =
true;
612 <<
"EcalBoundaryInfoCalculator: No valid next direction can be found: This should never happen!";
617 next = makeStepInDirection(currDirection, theEcalNav.get());
624 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"Path along gap reached starting position!";
629 <<
"Next step: " << (EcalDetId)next <<
" Status: " << status <<
" Start: " << (EcalDetId)
start;
631 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"End of gap cluster is found going left";
636 gapDetIds.push_back(next);
637 if (RecHits.
find(next) != RecHits.
end()) {
640 gapRecHits.push_back(nexthit);
641 gapEnergy += nexthit.
energy();
643 eta = cellGeom->getPosition().eta();
644 gapET += nexthit.
energy() / cosh(eta);
649 currDirection = turnRight(currDirection, reverseOrientation);
653 theEcalNav->setHome(start);
656 currDirection = startDirection;
657 currDirection = turnRight(currDirection, reverseOrientation);
663 while (!endIsFound) {
664 bool nextStepFound =
false;
667 while (!nextStepFound) {
668 next = makeStepInDirection(currDirection, theEcalNav.get());
669 theEcalNav->setHome(current);
672 status = (chit != ecalStatus.
end()) ? chit->getStatusCode() & 0x1F : -1;
677 }
else if (next == EcalDetId(0)) {
679 currDirection = turnRight(currDirection, reverseOrientation);
680 }
else if (status == 0) {
681 if (RecHits.
find(next) != RecHits.
end()) {
682 nextStepFound =
true;
691 <<
"EcalBoundaryInfoCalculator: No valid next direction can be found: This should never happen!";
696 next = makeStepInDirection(currDirection, theEcalNav.get());
701 <<
"Next step: " << (EcalDetId)next <<
" Status: " << status <<
" Start: " << (EcalDetId)
start;
703 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"End of gap cluster is found going right";
708 gapDetIds.push_back(next);
709 if (RecHits.
find(next) != RecHits.
end()) {
712 gapRecHits.push_back(nexthit);
713 gapEnergy += nexthit.
energy();
715 eta = cellGeom->getPosition().eta();
716 gapET += nexthit.
energy() / cosh(eta);
721 currDirection = turnLeft(currDirection, reverseOrientation);
726 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"<<<<<<<<<<<<<<< Final Gap object <<<<<<<<<<<<<<<";
727 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"No of RecHits along gap: " << gapRecHits.size();
728 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"No of DetIds along gap: " << gapDetIds.size();
729 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"Gap energy: " << gapEnergy;
730 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"Gap ET: " << gapET;
734 gapInfo.
subdet = hitdetid.subdet();
735 gapInfo.
detIds = gapDetIds;
740 std::vector<int> stati;
EcalBoundaryInfoCalculator()
BoundaryInformation gapRecHits(const EcalRecHitCollection &, const EcalRecHit *, const CaloTopology &theCaloTopology, const EcalChannelStatus &ecalStatus, const CaloGeometry &geometry) const
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
CdOrientation turnRight(CdOrientation currDirection, bool reverseOrientation) const
std::map< CdOrientation, CdOrientation > nextDirs
BoundaryInformation boundaryRecHits(const EcalRecHitCollection &, const EcalRecHit *, const CaloTopology &theCaloTopology, const EcalChannelStatus &ecalStatus, const CaloGeometry &geometry) const
EcalDetId makeStepInDirection(CdOrientation direction, const CaloNavigator< EcalDetId > *theNavi) const
static bool validDetId(int i, int j)
check if a valid index combination
int iphi() const
get the crystal iphi
void setHome(const T &startingPoint)
set the starting position
std::unique_ptr< CaloNavigator< EcalDetId > > initializeEcalNavigator(DetId startE, const CaloTopology &theCaloTopology, EcalSubdetector ecalSubDet) const
CdOrientation goBackOneCell(CdOrientation currDirection, EcalDetId prev, CaloNavigator< EcalDetId > *theEcalNav) const
T west() const
move the navigator west
T south() const
move the navigator south
int ieta() const
get the crystal ieta
T pos() const
get the current position
static const int ETAPHIMODE
bool checkRecHitHasDeadNeighbour(const EcalRecHit &hit, const EcalChannelStatus &ecalStatus, std::vector< int > &stati) const
const_iterator end() const
Log< level::Info, false > LogInfo
T east() const
move the navigator east
CdOrientation turnLeft(CdOrientation currDirection, bool reverseOrientation) const
DetId id() const
get the id
virtual std::shared_ptr< const CaloCellGeometry > getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
static bool validDetId(int crystal_ix, int crystal_iy, int iz)
std::vector< Item >::const_iterator const_iterator
std::map< CdOrientation, CdOrientation > prevDirs
const CaloSubdetectorTopology * getSubdetectorTopology(const DetId &id) const
access the subdetector Topology for the given subdetector directly
std::map< CdOrientation, CdOrientation > oppositeDirs
iterator find(key_type k)
const_iterator find(uint32_t rawId) const
T north() const
move the navigator north
const_iterator end() const
bool checkRecHitHasInvalidNeighbour(const EcalRecHit &hit, const EcalChannelStatus &ecalStatus) const
Log< level::Warning, false > LogWarning
~EcalBoundaryInfoCalculator()