Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-01-13 01:43:56

0001 #include "RecoLocalCalo/HGCalRecProducers/interface/HGCalImagingAlgo.h"
0002 
0003 // Geometry
0004 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0005 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0006 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0007 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0008 
0009 #include "RecoEcal/EgammaCoreTools/interface/PositionCalc.h"
0010 //
0011 #include "DataFormats/CaloRecHit/interface/CaloID.h"
0012 #include "oneapi/tbb/task_arena.h"
0013 #include "oneapi/tbb.h"
0014 
0015 using namespace hgcal_clustering;
0016 
0017 void HGCalImagingAlgo::getEventSetupPerAlgorithm(const edm::EventSetup &es) {
0018   points_.clear();
0019   minpos_.clear();
0020   maxpos_.clear();
0021   points_.resize(2 * (maxlayer_ + 1));
0022   minpos_.resize(2 * (maxlayer_ + 1), {{0.0f, 0.0f}});
0023   maxpos_.resize(2 * (maxlayer_ + 1), {{0.0f, 0.0f}});
0024 }
0025 
0026 void HGCalImagingAlgo::populate(const HGCRecHitCollection &hits) {
0027   // loop over all hits and create the Hexel structure, skip energies below ecut
0028 
0029   if (dependSensor_) {
0030     // for each layer and wafer calculate the thresholds (sigmaNoise and energy)
0031     // once
0032     computeThreshold();
0033   }
0034 
0035   std::vector<bool> firstHit(2 * (maxlayer_ + 1), true);
0036   for (unsigned int i = 0; i < hits.size(); ++i) {
0037     const HGCRecHit &hgrh = hits[i];
0038     DetId detid = hgrh.detid();
0039     unsigned int layer = rhtools_.getLayerWithOffset(detid);
0040     float thickness = 0.f;
0041     // set sigmaNoise default value 1 to use kappa value directly in case of
0042     // sensor-independent thresholds
0043     float sigmaNoise = 1.f;
0044     if (dependSensor_) {
0045       thickness = rhtools_.getSiThickness(detid);
0046       int thickness_index = rhtools_.getSiThickIndex(detid);
0047       if (thickness_index == -1)
0048         thickness_index = 3;
0049       double storedThreshold = thresholds_[layer - 1][thickness_index];
0050       sigmaNoise = sigmaNoise_[layer - 1][thickness_index];
0051 
0052       if (hgrh.energy() < storedThreshold)
0053         continue;  // this sets the ZS threshold at ecut times the sigma noise
0054                    // for the sensor
0055     }
0056     if (!dependSensor_ && hgrh.energy() < ecut_)
0057       continue;
0058 
0059     // map layers from positive endcap (z) to layer + maxlayer+1 to prevent
0060     // mixing up hits from different sides
0061     layer += int(rhtools_.zside(detid) > 0) * (maxlayer_ + 1);
0062 
0063     // determine whether this is a half-hexagon
0064     bool isHalf = rhtools_.isHalfCell(detid);
0065     const GlobalPoint position(rhtools_.getPosition(detid));
0066 
0067     // here's were the KDNode is passed its dims arguments - note that these are
0068     // *copied* from the Hexel
0069     points_[layer].emplace_back(
0070         Hexel(hgrh, detid, isHalf, sigmaNoise, thickness, &rhtools_), position.x(), position.y());
0071 
0072     // for each layer, store the minimum and maximum x and y coordinates for the
0073     // KDTreeBox boundaries
0074     if (firstHit[layer]) {
0075       minpos_[layer][0] = position.x();
0076       minpos_[layer][1] = position.y();
0077       maxpos_[layer][0] = position.x();
0078       maxpos_[layer][1] = position.y();
0079       firstHit[layer] = false;
0080     } else {
0081       minpos_[layer][0] = std::min((float)position.x(), minpos_[layer][0]);
0082       minpos_[layer][1] = std::min((float)position.y(), minpos_[layer][1]);
0083       maxpos_[layer][0] = std::max((float)position.x(), maxpos_[layer][0]);
0084       maxpos_[layer][1] = std::max((float)position.y(), maxpos_[layer][1]);
0085     }
0086 
0087   }  // end loop hits
0088 }
0089 // Create a vector of Hexels associated to one cluster from a collection of
0090 // HGCalRecHits - this can be used directly to make the final cluster list -
0091 // this method can be invoked multiple times for the same event with different
0092 // input (reset should be called between events)
0093 void HGCalImagingAlgo::makeClusters() {
0094   layerClustersPerLayer_.resize(2 * maxlayer_ + 2);
0095   // assign all hits in each layer to a cluster core or halo
0096   tbb::this_task_arena::isolate([&] {
0097     tbb::parallel_for(size_t(0), size_t(2 * maxlayer_ + 2), [&](size_t i) {
0098       KDTreeBox bounds(minpos_[i][0], maxpos_[i][0], minpos_[i][1], maxpos_[i][1]);
0099       KDTree hit_kdtree;
0100       hit_kdtree.build(points_[i], bounds);
0101 
0102       unsigned int actualLayer =
0103           i > maxlayer_ ? (i - (maxlayer_ + 1)) : i;  // maps back from index used for KD trees to actual layer
0104 
0105       double maxdensity = calculateLocalDensity(points_[i], hit_kdtree, actualLayer);  // also stores rho (energy
0106                                                                                        // density) for each point (node)
0107       // calculate distance to nearest point with higher density storing
0108       // distance (delta) and point's index
0109       calculateDistanceToHigher(points_[i]);
0110       findAndAssignClusters(points_[i], hit_kdtree, maxdensity, bounds, actualLayer, layerClustersPerLayer_[i]);
0111     });
0112   });
0113   //Now that we have the density per point we can store it
0114   for (auto const &p : points_) {
0115     setDensity(p);
0116   }
0117 }
0118 
0119 std::vector<reco::BasicCluster> HGCalImagingAlgo::getClusters(bool doSharing) {
0120   reco::CaloID caloID = reco::CaloID::DET_HGCAL_ENDCAP;
0121   std::vector<std::pair<DetId, float>> thisCluster;
0122   for (auto &clsOnLayer : layerClustersPerLayer_) {
0123     for (unsigned int i = 0; i < clsOnLayer.size(); ++i) {
0124       double energy = 0;
0125       Point position;
0126 
0127       //Will save the maximum density hit of the cluster
0128       size_t rsmax = max_index(clsOnLayer[i]);
0129 
0130       if (doSharing) {
0131         std::vector<unsigned> seeds = findLocalMaximaInCluster(clsOnLayer[i]);
0132         // sharing found seeds.size() sub-cluster seeds in cluster i
0133 
0134         std::vector<std::vector<double>> fractions;
0135         // first pass can have noise it in
0136         shareEnergy(clsOnLayer[i], seeds, fractions);
0137 
0138         // reset and run second pass after vetoing seeds
0139         // that result in trivial clusters (less than 2 effective cells)
0140 
0141         for (unsigned isub = 0; isub < fractions.size(); ++isub) {
0142           double effective_hits = 0.0;
0143           double energy = calculateEnergyWithFraction(clsOnLayer[i], fractions[isub]);
0144           Point position = calculatePositionWithFraction(clsOnLayer[i], fractions[isub]);
0145 
0146           for (unsigned ihit = 0; ihit < fractions[isub].size(); ++ihit) {
0147             const double fraction = fractions[isub][ihit];
0148             if (fraction > 1e-7) {
0149               effective_hits += fraction;
0150               thisCluster.emplace_back(clsOnLayer[i][ihit].data.detid, fraction);
0151             }
0152           }
0153 
0154           if (verbosity_ < pINFO) {
0155             std::cout << "\t******** NEW CLUSTER (SHARING) ********" << std::endl;
0156             std::cout << "\tEff. No. of cells = " << effective_hits << std::endl;
0157             std::cout << "\t     Energy       = " << energy << std::endl;
0158             std::cout << "\t     Phi          = " << position.phi() << std::endl;
0159             std::cout << "\t     Eta          = " << position.eta() << std::endl;
0160             std::cout << "\t*****************************" << std::endl;
0161           }
0162           clusters_v_.emplace_back(energy, position, caloID, thisCluster, algoId_);
0163           if (!clusters_v_.empty()) {
0164             clusters_v_.back().setSeed(clsOnLayer[i][rsmax].data.detid);
0165           }
0166           thisCluster.clear();
0167         }
0168       } else {
0169         position = calculatePosition(clsOnLayer[i]);  // energy-weighted position
0170         //   std::vector< KDNode >::iterator it;
0171         for (auto &it : clsOnLayer[i]) {
0172           energy += it.data.isHalo ? 0. : it.data.weight;
0173           // use fraction to store whether this is a Halo hit or not
0174           thisCluster.emplace_back(it.data.detid, (it.data.isHalo ? 0.f : 1.f));
0175         }
0176         if (verbosity_ < pINFO) {
0177           std::cout << "******** NEW CLUSTER (HGCIA) ********" << std::endl;
0178           std::cout << "Index          " << i << std::endl;
0179           std::cout << "No. of cells = " << clsOnLayer[i].size() << std::endl;
0180           std::cout << "     Energy     = " << energy << std::endl;
0181           std::cout << "     Phi        = " << position.phi() << std::endl;
0182           std::cout << "     Eta        = " << position.eta() << std::endl;
0183           std::cout << "*****************************" << std::endl;
0184         }
0185         clusters_v_.emplace_back(energy, position, caloID, thisCluster, algoId_);
0186         if (!clusters_v_.empty()) {
0187           clusters_v_.back().setSeed(clsOnLayer[i][rsmax].data.detid);
0188         }
0189         thisCluster.clear();
0190       }
0191     }
0192   }
0193   return clusters_v_;
0194 }
0195 
0196 math::XYZPoint HGCalImagingAlgo::calculatePosition(std::vector<KDNode> &v) const {
0197   float total_weight = 0.f;
0198   float x = 0.f;
0199   float y = 0.f;
0200 
0201   unsigned int v_size = v.size();
0202   unsigned int maxEnergyIndex = 0;
0203   float maxEnergyValue = 0;
0204 
0205   // loop over hits in cluster candidate
0206   // determining the maximum energy hit
0207   for (unsigned int i = 0; i < v_size; i++) {
0208     if (v[i].data.weight > maxEnergyValue) {
0209       maxEnergyValue = v[i].data.weight;
0210       maxEnergyIndex = i;
0211     }
0212   }
0213 
0214   // Si cell or Scintillator. Used to set approach and parameters
0215   int thick = rhtools_.getSiThickIndex(v[maxEnergyIndex].data.detid);
0216 
0217   // for hits within positionDeltaRho_c_ from maximum energy hit
0218   // build up weight for energy-weighted position
0219   // and save corresponding hits indices
0220   std::vector<unsigned int> innerIndices;
0221   for (unsigned int i = 0; i < v_size; i++) {
0222     if (thick == -1 || distance2(v[i].data, v[maxEnergyIndex].data) < positionDeltaRho_c_[thick]) {
0223       innerIndices.push_back(i);
0224 
0225       float rhEnergy = v[i].data.weight;
0226       total_weight += rhEnergy;
0227       // just fill x, y for scintillator
0228       // for Si it is overwritten later anyway
0229       if (thick == -1) {
0230         x += v[i].data.x * rhEnergy;
0231         y += v[i].data.y * rhEnergy;
0232       }
0233     }
0234   }
0235   // just loop on reduced vector of interesting indices
0236   // to compute log weighting
0237   if (thick != -1 && total_weight != 0.) {  // Silicon case
0238     float total_weight_log = 0.f;
0239     float x_log = 0.f;
0240     float y_log = 0.f;
0241     for (auto idx : innerIndices) {
0242       float rhEnergy = v[idx].data.weight;
0243       if (rhEnergy == 0.)
0244         continue;
0245       float Wi = std::max(thresholdW0_[thick] + log(rhEnergy / total_weight), 0.);
0246       x_log += v[idx].data.x * Wi;
0247       y_log += v[idx].data.y * Wi;
0248       total_weight_log += Wi;
0249     }
0250     total_weight = total_weight_log;
0251     x = x_log;
0252     y = y_log;
0253   }
0254 
0255   if (total_weight != 0.) {
0256     auto inv_tot_weight = 1. / total_weight;
0257     return math::XYZPoint(x * inv_tot_weight, y * inv_tot_weight, v[maxEnergyIndex].data.z);
0258   }
0259   return math::XYZPoint(0, 0, 0);
0260 }
0261 
0262 double HGCalImagingAlgo::calculateLocalDensity(std::vector<KDNode> &nd, KDTree &lp, const unsigned int layer) const {
0263   double maxdensity = 0.;
0264   float delta_c;  // maximum search distance (critical distance) for local
0265                   // density calculation
0266   if (layer <= lastLayerEE_)
0267     delta_c = vecDeltas_[0];
0268   else if (layer < firstLayerBH_)
0269     delta_c = vecDeltas_[1];
0270   else
0271     delta_c = vecDeltas_[2];
0272 
0273   // for each node calculate local density rho and store it
0274   for (unsigned int i = 0; i < nd.size(); ++i) {
0275     // speec up search by looking within +/- delta_c window only
0276     KDTreeBox search_box(
0277         nd[i].dims[0] - delta_c, nd[i].dims[0] + delta_c, nd[i].dims[1] - delta_c, nd[i].dims[1] + delta_c);
0278     std::vector<Hexel> found;
0279     lp.search(search_box, found);
0280     const unsigned int found_size = found.size();
0281     for (unsigned int j = 0; j < found_size; j++) {
0282       if (distance(nd[i].data, found[j]) < delta_c) {
0283         nd[i].data.rho += found[j].weight;
0284         maxdensity = std::max(maxdensity, nd[i].data.rho);
0285       }
0286     }  // end loop found
0287   }    // end loop nodes
0288   return maxdensity;
0289 }
0290 
0291 double HGCalImagingAlgo::calculateDistanceToHigher(std::vector<KDNode> &nd) const {
0292   // sort vector of Hexels by decreasing local density
0293   std::vector<size_t> rs = sorted_indices(nd);
0294 
0295   double maxdensity = 0.0;
0296   int nearestHigher = -1;
0297 
0298   if (!rs.empty())
0299     maxdensity = nd[rs[0]].data.rho;
0300   else
0301     return maxdensity;  // there are no hits
0302   double dist2 = 0.;
0303   // start by setting delta for the highest density hit to
0304   // the most distant hit - this is a convention
0305 
0306   for (auto &j : nd) {
0307     double tmp = distance2(nd[rs[0]].data, j.data);
0308     if (tmp > dist2)
0309       dist2 = tmp;
0310   }
0311   nd[rs[0]].data.delta = std::sqrt(dist2);
0312   nd[rs[0]].data.nearestHigher = nearestHigher;
0313 
0314   // now we save the largest distance as a starting point
0315   const double max_dist2 = dist2;
0316   const unsigned int nd_size = nd.size();
0317 
0318   for (unsigned int oi = 1; oi < nd_size; ++oi) {  // start from second-highest density
0319     dist2 = max_dist2;
0320     unsigned int i = rs[oi];
0321     // we only need to check up to oi since hits
0322     // are ordered by decreasing density
0323     // and all points coming BEFORE oi are guaranteed to have higher rho
0324     // and the ones AFTER to have lower rho
0325     for (unsigned int oj = 0; oj < oi; ++oj) {
0326       unsigned int j = rs[oj];
0327       // Limit the search box
0328       if ((nd[i].data.x - nd[j].data.x) * (nd[i].data.x - nd[j].data.x) > dist2)
0329         continue;
0330       if ((nd[i].data.y - nd[j].data.y) * (nd[i].data.y - nd[j].data.y) > dist2)
0331         continue;
0332       double tmp = distance2(nd[i].data, nd[j].data);
0333       if (tmp <= dist2) {  // this "<=" instead of "<" addresses the (rare) case
0334                            // when there are only two hits
0335         dist2 = tmp;
0336         nearestHigher = j;
0337       }
0338     }
0339     nd[i].data.delta = std::sqrt(dist2);
0340     nd[i].data.nearestHigher = nearestHigher;  // this uses the original unsorted hitlist
0341   }
0342   return maxdensity;
0343 }
0344 int HGCalImagingAlgo::findAndAssignClusters(std::vector<KDNode> &nd,
0345                                             KDTree &lp,
0346                                             double maxdensity,
0347                                             KDTreeBox<2> &bounds,
0348                                             const unsigned int layer,
0349                                             std::vector<std::vector<KDNode>> &clustersOnLayer) const {
0350   // this is called once per layer and endcap...
0351   // so when filling the cluster temporary vector of Hexels we resize each time
0352   // by the number  of clusters found. This is always equal to the number of
0353   // cluster centers...
0354 
0355   unsigned int nClustersOnLayer = 0;
0356   float delta_c;  // critical distance
0357   if (layer <= lastLayerEE_)
0358     delta_c = vecDeltas_[0];
0359   else if (layer < firstLayerBH_)
0360     delta_c = vecDeltas_[1];
0361   else
0362     delta_c = vecDeltas_[2];
0363 
0364   std::vector<size_t> rs = sorted_indices(nd);  // indices sorted by decreasing rho
0365   std::vector<size_t> ds = sort_by_delta(nd);   // sort in decreasing distance to higher
0366 
0367   const unsigned int nd_size = nd.size();
0368   for (unsigned int i = 0; i < nd_size; ++i) {
0369     if (nd[ds[i]].data.delta < delta_c)
0370       break;  // no more cluster centers to be looked at
0371     if (dependSensor_) {
0372       float rho_c = kappa_ * nd[ds[i]].data.sigmaNoise;
0373       if (nd[ds[i]].data.rho < rho_c)
0374         continue;  // set equal to kappa times noise threshold
0375 
0376     } else if (nd[ds[i]].data.rho * kappa_ < maxdensity)
0377       continue;
0378 
0379     nd[ds[i]].data.clusterIndex = nClustersOnLayer;
0380     if (verbosity_ < pINFO) {
0381       std::cout << "Adding new cluster with index " << nClustersOnLayer << std::endl;
0382       std::cout << "Cluster center is hit " << ds[i] << std::endl;
0383     }
0384     nClustersOnLayer++;
0385   }
0386 
0387   // at this point nClustersOnLayer is equal to the number of cluster centers -
0388   // if it is zero we are  done
0389   if (nClustersOnLayer == 0)
0390     return nClustersOnLayer;
0391 
0392   // assign remaining points to clusters, using the nearestHigher set from
0393   // previous step (always set except
0394   // for top density hit that is skipped...)
0395   for (unsigned int oi = 1; oi < nd_size; ++oi) {
0396     unsigned int i = rs[oi];
0397     int ci = nd[i].data.clusterIndex;
0398     if (ci == -1) {  // clusterIndex is initialised with -1 if not yet used in cluster
0399       nd[i].data.clusterIndex = nd[nd[i].data.nearestHigher].data.clusterIndex;
0400     }
0401   }
0402 
0403   // make room in the temporary cluster vector for the additional clusterIndex
0404   // clusters
0405   // from this layer
0406   if (verbosity_ < pINFO) {
0407     std::cout << "resizing cluster vector by " << nClustersOnLayer << std::endl;
0408   }
0409   clustersOnLayer.resize(nClustersOnLayer);
0410 
0411   // assign points closer than dc to other clusters to border region
0412   // and find critical border density
0413   std::vector<double> rho_b(nClustersOnLayer, 0.);
0414   lp.clear();
0415   lp.build(nd, bounds);
0416   // now loop on all hits again :( and check: if there are hits from another
0417   // cluster within d_c -> flag as border hit
0418   for (unsigned int i = 0; i < nd_size; ++i) {
0419     int ci = nd[i].data.clusterIndex;
0420     bool flag_isolated = true;
0421     if (ci != -1) {
0422       KDTreeBox search_box(
0423           nd[i].dims[0] - delta_c, nd[i].dims[0] + delta_c, nd[i].dims[1] - delta_c, nd[i].dims[1] + delta_c);
0424       std::vector<Hexel> found;
0425       lp.search(search_box, found);
0426 
0427       const unsigned int found_size = found.size();
0428       for (unsigned int j = 0; j < found_size; j++) {  // start from 0 here instead of 1
0429         // check if the hit is not within d_c of another cluster
0430         if (found[j].clusterIndex != -1) {
0431           float dist = distance(found[j], nd[i].data);
0432           if (dist < delta_c && found[j].clusterIndex != ci) {
0433             // in which case we assign it to the border
0434             nd[i].data.isBorder = true;
0435             break;
0436           }
0437           // because we are using two different containers, we have to make sure
0438           // that we don't unflag the
0439           // hit when it finds *itself* closer than delta_c
0440           if (dist < delta_c && dist != 0. && found[j].clusterIndex == ci) {
0441             // in this case it is not an isolated hit
0442             // the dist!=0 is because the hit being looked at is also inside the
0443             // search box and at dist==0
0444             flag_isolated = false;
0445           }
0446         }
0447       }
0448       if (flag_isolated)
0449         nd[i].data.isBorder = true;  // the hit is more than delta_c from any of its brethren
0450     }
0451     // check if this border hit has density larger than the current rho_b and
0452     // update
0453     if (nd[i].data.isBorder && rho_b[ci] < nd[i].data.rho)
0454       rho_b[ci] = nd[i].data.rho;
0455   }  // end loop all hits
0456 
0457   // flag points in cluster with density < rho_b as halo points, then fill the
0458   // cluster vector
0459   for (unsigned int i = 0; i < nd_size; ++i) {
0460     int ci = nd[i].data.clusterIndex;
0461     if (ci != -1) {
0462       if (nd[i].data.rho <= rho_b[ci])
0463         nd[i].data.isHalo = true;
0464       clustersOnLayer[ci].push_back(nd[i]);
0465       if (verbosity_ < pINFO) {
0466         std::cout << "Pushing hit " << i << " into cluster with index " << ci << std::endl;
0467       }
0468     }
0469   }
0470 
0471   // prepare the offset for the next layer if there is one
0472   if (verbosity_ < pINFO) {
0473     std::cout << "moving cluster offset by " << nClustersOnLayer << std::endl;
0474   }
0475   return nClustersOnLayer;
0476 }
0477 
0478 // find local maxima within delta_c, marking the indices in the cluster
0479 std::vector<unsigned> HGCalImagingAlgo::findLocalMaximaInCluster(const std::vector<KDNode> &cluster) {
0480   std::vector<unsigned> result;
0481   std::vector<bool> seed(cluster.size(), true);
0482   float delta_c = 2.;
0483 
0484   for (unsigned i = 0; i < cluster.size(); ++i) {
0485     for (unsigned j = 0; j < cluster.size(); ++j) {
0486       if (i != j and distance(cluster[i].data, cluster[j].data) < delta_c) {
0487         if (cluster[i].data.weight < cluster[j].data.weight) {
0488           seed[i] = false;
0489           break;
0490         }
0491       }
0492     }
0493   }
0494 
0495   for (unsigned i = 0; i < cluster.size(); ++i) {
0496     if (seed[i] && cluster[i].data.weight > 5e-4) {
0497       // seed at i with energy cluster[i].weight
0498       result.push_back(i);
0499     }
0500   }
0501 
0502   // Found result.size() sub-clusters in input cluster of length cluster.size()
0503 
0504   return result;
0505 }
0506 
0507 math::XYZPoint HGCalImagingAlgo::calculatePositionWithFraction(const std::vector<KDNode> &hits,
0508                                                                const std::vector<double> &fractions) {
0509   double norm(0.0), x(0.0), y(0.0), z(0.0);
0510   for (unsigned i = 0; i < hits.size(); ++i) {
0511     const double weight = fractions[i] * hits[i].data.weight;
0512     norm += weight;
0513     x += weight * hits[i].data.x;
0514     y += weight * hits[i].data.y;
0515     z += weight * hits[i].data.z;
0516   }
0517   math::XYZPoint result(x, y, z);
0518   result /= norm;
0519   return result;
0520 }
0521 
0522 double HGCalImagingAlgo::calculateEnergyWithFraction(const std::vector<KDNode> &hits,
0523                                                      const std::vector<double> &fractions) {
0524   double result = 0.0;
0525   for (unsigned i = 0; i < hits.size(); ++i) {
0526     result += fractions[i] * hits[i].data.weight;
0527   }
0528   return result;
0529 }
0530 
0531 void HGCalImagingAlgo::shareEnergy(const std::vector<KDNode> &incluster,
0532                                    const std::vector<unsigned> &seeds,
0533                                    std::vector<std::vector<double>> &outclusters) {
0534   std::vector<bool> isaseed(incluster.size(), false);
0535   outclusters.clear();
0536   outclusters.resize(seeds.size());
0537   std::vector<Point> centroids(seeds.size());
0538   std::vector<double> energies(seeds.size());
0539 
0540   if (seeds.size() == 1) {  // short circuit the case of a lone cluster
0541     outclusters[0].clear();
0542     outclusters[0].resize(incluster.size(), 1.0);
0543     return;
0544   }
0545 
0546   // saving seeds
0547 
0548   // create quick seed lookup
0549   for (unsigned i = 0; i < seeds.size(); ++i) {
0550     isaseed[seeds[i]] = true;
0551   }
0552 
0553   // initialize clusters to be shared
0554   // centroids start off at seed positions
0555   // seeds always have fraction 1.0, to stabilize fit
0556   // initializing fit
0557   for (unsigned i = 0; i < seeds.size(); ++i) {
0558     outclusters[i].resize(incluster.size(), 0.0);
0559     for (unsigned j = 0; j < incluster.size(); ++j) {
0560       if (j == seeds[i]) {
0561         outclusters[i][j] = 1.0;
0562         centroids[i] = math::XYZPoint(incluster[j].data.x, incluster[j].data.y, incluster[j].data.z);
0563         energies[i] = incluster[j].data.weight;
0564       }
0565     }
0566   }
0567 
0568   // run the fit while we are less than max iterations, and clusters are still
0569   // moving
0570   const double minFracTot = 1e-20;
0571   unsigned iter = 0;
0572   const unsigned iterMax = 50;
0573   double diff = std::numeric_limits<double>::max();
0574   const double stoppingTolerance = 1e-8;
0575   const auto numberOfSeeds = seeds.size();
0576   auto toleranceScaling = numberOfSeeds > 2 ? (numberOfSeeds - 1) * (numberOfSeeds - 1) : 1;
0577   std::vector<Point> prevCentroids;
0578   std::vector<double> frac(numberOfSeeds), dist2(numberOfSeeds);
0579   while (iter++ < iterMax && diff > stoppingTolerance * toleranceScaling) {
0580     for (unsigned i = 0; i < incluster.size(); ++i) {
0581       const Hexel &ihit = incluster[i].data;
0582       double fracTot(0.0);
0583       for (unsigned j = 0; j < numberOfSeeds; ++j) {
0584         double fraction = 0.0;
0585         double d2 = (std::pow(ihit.x - centroids[j].x(), 2.0) + std::pow(ihit.y - centroids[j].y(), 2.0) +
0586                      std::pow(ihit.z - centroids[j].z(), 2.0)) /
0587                     sigma2_;
0588         dist2[j] = d2;
0589         // now we set the fractions up based on hit type
0590         if (i == seeds[j]) {  // this cluster's seed
0591           fraction = 1.0;
0592         } else if (isaseed[i]) {
0593           fraction = 0.0;
0594         } else {
0595           fraction = energies[j] * std::exp(-0.5 * d2);
0596         }
0597         fracTot += fraction;
0598         frac[j] = fraction;
0599       }
0600       // now that we have calculated all fractions for all hits
0601       // assign the new fractions
0602       for (unsigned j = 0; j < numberOfSeeds; ++j) {
0603         if (fracTot > minFracTot || (i == seeds[j] && fracTot > 0.0)) {
0604           outclusters[j][i] = frac[j] / fracTot;
0605         } else {
0606           outclusters[j][i] = 0.0;
0607         }
0608       }
0609     }
0610 
0611     // save previous centroids
0612     prevCentroids = std::move(centroids);
0613     // finally update the position of the centroids from the last iteration
0614     centroids.resize(numberOfSeeds);
0615     double diff2 = 0.0;
0616     for (unsigned i = 0; i < numberOfSeeds; ++i) {
0617       centroids[i] = calculatePositionWithFraction(incluster, outclusters[i]);
0618       energies[i] = calculateEnergyWithFraction(incluster, outclusters[i]);
0619       // calculate convergence parameters
0620       const double delta2 = (prevCentroids[i] - centroids[i]).perp2();
0621       diff2 = std::max(delta2, diff2);
0622     }
0623     // update convergance parameter outside loop
0624     diff = std::sqrt(diff2);
0625   }
0626 }
0627 
0628 void HGCalImagingAlgo::computeThreshold() {
0629   // To support the TDR geometry and also the post-TDR one (v9 onwards), we
0630   // need to change the logic of the vectors containing signal to noise and
0631   // thresholds. The first 3 indices will keep on addressing the different
0632   // thicknesses of the Silicon detectors, while the last one, number 3 (the
0633   // fourth) will address the Scintillators. This change will support both
0634   // geometries at the same time.
0635 
0636   if (initialized_)
0637     return;  // only need to calculate thresholds once
0638 
0639   initialized_ = true;
0640 
0641   std::vector<double> dummy;
0642   const unsigned maxNumberOfThickIndices = 3;
0643   dummy.resize(maxNumberOfThickIndices + 1, 0);  // +1 to accomodate for the Scintillators
0644   thresholds_.resize(maxlayer_, dummy);
0645   sigmaNoise_.resize(maxlayer_, dummy);
0646 
0647   for (unsigned ilayer = 1; ilayer <= maxlayer_; ++ilayer) {
0648     for (unsigned ithick = 0; ithick < maxNumberOfThickIndices; ++ithick) {
0649       float sigmaNoise = 0.001f * fcPerEle_ * nonAgedNoises_[ithick] * dEdXweights_[ilayer] /
0650                          (fcPerMip_[ithick] * thicknessCorrection_[ithick]);
0651       thresholds_[ilayer - 1][ithick] = sigmaNoise * ecut_;
0652       sigmaNoise_[ilayer - 1][ithick] = sigmaNoise;
0653     }
0654     float scintillators_sigmaNoise = 0.001f * noiseMip_ * dEdXweights_[ilayer];
0655     thresholds_[ilayer - 1][maxNumberOfThickIndices] = ecut_ * scintillators_sigmaNoise;
0656     sigmaNoise_[ilayer - 1][maxNumberOfThickIndices] = scintillators_sigmaNoise;
0657   }
0658 }
0659 
0660 void HGCalImagingAlgo::setDensity(const std::vector<KDNode> &nd) {
0661   // for each node calculate local density rho and store it
0662   for (auto &i : nd) {
0663     density_[i.data.detid] = i.data.rho;
0664   }  // end loop nodes
0665 }
0666 
0667 //Density
0668 Density HGCalImagingAlgo::getDensity() { return density_; }