DBScanClusterBuilder.cxx

// #include "DetDescr/EcalGeometry.h"
// #include "Recon/Event/HgcrocDigiCollection.h"
#include "Recon/DBScanClusterBuilder.h"
 
#include <set>
 
#include "TFitResult.h"
#include "TGraph.h"
 
namespace recon {
 
DBScanClusterBuilder::DBScanClusterBuilder() {
  min_hit_energy_ = 0;
  cluster_hit_dist_ = 100;
  cluster_z_bias_ = 1;  // defaults to 1
  min_cluster_hit_mult_ = 2;
}
 
DBScanClusterBuilder::DBScanClusterBuilder(float minHitEnergy,
                                           float clusterHitDist,
                                           float clusterZBias,
                                           float minClusterHitMult) {
  min_hit_energy_ = minHitEnergy;
  cluster_hit_dist_ = clusterHitDist;
  cluster_z_bias_ = clusterZBias;  // clustering bias in the z_ direction
  min_cluster_hit_mult_ = minClusterHitMult;
}
 
std::vector<std::vector<const ldmx::CalorimeterHit*> >
DBScanClusterBuilder::runDBSCAN(
    const std::vector<const ldmx::CalorimeterHit*>& hits_) {
  const int n = hits_.size();
  std::vector<std::vector<const ldmx::CalorimeterHit*> > idx_clusters;
  std::vector<unsigned int> tried;
  tried.reserve(n);
  std::vector<unsigned int> used;
  used.reserve(n);
  for (unsigned int i = 0; i < n; i++) {
    if (isIn(i, tried)) continue;
    tried.push_back(i);
    ldmx_log(debug) << "trying " << i;
    if (hits_[i]->getEnergy() < min_hit_energy_) continue;
    std::set<unsigned int> neighbors;
    unsigned int n_nearby = 1;
    // find neighbors
    for (unsigned int j = 0; j < n; j++) {
      if (i != j &&
          dist(hits_[i], hits_[j]) < cluster_hit_dist_) {  // pair-wise distance
        neighbors.insert(j);
        if (hits_[j]->getEnergy() >= min_hit_energy_) n_nearby++;
      }
    }
    if (n_nearby >= min_cluster_hit_mult_) {
      std::vector<const ldmx::CalorimeterHit*> idx_cluster{
          hits_[i]};  // start a cluster
      used.push_back(i);
      ldmx_log(debug) << "- starting a cluster from " << i;
      for (unsigned int j : neighbors) {
        if (!isIn(j, tried)) {
          tried.push_back(j);
          ldmx_log(debug) << "== tried " << j;
          std::vector<unsigned int> neighbors2;
          for (unsigned int k = 0; k < n; k++) {
            if (dist(hits_[k], hits_[j]) < cluster_hit_dist_) {
              neighbors2.push_back(k);
            }
          }
          for (unsigned int k : neighbors2) neighbors.insert(k);
        }
        if (!isIn(j, used)) {
          ldmx_log(debug) << "== used " << j;
          used.push_back(j);
          idx_cluster.push_back(hits_[j]);
        }
      }
      idx_clusters.push_back(idx_cluster);
    }
  }
  ldmx_log(debug) << "done. writing this many clusters out: "
                  << idx_clusters.size();
  return idx_clusters;
}
 
void DBScanClusterBuilder::fillClusterInfoFromHits(
    ldmx::CaloCluster* cl, std::vector<const ldmx::CalorimeterHit*> hits_,
    bool logEnergyWeight, bool saveHitContribs) {
  float e(0), x(0), y(0), z(0), xx(0), yy(0), zz(0), n(0);
  float w = 1;  // weight
  float sumw = 0;
  std::vector<float> raw_xvals{};
  std::vector<float> raw_yvals{};
  std::vector<float> raw_zvals{};
  std::vector<float> raw_evals{};
  std::vector<const ldmx::CalorimeterHit*> constituent_hits;
 
  for (const ldmx::CalorimeterHit* h : hits_) {
    if (h->getEnergy() < min_hit_energy_) continue;
    if (logEnergyWeight) w = log(h->getEnergy()) - log(min_hit_energy_);
    e += h->getEnergy();
    x += w * h->getXPos();
    y += w * h->getYPos();
    z += w * h->getZPos();
    xx += w * h->getXPos() * h->getXPos();
    yy += w * h->getYPos() * h->getYPos();
    zz += w * h->getZPos() * h->getZPos();
    n += 1;
    sumw += w;
    if (saveHitContribs) {
      raw_xvals.push_back(h->getXPos());
      raw_yvals.push_back(h->getYPos());
      raw_zvals.push_back(h->getZPos());
      raw_evals.push_back(h->getEnergy());
      constituent_hits.emplace_back(h);
    }
  }  // over hits_
  x /= sumw;  // now is <x_>
  y /= sumw;
  z /= sumw;
  xx /= sumw;  // now is <x_^2>
  yy /= sumw;
  zz /= sumw;
  xx = sqrt(xx - x * x);  // now is sqrt(<x_^2>-<x_>^2)
  yy = sqrt(yy - y * y);
  zz = sqrt(zz - z * z);
  cl->setEnergy(e);
  cl->setNHits(n);
  cl->setCentroidXYZ(x, y, z);
  cl->setRMSXYZ(xx, yy, zz);
  if (saveHitContribs) {
    cl->setHitValsX(raw_xvals);
    cl->setHitValsY(raw_yvals);
    cl->setHitValsZ(raw_zvals);
    cl->setHitValsE(raw_evals);
    cl->addHits(constituent_hits);  // associate used hits_ to cluster
  }
 
  if (raw_xvals.size() > 2) {
    // skip fits for 'vertical' clusters
    std::vector<float> sorted_z = raw_zvals;
    std::sort(sorted_z.begin(), sorted_z.end());
    if ((sorted_z.size() > 2) and (sorted_z.back() - sorted_z.front() > 1e3)) {
      for (int i = 0; i < raw_xvals.size(); i++) {  // mean subtract
        raw_xvals[i] = raw_xvals[i] - x;
        raw_yvals[i] = raw_yvals[i] - y;
        raw_zvals[i] = raw_zvals[i] - z;
      }
 
      TGraph gxz(raw_zvals.size(), raw_zvals.data(), raw_xvals.data());
      auto r_xz = gxz.Fit("pol1", "SQ");  // p0 + x_*p1
      cl->setDXDZ(r_xz->Value(1));
      cl->setEDXDZ(r_xz->ParError(1));
 
      TGraph gyz(raw_zvals.size(), raw_zvals.data(), raw_yvals.data());
      auto r_yz = gyz.Fit("pol1", "SQ");  // p0 + x_*p1
      cl->setDYDZ(r_yz->Value(1));
      cl->setEDYDZ(r_yz->ParError(1));
    }
  }
  return;
}
 
}  // namespace recon