Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:24:32

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   unsigned int used = 0;
0061 
0062   unsigned int es_size = es.size();
0063   for (unsigned int i = 0; i < es_size; ++i) {
0064     if (vused[i] == 0) {
0065       reco::HGCalMultiCluster temp;
0066       temp.push_back(thecls[es[i]]);
0067       vused[i] = (thecls[es[i]]->z() > 0) ? 1 : -1;
0068       ++used;
0069       // Starting from cluster es[i] at from[0] - from[1] - from[2]
0070       std::array<double, 3> from{{thecls[es[i]]->x(), thecls[es[i]]->y(), thecls[es[i]]->z()}};
0071       unsigned int firstlayer = int(thecls[es[i]]->z() > 0) * (maxlayer + 1);
0072       unsigned int lastlayer = firstlayer + maxlayer + 1;
0073       for (unsigned int j = firstlayer; j < lastlayer; ++j) {
0074         if (zees[j] == 0.) {
0075           // layer j not yet ever reached?
0076           continue;
0077         }
0078         std::array<double, 3> to{{0., 0., zees[j]}};
0079         layerIntersection(to, from);
0080         unsigned int layer =
0081             j > maxlayer ? (j - (maxlayer + 1)) : j;  //maps back from index used for KD trees to actual layer
0082         float radius = radii[2];
0083         if (layer <= rhtools_.lastLayerEE())
0084           radius = radii[0];
0085         else if (layer < rhtools_.firstLayerBH())
0086           radius = radii[1];
0087         float radius2 = radius * radius;
0088         KDTreeBox search_box(
0089             float(to[0]) - radius, float(to[0]) + radius, float(to[1]) - radius, float(to[1]) + radius);
0090         std::vector<ClusterRef> found;
0091         // at layer j in box float(to[0])+/-radius - float(to[1])+/-radius
0092         hit_kdtree[j].search(search_box, found);
0093         // found found.size() clusters within box
0094         for (unsigned int k = 0; k < found.size(); k++) {
0095           if (vused[found[k].ind] == 0 && distReal2(thecls[es[found[k].ind]], to) < radius2) {
0096             temp.push_back(thecls[es[found[k].ind]]);
0097             vused[found[k].ind] = vused[i];
0098             ++used;
0099           }
0100         }
0101       }
0102       if (temp.size() > minClusters) {
0103         math::XYZPoint position = clusterTools->getMultiClusterPosition(temp);
0104         if (std::abs(position.z()) <= 0.)
0105           continue;
0106         // only store multiclusters that pass the energy threshold in getMultiClusterPosition
0107         // giving them a position inside the HGCal
0108         thePreClusters.push_back(temp);
0109         auto &back = thePreClusters.back();
0110         back.setPosition(position);
0111         back.setEnergy(clusterTools->getMultiClusterEnergy(back));
0112       }
0113     }
0114   }
0115 
0116   return thePreClusters;
0117 }
0118 
0119 void HGCal3DClustering::layerIntersection(std::array<double, 3> &to, const std::array<double, 3> &from) const {
0120   if (from[2] != 0) {
0121     to[0] = from[0] / from[2] * to[2];
0122     to[1] = from[1] / from[2] * to[2];
0123   } else {
0124     to[0] = 0;
0125     to[1] = 0;
0126   }
0127 }