File indexing completed on 2024-04-06 12:24:40
0001 #include <vector>
0002 #include <map>
0003 #include <iostream>
0004
0005 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0006 #include "RecoEcal/EgammaClusterAlgos/interface/PreshowerPhiClusterAlgo.h"
0007 #include "RecoCaloTools/Navigation/interface/EcalPreshowerNavigator.h"
0008 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0009 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011 #include "DataFormats/Math/interface/deltaPhi.h"
0012
0013 #include "TH1.h"
0014
0015 reco::PreshowerCluster PreshowerPhiClusterAlgo::makeOneCluster(ESDetId strip,
0016 HitsID* used_strips,
0017 RecHitsMap* the_rechitsMap_p,
0018 const CaloSubdetectorGeometry* geometry_p,
0019 double deltaEta,
0020 double minDeltaPhi,
0021 double maxDeltaPhi) {
0022 rechits_map = the_rechitsMap_p;
0023 used_s = used_strips;
0024
0025 int plane = strip.plane();
0026
0027
0028 std::vector<std::pair<DetId, float> > dummy;
0029 Point posi(0, 0, 0);
0030 reco::PreshowerCluster nullcluster = reco::PreshowerCluster(0., posi, dummy, plane);
0031 if (strip == ESDetId(0))
0032 return nullcluster;
0033
0034 auto refCell = geometry_p->getGeometry(strip);
0035 const GlobalPoint& refpos = refCell->getPosition();
0036 double refEta = refpos.eta();
0037 double refPhi = refpos.phi();
0038
0039
0040 EcalRecHitCollection clusterRecHits;
0041
0042 double x_pos = 0;
0043 double y_pos = 0;
0044 double z_pos = 0;
0045
0046 RecHitsMap::iterator strip_it;
0047 for (strip_it = rechits_map->begin(); strip_it != rechits_map->end(); ++strip_it) {
0048 if (!goodStrip(strip_it))
0049 continue;
0050
0051 ESDetId mystrip = (strip_it->first == DetId(0)) ? ESDetId(0) : ESDetId(strip_it->first);
0052 if (mystrip.plane() != plane)
0053 continue;
0054
0055 auto thisCell = geometry_p->getGeometry(strip_it->first);
0056 const GlobalPoint& position = thisCell->getPosition();
0057
0058 if (fabs(position.eta() - refEta) < deltaEta) {
0059
0060
0061 if (reco::deltaPhi(position.phi(), refPhi) > 0 && reco::deltaPhi(position.phi(), refPhi) > maxDeltaPhi)
0062 continue;
0063 if (reco::deltaPhi(position.phi(), refPhi) < 0 && reco::deltaPhi(position.phi(), refPhi) < minDeltaPhi)
0064 continue;
0065
0066
0067
0068 clusterRecHits.push_back(strip_it->second);
0069 used_s->insert(strip_it->first);
0070
0071 x_pos += strip_it->second.energy() * position.x();
0072 y_pos += strip_it->second.energy() * position.y();
0073 z_pos += strip_it->second.energy() * position.z();
0074 }
0075 }
0076
0077 EcalRecHitCollection::iterator it;
0078 double Eclust = 0;
0079
0080 std::vector<std::pair<DetId, float> > usedHits;
0081 for (it = clusterRecHits.begin(); it != clusterRecHits.end(); it++) {
0082 Eclust += it->energy();
0083 usedHits.push_back(std::pair<DetId, float>(it->id(), 1.));
0084 }
0085
0086 if (Eclust > 0.) {
0087 x_pos /= Eclust;
0088 y_pos /= Eclust;
0089 z_pos /= Eclust;
0090 }
0091 Point pos(x_pos, y_pos, z_pos);
0092
0093 reco::PreshowerCluster cluster = reco::PreshowerCluster(Eclust, pos, usedHits, plane);
0094
0095 return cluster;
0096 }
0097
0098
0099 bool PreshowerPhiClusterAlgo::goodStrip(RecHitsMap::iterator candidate_it) {
0100 if (used_s->find(candidate_it->first) != used_s->end())
0101 LogTrace("PreShowerPhiClusterAlgo") << " This strip is in use";
0102 if (candidate_it == rechits_map->end())
0103 LogTrace("PreShowerPhiClusterAlgo") << " No such a strip in rechits_map";
0104 if (candidate_it->second.energy() <= esStripEnergyCut_)
0105 LogTrace("PreShowerPhiClusterAlgo") << "Strip energy" << candidate_it->second.energy() << "is below threshold";
0106
0107
0108 if ((used_s->find(candidate_it->first) != used_s->end()) ||
0109 (candidate_it == rechits_map->end()) ||
0110 (candidate_it->second.energy() <= esStripEnergyCut_))
0111 {
0112 return false;
0113 }
0114 return true;
0115 }