HGCalRecAlgos/src/HGCal3DClustering.cc

0001 #include "RecoLocalCalo/HGCalRecAlgos/interface/HGCal3DClustering.h"
0002 #include "DataFormats/Math/interface/deltaR.h"
0003
0004 namespace {
0005   std::vector<size_t> sorted_indices(const reco::HGCalMultiCluster::ClusterCollection &v) {
0006     // initialize original index locations
0007     std::vector<size_t> idx(v.size());
0008     std::iota(std::begin(idx), std::end(idx), 0);
0009
0010     // sort indices based on comparing values in v
0011     std::sort(idx.begin(), idx.end(), [&v](size_t i1, size_t i2) { return (*v[i1]) > (*v[i2]); });
0012
0013     return idx;
0014   }
0015
0016   float distReal2(const edm::Ptr<reco::BasicCluster> &a, const std::array<double, 3> &to) {
0017     return (a->x() - to[0]) * (a->x() - to[0]) + (a->y() - to[1]) * (a->y() - to[1]);
0018   }
0019 }  // namespace
0020
0021 void HGCal3DClustering::organizeByLayer(const reco::HGCalMultiCluster::ClusterCollection &thecls) {
0022   es = sorted_indices(thecls);
0023   unsigned int es_size = es.size();
0024   for (unsigned int i = 0; i < es_size; ++i) {
0025     int layer = rhtools_.getLayerWithOffset(thecls[es[i]]->hitsAndFractions()[0].first);
0026     layer += int(thecls[es[i]]->z() > 0) * (maxlayer + 1);
0027     float x = thecls[es[i]]->x();
0028     float y = thecls[es[i]]->y();
0029     float z = thecls[es[i]]->z();
0030     points[layer].emplace_back(ClusterRef(i, z), x, y);
0031     if (zees[layer] == 0.) {
0032       // At least one cluster for layer at z
0033       zees[layer] = z;
0034     }
0035     if (points[layer].empty()) {
0036       minpos[layer][0] = x;
0037       minpos[layer][1] = y;
0038       maxpos[layer][0] = x;
0039       maxpos[layer][1] = y;
0040     } else {
0041       minpos[layer][0] = std::min(x, minpos[layer][0]);
0042       minpos[layer][1] = std::min(y, minpos[layer][1]);
0043       maxpos[layer][0] = std::max(x, maxpos[layer][0]);
0044       maxpos[layer][1] = std::max(y, maxpos[layer][1]);
0045     }
0046   }
0047 }
0048 std::vector<reco::HGCalMultiCluster> HGCal3DClustering::makeClusters(
0049     const reco::HGCalMultiCluster::ClusterCollection &thecls) {
0050   reset();
0051   organizeByLayer(thecls);
0052   std::vector<reco::HGCalMultiCluster> thePreClusters;
0053
0054   std::vector<KDTree> hit_kdtree(2 * (maxlayer + 1));
0055   for (unsigned int i = 0; i <= 2 * maxlayer + 1; ++i) {
0056     KDTreeBox bounds(minpos[i][0], maxpos[i][0], minpos[i][1], maxpos[i][1]);
0057     hit_kdtree[i].build(points[i], bounds);
0058   }
0059   std::vector<int> vused(es.size(), 0);
0060
0061   unsigned int es_size = es.size();
0062   for (unsigned int i = 0; i < es_size; ++i) {
0063     if (vused[i] == 0) {
0064       reco::HGCalMultiCluster temp;
0065       temp.push_back(thecls[es[i]]);
0066       vused[i] = (thecls[es[i]]->z() > 0) ? 1 : -1;
0067       // Starting from cluster es[i] at from[0] - from[1] - from[2]
0068       std::array<double, 3> from{{thecls[es[i]]->x(), thecls[es[i]]->y(), thecls[es[i]]->z()}};
0069       unsigned int firstlayer = int(thecls[es[i]]->z() > 0) * (maxlayer + 1);
0070       unsigned int lastlayer = firstlayer + maxlayer + 1;
0071       for (unsigned int j = firstlayer; j < lastlayer; ++j) {
0072         if (zees[j] == 0.) {
0073           // layer j not yet ever reached?
0074           continue;
0075         }
0076         std::array<double, 3> to{{0., 0., zees[j]}};
0077         layerIntersection(to, from);
0078         unsigned int layer =
0079             j > maxlayer ? (j - (maxlayer + 1)) : j;  //maps back from index used for KD trees to actual layer
0080         float radius = radii[2];
0081         if (layer <= rhtools_.lastLayerEE())
0082           radius = radii[0];
0083         else if (layer < rhtools_.firstLayerBH())
0084           radius = radii[1];
0085         float radius2 = radius * radius;
0086         KDTreeBox search_box(
0087             float(to[0]) - radius, float(to[0]) + radius, float(to[1]) - radius, float(to[1]) + radius);
0088         std::vector<ClusterRef> found;
0089         // at layer j in box float(to[0])+/-radius - float(to[1])+/-radius
0090         hit_kdtree[j].search(search_box, found);
0091         // found found.size() clusters within box
0092         for (unsigned int k = 0; k < found.size(); k++) {
0093           if (vused[found[k].ind] == 0 && distReal2(thecls[es[found[k].ind]], to) < radius2) {
0094             temp.push_back(thecls[es[found[k].ind]]);
0095             vused[found[k].ind] = vused[i];
0096           }
0097         }
0098       }
0099       if (temp.size() > minClusters) {
0100         math::XYZPoint position = clusterTools->getMultiClusterPosition(temp);
0101         if (std::abs(position.z()) <= 0.)
0102           continue;
0103         // only store multiclusters that pass the energy threshold in getMultiClusterPosition
0104         // giving them a position inside the HGCal
0105         thePreClusters.push_back(temp);
0106         auto &back = thePreClusters.back();
0107         back.setPosition(position);
0108         back.setEnergy(clusterTools->getMultiClusterEnergy(back));
0109       }
0110     }
0111   }
0112
0113   return thePreClusters;
0114 }
0115
0116 void HGCal3DClustering::layerIntersection(std::array<double, 3> &to, const std::array<double, 3> &from) const {
0117   if (from[2] != 0) {
0118     to[0] = from[0] / from[2] * to[2];
0119     to[1] = from[1] / from[2] * to[2];
0120   } else {
0121     to[0] = 0;
0122     to[1] = 0;
0123   }
0124 }