File indexing completed on 2023-05-08 23:19:15
0001 #include "RecoLocalCalo/HGCalRecAlgos/interface/HGCalDepthPreClusterer.h"
0002 #include "DataFormats/Math/interface/deltaR.h"
0003
0004 #include <list>
0005
0006 namespace {
0007 std::vector<size_t> sorted_indices(const reco::HGCalMultiCluster::ClusterCollection &v) {
0008
0009 std::vector<size_t> idx(v.size());
0010 for (size_t i = 0; i != idx.size(); ++i)
0011 idx[i] = i;
0012
0013
0014 std::sort(idx.begin(), idx.end(), [&v](size_t i1, size_t i2) { return (*v[i1]) > (*v[i2]); });
0015
0016 return idx;
0017 }
0018
0019 float dist2(const edm::Ptr<reco::BasicCluster> &a, const edm::Ptr<reco::BasicCluster> &b) {
0020 return reco::deltaR2(*a, *b);
0021 }
0022
0023
0024
0025 float distAxisCluster2(const edm::Ptr<reco::BasicCluster> &a, const edm::Ptr<reco::BasicCluster> &b) {
0026 float tanTheta = tan(2 * atan(exp(-1 * a->eta())));
0027 float ax = b->z() * tanTheta * cos(a->phi());
0028 float ay = b->z() * tanTheta * sin(a->phi());
0029 return (ax - b->x()) * (ax - b->x()) + (ay - b->y()) * (ay - b->y());
0030 }
0031 }
0032
0033 std::vector<reco::HGCalMultiCluster> HGCalDepthPreClusterer::makePreClusters(
0034 const reco::HGCalMultiCluster::ClusterCollection &thecls) const {
0035 std::vector<reco::HGCalMultiCluster> thePreClusters;
0036 std::vector<size_t> es = sorted_indices(thecls);
0037 std::vector<int> vused(es.size(), 0);
0038
0039 for (unsigned int i = 0; i < es.size(); ++i) {
0040 if (vused[i] == 0) {
0041 reco::HGCalMultiCluster temp;
0042 temp.push_back(thecls[es[i]]);
0043 vused[i] = (thecls[es[i]]->z() > 0) ? 1 : -1;
0044 for (unsigned int j = i + 1; j < es.size(); ++j) {
0045 if (vused[j] == 0) {
0046 float distanceCheck = 9999.;
0047 if (realSpaceCone)
0048 distanceCheck = distAxisCluster2(thecls[es[i]], thecls[es[j]]);
0049 else
0050 distanceCheck = dist2(thecls[es[i]], thecls[es[j]]);
0051 DetId detid = thecls[es[j]]->hitsAndFractions()[0].first();
0052 unsigned int layer = clusterTools->getLayer(detid);
0053 float radius = radii[2];
0054 if (layer <= rhtools_.lastLayerEE())
0055 radius = radii[0];
0056 else if (layer < rhtools_.firstLayerBH())
0057 radius = radii[1];
0058 float radius2 = radius * radius;
0059 if (distanceCheck<radius2 &&int(thecls[es[j]]->z() * vused[i])> 0) {
0060 temp.push_back(thecls[es[j]]);
0061 vused[j] = vused[i];
0062 }
0063 }
0064 }
0065 if (temp.size() > minClusters) {
0066 thePreClusters.push_back(temp);
0067 auto &back = thePreClusters.back();
0068 back.setPosition(clusterTools->getMultiClusterPosition(back));
0069 back.setEnergy(clusterTools->getMultiClusterEnergy(back));
0070 }
0071 }
0072 }
0073
0074 return thePreClusters;
0075 }