HGCalRecAlgos/src/HGCalDepthPreClusterer.cc

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     // initialize original index locations
0009     std::vector<size_t> idx(v.size());
0010     for (size_t i = 0; i != idx.size(); ++i)
0011       idx[i] = i;
0012
0013     // sort indices based on comparing values in v
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   //get distance between cluster and multicluster axis (defined by remaning cluster with highest energy)
0024   // N.B. the order of the clusters matters
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 }  // namespace
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 }