Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-29 11:58:16

0001 #include <cfloat>
0002 
0003 #include "TSToSimTSHitLCAssociatorByEnergyScoreImpl.h"
0004 #include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 
0007 TSToSimTSHitLCAssociatorByEnergyScoreImpl::TSToSimTSHitLCAssociatorByEnergyScoreImpl(
0008     edm::EDProductGetter const& productGetter,
0009     bool hardScatterOnly,
0010     std::shared_ptr<hgcal::RecHitTools> recHitTools,
0011     const std::unordered_map<DetId, const HGCRecHit*>* hitMap)
0012     : hardScatterOnly_(hardScatterOnly), recHitTools_(recHitTools), hitMap_(hitMap), productGetter_(&productGetter) {
0013   layers_ = recHitTools_->lastLayerBH();
0014 }
0015 
0016 hgcal::association_t TSToSimTSHitLCAssociatorByEnergyScoreImpl::makeConnections(
0017     const edm::Handle<ticl::TracksterCollection>& tCH,
0018     const edm::Handle<reco::CaloClusterCollection>& lCCH,
0019     const edm::Handle<SimClusterCollection>& sCCH,
0020     const edm::Handle<CaloParticleCollection>& cPCH,
0021     const edm::Handle<ticl::TracksterCollection>& sTCH) const {
0022   // Get collections
0023   const auto& tracksters = *tCH.product();
0024   const auto& layerClusters = *lCCH.product();
0025   const auto& sC = *sCCH.product();
0026   const auto& cP = *cPCH.product();
0027   const auto cPHandle_id = cPCH.id();
0028   const auto& simTSs = *sTCH.product();
0029   const auto nTracksters = tracksters.size();
0030   const auto nSimTracksters = simTSs.size();
0031 
0032   std::unordered_map<DetId, std::vector<std::pair<int, float>>> detIdSimTSId_Map;
0033   std::unordered_map<DetId, std::vector<std::pair<int, float>>> detIdSimClusterId_Map;
0034   std::unordered_map<DetId, std::vector<std::pair<int, float>>> detIdCaloParticleId_Map;
0035   std::unordered_map<DetId, std::vector<std::pair<int, float>>> detIdToRecoTSId_Map;
0036 
0037   hgcal::sharedEnergyAndScore_t recoToSim_sharedEnergyAndScore;
0038   hgcal::sharedEnergyAndScore_t simToReco_sharedEnergyAndScore;
0039 
0040   recoToSim_sharedEnergyAndScore.resize(nTracksters);
0041   simToReco_sharedEnergyAndScore.resize(nSimTracksters);
0042 
0043   for (size_t i = 0; i < nTracksters; ++i)
0044     recoToSim_sharedEnergyAndScore[i].resize(nSimTracksters);
0045   for (size_t i = 0; i < nSimTracksters; ++i)
0046     simToReco_sharedEnergyAndScore[i].resize(nTracksters);
0047 
0048   // fill sim maps
0049 
0050   for (size_t i = 0; i < sC.size(); ++i) {
0051     for (const auto& haf : sC[i].hits_and_fractions()) {
0052       detIdSimClusterId_Map[haf.first].emplace_back(i, haf.second);
0053     }
0054   }
0055 
0056   for (size_t i = 0; i < cP.size(); ++i) {
0057     for (const auto& sc : cP[i].simClusters()) {
0058       for (const auto& haf : sc->hits_and_fractions()) {
0059         auto hitId = haf.first;
0060         auto found = std::find_if(detIdCaloParticleId_Map[hitId].begin(),
0061                                   detIdCaloParticleId_Map[hitId].end(),
0062                                   [=](const std::pair<int, float>& v) { return v.first == static_cast<int>(i); });
0063         if (found == detIdCaloParticleId_Map[hitId].end())
0064           detIdCaloParticleId_Map[haf.first].emplace_back(i, haf.second);
0065         else
0066           found->second += haf.second;
0067       }
0068     }
0069   }
0070 
0071   for (size_t i = 0; i < nSimTracksters; ++i) {
0072     const auto& lcsInSimTrackster = simTSs[i].vertices();
0073     const auto& multiplicities = simTSs[i].vertex_multiplicity();
0074     for (size_t j = 0; j < lcsInSimTrackster.size(); ++j) {
0075       assert(multiplicities[j] > 0.f);
0076       const auto& v = lcsInSimTrackster[j];
0077       float fraction = 1.f / multiplicities[j];
0078       for (const auto& haf : layerClusters[v].hitsAndFractions()) {
0079         detIdSimTSId_Map[haf.first].emplace_back(i, haf.second * fraction);
0080       }
0081     }
0082   }
0083 
0084   // fill reco map
0085 
0086   for (size_t i = 0; i < nTracksters; ++i) {
0087     const auto& lcsInSimTrackster = tracksters[i].vertices();
0088     const auto& multiplicities = tracksters[i].vertex_multiplicity();
0089     for (size_t j = 0; j < lcsInSimTrackster.size(); ++j) {
0090       assert(multiplicities[j] > 0.f);
0091       const auto& v = lcsInSimTrackster[j];
0092       float fraction = 1.f / multiplicities[j];
0093       for (const auto& haf : layerClusters[v].hitsAndFractions()) {
0094         detIdToRecoTSId_Map[haf.first].emplace_back(i, haf.second * fraction);
0095       }
0096     }
0097   }
0098 
0099   std::vector<float> denominator_simToReco(nSimTracksters, 0.f);
0100   std::vector<std::vector<float>> numerator_simToReco(nSimTracksters);
0101   std::vector<std::vector<float>> sharedEnergy(nSimTracksters);
0102 
0103   for (size_t i = 0; i < nSimTracksters; ++i) {
0104     numerator_simToReco[i].resize(nTracksters, 0.f);
0105     sharedEnergy[i].resize(nTracksters, 0.f);
0106 
0107     const auto seedIndex = simTSs[i].seedIndex();
0108     const auto& lcsInSimTrackster = simTSs[i].vertices();
0109 
0110     for (const auto& v : lcsInSimTrackster) {
0111       for (const auto& haf : layerClusters[v].hitsAndFractions()) {
0112         const auto hitId = haf.first;
0113         float simFraction = 0.f;
0114 
0115         std::vector<std::pair<int, float>>::iterator found;
0116         if (simTSs[i].seedID() == cPHandle_id) {
0117           found = std::find_if(detIdSimTSId_Map[hitId].begin(),
0118                                detIdSimTSId_Map[hitId].end(),
0119                                [=](const std::pair<int, float>& v) { return v.first == seedIndex; });
0120           if (found != detIdSimTSId_Map[hitId].end()) {
0121             const auto iLC = std::find(simTSs[i].vertices().begin(), simTSs[i].vertices().end(), v);
0122             const auto lcFraction =
0123                 1.f / simTSs[i].vertex_multiplicity(std::distance(std::begin(simTSs[i].vertices()), iLC));
0124             simFraction = found->second * lcFraction;
0125           }
0126         } else {
0127           found = std::find_if(detIdSimClusterId_Map[hitId].begin(),
0128                                detIdSimClusterId_Map[hitId].end(),
0129                                [=](const std::pair<int, float>& v) { return v.first == seedIndex; });
0130           if (found != detIdSimClusterId_Map[hitId].end()) {
0131             simFraction = found->second;
0132           }
0133         }
0134 
0135         float hitEnergy = hitMap_->find(hitId)->second->energy();
0136         float hitEnergySquared = hitEnergy * hitEnergy;
0137         float simFractionSquared = simFraction * simFraction;
0138         denominator_simToReco[i] += simFractionSquared * hitEnergySquared;
0139         for (size_t j = 0; j < nTracksters; ++j) {
0140           float recoFraction = 0.f;
0141 
0142           auto found_reco =
0143               std::find_if(detIdToRecoTSId_Map[hitId].begin(),
0144                            detIdToRecoTSId_Map[hitId].end(),
0145                            [=](const std::pair<int, float>& v) { return v.first == static_cast<int>(j); });
0146           if (found_reco != detIdToRecoTSId_Map[hitId].end())
0147             recoFraction = found_reco->second;
0148           numerator_simToReco[i][j] +=
0149               std::min(simFractionSquared, (simFraction - recoFraction) * (simFraction - recoFraction)) *
0150               hitEnergySquared;
0151           sharedEnergy[i][j] += std::min(simFraction, recoFraction) * hitEnergy;
0152         }
0153       }
0154     }
0155   }
0156 
0157   std::vector<float> denominator_recoToSim(nTracksters, 0.f);
0158   std::vector<std::vector<float>> numerator_recoToSim(nTracksters);
0159 
0160   for (unsigned int i = 0; i < nTracksters; ++i) {
0161     numerator_recoToSim[i].resize(nSimTracksters, 0.f);
0162     const auto& lcsInTrackster = tracksters[i].vertices();
0163     for (const auto& v : lcsInTrackster) {
0164       for (const auto& haf : layerClusters[v].hitsAndFractions()) {
0165         const auto hitId = haf.first;
0166         float recoFraction = 0.f;
0167 
0168         auto found = std::find_if(detIdToRecoTSId_Map[hitId].begin(),
0169                                   detIdToRecoTSId_Map[hitId].end(),
0170                                   [=](const std::pair<int, float>& v) { return v.first == static_cast<int>(i); });
0171         if (found != detIdToRecoTSId_Map[hitId].end())
0172           recoFraction = found->second;
0173 
0174         float hitEnergy = hitMap_->find(hitId)->second->energy();
0175         float hitEnergySquared = hitEnergy * hitEnergy;
0176         float recoFractionSquared = recoFraction * recoFraction;
0177         denominator_recoToSim[i] += recoFractionSquared * hitEnergySquared;
0178 
0179         for (size_t j = 0; j < nSimTracksters; ++j) {
0180           float simFraction = 0.f;
0181 
0182           auto found_sim = std::find_if(detIdSimTSId_Map[hitId].begin(),
0183                                         detIdSimTSId_Map[hitId].end(),
0184                                         [=](const std::pair<int, float>& v) { return v.first == static_cast<int>(j); });
0185           if (found_sim != detIdSimTSId_Map[hitId].end())
0186             simFraction = found_sim->second;
0187           numerator_recoToSim[i][j] +=
0188               std::min(recoFractionSquared, (simFraction - recoFraction) * (simFraction - recoFraction)) *
0189               hitEnergySquared;
0190         }
0191       }
0192     }
0193   }
0194 
0195   // compute score
0196 
0197   for (size_t i = 0; i < nSimTracksters; ++i) {
0198     for (size_t j = 0; j < nTracksters; ++j) {
0199       simToReco_sharedEnergyAndScore[i][j].first = sharedEnergy[i][j];
0200       simToReco_sharedEnergyAndScore[i][j].second = numerator_simToReco[i][j] / denominator_simToReco[i];
0201       recoToSim_sharedEnergyAndScore[j][i].first = sharedEnergy[i][j];
0202       recoToSim_sharedEnergyAndScore[j][i].second = numerator_recoToSim[j][i] / denominator_recoToSim[j];
0203     }
0204   }
0205 
0206   return {recoToSim_sharedEnergyAndScore, simToReco_sharedEnergyAndScore};
0207 }
0208 
0209 hgcal::RecoToSimCollectionSimTracksters TSToSimTSHitLCAssociatorByEnergyScoreImpl::associateRecoToSim(
0210     const edm::Handle<ticl::TracksterCollection>& tCH,
0211     const edm::Handle<reco::CaloClusterCollection>& lCCH,
0212     const edm::Handle<SimClusterCollection>& sCCH,
0213     const edm::Handle<CaloParticleCollection>& cPCH,
0214     const edm::Handle<ticl::TracksterCollection>& sTCH) const {
0215   hgcal::RecoToSimCollectionSimTracksters returnValue(productGetter_);
0216   const auto& links = makeConnections(tCH, lCCH, sCCH, cPCH, sTCH);
0217   const auto& recoToSim_sharedEnergyAndScore = std::get<0>(links);
0218   for (std::size_t tsId = 0; tsId < recoToSim_sharedEnergyAndScore.size(); ++tsId) {
0219     std::size_t numSimTracksters = recoToSim_sharedEnergyAndScore[tsId].size();
0220     for (std::size_t simTsId = 0; simTsId < numSimTracksters; ++simTsId) {
0221       LogDebug("TSToSimTSHitLCAssociatorByEnergyScoreImpl")
0222           << " Trackster Id:\t" << tsId << "\tSimTrackster id:\t" << recoToSim_sharedEnergyAndScore[tsId][simTsId].first
0223           << "\tscore:\t" << recoToSim_sharedEnergyAndScore[tsId][simTsId].second << "\n";
0224       // Fill AssociationMap
0225       returnValue.insert(
0226           edm::Ref<ticl::TracksterCollection>(tCH, tsId),  // Ref to TS
0227           std::make_pair(
0228               edm::Ref<ticl::TracksterCollection>(sTCH, simTsId),
0229               std::make_pair(recoToSim_sharedEnergyAndScore[tsId][simTsId].first,
0230                              recoToSim_sharedEnergyAndScore[tsId][simTsId].second))  // Pair <Ref to ST, score>
0231       );
0232     }
0233   }
0234   return returnValue;
0235 }
0236 
0237 hgcal::SimToRecoCollectionSimTracksters TSToSimTSHitLCAssociatorByEnergyScoreImpl::associateSimToReco(
0238     const edm::Handle<ticl::TracksterCollection>& tCH,
0239     const edm::Handle<reco::CaloClusterCollection>& lCCH,
0240     const edm::Handle<SimClusterCollection>& sCCH,
0241     const edm::Handle<CaloParticleCollection>& cPCH,
0242     const edm::Handle<ticl::TracksterCollection>& sTCH) const {
0243   hgcal::SimToRecoCollectionSimTracksters returnValue(productGetter_);
0244   const auto& links = makeConnections(tCH, lCCH, sCCH, cPCH, sTCH);
0245   const auto& simToReco_sharedEnergyAndScore = std::get<1>(links);
0246   for (std::size_t simTsId = 0; simTsId < simToReco_sharedEnergyAndScore.size(); ++simTsId) {
0247     std::size_t numTracksters = simToReco_sharedEnergyAndScore[simTsId].size();
0248     for (std::size_t tsId = 0; tsId < numTracksters; ++tsId) {
0249       LogDebug("TSToSimTSHitLCAssociatorByEnergyScoreImpl")
0250           << "Trackster Id:\t" << tsId << "\tSimTrackster id:\t" << simTsId << " Shared energy "
0251           << simToReco_sharedEnergyAndScore[simTsId][tsId].first << "\tscore:\t"
0252           << simToReco_sharedEnergyAndScore[simTsId][tsId].second << "\n";
0253       // Fill AssociationMap
0254       returnValue.insert(edm::Ref<ticl::TracksterCollection>(sTCH, simTsId),
0255                          std::make_pair(edm::Ref<ticl::TracksterCollection>(tCH, tsId),
0256                                         std::make_pair(simToReco_sharedEnergyAndScore[simTsId][tsId].first,
0257                                                        simToReco_sharedEnergyAndScore[simTsId][tsId].second)));
0258     }
0259   }
0260   return returnValue;
0261 }