File indexing completed on 2023-03-17 11:19:00
0001 #include "RecoLocalCalo/HGCalRecProducers/interface/HGCalImagingAlgo.h"
0002
0003
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
0028
0029 if (dependSensor_) {
0030
0031
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
0042
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;
0054
0055 }
0056 if (!dependSensor_ && hgrh.energy() < ecut_)
0057 continue;
0058
0059
0060
0061 layer += int(rhtools_.zside(detid) > 0) * (maxlayer_ + 1);
0062
0063
0064 bool isHalf = rhtools_.isHalfCell(detid);
0065 const GlobalPoint position(rhtools_.getPosition(detid));
0066
0067
0068
0069 points_[layer].emplace_back(
0070 Hexel(hgrh, detid, isHalf, sigmaNoise, thickness, &rhtools_), position.x(), position.y());
0071
0072
0073
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 }
0088 }
0089
0090
0091
0092
0093 void HGCalImagingAlgo::makeClusters() {
0094 layerClustersPerLayer_.resize(2 * maxlayer_ + 2);
0095
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;
0104
0105 double maxdensity = calculateLocalDensity(points_[i], hit_kdtree, actualLayer);
0106
0107
0108
0109 calculateDistanceToHigher(points_[i]);
0110 findAndAssignClusters(points_[i], hit_kdtree, maxdensity, bounds, actualLayer, layerClustersPerLayer_[i]);
0111 });
0112 });
0113
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
0128 size_t rsmax = max_index(clsOnLayer[i]);
0129
0130 if (doSharing) {
0131 std::vector<unsigned> seeds = findLocalMaximaInCluster(clsOnLayer[i]);
0132
0133
0134 std::vector<std::vector<double>> fractions;
0135
0136 shareEnergy(clsOnLayer[i], seeds, fractions);
0137
0138
0139
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]);
0170
0171 for (auto &it : clsOnLayer[i]) {
0172 energy += it.data.isHalo ? 0. : it.data.weight;
0173
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
0206
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
0215 int thick = rhtools_.getSiThickIndex(v[maxEnergyIndex].data.detid);
0216
0217
0218
0219
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
0228
0229 if (thick == -1) {
0230 x += v[i].data.x * rhEnergy;
0231 y += v[i].data.y * rhEnergy;
0232 }
0233 }
0234 }
0235
0236
0237 if (thick != -1 && total_weight != 0.) {
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;
0265
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
0274 for (unsigned int i = 0; i < nd.size(); ++i) {
0275
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 }
0287 }
0288 return maxdensity;
0289 }
0290
0291 double HGCalImagingAlgo::calculateDistanceToHigher(std::vector<KDNode> &nd) const {
0292
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;
0302 double dist2 = 0.;
0303
0304
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
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) {
0319 dist2 = max_dist2;
0320 unsigned int i = rs[oi];
0321
0322
0323
0324
0325 for (unsigned int oj = 0; oj < oi; ++oj) {
0326 unsigned int j = rs[oj];
0327
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) {
0334
0335 dist2 = tmp;
0336 nearestHigher = j;
0337 }
0338 }
0339 nd[i].data.delta = std::sqrt(dist2);
0340 nd[i].data.nearestHigher = nearestHigher;
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
0351
0352
0353
0354
0355 unsigned int nClustersOnLayer = 0;
0356 float delta_c;
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);
0365 std::vector<size_t> ds = sort_by_delta(nd);
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;
0371 if (dependSensor_) {
0372 float rho_c = kappa_ * nd[ds[i]].data.sigmaNoise;
0373 if (nd[ds[i]].data.rho < rho_c)
0374 continue;
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
0388
0389 if (nClustersOnLayer == 0)
0390 return nClustersOnLayer;
0391
0392
0393
0394
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) {
0399 nd[i].data.clusterIndex = nd[nd[i].data.nearestHigher].data.clusterIndex;
0400 }
0401 }
0402
0403
0404
0405
0406 if (verbosity_ < pINFO) {
0407 std::cout << "resizing cluster vector by " << nClustersOnLayer << std::endl;
0408 }
0409 clustersOnLayer.resize(nClustersOnLayer);
0410
0411
0412
0413 std::vector<double> rho_b(nClustersOnLayer, 0.);
0414 lp.clear();
0415 lp.build(nd, bounds);
0416
0417
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++) {
0429
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
0434 nd[i].data.isBorder = true;
0435 break;
0436 }
0437
0438
0439
0440 if (dist < delta_c && dist != 0. && found[j].clusterIndex == ci) {
0441
0442
0443
0444 flag_isolated = false;
0445 }
0446 }
0447 }
0448 if (flag_isolated)
0449 nd[i].data.isBorder = true;
0450 }
0451
0452
0453 if (nd[i].data.isBorder && rho_b[ci] < nd[i].data.rho)
0454 rho_b[ci] = nd[i].data.rho;
0455 }
0456
0457
0458
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
0472 if (verbosity_ < pINFO) {
0473 std::cout << "moving cluster offset by " << nClustersOnLayer << std::endl;
0474 }
0475 return nClustersOnLayer;
0476 }
0477
0478
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
0498 result.push_back(i);
0499 }
0500 }
0501
0502
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) {
0541 outclusters[0].clear();
0542 outclusters[0].resize(incluster.size(), 1.0);
0543 return;
0544 }
0545
0546
0547
0548
0549 for (unsigned i = 0; i < seeds.size(); ++i) {
0550 isaseed[seeds[i]] = true;
0551 }
0552
0553
0554
0555
0556
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
0569
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
0590 if (i == seeds[j]) {
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
0601
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
0612 prevCentroids = std::move(centroids);
0613
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
0620 const double delta2 = (prevCentroids[i] - centroids[i]).perp2();
0621 diff2 = std::max(delta2, diff2);
0622 }
0623
0624 diff = std::sqrt(diff2);
0625 }
0626 }
0627
0628 void HGCalImagingAlgo::computeThreshold() {
0629
0630
0631
0632
0633
0634
0635
0636 if (initialized_)
0637 return;
0638
0639 initialized_ = true;
0640
0641 std::vector<double> dummy;
0642 const unsigned maxNumberOfThickIndices = 3;
0643 dummy.resize(maxNumberOfThickIndices + 1, 0);
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
0662 for (auto &i : nd) {
0663 density_[i.data.detid] = i.data.rho;
0664 }
0665 }
0666
0667
0668 Density HGCalImagingAlgo::getDensity() { return density_; }