HGCalDepthPreClusterer.cc

Go to the documentation of this file.

#include "RecoLocalCalo/HGCalRecAlgos/interface/HGCalDepthPreClusterer.h"
#include "DataFormats/Math/interface/deltaR.h"

#include <list>

namespace {
  std::vector<size_t> sorted_indices(const reco::HGCalMultiCluster::ClusterCollection& v) {

    // initialize original index locations
    std::vector<size_t> idx(v.size());
    for (size_t i = 0; i != idx.size(); ++i) idx[i] = i;

    // sort indices based on comparing values in v
    std::sort(idx.begin(), idx.end(),
         [&v](size_t i1, size_t i2) {return (*v[i1]) > (*v[i2]);});

    return idx;
  }

  float dist2(const edm::Ptr<reco::BasicCluster> &a,
          const edm::Ptr<reco::BasicCluster> &b) {
    return reco::deltaR2(*a,*b);
  }

  //get distance between cluster and multicluster axis (defined by remaning cluster with highest energy)
  // N.B. the order of the clusters matters
  float distAxisCluster2(const edm::Ptr<reco::BasicCluster> &a,
             const edm::Ptr<reco::BasicCluster> &b) {
    float tanTheta = tan(2*atan(exp(-1*a->eta())));
    float ax = b->z()*tanTheta*cos(a->phi());
    float ay = b->z()*tanTheta*sin(a->phi());
    return (ax-b->x())*(ax-b->x()) + (ay-b->y())*(ay-b->y());
  }
}

std::vector<reco::HGCalMultiCluster> HGCalDepthPreClusterer::makePreClusters(const reco::HGCalMultiCluster::ClusterCollection &thecls) const {

  std::vector<reco::HGCalMultiCluster> thePreClusters;
  std::vector<size_t> es = sorted_indices(thecls);
  std::vector<int> vused(es.size(),0);
  unsigned int used = 0;

  for(unsigned int i = 0; i < es.size(); ++i) {
    if(vused[i]==0) {
      reco::HGCalMultiCluster temp;
      temp.push_back(thecls[es[i]]);
      vused[i]=(thecls[es[i]]->z()>0)? 1 : -1;
      ++used;
      for(unsigned int j = i+1; j < es.size(); ++j) {
    if(vused[j]==0) {
          float distanceCheck = 9999.;
          if( realSpaceCone ) distanceCheck = distAxisCluster2(thecls[es[i]],thecls[es[j]]);
          else distanceCheck = dist2(thecls[es[i]],thecls[es[j]]);
          DetId detid = thecls[es[j]]->hitsAndFractions()[0].first();
          unsigned int layer = clusterTools->getLayer(detid);
          float radius2 = 9999.;
          if(layer <= lastLayerEE) radius2 = radii[0]*radii[0];
          else if(layer <= lastLayerFH) radius2 = radii[1]*radii[1];
          else if(layer <= lastLayerBH) radius2 = radii[2]*radii[2];
          else assert(radius2<100. && "nonsense layer value - cannot assign multicluster radius");
      if( distanceCheck<radius2 && int(thecls[es[j]]->z()*vused[i])>0 ) {
        temp.push_back(thecls[es[j]]);
        vused[j]=vused[i];
        ++used;
      }
    }
      }
      if( temp.size() > minClusters ) {
        thePreClusters.push_back(temp);
        auto& back = thePreClusters.back();
        back.setPosition(clusterTools->getMultiClusterPosition(back));
        back.setEnergy(clusterTools->getMultiClusterEnergy(back));
      }
    }
  }



  return thePreClusters;
}