Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-07-09 05:00:31 

0001 // Original Author: Marco Rovere
0002 //
0003 
0004 #include "LCToCPAssociatorByEnergyScoreImpl.h"
0005 
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
0008 #include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
0009 
0010 #include "SimCalorimetry/HGCalAssociatorProducers/interface/AssociatorTools.h"
0011 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0012 
0013 template <typename HIT>
0014 LCToCPAssociatorByEnergyScoreImpl<HIT>::LCToCPAssociatorByEnergyScoreImpl(
0015     edm::EDProductGetter const& productGetter,
0016     bool hardScatterOnly,
0017     std::shared_ptr<hgcal::RecHitTools> recHitTools,
0018     const std::unordered_map<DetId, const unsigned int>* hitMap,
0019     const std::vector<const HIT*>& hits)
0020     : hardScatterOnly_(hardScatterOnly),
0021       recHitTools_(recHitTools),
0022       hitMap_(hitMap),
0023       productGetter_(&productGetter),
0024       hits_(hits) {
0025   if constexpr (std::is_same_v<HIT, HGCRecHit>)
0026     layers_ = recHitTools_->lastLayerBH();
0027   else
0028     layers_ = recHitTools_->lastLayerBarrel() + 1;
0029 }
0030 
0031 template <typename HIT>
0032 ticl::association LCToCPAssociatorByEnergyScoreImpl<HIT>::makeConnections(
0033     const edm::Handle<reco::CaloClusterCollection>& cCCH, const edm::Handle<CaloParticleCollection>& cPCH) const {
0034   // Get collections
0035   const auto& clusters = *cCCH.product();
0036   const auto& caloParticles = *cPCH.product();
0037   auto nLayerClusters = clusters.size();
0038 
0039   //Consider CaloParticles coming from the hard scatterer, excluding the PU contribution and save the indices.
0040   std::vector<size_t> cPIndices;
0041   removeCPFromPU(caloParticles, cPIndices, hardScatterOnly_);
0042   auto nCaloParticles = cPIndices.size();
0043 
0044   // Initialize cPOnLayer. It contains the caloParticleOnLayer structure for all CaloParticles in each layer and
0045   // among other the information to compute the CaloParticle-To-LayerCluster score. It is one of the two objects that
0046   // build the output of the makeConnections function.
0047   // cPOnLayer[cpId][layerId]
0048   ticl::caloParticleToLayerCluster cPOnLayer;
0049   cPOnLayer.resize(nCaloParticles);
0050   for (unsigned int i = 0; i < nCaloParticles; ++i) {
0051     unsigned int nLayers = layers_ * 2;
0052     if constexpr (std::is_same_v<HIT, reco::PFRecHit>)
0053       nLayers = layers_;
0054     cPOnLayer[i].resize(nLayers);
0055     for (unsigned int j = 0; j < nLayers; ++j) {
0056       cPOnLayer[i][j].caloParticleId = i;
0057       cPOnLayer[i][j].energy = 0.f;
0058       cPOnLayer[i][j].hits_and_fractions.clear();
0059       //cPOnLayer[i][j].layerClusterIdToEnergyAndScore.reserve(nLayerClusters); // Not necessary but may improve performance
0060     }
0061   }
0062 
0063   // Fill detIdToCaloParticleId_Map and update cPOnLayer
0064   // The detIdToCaloParticleId_Map is used to connect a hit Detid (key) with all the CaloParticles that
0065   // contributed to that hit by storing the CaloParticle id and the fraction of the hit. Observe here
0066   // that all the different contributions of the same CaloParticle to a single hit (coming from their
0067   // internal SimClusters) are merged into a single entry with the fractions properly summed.
0068   std::unordered_map<DetId, std::vector<ticl::detIdInfoInCluster>> detIdToCaloParticleId_Map;
0069   for (const auto& cpId : cPIndices) {
0070     const SimClusterRefVector& simClusterRefVector = caloParticles[cpId].simClusters();
0071     for (const auto& it_sc : simClusterRefVector) {
0072       const SimCluster& simCluster = (*(it_sc));
0073       std::vector<std::pair<uint32_t, float>> hits_and_fractions;
0074       if constexpr (std::is_same_v<HIT, HGCRecHit>)
0075         hits_and_fractions = simCluster.endcap_hits_and_fractions();
0076       else
0077         hits_and_fractions = simCluster.barrel_hits_and_fractions();
0078       for (const auto& it_haf : hits_and_fractions) {
0079         const auto hitid = (it_haf.first);
0080         unsigned int cpLayerId = recHitTools_->getLayerWithOffset(hitid);
0081         if constexpr (std::is_same_v<HIT, HGCRecHit>)
0082           cpLayerId += layers_ * ((recHitTools_->zside(hitid) + 1) >> 1) - 1;
0083 
0084         const auto itcheck = hitMap_->find(hitid);
0085 
0086         if (itcheck != hitMap_->end()) {
0087           auto hit_find_it = detIdToCaloParticleId_Map.find(hitid);
0088           if (hit_find_it == detIdToCaloParticleId_Map.end()) {
0089             detIdToCaloParticleId_Map[hitid] = std::vector<ticl::detIdInfoInCluster>();
0090             detIdToCaloParticleId_Map[hitid].emplace_back(cpId, it_haf.second);
0091           } else {
0092             auto findHitIt = std::find(detIdToCaloParticleId_Map[hitid].begin(),
0093                                        detIdToCaloParticleId_Map[hitid].end(),
0094                                        ticl::detIdInfoInCluster{cpId, it_haf.second});
0095             if (findHitIt != detIdToCaloParticleId_Map[hitid].end()) {
0096               findHitIt->fraction += it_haf.second;
0097             } else {
0098               detIdToCaloParticleId_Map[hitid].emplace_back(cpId, it_haf.second);
0099             }
0100           }
0101           const HIT* hit = hits_[itcheck->second];
0102           cPOnLayer[cpId][cpLayerId].energy += it_haf.second * hit->energy();
0103           // We need to compress the hits and fractions in order to have a
0104           // reasonable score between CP and LC. Imagine, for example, that a
0105           // CP has detID X used by 2 SimClusters with different fractions. If
0106           // a single LC uses X with fraction 1 and is compared to the 2
0107           // contributions separately, it will be assigned a score != 0, which
0108           // is wrong.
0109           auto& haf = cPOnLayer[cpId][cpLayerId].hits_and_fractions;
0110           auto found = std::find_if(
0111               std::begin(haf), std::end(haf), [&hitid](const std::pair<DetId, float>& v) { return v.first == hitid; });
0112           if (found != haf.end()) {
0113             found->second += it_haf.second;
0114           } else {
0115             cPOnLayer[cpId][cpLayerId].hits_and_fractions.emplace_back(hitid, it_haf.second);
0116           }
0117         }
0118       }
0119     }
0120   }
0121 
0122 #ifdef EDM_ML_DEBUG
0123   LogDebug("LCToCPAssociatorByEnergyScoreImpl") << "cPOnLayer INFO" << std::endl;
0124   for (size_t cp = 0; cp < cPOnLayer.size(); ++cp) {
0125     LogDebug("LCToCPAssociatorByEnergyScoreImpl") << "For CaloParticle Idx: " << cp << " we have: " << std::endl;
0126     for (size_t cpp = 0; cpp < cPOnLayer[cp].size(); ++cpp) {
0127       LogDebug("LCToCPAssociatorByEnergyScoreImpl") << "  On Layer: " << cpp << " we have:" << std::endl;
0128       LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0129           << "    CaloParticleIdx: " << cPOnLayer[cp][cpp].caloParticleId << std::endl;
0130       LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0131           << "    Energy:          " << cPOnLayer[cp][cpp].energy << std::endl;
0132       double tot_energy = 0.;
0133       for (auto const& haf : cPOnLayer[cp][cpp].hits_and_fractions) {
0134         //const HIT* hit = &(hits_[hitMap_->at(DetId(haf.first))]);
0135         const HIT* hit = hits_[hitMap_->at(haf.first)];
0136         LogDebug("LCToCPAssociatorByEnergyScoreImpl") << "      Hits/fraction/energy: " << (uint32_t)haf.first << "/"
0137                                                       << haf.second << "/" << haf.second * hit->energy() << std::endl;
0138         tot_energy += haf.second * hit->energy();
0139       }
0140       LogDebug("LCToCPAssociatorByEnergyScoreImpl") << "    Tot Sum haf: " << tot_energy << std::endl;
0141       for (auto const& lc : cPOnLayer[cp][cpp].layerClusterIdToEnergyAndScore) {
0142         LogDebug("LCToCPAssociatorByEnergyScoreImpl") << "      lcIdx/energy/score: " << lc.first << "/"
0143                                                       << lc.second.first << "/" << lc.second.second << std::endl;
0144       }
0145     }
0146   }
0147 
0148   LogDebug("LCToCPAssociatorByEnergyScoreImpl") << "detIdToCaloParticleId_Map INFO" << std::endl;
0149   for (auto const& cp : detIdToCaloParticleId_Map) {
0150     LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0151         << "For detId: " << (uint32_t)cp.first
0152         << " we have found the following connections with CaloParticles:" << std::endl;
0153     const HIT* hit = hits_[hitMap_->at(cp.first)];
0154     for (auto const& cpp : cp.second) {
0155       LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0156           << "  CaloParticle Id: " << cpp.clusterId << " with fraction: " << cpp.fraction
0157           << " and energy: " << cpp.fraction * hit->energy() << std::endl;
0158     }
0159   }
0160 #endif
0161 
0162   // Fill detIdToLayerClusterId_Map and cpsInLayerCluster; update cPOnLayer
0163   std::unordered_map<DetId, std::vector<ticl::detIdInfoInCluster>> detIdToLayerClusterId_Map;
0164   // this contains the ids of the caloparticles contributing with at least one
0165   // hit to the layer cluster and the reconstruction error. To be returned
0166   // since this contains the information to compute the
0167   // LayerCluster-To-CaloParticle score.
0168   ticl::layerClusterToCaloParticle cpsInLayerCluster;
0169   cpsInLayerCluster.resize(nLayerClusters);
0170 
0171   for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) {
0172     const std::vector<std::pair<DetId, float>>& hits_and_fractions = clusters[lcId].hitsAndFractions();
0173     unsigned int numberOfHitsInLC = hits_and_fractions.size();
0174     const auto firstHitDetId = hits_and_fractions[0].first;
0175     unsigned int lcLayerId = recHitTools_->getLayerWithOffset(firstHitDetId);
0176     if constexpr (std::is_same_v<HIT, HGCRecHit>)
0177       lcLayerId += layers_ * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
0178     for (unsigned int hitId = 0; hitId < numberOfHitsInLC; hitId++) {
0179       const auto rh_detid = hits_and_fractions[hitId].first;
0180       const auto rhFraction = hits_and_fractions[hitId].second;
0181 
0182       auto hit_find_in_LC = detIdToLayerClusterId_Map.find(rh_detid);
0183       if (hit_find_in_LC == detIdToLayerClusterId_Map.end()) {
0184         detIdToLayerClusterId_Map[rh_detid] = std::vector<ticl::detIdInfoInCluster>();
0185       }
0186       detIdToLayerClusterId_Map[rh_detid].emplace_back(lcId, rhFraction);
0187 
0188       auto hit_find_in_CP = detIdToCaloParticleId_Map.find(rh_detid);
0189 
0190       if (hit_find_in_CP != detIdToCaloParticleId_Map.end()) {
0191         const auto itcheck = hitMap_->find(rh_detid);
0192         const HIT* hit = hits_[itcheck->second];
0193         for (auto& h : hit_find_in_CP->second) {
0194           cPOnLayer[h.clusterId][lcLayerId].layerClusterIdToEnergyAndScore[lcId].first += h.fraction * hit->energy();
0195           cpsInLayerCluster[lcId].emplace_back(h.clusterId, 0.f);
0196         }
0197       }
0198     }  // End loop over hits on a LayerCluster
0199   }  // End of loop over LayerClusters
0200 
0201 #ifdef EDM_ML_DEBUG
0202   for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) {
0203     const auto& hits_and_fractions = clusters[lcId].hitsAndFractions();
0204     unsigned int numberOfHitsInLC = hits_and_fractions.size();
0205     const auto firstHitDetId = hits_and_fractions[0].first;
0206     int lcLayerId = recHitTools_->getLayerWithOffset(firstHitDetId);
0207     if constexpr (std::is_same_v<HIT, HGCRecHit>)
0208       lcLayerId += layers_ * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
0209     // This vector will store, for each hit in the Layercluster, the index of
0210     // the CaloParticle that contributed the most, in terms of energy, to it.
0211     // Special values are:
0212     //
0213     // -2  --> the reconstruction fraction of the RecHit is 0 (used in the past to monitor Halo Hits)
0214     // -3  --> same as before with the added condition that no CaloParticle has been linked to this RecHit
0215     // -1  --> the reco fraction is >0, but no CaloParticle has been linked to it
0216     // >=0 --> index of the linked CaloParticle
0217     std::vector<int> hitsToCaloParticleId(numberOfHitsInLC);
0218     // This will store the index of the CaloParticle linked to the LayerCluster that has the most number of hits in common.
0219     int maxCPId_byNumberOfHits = -1;
0220     // This will store the maximum number of shared hits between a Layercluster andd a CaloParticle
0221     unsigned int maxCPNumberOfHitsInLC = 0;
0222     // This will store the index of the CaloParticle linked to the LayerCluster that has the most energy in common.
0223     int maxCPId_byEnergy = -1;
0224     // This will store the maximum number of shared energy between a Layercluster and a CaloParticle
0225     float maxEnergySharedLCandCP = 0.f;
0226     // This will store the fraction of the LayerCluster energy shared with the best(energy) CaloParticle: e_shared/lc_energy
0227     float energyFractionOfLCinCP = 0.f;
0228     // This will store the fraction of the CaloParticle energy shared with the LayerCluster: e_shared/cp_energy
0229     float energyFractionOfCPinLC = 0.f;
0230     std::unordered_map<unsigned, unsigned> occurrencesCPinLC;
0231     unsigned int numberOfNoiseHitsInLC = 0;
0232     std::unordered_map<unsigned, float> CPEnergyInLC;
0233     for (unsigned int hitId = 0; hitId < numberOfHitsInLC; hitId++) {
0234       const auto rh_detid = hits_and_fractions[hitId].first;
0235       const auto rhFraction = hits_and_fractions[hitId].second;
0236 
0237       auto hit_find_in_CP = detIdToCaloParticleId_Map.find(rh_detid);
0238 
0239       // if the fraction is zero or the hit does not belong to any calo
0240       // particle, set the caloparticleId for the hit to -1 this will
0241       // contribute to the number of noise hits
0242 
0243       // MR Remove the case in which the fraction is 0, since this could be a
0244       // real hit that has been marked as halo.
0245       if (rhFraction == 0.) {
0246         hitsToCaloParticleId[hitId] = -2;
0247       }
0248       if (hit_find_in_CP == detIdToCaloParticleId_Map.end()) {
0249         hitsToCaloParticleId[hitId] -= 1;
0250       } else {
0251         const HIT* hit = hits_[hitMap_->at(rh_detid)];
0252         auto maxCPEnergyInLC = 0.f;
0253         auto maxCPId = -1;
0254         for (auto& h : hit_find_in_CP->second) {
0255           CPEnergyInLC[h.clusterId] += h.fraction * hit->energy();
0256           // Keep track of which CaloParticle ccontributed the most, in terms
0257           // of energy, to this specific LayerCluster.
0258           if (CPEnergyInLC[h.clusterId] > maxCPEnergyInLC) {
0259             maxCPEnergyInLC = CPEnergyInLC[h.clusterId];
0260             maxCPId = h.clusterId;
0261           }
0262         }
0263         hitsToCaloParticleId[hitId] = maxCPId;
0264       }
0265     }  // End loop over hits on a LayerCluster
0266     for (const auto& c : hitsToCaloParticleId) {
0267       if (c < 0) {
0268         numberOfNoiseHitsInLC++;
0269       } else {
0270         occurrencesCPinLC[c]++;
0271       }
0272     }
0273 
0274     for (const auto& c : occurrencesCPinLC) {
0275       if (c.second > maxCPNumberOfHitsInLC) {
0276         maxCPId_byNumberOfHits = c.first;
0277         maxCPNumberOfHitsInLC = c.second;
0278       }
0279     }
0280 
0281     for (const auto& c : CPEnergyInLC) {
0282       if (c.second > maxEnergySharedLCandCP) {
0283         maxCPId_byEnergy = c.first;
0284         maxEnergySharedLCandCP = c.second;
0285       }
0286     }
0287 
0288     float totalCPEnergyOnLayer = 0.f;
0289     if (maxCPId_byEnergy >= 0) {
0290       totalCPEnergyOnLayer = cPOnLayer[maxCPId_byEnergy][lcLayerId].energy;
0291       energyFractionOfCPinLC = maxEnergySharedLCandCP / totalCPEnergyOnLayer;
0292       if (clusters[lcId].energy() > 0.f) {
0293         energyFractionOfLCinCP = maxEnergySharedLCandCP / clusters[lcId].energy();
0294       }
0295     }
0296 
0297     LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0298         << std::setw(10) << "LayerId:\t" << std::setw(12) << "layerCluster\t" << std::setw(10) << "lc energy\t"
0299         << std::setw(5) << "nhits\t" << std::setw(12) << "noise hits\t" << std::setw(22) << "maxCPId_byNumberOfHits\t"
0300         << std::setw(8) << "nhitsCP\t" << std::setw(13) << "maxCPId_byEnergy\t" << std::setw(20)
0301         << "maxEnergySharedLCandCP\t" << std::setw(22) << "totalCPEnergyOnLayer\t" << std::setw(22)
0302         << "energyFractionOfLCinCP\t" << std::setw(25) << "energyFractionOfCPinLC\t"
0303         << "\n";
0304     LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0305         << std::setw(10) << lcLayerId << "\t" << std::setw(12) << lcId << "\t" << std::setw(10)
0306         << clusters[lcId].energy() << "\t" << std::setw(5) << numberOfHitsInLC << "\t" << std::setw(12)
0307         << numberOfNoiseHitsInLC << "\t" << std::setw(22) << maxCPId_byNumberOfHits << "\t" << std::setw(8)
0308         << maxCPNumberOfHitsInLC << "\t" << std::setw(13) << maxCPId_byEnergy << "\t" << std::setw(20)
0309         << maxEnergySharedLCandCP << "\t" << std::setw(22) << totalCPEnergyOnLayer << "\t" << std::setw(22)
0310         << energyFractionOfLCinCP << "\t" << std::setw(25) << energyFractionOfCPinLC << "\n";
0311   }  // End of loop over LayerClusters
0312 
0313   LogDebug("LCToCPAssociatorByEnergyScoreImpl") << "Improved cPOnLayer INFO" << std::endl;
0314   for (size_t cp = 0; cp < cPOnLayer.size(); ++cp) {
0315     LogDebug("LCToCPAssociatorByEnergyScoreImpl") << "For CaloParticle Idx: " << cp << " we have: " << std::endl;
0316     for (size_t cpp = 0; cpp < cPOnLayer[cp].size(); ++cpp) {
0317       LogDebug("LCToCPAssociatorByEnergyScoreImpl") << "  On Layer: " << cpp << " we have:" << std::endl;
0318       LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0319           << "    CaloParticleIdx: " << cPOnLayer[cp][cpp].caloParticleId << std::endl;
0320       LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0321           << "    Energy:          " << cPOnLayer[cp][cpp].energy << std::endl;
0322       double tot_energy = 0.;
0323       for (auto const& haf : cPOnLayer[cp][cpp].hits_and_fractions) {
0324         const HIT* hit = hits_[hitMap_->at(haf.first)];
0325         LogDebug("LCToCPAssociatorByEnergyScoreImpl") << "      Hits/fraction/energy: " << (uint32_t)haf.first << "/"
0326                                                       << haf.second << "/" << haf.second * hit->energy() << std::endl;
0327         tot_energy += haf.second * hit->energy();
0328       }
0329       LogDebug("LCToCPAssociatorByEnergyScoreImpl") << "    Tot Sum haf: " << tot_energy << std::endl;
0330       for (auto const& lc : cPOnLayer[cp][cpp].layerClusterIdToEnergyAndScore) {
0331         LogDebug("LCToCPAssociatorByEnergyScoreImpl") << "      lcIdx/energy/score: " << lc.first << "/"
0332                                                       << lc.second.first << "/" << lc.second.second << std::endl;
0333       }
0334     }
0335   }
0336 
0337   LogDebug("LCToCPAssociatorByEnergyScoreImpl") << "Improved detIdToCaloParticleId_Map INFO" << std::endl;
0338   for (auto const& cp : detIdToCaloParticleId_Map) {
0339     LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0340         << "For detId: " << (uint32_t)cp.first
0341         << " we have found the following connections with CaloParticles:" << std::endl;
0342     const HIT* hit = hits_[hitMap_->at(cp.first)];
0343     for (auto const& cpp : cp.second) {
0344       LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0345           << "  CaloParticle Id: " << cpp.clusterId << " with fraction: " << cpp.fraction
0346           << " and energy: " << cpp.fraction * hit->energy() << std::endl;
0347     }
0348   }
0349 #endif
0350 
0351   // Update cpsInLayerCluster; compute the score LayerCluster-to-CaloParticle,
0352   // together with the returned AssociationMap
0353   for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) {
0354     // find the unique caloparticles id contributing to the layer clusters
0355     if constexpr (std::is_same_v<HIT, HGCRecHit>) {
0356       if (recHitTools_->isBarrel(clusters[lcId].seed()))
0357         continue;
0358     } else {
0359       if (!recHitTools_->isBarrel(clusters[lcId].seed()))
0360         continue;
0361     }
0362     std::sort(cpsInLayerCluster[lcId].begin(), cpsInLayerCluster[lcId].end());
0363     auto last = std::unique(cpsInLayerCluster[lcId].begin(), cpsInLayerCluster[lcId].end());
0364     cpsInLayerCluster[lcId].erase(last, cpsInLayerCluster[lcId].end());
0365     const auto& hits_and_fractions = clusters[lcId].hitsAndFractions();
0366     unsigned int numberOfHitsInLC = hits_and_fractions.size();
0367     // If a reconstructed LayerCluster has energy 0 but is linked to a
0368     // CaloParticle, assigned score 1
0369     if (clusters[lcId].energy() == 0. && !cpsInLayerCluster[lcId].empty()) {
0370       for (auto& cpPair : cpsInLayerCluster[lcId]) {
0371         cpPair.second = 1.;
0372         LogDebug("LCToCPAssociatorByEnergyScoreImpl") << "layerClusterId : \t " << lcId << "\t CP id : \t"
0373                                                       << cpPair.first << "\t score \t " << cpPair.second << "\n";
0374       }
0375       continue;
0376     }
0377 
0378     // Compute the correct normalization
0379     // It is the inverse of the denominator of the LCToCP score formula. Observe that this is the sum of the squares.
0380     float invLayerClusterEnergyWeight = 0.f;
0381     for (auto const& haf : hits_and_fractions) {
0382       const HIT* hit = hits_[hitMap_->at(haf.first)];
0383       invLayerClusterEnergyWeight += (haf.second * hit->energy()) * (haf.second * hit->energy());
0384     }
0385     invLayerClusterEnergyWeight = 1.f / invLayerClusterEnergyWeight;
0386     for (unsigned int i = 0; i < numberOfHitsInLC; ++i) {
0387       DetId rh_detid = hits_and_fractions[i].first;
0388       float rhFraction = hits_and_fractions[i].second;
0389 
0390       bool hitWithNoCP = (detIdToCaloParticleId_Map.find(rh_detid) == detIdToCaloParticleId_Map.end());
0391 
0392       auto itcheck = hitMap_->find(rh_detid);
0393       if (itcheck == hitMap_->end())
0394         continue;
0395       const HIT* hit = hits_[itcheck->second];
0396       float hitEnergyWeight = hit->energy() * hit->energy();
0397 
0398       for (auto& cpPair : cpsInLayerCluster[lcId]) {
0399         float cpFraction = 0.f;
0400         if (!hitWithNoCP) {
0401           auto findHitIt = std::find(detIdToCaloParticleId_Map[rh_detid].begin(),
0402                                      detIdToCaloParticleId_Map[rh_detid].end(),
0403                                      ticl::detIdInfoInCluster{cpPair.first, 0.f});
0404           if (findHitIt != detIdToCaloParticleId_Map[rh_detid].end())
0405             cpFraction = findHitIt->fraction;
0406         }
0407         cpPair.second += std::min(std::pow(rhFraction - cpFraction, 2), std::pow(rhFraction, 2)) * hitEnergyWeight *
0408                          invLayerClusterEnergyWeight;
0409       }  //End of loop over CaloParticles related the this LayerCluster.
0410     }  // End of loop over Hits within a LayerCluster
0411 #ifdef EDM_ML_DEBUG
0412     if (cpsInLayerCluster[lcId].empty())
0413       LogDebug("LCToCPAssociatorByEnergyScoreImpl") << "layerCluster Id: \t" << lcId << "\tCP id:\t-1 "
0414                                                     << "\t score \t-1\n";
0415 #endif
0416   }  // End of loop over LayerClusters
0417 
0418   // Compute the CaloParticle-To-LayerCluster score
0419   for (const auto& cpId : cPIndices) {
0420     unsigned int nLayers = layers_ * 2;
0421     if constexpr (std::is_same_v<HIT, reco::PFRecHit>)
0422       nLayers = layers_;
0423     for (unsigned int layerId = 0; layerId < nLayers; ++layerId) {
0424       unsigned int CPNumberOfHits = cPOnLayer[cpId][layerId].hits_and_fractions.size();
0425       if (CPNumberOfHits == 0)
0426         continue;
0427 #ifdef EDM_ML_DEBUG
0428       int lcWithMaxEnergyInCP = -1;
0429       float maxEnergyLCinCP = 0.f;
0430       float CPenergy = cPOnLayer[cpId][layerId].energy;
0431       float CPEnergyFractionInLC = 0.f;
0432       for (auto& lc : cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore) {
0433         if (lc.second.first > maxEnergyLCinCP) {
0434           maxEnergyLCinCP = lc.second.first;
0435           lcWithMaxEnergyInCP = lc.first;
0436         }
0437       }
0438       if (CPenergy > 0.f)
0439         CPEnergyFractionInLC = maxEnergyLCinCP / CPenergy;
0440 
0441       LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0442           << std::setw(8) << "LayerId:\t" << std::setw(12) << "caloparticle\t" << std::setw(15) << "cp total energy\t"
0443           << std::setw(15) << "cpEnergyOnLayer\t" << std::setw(14) << "CPNhitsOnLayer\t" << std::setw(18)
0444           << "lcWithMaxEnergyInCP\t" << std::setw(15) << "maxEnergyLCinCP\t" << std::setw(20) << "CPEnergyFractionInLC"
0445           << "\n";
0446       LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0447           << std::setw(8) << layerId << "\t" << std::setw(12) << cpId << "\t" << std::setw(15)
0448           << caloParticles[cpId].energy() << "\t" << std::setw(15) << CPenergy << "\t" << std::setw(14)
0449           << CPNumberOfHits << "\t" << std::setw(18) << lcWithMaxEnergyInCP << "\t" << std::setw(15) << maxEnergyLCinCP
0450           << "\t" << std::setw(20) << CPEnergyFractionInLC << "\n";
0451 #endif
0452       // Compute the correct normalization. Observe that this is the sum of the squares.
0453       float invCPEnergyWeight = 0.f;
0454       for (auto const& haf : cPOnLayer[cpId][layerId].hits_and_fractions) {
0455         const HIT* hit = hits_[hitMap_->at(haf.first)];
0456         invCPEnergyWeight += std::pow(haf.second * hit->energy(), 2);
0457       }
0458       invCPEnergyWeight = 1.f / invCPEnergyWeight;
0459       for (unsigned int i = 0; i < CPNumberOfHits; ++i) {
0460         auto& cp_hitDetId = cPOnLayer[cpId][layerId].hits_and_fractions[i].first;
0461         auto& cpFraction = cPOnLayer[cpId][layerId].hits_and_fractions[i].second;
0462 
0463         bool hitWithNoLC = false;
0464         if (cpFraction == 0.f)
0465           continue;  //hopefully this should never happen
0466         auto hit_find_in_LC = detIdToLayerClusterId_Map.find(cp_hitDetId);
0467         if (hit_find_in_LC == detIdToLayerClusterId_Map.end())
0468           hitWithNoLC = true;
0469         auto itcheck = hitMap_->find(cp_hitDetId);
0470         const HIT* hit = hits_[itcheck->second];
0471         float hitEnergyWeight = hit->energy() * hit->energy();
0472         for (auto& lcPair : cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore) {
0473           unsigned int layerClusterId = lcPair.first;
0474           float lcFraction = 0.f;
0475 
0476           if (!hitWithNoLC) {
0477             auto findHitIt = std::find(detIdToLayerClusterId_Map[cp_hitDetId].begin(),
0478                                        detIdToLayerClusterId_Map[cp_hitDetId].end(),
0479                                        ticl::detIdInfoInCluster{layerClusterId, 0.f});
0480             if (findHitIt != detIdToLayerClusterId_Map[cp_hitDetId].end())
0481               lcFraction = findHitIt->fraction;
0482           }
0483           lcPair.second.second += std::min(std::pow(lcFraction - cpFraction, 2), std::pow(cpFraction, 2)) *
0484                                   hitEnergyWeight * invCPEnergyWeight;
0485 #ifdef EDM_ML_DEBUG
0486           LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0487               << "cpDetId:\t" << (uint32_t)cp_hitDetId << "\tlayerClusterId:\t" << layerClusterId << "\t"
0488               << "lcfraction,cpfraction:\t" << lcFraction << ", " << cpFraction << "\t"
0489               << "hitEnergyWeight:\t" << hitEnergyWeight << "\t"
0490               << "current score:\t" << lcPair.second.second << "\t"
0491               << "invCPEnergyWeight:\t" << invCPEnergyWeight << "\n";
0492 #endif
0493         }  // End of loop over LayerClusters linked to hits of this CaloParticle
0494       }  // End of loop over hits of CaloParticle on a Layer
0495 #ifdef EDM_ML_DEBUG
0496       if (cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore.empty())
0497         LogDebug("LCToCPAssociatorByEnergyScoreImpl") << "CP Id: \t" << cpId << "\tLC id:\t-1 "
0498                                                       << "\t score \t-1\n";
0499 
0500       for (const auto& lcPair : cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore) {
0501         LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0502             << "CP Id: \t" << cpId << "\t LC id: \t" << lcPair.first << "\t score \t" << lcPair.second.second
0503             << "\t shared energy:\t" << lcPair.second.first << "\t shared energy fraction:\t"
0504             << (lcPair.second.first / CPenergy) << "\n";
0505       }
0506 #endif
0507     }  // End of loop over layers
0508   }  // End of loop over CaloParticles
0509 
0510   return {cpsInLayerCluster, cPOnLayer};
0511 }
0512 
0513 template <typename HIT>
0514 ticl::RecoToSimCollection LCToCPAssociatorByEnergyScoreImpl<HIT>::associateRecoToSim(
0515     const edm::Handle<reco::CaloClusterCollection>& cCCH, const edm::Handle<CaloParticleCollection>& cPCH) const {
0516   ticl::RecoToSimCollection returnValue(productGetter_);
0517   const auto& links = makeConnections(cCCH, cPCH);
0518 
0519   const auto& cpsInLayerCluster = std::get<0>(links);
0520   for (size_t lcId = 0; lcId < cpsInLayerCluster.size(); ++lcId) {
0521     for (auto& cpPair : cpsInLayerCluster[lcId]) {
0522       LogDebug("LCToCPAssociatorByEnergyScoreImpl")
0523           << "layerCluster Id: \t" << lcId << "\t CP id: \t" << cpPair.first << "\t score \t" << cpPair.second << "\n";
0524       // Fill AssociationMap
0525       returnValue.insert(edm::Ref<reco::CaloClusterCollection>(cCCH, lcId),  // Ref to LC
0526                          std::make_pair(edm::Ref<CaloParticleCollection>(cPCH, cpPair.first),
0527                                         cpPair.second)  // Pair <Ref to CP, score>
0528       );
0529     }
0530   }
0531   return returnValue;
0532 }
0533 
0534 template <typename HIT>
0535 ticl::SimToRecoCollection LCToCPAssociatorByEnergyScoreImpl<HIT>::associateSimToReco(
0536     const edm::Handle<reco::CaloClusterCollection>& cCCH, const edm::Handle<CaloParticleCollection>& cPCH) const {
0537   ticl::SimToRecoCollection returnValue(productGetter_);
0538   const auto& links = makeConnections(cCCH, cPCH);
0539   const auto& cPOnLayer = std::get<1>(links);
0540   for (size_t cpId = 0; cpId < cPOnLayer.size(); ++cpId) {
0541     for (size_t layerId = 0; layerId < cPOnLayer[cpId].size(); ++layerId) {
0542       for (auto& lcPair : cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore) {
0543         returnValue.insert(
0544             edm::Ref<CaloParticleCollection>(cPCH, cpId),                              // Ref to CP
0545             std::make_pair(edm::Ref<reco::CaloClusterCollection>(cCCH, lcPair.first),  // Pair <Ref to LC,
0546                            std::make_pair(lcPair.second.first, lcPair.second.second))  // pair <energy, score> >
0547         );
0548       }
0549     }
0550   }
0551   return returnValue;
0552 }
0553 
0554 template class LCToCPAssociatorByEnergyScoreImpl<HGCRecHit>;
0555 template class LCToCPAssociatorByEnergyScoreImpl<reco::PFRecHit>;
[ 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!