Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-07-03 04:18:12 

0001 #include "RecoLocalCalo/HGCalRecProducers/plugins/HGCalCLUEAlgo.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 #include <limits>
0015 #include "DataFormats/DetId/interface/DetId.h"
0016 
0017 using namespace hgcal_clustering;
0018 
0019 template <typename T, typename STRATEGY>
0020 void HGCalCLUEAlgoT<T, STRATEGY>::getEventSetupPerAlgorithm(const edm::EventSetup& es) {
0021   cells_.clear();
0022   numberOfClustersPerLayer_.clear();
0023   cells_.resize(2 * (maxlayer_ + 1));
0024   numberOfClustersPerLayer_.resize(2 * (maxlayer_ + 1), 0);
0025 }
0026 
0027 template <typename T, typename STRATEGY>
0028 void HGCalCLUEAlgoT<T, STRATEGY>::populate(const HGCRecHitCollection& hits) {
0029   // loop over all hits and create the Hexel structure, skip energies below ecut
0030   if (dependSensor_) {
0031     // for each layer and wafer calculate the thresholds (sigmaNoise and energy)
0032     // once
0033     computeThreshold();
0034   }
0035 
0036   for (unsigned int i = 0; i < hits.size(); ++i) {
0037     const HGCRecHit& hgrh = hits[i];
0038     DetId detid = hgrh.detid();
0039     unsigned int layerOnSide = (rhtools_.getLayerWithOffset(detid) - 1);
0040 
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       int thickness_index = rhtools_.getSiThickIndex(detid);
0046       if (thickness_index == -1)
0047         thickness_index = maxNumberOfThickIndices_;
0048 
0049       double storedThreshold = thresholds_[layerOnSide][thickness_index];
0050       if (detid.det() == DetId::HGCalHSi || detid.subdetId() == HGCHEF) {
0051         storedThreshold = thresholds_[layerOnSide][thickness_index + deltasi_index_regemfac_];
0052       }
0053       sigmaNoise = v_sigmaNoise_[layerOnSide][thickness_index];
0054 
0055       if (hgrh.energy() < storedThreshold)
0056         continue;  // this sets the ZS threshold at ecut times the sigma noise
0057                    // for the sensor
0058     }
0059     if (!dependSensor_ && hgrh.energy() < ecut_)
0060       continue;
0061     const GlobalPoint position(rhtools_.getPosition(detid));
0062     int offset = ((rhtools_.zside(detid) + 1) >> 1) * maxlayer_;
0063     int layer = layerOnSide + offset;
0064     // setting the layer position only once per layer
0065     if (cells_[layer].layerDim3 == std::numeric_limits<float>::infinity())
0066       cells_[layer].layerDim3 = position.z();
0067 
0068     cells_[layer].detid.emplace_back(detid);
0069     if constexpr (std::is_same_v<STRATEGY, HGCalScintillatorStrategy>) {
0070       cells_[layer].dim1.emplace_back(position.eta());
0071       cells_[layer].dim2.emplace_back(position.phi());
0072     }  // else, isSilicon == true and eta phi values will not be used
0073     else {
0074       cells_[layer].dim1.emplace_back(position.x());
0075       cells_[layer].dim2.emplace_back(position.y());
0076     }
0077     cells_[layer].weight.emplace_back(hgrh.energy());
0078     cells_[layer].sigmaNoise.emplace_back(sigmaNoise);
0079   }
0080 }
0081 
0082 template <typename T, typename STRATEGY>
0083 void HGCalCLUEAlgoT<T, STRATEGY>::prepareDataStructures(unsigned int l) {
0084   auto cellsSize = cells_[l].detid.size();
0085   cells_[l].rho.resize(cellsSize, 0.f);
0086   cells_[l].delta.resize(cellsSize, 9999999);
0087   cells_[l].nearestHigher.resize(cellsSize, -1);
0088   cells_[l].clusterIndex.resize(cellsSize, -1);
0089   cells_[l].followers.resize(cellsSize);
0090   cells_[l].isSeed.resize(cellsSize, false);
0091 }
0092 
0093 // Create a vector of Hexels associated to one cluster from a collection of
0094 // HGCalRecHits - this can be used directly to make the final cluster list -
0095 // this method can be invoked multiple times for the same event with different
0096 // input (reset should be called between events)
0097 template <typename T, typename STRATEGY>
0098 void HGCalCLUEAlgoT<T, STRATEGY>::makeClusters() {
0099   // assign all hits in each layer to a cluster core
0100   tbb::this_task_arena::isolate([&] {
0101     tbb::parallel_for(size_t(0), size_t(2 * maxlayer_ + 2), [&](size_t i) {
0102       prepareDataStructures(i);
0103       T lt;
0104       lt.clear();
0105       lt.fill(cells_[i].dim1, cells_[i].dim2);
0106 
0107       float delta;
0108       if constexpr (std::is_same_v<STRATEGY, HGCalSiliconStrategy>) {
0109         // maximum search distance (critical distance) for local density calculation
0110         float delta_c;
0111         if (i % maxlayer_ < lastLayerEE_)
0112           delta_c = vecDeltas_[0];
0113         else if (i % maxlayer_ < (firstLayerBH_ - 1))
0114           delta_c = vecDeltas_[1];
0115         else
0116           delta_c = vecDeltas_[2];
0117         delta = delta_c;
0118       } else {
0119         float delta_r = vecDeltas_[3];
0120         delta = delta_r;
0121       }
0122       LogDebug("HGCalCLUEAlgo") << "maxlayer: " << maxlayer_ << " lastLayerEE: " << lastLayerEE_
0123                                 << " firstLayerBH: " << firstLayerBH_ << "\n";
0124 
0125       calculateLocalDensity(lt, i, delta);
0126       calculateDistanceToHigher(lt, i, delta);
0127       numberOfClustersPerLayer_[i] = findAndAssignClusters(i, delta);
0128     });
0129   });
0130 #if DEBUG_CLUSTERS_ALPAKA
0131   hgcalUtils::DumpLegacySoA dumperLegacySoA;
0132   dumperLegacySoA.dumpInfos(cells_, moduleType_);
0133 #endif
0134 }
0135 
0136 template <typename T, typename STRATEGY>
0137 std::vector<reco::BasicCluster> HGCalCLUEAlgoT<T, STRATEGY>::getClusters(bool) {
0138   std::vector<int> offsets(numberOfClustersPerLayer_.size(), 0);
0139 
0140   int maxClustersOnLayer = numberOfClustersPerLayer_[0];
0141 
0142   for (unsigned layerId = 1; layerId < offsets.size(); ++layerId) {
0143     offsets[layerId] = offsets[layerId - 1] + numberOfClustersPerLayer_[layerId - 1];
0144 
0145     maxClustersOnLayer = std::max(maxClustersOnLayer, numberOfClustersPerLayer_[layerId]);
0146   }
0147 
0148   auto totalNumberOfClusters = offsets.back() + numberOfClustersPerLayer_.back();
0149   clusters_v_.resize(totalNumberOfClusters);
0150   std::vector<std::vector<int>> cellsIdInCluster;
0151   cellsIdInCluster.reserve(maxClustersOnLayer);
0152 
0153   for (unsigned int layerId = 0; layerId < 2 * maxlayer_ + 2; ++layerId) {
0154     cellsIdInCluster.resize(numberOfClustersPerLayer_[layerId]);
0155     auto& cellsOnLayer = cells_[layerId];
0156     unsigned int numberOfCells = cellsOnLayer.detid.size();
0157     auto firstClusterIdx = offsets[layerId];
0158 
0159     for (unsigned int i = 0; i < numberOfCells; ++i) {
0160       auto clusterIndex = cellsOnLayer.clusterIndex[i];
0161       if (clusterIndex != -1)
0162         cellsIdInCluster[clusterIndex].push_back(i);
0163     }
0164 
0165     std::vector<std::pair<DetId, float>> thisCluster;
0166 
0167     for (auto& cl : cellsIdInCluster) {
0168       float maxEnergyValue = std::numeric_limits<float>::min();
0169       int maxEnergyCellIndex = -1;
0170       DetId maxEnergyDetId;
0171       float energy = 0.f;
0172       int seedDetId = -1;
0173       float x = 0.f;
0174       float y = 0.f;
0175       float z = cellsOnLayer.layerDim3;
0176       // TODO Felice: maybe use the seed for the position calculation
0177       for (auto cellIdx : cl) {
0178         energy += cellsOnLayer.weight[cellIdx];
0179         if (cellsOnLayer.weight[cellIdx] > maxEnergyValue) {
0180           maxEnergyValue = cellsOnLayer.weight[cellIdx];
0181           maxEnergyCellIndex = cellIdx;
0182           maxEnergyDetId = cellsOnLayer.detid[cellIdx];
0183         }
0184         thisCluster.emplace_back(cellsOnLayer.detid[cellIdx], 1.f);
0185         if (cellsOnLayer.isSeed[cellIdx]) {
0186           seedDetId = cellsOnLayer.detid[cellIdx];
0187         }
0188       }
0189 
0190       float total_weight_log = 0.f;
0191       float total_weight = energy;
0192 
0193       if constexpr (std::is_same_v<STRATEGY, HGCalSiliconStrategy>) {
0194         auto thick = rhtools_.getSiThickIndex(maxEnergyDetId);
0195         for (auto cellIdx : cl) {
0196           const float d1 = cellsOnLayer.dim1[cellIdx] - cellsOnLayer.dim1[maxEnergyCellIndex];
0197           const float d2 = cellsOnLayer.dim2[cellIdx] - cellsOnLayer.dim2[maxEnergyCellIndex];
0198           if ((d1 * d1 + d2 * d2) < positionDeltaRho2_) {
0199             float Wi = std::max(thresholdW0_[thick] + std::log(cellsOnLayer.weight[cellIdx] / energy), 0.);
0200             x += cellsOnLayer.dim1[cellIdx] * Wi;
0201             y += cellsOnLayer.dim2[cellIdx] * Wi;
0202             total_weight_log += Wi;
0203           }
0204         }
0205       } else {
0206         for (auto cellIdx : cl) {
0207           auto position = rhtools_.getPosition(cellsOnLayer.detid[cellIdx]);
0208           x += position.x() * cellsOnLayer.weight[cellIdx];
0209           y += position.y() * cellsOnLayer.weight[cellIdx];
0210         }
0211       }
0212 
0213       if constexpr (std::is_same_v<STRATEGY, HGCalSiliconStrategy>) {
0214         total_weight = total_weight_log;
0215       }
0216 
0217       if (total_weight != 0.) {
0218         float inv_tot_weight = 1.f / total_weight;
0219         x *= inv_tot_weight;
0220         y *= inv_tot_weight;
0221       } else {
0222         x = cellsOnLayer.dim1[maxEnergyCellIndex];
0223         y = cellsOnLayer.dim2[maxEnergyCellIndex];
0224       }
0225       math::XYZPoint position = math::XYZPoint(x, y, z);
0226 
0227       auto globalClusterIndex = cellsOnLayer.clusterIndex[cl[0]] + firstClusterIdx;
0228 
0229       clusters_v_[globalClusterIndex] =
0230           reco::BasicCluster(energy, position, reco::CaloID::DET_HGCAL_ENDCAP, thisCluster, algoId_);
0231       clusters_v_[globalClusterIndex].setSeed(seedDetId);
0232       thisCluster.clear();
0233     }
0234 
0235     cellsIdInCluster.clear();
0236   }
0237   return clusters_v_;
0238 }
0239 template <typename T, typename STRATEGY>
0240 void HGCalCLUEAlgoT<T, STRATEGY>::calculateLocalDensity(const T& lt,
0241                                                         const unsigned int layerId,
0242                                                         float delta,
0243                                                         HGCalSiliconStrategy strategy) {
0244   auto& cellsOnLayer = cells_[layerId];
0245   unsigned int numberOfCells = cellsOnLayer.detid.size();
0246   for (unsigned int i = 0; i < numberOfCells; i++) {
0247     std::array<int, 4> search_box = lt.searchBox(cellsOnLayer.dim1[i] - delta,
0248                                                  cellsOnLayer.dim1[i] + delta,
0249                                                  cellsOnLayer.dim2[i] - delta,
0250                                                  cellsOnLayer.dim2[i] + delta);
0251 
0252     for (int xBin = search_box[0]; xBin < search_box[1] + 1; ++xBin) {
0253       for (int yBin = search_box[2]; yBin < search_box[3] + 1; ++yBin) {
0254         int binId = lt.getGlobalBinByBin(xBin, yBin);
0255         size_t binSize = lt[binId].size();
0256 
0257         for (unsigned int j = 0; j < binSize; j++) {
0258           unsigned int otherId = lt[binId][j];
0259           if (distance(lt, i, otherId, layerId) < delta) {
0260             cellsOnLayer.rho[i] += (i == otherId ? 1.f : 0.5f) * cellsOnLayer.weight[otherId];
0261           }
0262         }
0263       }
0264     }
0265     LogDebug("HGCalCLUEAlgo") << "Debugging calculateLocalDensity: \n"
0266                               << "  cell: " << i << " eta: " << cellsOnLayer.dim1[i] << " phi: " << cellsOnLayer.dim2[i]
0267                               << " energy: " << cellsOnLayer.weight[i] << " density: " << cellsOnLayer.rho[i] << "\n";
0268   }
0269 }
0270 template <typename T, typename STRATEGY>
0271 void HGCalCLUEAlgoT<T, STRATEGY>::calculateLocalDensity(const T& lt,
0272                                                         const unsigned int layerId,
0273                                                         float delta,
0274                                                         HGCalScintillatorStrategy strategy) {
0275   auto& cellsOnLayer = cells_[layerId];
0276   unsigned int numberOfCells = cellsOnLayer.detid.size();
0277   for (unsigned int i = 0; i < numberOfCells; i++) {
0278     std::array<int, 4> search_box = lt.searchBox(cellsOnLayer.dim1[i] - delta,
0279                                                  cellsOnLayer.dim1[i] + delta,
0280                                                  cellsOnLayer.dim2[i] - delta,
0281                                                  cellsOnLayer.dim2[i] + delta);
0282     cellsOnLayer.rho[i] += cellsOnLayer.weight[i];
0283     float northeast(0), northwest(0), southeast(0), southwest(0), all(0);
0284     for (int etaBin = search_box[0]; etaBin < search_box[1] + 1; ++etaBin) {
0285       for (int phiBin = search_box[2]; phiBin < search_box[3] + 1; ++phiBin) {
0286         int phi = (phiBin % T::type::nRows);
0287         int binId = lt.getGlobalBinByBin(etaBin, phi);
0288         size_t binSize = lt[binId].size();
0289         for (unsigned int j = 0; j < binSize; j++) {
0290           unsigned int otherId = lt[binId][j];
0291           if (distance(lt, i, otherId, layerId) < delta) {
0292             int iPhi = HGCScintillatorDetId(cellsOnLayer.detid[i]).iphi();
0293             int otherIPhi = HGCScintillatorDetId(cellsOnLayer.detid[otherId]).iphi();
0294             int iEta = HGCScintillatorDetId(cellsOnLayer.detid[i]).ieta();
0295             int otherIEta = HGCScintillatorDetId(cellsOnLayer.detid[otherId]).ieta();
0296             int dIPhi = otherIPhi - iPhi;
0297             dIPhi += abs(dIPhi) < 2 ? 0
0298                      : dIPhi < 0    ? scintMaxIphi_
0299                                     : -scintMaxIphi_;  // cells with iPhi=288 and iPhi=1 should be neiboring cells
0300             int dIEta = otherIEta - iEta;
0301             LogDebug("HGCalCLUEAlgo") << "  Debugging calculateLocalDensity for Scintillator: \n"
0302                                       << "    cell: " << otherId << " energy: " << cellsOnLayer.weight[otherId]
0303                                       << " otherIPhi: " << otherIPhi << " iPhi: " << iPhi << " otherIEta: " << otherIEta
0304                                       << " iEta: " << iEta << "\n";
0305 
0306             if (otherId != i) {
0307               auto neighborCellContribution = 0.5f * cellsOnLayer.weight[otherId];
0308               all += neighborCellContribution;
0309               if (dIPhi >= 0 && dIEta >= 0)
0310                 northeast += neighborCellContribution;
0311               if (dIPhi <= 0 && dIEta >= 0)
0312                 southeast += neighborCellContribution;
0313               if (dIPhi >= 0 && dIEta <= 0)
0314                 northwest += neighborCellContribution;
0315               if (dIPhi <= 0 && dIEta <= 0)
0316                 southwest += neighborCellContribution;
0317             }
0318             LogDebug("HGCalCLUEAlgo") << "  Debugging calculateLocalDensity for Scintillator: \n"
0319                                       << "    northeast: " << northeast << " southeast: " << southeast
0320                                       << " northwest: " << northwest << " southwest: " << southwest << "\n";
0321           }
0322         }
0323       }
0324     }
0325     float neighborsval = (std::max(northeast, northwest) > std::max(southeast, southwest))
0326                              ? std::max(northeast, northwest)
0327                              : std::max(southeast, southwest);
0328     if (use2x2_)
0329       cellsOnLayer.rho[i] += neighborsval;
0330     else
0331       cellsOnLayer.rho[i] += all;
0332     LogDebug("HGCalCLUEAlgo") << "Debugging calculateLocalDensity: \n"
0333                               << "  cell: " << i << " eta: " << cellsOnLayer.dim1[i] << " phi: " << cellsOnLayer.dim2[i]
0334                               << " energy: " << cellsOnLayer.weight[i] << " density: " << cellsOnLayer.rho[i] << "\n";
0335   }
0336 }
0337 template <typename T, typename STRATEGY>
0338 void HGCalCLUEAlgoT<T, STRATEGY>::calculateLocalDensity(const T& lt, const unsigned int layerId, float delta) {
0339   if constexpr (std::is_same_v<STRATEGY, HGCalSiliconStrategy>) {
0340     calculateLocalDensity(lt, layerId, delta, HGCalSiliconStrategy());
0341   } else {
0342     calculateLocalDensity(lt, layerId, delta, HGCalScintillatorStrategy());
0343   }
0344 }
0345 
0346 template <typename T, typename STRATEGY>
0347 void HGCalCLUEAlgoT<T, STRATEGY>::calculateDistanceToHigher(const T& lt, const unsigned int layerId, float delta) {
0348   auto& cellsOnLayer = cells_[layerId];
0349   unsigned int numberOfCells = cellsOnLayer.detid.size();
0350 
0351   for (unsigned int i = 0; i < numberOfCells; i++) {
0352     // initialize delta and nearest higher for i
0353     float maxDelta = std::numeric_limits<float>::max();
0354     float i_delta = maxDelta;
0355     int i_nearestHigher = -1;
0356     float rho_max = 0.f;
0357     auto range = outlierDeltaFactor_ * delta;
0358     std::array<int, 4> search_box = lt.searchBox(cellsOnLayer.dim1[i] - range,
0359                                                  cellsOnLayer.dim1[i] + range,
0360                                                  cellsOnLayer.dim2[i] - range,
0361                                                  cellsOnLayer.dim2[i] + range);
0362     // loop over all bins in the search box
0363     for (int dim1Bin = search_box[0]; dim1Bin < search_box[1] + 1; ++dim1Bin) {
0364       for (int dim2Bin = search_box[2]; dim2Bin < search_box[3] + 1; ++dim2Bin) {
0365         // get the id of this bin
0366         size_t binId = lt.getGlobalBinByBin(dim1Bin, dim2Bin);
0367         if constexpr (std::is_same_v<STRATEGY, HGCalScintillatorStrategy>)
0368           binId = lt.getGlobalBinByBin(dim1Bin, (dim2Bin % T::type::nRows));
0369         // get the size of this bin
0370         size_t binSize = lt[binId].size();
0371 
0372         // loop over all hits in this bin
0373         for (unsigned int j = 0; j < binSize; j++) {
0374           unsigned int otherId = lt[binId][j];
0375           float dist = distance2(lt, i, otherId, layerId);
0376           bool foundHigher =
0377               (cellsOnLayer.rho[otherId] > cellsOnLayer.rho[i]) ||
0378               (cellsOnLayer.rho[otherId] == cellsOnLayer.rho[i] && cellsOnLayer.detid[otherId] > cellsOnLayer.detid[i]);
0379           if (!foundHigher) {
0380             continue;
0381           }
0382           if ((dist < i_delta) || ((dist == i_delta) && (cellsOnLayer.rho[otherId] > rho_max)) ||
0383               ((dist == i_delta) && (cellsOnLayer.rho[otherId] == rho_max) &&
0384                (cellsOnLayer.detid[otherId] > cellsOnLayer.detid[i]))) {
0385             rho_max = cellsOnLayer.rho[otherId];
0386             i_delta = dist;
0387             i_nearestHigher = otherId;
0388           }
0389         }
0390       }
0391     }
0392     bool foundNearestHigherInSearchBox = (i_delta != maxDelta);
0393     if (foundNearestHigherInSearchBox) {
0394       cellsOnLayer.delta[i] = std::sqrt(i_delta);
0395       cellsOnLayer.nearestHigher[i] = i_nearestHigher;
0396     } else {
0397       // otherwise delta is guaranteed to be larger outlierDeltaFactor_*delta_c
0398       // we can safely maximize delta to be maxDelta
0399       cellsOnLayer.delta[i] = maxDelta;
0400       cellsOnLayer.nearestHigher[i] = -1;
0401     }
0402 
0403     LogDebug("HGCalCLUEAlgo") << "Debugging calculateDistanceToHigher: \n"
0404                               << "  cell: " << i << " eta: " << cellsOnLayer.dim1[i] << " phi: " << cellsOnLayer.dim2[i]
0405                               << " energy: " << cellsOnLayer.weight[i] << " density: " << cellsOnLayer.rho[i]
0406                               << " nearest higher: " << cellsOnLayer.nearestHigher[i]
0407                               << " distance: " << cellsOnLayer.delta[i] << "\n";
0408   }
0409 }
0410 
0411 template <typename T, typename STRATEGY>
0412 int HGCalCLUEAlgoT<T, STRATEGY>::findAndAssignClusters(const unsigned int layerId, float delta) {
0413   // this is called once per layer and endcap...
0414   // so when filling the cluster temporary vector of Hexels we resize each time
0415   // by the number  of clusters found. This is always equal to the number of
0416   // cluster centers...
0417   unsigned int nClustersOnLayer = 0;
0418   auto& cellsOnLayer = cells_[layerId];
0419   unsigned int numberOfCells = cellsOnLayer.detid.size();
0420   std::vector<int> localStack;
0421   // find cluster seeds and outlier
0422   for (unsigned int i = 0; i < numberOfCells; i++) {
0423     float rho_c = kappa_ * cellsOnLayer.sigmaNoise[i];
0424     // initialize clusterIndex
0425     cellsOnLayer.clusterIndex[i] = -1;
0426     bool isSeed = (cellsOnLayer.delta[i] > delta) && (cellsOnLayer.rho[i] >= rho_c);
0427     bool isOutlier = (cellsOnLayer.delta[i] > outlierDeltaFactor_ * delta) && (cellsOnLayer.rho[i] < rho_c);
0428     if (isSeed) {
0429       cellsOnLayer.clusterIndex[i] = nClustersOnLayer;
0430       cellsOnLayer.isSeed[i] = true;
0431       nClustersOnLayer++;
0432       localStack.push_back(i);
0433 
0434     } else if (!isOutlier) {
0435       cellsOnLayer.followers[cellsOnLayer.nearestHigher[i]].push_back(i);
0436     }
0437   }
0438 
0439   // need to pass clusterIndex to their followers
0440   while (!localStack.empty()) {
0441     int endStack = localStack.back();
0442     auto& thisSeed = cellsOnLayer.followers[endStack];
0443     localStack.pop_back();
0444 
0445     // loop over followers
0446     for (int j : thisSeed) {
0447       // pass id to a follower
0448       cellsOnLayer.clusterIndex[j] = cellsOnLayer.clusterIndex[endStack];
0449       // push this follower to localStack
0450       localStack.push_back(j);
0451     }
0452   }
0453   return nClustersOnLayer;
0454 }
0455 
0456 template <typename T, typename STRATEGY>
0457 void HGCalCLUEAlgoT<T, STRATEGY>::computeThreshold() {
0458   // To support the TDR geometry and also the post-TDR one (v9 onwards), we
0459   // need to change the logic of the vectors containing signal to noise and
0460   // thresholds. The first 3 indices will keep on addressing the different
0461   // thicknesses of the Silicon detectors in CE_E , the next 3 indices will address
0462   // the thicknesses of the Silicon detectors in CE_H, while the last one, number 6 (the
0463   // seventh) will address the Scintillators. This change will support both
0464   // geometries at the same time.
0465 
0466   if (initialized_)
0467     return;  // only need to calculate thresholds once
0468 
0469   initialized_ = true;
0470 
0471   std::vector<double> dummy;
0472 
0473   dummy.resize(maxNumberOfThickIndices_ + !isNose_, 0);  // +1 to accomodate for the Scintillators
0474   thresholds_.resize(maxlayer_, dummy);
0475   v_sigmaNoise_.resize(maxlayer_, dummy);
0476 
0477   for (unsigned ilayer = 1; ilayer <= maxlayer_; ++ilayer) {
0478     for (unsigned ithick = 0; ithick < maxNumberOfThickIndices_; ++ithick) {
0479       float sigmaNoise = 0.001f * fcPerEle_ * nonAgedNoises_[ithick] * dEdXweights_[ilayer] /
0480                          (fcPerMip_[ithick] * thicknessCorrection_[ithick]);
0481       thresholds_[ilayer - 1][ithick] = sigmaNoise * ecut_;
0482       v_sigmaNoise_[ilayer - 1][ithick] = sigmaNoise;
0483       LogDebug("HGCalCLUEAlgo") << "ilayer: " << ilayer << " nonAgedNoises: " << nonAgedNoises_[ithick]
0484                                 << " fcPerEle: " << fcPerEle_ << " fcPerMip: " << fcPerMip_[ithick]
0485                                 << " noiseMip: " << fcPerEle_ * nonAgedNoises_[ithick] / fcPerMip_[ithick]
0486                                 << " sigmaNoise: " << sigmaNoise << "\n";
0487     }
0488 
0489     if (!isNose_) {
0490       float scintillators_sigmaNoise = 0.001f * noiseMip_ * dEdXweights_[ilayer] / sciThicknessCorrection_;
0491       thresholds_[ilayer - 1][maxNumberOfThickIndices_] = ecut_ * scintillators_sigmaNoise;
0492       v_sigmaNoise_[ilayer - 1][maxNumberOfThickIndices_] = scintillators_sigmaNoise;
0493       LogDebug("HGCalCLUEAlgo") << "ilayer: " << ilayer << " noiseMip: " << noiseMip_
0494                                 << " scintillators_sigmaNoise: " << scintillators_sigmaNoise << "\n";
0495     }
0496   }
0497 }
0498 
0499 // explicit template instantiation
0500 template class HGCalCLUEAlgoT<HGCalSiliconLayerTiles, HGCalSiliconStrategy>;
0501 template class HGCalCLUEAlgoT<HGCalScintillatorLayerTiles, HGCalScintillatorStrategy>;
0502 template class HGCalCLUEAlgoT<HFNoseLayerTiles, HGCalSiliconStrategy>;
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!