HGCalAssociatorProducers/plugins/TSToSimTSHitLCAssociatorByEnergyScoreImpl.cc

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