Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:38:06 

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