Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 #include "LCToSCAssociatorByEnergyScoreImpl.h"
0002 
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
0005 
0006 template <typename HIT>
0007 LCToSCAssociatorByEnergyScoreImpl<HIT>::LCToSCAssociatorByEnergyScoreImpl(
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 HIT*>& hits)
0013     : hardScatterOnly_(hardScatterOnly),
0014       recHitTools_(recHitTools),
0015       hitMap_(hitMap),
0016       productGetter_(&productGetter),
0017       hits_(hits) {
0018   if constexpr (std::is_same_v<HIT, HGCRecHit>)
0019     layers_ = recHitTools_->lastLayerBH();
0020   else
0021     layers_ = recHitTools_->lastLayerBarrel() + 1;  //EB + 4 HB
0022 }
0023 
0024 template <typename HIT>
0025 ticl::association LCToSCAssociatorByEnergyScoreImpl<HIT>::makeConnections(
0026     const edm::Handle<reco::CaloClusterCollection>& cCCH, const edm::Handle<SimClusterCollection>& sCCH) const {
0027   // Get collections
0028   const auto& clusters = *cCCH.product();
0029   const auto& simClusters = *sCCH.product();
0030   auto nLayerClusters = clusters.size();
0031 
0032   //There shouldn't be any SimTracks from different crossings, but maybe they will be added later.
0033   //At the moment there should be one SimTrack in each SimCluster.
0034   auto nSimClusters = simClusters.size();
0035   std::vector<size_t> sCIndices;
0036   for (unsigned int scId = 0; scId < nSimClusters; ++scId) {
0037     if (hardScatterOnly_ && (simClusters[scId].g4Tracks()[0].eventId().event() != 0 or
0038                              simClusters[scId].g4Tracks()[0].eventId().bunchCrossing() != 0)) {
0039       LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0040           << "Excluding SimCluster from event: " << simClusters[scId].g4Tracks()[0].eventId().event()
0041           << " with BX: " << simClusters[scId].g4Tracks()[0].eventId().bunchCrossing() << std::endl;
0042       continue;
0043     }
0044     sCIndices.emplace_back(scId);
0045   }
0046   nSimClusters = sCIndices.size();
0047 
0048   // Initialize lcsInSimCluster. It contains the simClusterOnLayer structure for all simClusters in each layer and
0049   // among other the information to compute the SimCluster-To-LayerCluster score. It is one of the two objects that
0050   // build the output of the makeConnections function.
0051   // lcsInSimCluster[scId][layerId]
0052   ticl::simClusterToLayerCluster lcsInSimCluster;
0053   lcsInSimCluster.resize(nSimClusters);
0054   for (unsigned int i = 0; i < nSimClusters; ++i) {
0055     unsigned int nLayers = 2 * layers_;
0056     if constexpr (std::is_same_v<HIT, reco::PFRecHit>)
0057       nLayers = layers_;
0058     lcsInSimCluster[i].resize(nLayers);
0059     for (unsigned int j = 0; j < nLayers; ++j) {
0060       lcsInSimCluster[i][j].simClusterId = i;
0061       lcsInSimCluster[i][j].energy = 0.f;
0062       lcsInSimCluster[i][j].hits_and_fractions.clear();
0063     }
0064   }
0065 
0066   // Fill detIdToSimClusterId_Map and update lcsInSimCluster
0067   // The detIdToSimClusterId_Map is used to connect a hit Detid (key) with all the SimClusters that
0068   // contributed to that hit by storing the SimCluster id and the fraction of the hit. Observe here
0069   // that in contrast to the CaloParticle case there is no merging and summing of the fractions, which
0070   // in the CaloParticle's case was necessary due to the multiple SimClusters of a single CaloParticle.
0071   std::unordered_map<DetId, std::vector<ticl::detIdInfoInCluster>> detIdToSimClusterId_Map;
0072   for (const auto& scId : sCIndices) {
0073     std::vector<std::pair<uint32_t, float>> hits_and_fractions = simClusters[scId].hits_and_fractions();
0074     if constexpr (std::is_same_v<HIT, HGCRecHit>)
0075       hits_and_fractions = simClusters[scId].endcap_hits_and_fractions();
0076     else
0077       hits_and_fractions = simClusters[scId].barrel_hits_and_fractions();
0078     for (const auto& it_haf : hits_and_fractions) {
0079       const auto hitid = (it_haf.first);
0080       unsigned int scLayerId = recHitTools_->getLayer(hitid);
0081       if constexpr (std::is_same_v<HIT, HGCRecHit>)
0082         scLayerId += layers_ * ((recHitTools_->zside(hitid) + 1) >> 1) - 1;
0083       const auto itcheck = hitMap_->find(hitid);
0084       if (itcheck != hitMap_->end()) {
0085         auto hit_find_it = detIdToSimClusterId_Map.find(hitid);
0086         if (hit_find_it == detIdToSimClusterId_Map.end()) {
0087           detIdToSimClusterId_Map[hitid] = std::vector<ticl::detIdInfoInCluster>();
0088         }
0089         detIdToSimClusterId_Map[hitid].emplace_back(scId, it_haf.second);
0090         const HIT* hit = hits_[itcheck->second];
0091         lcsInSimCluster[scId][scLayerId].energy += it_haf.second * hit->energy();
0092         lcsInSimCluster[scId][scLayerId].hits_and_fractions.emplace_back(hitid, it_haf.second);
0093       }
0094     }
0095   }  // end of loop over SimClusters
0096 
0097 #ifdef EDM_ML_DEBUG
0098   LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0099       << "lcsInSimCluster INFO (Only SimCluster filled at the moment)" << std::endl;
0100   LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "    # of clusters :          " << nLayerClusters << std::endl;
0101   for (size_t sc = 0; sc < lcsInSimCluster.size(); ++sc) {
0102     LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "For SimCluster Idx: " << sc << " we have: " << std::endl;
0103     for (size_t sclay = 0; sclay < lcsInSimCluster[sc].size(); ++sclay) {
0104       LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "  On Layer: " << sclay << " we have:" << std::endl;
0105       LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0106           << "    SimClusterIdx: " << lcsInSimCluster[sc][sclay].simClusterId << std::endl;
0107       LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0108           << "    Energy:          " << lcsInSimCluster[sc][sclay].energy << std::endl;
0109       double tot_energy = 0.;
0110       for (auto const& haf : lcsInSimCluster[sc][sclay].hits_and_fractions) {
0111         const HIT* hit = hits_[hitMap_->at(haf.first)];
0112         LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "      Hits/fraction/energy: " << (uint32_t)haf.first << "/"
0113                                                       << haf.second << "/" << haf.second * hit->energy() << std::endl;
0114         tot_energy += haf.second * hit->energy();
0115       }
0116       LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "    Tot Sum haf: " << tot_energy << std::endl;
0117       for (auto const& lc : lcsInSimCluster[sc][sclay].layerClusterIdToEnergyAndScore) {
0118         LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "      lcIdx/energy/score: " << lc.first << "/"
0119                                                       << lc.second.first << "/" << lc.second.second << std::endl;
0120       }
0121     }
0122   }
0123 
0124   LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "detIdToSimClusterId_Map INFO" << std::endl;
0125   for (auto const& sc : detIdToSimClusterId_Map) {
0126     LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0127         << "For detId: " << (uint32_t)sc.first
0128         << " we have found the following connections with SimClusters:" << std::endl;
0129     // At this point here if you activate the printing you will notice cases where in a
0130     // specific detId there are more that one SimClusters contributing with fractions less than 1.
0131     // This is important since it effects the score computation, since the fraction is also in the
0132     // denominator of the score formula.
0133     const HIT* hit = hits_[hitMap_->at(sc.first)];
0134     for (auto const& sclu : sc.second) {
0135       LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0136           << "  SimCluster Id: " << sclu.clusterId << " with fraction: " << sclu.fraction
0137           << " and energy: " << sclu.fraction * hit->energy() << std::endl;
0138     }
0139   }
0140 #endif
0141 
0142   // Fill detIdToLayerClusterId_Map and scsInLayerCluster; update lcsInSimCluster
0143   // The detIdToLayerClusterId_Map is used to connect a hit Detid (key) with all the LayerClusters that
0144   // contributed to that hit by storing the LayerCluster id and the fraction of the corresponding hit.
0145   std::unordered_map<DetId, std::vector<ticl::detIdInfoInCluster>> detIdToLayerClusterId_Map;
0146   // scsInLayerCluster together with lcsInSimCluster are the two objects that are used to build the
0147   // output of the makeConnections function. scsInLayerCluster connects a LayerCluster with
0148   // all the SimClusters that share at least one cell with the LayerCluster and for each pair (LC,SC)
0149   // it stores the score.
0150   ticl::layerClusterToSimCluster scsInLayerCluster;  //[lcId][scId]->(score)
0151   scsInLayerCluster.resize(nLayerClusters);
0152 
0153   for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) {
0154     const std::vector<std::pair<DetId, float>>& hits_and_fractions = clusters[lcId].hitsAndFractions();
0155     unsigned int numberOfHitsInLC = hits_and_fractions.size();
0156     const auto firstHitDetId = hits_and_fractions[0].first;
0157     int lcLayerId = recHitTools_->getLayer(firstHitDetId);
0158     if constexpr (std::is_same_v<HIT, HGCRecHit>)
0159       lcLayerId += layers_ * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
0160     for (unsigned int hitId = 0; hitId < numberOfHitsInLC; hitId++) {
0161       const auto rh_detid = hits_and_fractions[hitId].first;
0162       const auto rhFraction = hits_and_fractions[hitId].second;
0163 
0164       auto hit_find_in_LC = detIdToLayerClusterId_Map.find(rh_detid);
0165       if (hit_find_in_LC == detIdToLayerClusterId_Map.end()) {
0166         detIdToLayerClusterId_Map[rh_detid] = std::vector<ticl::detIdInfoInCluster>();
0167       }
0168       detIdToLayerClusterId_Map[rh_detid].emplace_back(lcId, rhFraction);
0169 
0170       auto hit_find_in_SC = detIdToSimClusterId_Map.find(rh_detid);
0171 
0172       if (hit_find_in_SC != detIdToSimClusterId_Map.end()) {
0173         const auto itcheck = hitMap_->find(rh_detid);
0174         const HIT* hit = hits_[itcheck->second];
0175         //Loops through all the simclusters that have the layer cluster rechit under study
0176         //Here is time to update the lcsInSimCluster and connect the SimCluster with all
0177         //the layer clusters that have the current rechit detid matched.
0178         for (auto& h : hit_find_in_SC->second) {
0179           //lcsInSimCluster[simclusterId][layerId][layerclusterId]-> (energy,score)
0180           //SC_i - > LC_j, LC_k, ...
0181           lcsInSimCluster[h.clusterId][lcLayerId].layerClusterIdToEnergyAndScore[lcId].first +=
0182               h.fraction * hit->energy();
0183           //LC_i -> SC_j, SC_k, ...
0184           scsInLayerCluster[lcId].emplace_back(h.clusterId, 0.f);
0185         }
0186       }
0187     }  // End loop over hits on a LayerCluster
0188   }  // End of loop over LayerClusters
0189 
0190 #ifdef EDM_ML_DEBUG
0191   for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) {
0192     const auto& hits_and_fractions = clusters[lcId].hitsAndFractions();
0193     unsigned int numberOfHitsInLC = hits_and_fractions.size();
0194     const auto firstHitDetId = hits_and_fractions[0].first;
0195     int lcLayerId = recHitTools_->getLayer(firstHitDetId);
0196     if constexpr (std::is_same_v<HIT, HGCRecHit>)
0197       lcLayerId += layers_ * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1;
0198     // This vector will store, for each hit in the Layercluster, the index of
0199     // the SimCluster that contributed the most, in terms of energy, to it.
0200     // Special values are:
0201     //
0202     // -2  --> the reconstruction fraction of the RecHit is 0 (used in the past to monitor Halo Hits)
0203     // -3  --> same as before with the added condition that no SimCluster has been linked to this RecHit
0204     // -1  --> the reco fraction is >0, but no SimCluster has been linked to it
0205     // >=0 --> index of the linked SimCluster
0206     std::vector<int> hitsToSimClusterId(numberOfHitsInLC);
0207     // This will store the index of the SimCluster linked to the LayerCluster that has the most number of hits in common.
0208     int maxSCId_byNumberOfHits = -1;
0209     // This will store the maximum number of shared hits between a Layercluster and a SimCluster
0210     unsigned int maxSCNumberOfHitsInLC = 0;
0211     // This will store the index of the SimCluster linked to the LayerCluster that has the most energy in common.
0212     int maxSCId_byEnergy = -1;
0213     // This will store the maximum number of shared energy between a Layercluster and a SimCluster
0214     float maxEnergySharedLCandSC = 0.f;
0215     // This will store the fraction of the LayerCluster energy shared with the best(energy) SimCluster: e_shared/lc_energy
0216     float energyFractionOfLCinSC = 0.f;
0217     // This will store the fraction of the SimCluster energy shared with the LayerCluster: e_shared/sc_energy
0218     float energyFractionOfSCinLC = 0.f;
0219     std::unordered_map<unsigned, unsigned> occurrencesSCinLC;
0220     unsigned int numberOfNoiseHitsInLC = 0;
0221     std::unordered_map<unsigned, float> SCEnergyInLC;
0222 
0223     for (unsigned int hitId = 0; hitId < numberOfHitsInLC; hitId++) {
0224       const auto rh_detid = hits_and_fractions[hitId].first;
0225       const auto rhFraction = hits_and_fractions[hitId].second;
0226 
0227       auto hit_find_in_SC = detIdToSimClusterId_Map.find(rh_detid);
0228 
0229       // if the fraction is zero or the hit does not belong to any simcluster
0230       // set the caloparticleId for the hit to -1 this will
0231       // contribute to the number of noise hits
0232 
0233       // MR Remove the case in which the fraction is 0, since this could be a
0234       // real hit that has been marked as halo.
0235       if (rhFraction == 0.) {
0236         hitsToSimClusterId[hitId] = -2;
0237       }
0238       //Now check if there are SimClusters linked to this rechit of the layercluster
0239       if (hit_find_in_SC == detIdToSimClusterId_Map.end()) {
0240         hitsToSimClusterId[hitId] -= 1;
0241       } else {
0242         const auto itcheck = hitMap_->find(rh_detid);
0243         const HIT* hit = hits_[itcheck->second];
0244         auto maxSCEnergyInLC = 0.f;
0245         auto maxSCId = -1;
0246         //Loop through all the linked SimClusters
0247         for (auto& h : hit_find_in_SC->second) {
0248           SCEnergyInLC[h.clusterId] += h.fraction * hit->energy();
0249           // Keep track of which SimCluster ccontributed the most, in terms
0250           // of energy, to this specific LayerCluster.
0251           if (SCEnergyInLC[h.clusterId] > maxSCEnergyInLC) {
0252             maxSCEnergyInLC = SCEnergyInLC[h.clusterId];
0253             maxSCId = h.clusterId;
0254           }
0255         }
0256         hitsToSimClusterId[hitId] = maxSCId;
0257       }
0258     }  // End loop over hits on a LayerCluster
0259 
0260     for (const auto& c : hitsToSimClusterId) {
0261       if (c < 0) {
0262         numberOfNoiseHitsInLC++;
0263       } else {
0264         occurrencesSCinLC[c]++;
0265       }
0266     }
0267 
0268     for (const auto& c : occurrencesSCinLC) {
0269       if (c.second > maxSCNumberOfHitsInLC) {
0270         maxSCId_byNumberOfHits = c.first;
0271         maxSCNumberOfHitsInLC = c.second;
0272       }
0273     }
0274 
0275     for (const auto& c : SCEnergyInLC) {
0276       if (c.second > maxEnergySharedLCandSC) {
0277         maxSCId_byEnergy = c.first;
0278         maxEnergySharedLCandSC = c.second;
0279       }
0280     }
0281 
0282     float totalSCEnergyOnLayer = 0.f;
0283     if (maxSCId_byEnergy >= 0) {
0284       totalSCEnergyOnLayer = lcsInSimCluster[maxSCId_byEnergy][lcLayerId].energy;
0285       energyFractionOfSCinLC = maxEnergySharedLCandSC / totalSCEnergyOnLayer;
0286       if (clusters[lcId].energy() > 0.f) {
0287         energyFractionOfLCinSC = maxEnergySharedLCandSC / clusters[lcId].energy();
0288       }
0289     }
0290 
0291     LogDebug("LCToSCAssociatorByEnergyScoreImpl") << std::setw(10) << "LayerId:"
0292                                                   << "\t" << std::setw(12) << "layerCluster"
0293                                                   << "\t" << std::setw(10) << "lc energy"
0294                                                   << "\t" << std::setw(5) << "nhits"
0295                                                   << "\t" << std::setw(12) << "noise hits"
0296                                                   << "\t" << std::setw(22) << "maxSCId_byNumberOfHits"
0297                                                   << "\t" << std::setw(8) << "nhitsSC"
0298                                                   << "\t" << std::setw(13) << "maxSCId_byEnergy"
0299                                                   << "\t" << std::setw(20) << "maxEnergySharedLCandSC"
0300                                                   << "\t" << std::setw(22) << "totalSCEnergyOnLayer"
0301                                                   << "\t" << std::setw(22) << "energyFractionOfLCinSC"
0302                                                   << "\t" << std::setw(25) << "energyFractionOfSCinLC"
0303                                                   << "\t"
0304                                                   << "\n";
0305     LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0306         << std::setw(10) << lcLayerId << "\t" << std::setw(12) << lcId << "\t" << std::setw(10)
0307         << clusters[lcId].energy() << "\t" << std::setw(5) << numberOfHitsInLC << "\t" << std::setw(12)
0308         << numberOfNoiseHitsInLC << "\t" << std::setw(22) << maxSCId_byNumberOfHits << "\t" << std::setw(8)
0309         << maxSCNumberOfHitsInLC << "\t" << std::setw(13) << maxSCId_byEnergy << "\t" << std::setw(20)
0310         << maxEnergySharedLCandSC << "\t" << std::setw(22) << totalSCEnergyOnLayer << "\t" << std::setw(22)
0311         << energyFractionOfLCinSC << "\t" << std::setw(25) << energyFractionOfSCinLC << "\n";
0312   }  // End of loop over LayerClusters
0313 
0314   LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0315       << "Improved lcsInSimCluster INFO (Now containing the linked layer clusters id and energy - score still empty)"
0316       << std::endl;
0317   for (size_t sc = 0; sc < lcsInSimCluster.size(); ++sc) {
0318     LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "For SimCluster Idx: " << sc << " we have: " << std::endl;
0319     for (size_t sclay = 0; sclay < lcsInSimCluster[sc].size(); ++sclay) {
0320       LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "  On Layer: " << sclay << " we have:" << std::endl;
0321       LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0322           << "    SimClusterIdx: " << lcsInSimCluster[sc][sclay].simClusterId << std::endl;
0323       LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0324           << "    Energy:          " << lcsInSimCluster[sc][sclay].energy << std::endl;
0325       double tot_energy = 0.;
0326       for (auto const& haf : lcsInSimCluster[sc][sclay].hits_and_fractions) {
0327         const HIT* hit = hits_[hitMap_->at(haf.first)];
0328         LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "      Hits/fraction/energy: " << (uint32_t)haf.first << "/"
0329                                                       << haf.second << "/" << haf.second * hit->energy() << std::endl;
0330         tot_energy += haf.second * hit->energy();
0331       }
0332       LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "    Tot Sum haf: " << tot_energy << std::endl;
0333       for (auto const& lc : lcsInSimCluster[sc][sclay].layerClusterIdToEnergyAndScore) {
0334         LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "      lcIdx/energy/score: " << lc.first << "/"
0335                                                       << lc.second.first << "/" << lc.second.second << std::endl;
0336       }
0337     }
0338   }
0339 
0340   LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "Improved detIdToSimClusterId_Map INFO" << std::endl;
0341   for (auto const& sc : detIdToSimClusterId_Map) {
0342     const HIT* hit = hits_[hitMap_->at(sc.first)];
0343     LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0344         << "For detId: " << (uint32_t)sc.first
0345         << " we have found the following connections with SimClusters:" << std::endl;
0346     for (auto const& sclu : sc.second) {
0347       LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0348           << "  SimCluster Id: " << sclu.clusterId << " with fraction: " << sclu.fraction
0349           << " and energy: " << sclu.fraction * hit->energy() << std::endl;
0350     }
0351   }
0352 #endif
0353 
0354   // Update scsInLayerCluster; compute the score LayerCluster-to-SimCluster,
0355   // together with the returned AssociationMap
0356   for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) {
0357     // The simclusters contributing to the layer clusters should already be unique.
0358     // find the unique simclusters id contributing to the layer clusters
0359     if constexpr (std::is_same_v<HIT, HGCRecHit>) {
0360       if (recHitTools_->isBarrel(clusters[lcId].seed()))
0361         continue;
0362     } else {
0363       if (!recHitTools_->isBarrel(clusters[lcId].seed()))
0364         continue;
0365     }
0366     std::sort(scsInLayerCluster[lcId].begin(), scsInLayerCluster[lcId].end());
0367     auto last = std::unique(scsInLayerCluster[lcId].begin(), scsInLayerCluster[lcId].end());
0368     scsInLayerCluster[lcId].erase(last, scsInLayerCluster[lcId].end());
0369     const auto& hits_and_fractions = clusters[lcId].hitsAndFractions();
0370     unsigned int numberOfHitsInLC = hits_and_fractions.size();
0371     // If a reconstructed LayerCluster has energy 0 but is linked to a
0372     // SimCluster, assigned score 1
0373     if (clusters[lcId].energy() == 0. && !scsInLayerCluster[lcId].empty()) {
0374       for (auto& scPair : scsInLayerCluster[lcId]) {
0375         scPair.second = 1.;
0376         LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "layerClusterId : \t " << lcId << "\t SC id : \t"
0377                                                       << scPair.first << "\t score \t " << scPair.second << "\n";
0378       }
0379       continue;
0380     }
0381 
0382     // Compute the correct normalization.
0383     // It is the inverse of the denominator of the LCToSC score formula. Observe that this is the sum of the squares.
0384     float invLayerClusterEnergyWeight = 0.f;
0385     for (auto const& haf : hits_and_fractions) {
0386       const HIT* hit = hits_[hitMap_->at(haf.first)];
0387       invLayerClusterEnergyWeight += (haf.second * hit->energy()) * (haf.second * hit->energy());
0388     }
0389     invLayerClusterEnergyWeight = 1.f / invLayerClusterEnergyWeight;
0390     for (unsigned int i = 0; i < numberOfHitsInLC; ++i) {
0391       DetId rh_detid = hits_and_fractions[i].first;
0392       float rhFraction = hits_and_fractions[i].second;
0393 
0394       bool hitWithSC = (detIdToSimClusterId_Map.find(rh_detid) != detIdToSimClusterId_Map.end());
0395 
0396       auto itcheck = hitMap_->find(rh_detid);
0397       const HIT* hit = hits_[itcheck->second];
0398       float hitEnergyWeight = hit->energy() * hit->energy();
0399 
0400       for (auto& scPair : scsInLayerCluster[lcId]) {
0401         float scFraction = 0.f;
0402         if (hitWithSC) {
0403           auto findHitIt = std::find(detIdToSimClusterId_Map[rh_detid].begin(),
0404                                      detIdToSimClusterId_Map[rh_detid].end(),
0405                                      ticl::detIdInfoInCluster{scPair.first, 0.f});
0406           if (findHitIt != detIdToSimClusterId_Map[rh_detid].end())
0407             scFraction = findHitIt->fraction;
0408         }
0409         scPair.second += std::min(std::pow(rhFraction - scFraction, 2), std::pow(rhFraction, 2)) * hitEnergyWeight *
0410                          invLayerClusterEnergyWeight;
0411 #ifdef EDM_ML_DEBUG
0412         LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0413             << "rh_detid:\t" << (uint32_t)rh_detid << "\tlayerClusterId:\t" << lcId << "\t"
0414             << "rhfraction,scfraction:\t" << rhFraction << ", " << scFraction << "\t"
0415             << "hitEnergyWeight:\t" << hitEnergyWeight << "\t"
0416             << "current score:\t" << scPair.second << "\t"
0417             << "invLayerClusterEnergyWeight:\t" << invLayerClusterEnergyWeight << "\n";
0418 #endif
0419       }
0420     }  // End of loop over Hits within a LayerCluster
0421 #ifdef EDM_ML_DEBUG
0422     if (scsInLayerCluster[lcId].empty())
0423       LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "layerCluster Id: \t" << lcId << "\tSC id:\t-1 "
0424                                                     << "\t score \t-1"
0425                                                     << "\n";
0426 #endif
0427   }  // End of loop over LayerClusters
0428 
0429   // Compute the SimCluster-To-LayerCluster score
0430   for (const auto& scId : sCIndices) {
0431     unsigned int nLayers = layers_ * 2;
0432     if constexpr (std::is_same_v<HIT, reco::PFRecHit>)
0433       nLayers = layers_;
0434     for (unsigned int layerId = 0; layerId < nLayers; ++layerId) {
0435       unsigned int SCNumberOfHits = lcsInSimCluster[scId][layerId].hits_and_fractions.size();
0436       if (SCNumberOfHits == 0)
0437         continue;
0438 #ifdef EDM_ML_DEBUG
0439       int lcWithMaxEnergyInSC = -1;
0440       //energy of the most energetic LC from all that were linked to SC
0441       float maxEnergyLCinSC = 0.f;
0442       //Energy of the SC scId on layer layerId that was reconstructed.
0443       float SCenergy = lcsInSimCluster[scId][layerId].energy;
0444       //most energetic LC from all LCs linked to SC over SC energy.
0445       float SCEnergyFractionInLC = 0.f;
0446       for (auto& lc : lcsInSimCluster[scId][layerId].layerClusterIdToEnergyAndScore) {
0447         if (lc.second.first > maxEnergyLCinSC) {
0448           maxEnergyLCinSC = lc.second.first;
0449           lcWithMaxEnergyInSC = lc.first;
0450         }
0451       }
0452       if (SCenergy > 0.f)
0453         SCEnergyFractionInLC = maxEnergyLCinSC / SCenergy;
0454 
0455       LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0456           << std::setw(8) << "LayerId:\t" << std::setw(12) << "simcluster\t" << std::setw(15) << "sc total energy\t"
0457           << std::setw(15) << "scEnergyOnLayer\t" << std::setw(14) << "SCNhitsOnLayer\t" << std::setw(18)
0458           << "lcWithMaxEnergyInSC\t" << std::setw(15) << "maxEnergyLCinSC\t" << std::setw(20) << "SCEnergyFractionInLC"
0459           << "\n";
0460       LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0461           << std::setw(8) << layerId << "\t" << std::setw(12) << scId << "\t" << std::setw(15)
0462           << simClusters[scId].energy() << "\t" << std::setw(15) << SCenergy << "\t" << std::setw(14) << SCNumberOfHits
0463           << "\t" << std::setw(18) << lcWithMaxEnergyInSC << "\t" << std::setw(15) << maxEnergyLCinSC << "\t"
0464           << std::setw(20) << SCEnergyFractionInLC << "\n";
0465 #endif
0466       // Compute the correct normalization. Observe that this is the sum of the squares.
0467       float invSCEnergyWeight = 0.f;
0468       for (auto const& haf : lcsInSimCluster[scId][layerId].hits_and_fractions) {
0469         const HIT* hit = hits_[hitMap_->at(haf.first)];
0470         invSCEnergyWeight += std::pow(haf.second * hit->energy(), 2);
0471       }
0472       invSCEnergyWeight = 1.f / invSCEnergyWeight;
0473       for (unsigned int i = 0; i < SCNumberOfHits; ++i) {
0474         auto& sc_hitDetId = lcsInSimCluster[scId][layerId].hits_and_fractions[i].first;
0475         auto& scFraction = lcsInSimCluster[scId][layerId].hits_and_fractions[i].second;
0476 
0477         bool hitWithLC = false;
0478         if (scFraction == 0.f)
0479           continue;  //hopefully this should never happen
0480         auto hit_find_in_LC = detIdToLayerClusterId_Map.find(sc_hitDetId);
0481         if (hit_find_in_LC != detIdToLayerClusterId_Map.end())
0482           hitWithLC = true;
0483         auto itcheck = hitMap_->find(sc_hitDetId);
0484         const HIT* hit = hits_[itcheck->second];
0485         float hitEnergyWeight = hit->energy() * hit->energy();
0486         for (auto& lcPair : lcsInSimCluster[scId][layerId].layerClusterIdToEnergyAndScore) {
0487           unsigned int layerClusterId = lcPair.first;
0488           float lcFraction = 0.f;
0489 
0490           if (hitWithLC) {
0491             auto findHitIt = std::find(detIdToLayerClusterId_Map[sc_hitDetId].begin(),
0492                                        detIdToLayerClusterId_Map[sc_hitDetId].end(),
0493                                        ticl::detIdInfoInCluster{layerClusterId, 0.f});
0494             if (findHitIt != detIdToLayerClusterId_Map[sc_hitDetId].end())
0495               lcFraction = findHitIt->fraction;
0496           }
0497           lcPair.second.second += std::min(std::pow(lcFraction - scFraction, 2), std::pow(scFraction, 2)) *
0498                                   hitEnergyWeight * invSCEnergyWeight;
0499 #ifdef EDM_ML_DEBUG
0500           LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0501               << "scDetId:\t" << (uint32_t)sc_hitDetId << "\tlayerClusterId:\t" << layerClusterId << "\t"
0502               << "lcfraction,scfraction:\t" << lcFraction << ", " << scFraction << "\t"
0503               << "hitEnergyWeight:\t" << hitEnergyWeight << "\t"
0504               << "current score:\t" << lcPair.second.second << "\t"
0505               << "invSCEnergyWeight:\t" << invSCEnergyWeight << "\n";
0506 #endif
0507         }  // End of loop over LayerClusters linked to hits of this SimCluster
0508       }  // End of loop over hits of SimCluster on a Layer
0509 #ifdef EDM_ML_DEBUG
0510       if (lcsInSimCluster[scId][layerId].layerClusterIdToEnergyAndScore.empty())
0511         LogDebug("LCToSCAssociatorByEnergyScoreImpl") << "SC Id: \t" << scId << "\tLC id:\t-1 "
0512                                                       << "\t score \t-1"
0513                                                       << "\n";
0514 
0515       for (const auto& lcPair : lcsInSimCluster[scId][layerId].layerClusterIdToEnergyAndScore) {
0516         LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0517             << "SC Id: \t" << scId << "\t LC id: \t" << lcPair.first << "\t score \t" << lcPair.second.second
0518             << "\t shared energy:\t" << lcPair.second.first << "\t shared energy fraction:\t"
0519             << (lcPair.second.first / SCenergy) << "\n";
0520       }
0521 #endif
0522     }  // End of loop over layers
0523   }  // End of loop over SimClusters
0524 
0525   return {scsInLayerCluster, lcsInSimCluster};
0526 }
0527 
0528 template <typename HIT>
0529 ticl::RecoToSimCollectionWithSimClusters LCToSCAssociatorByEnergyScoreImpl<HIT>::associateRecoToSim(
0530     const edm::Handle<reco::CaloClusterCollection>& cCCH, const edm::Handle<SimClusterCollection>& sCCH) const {
0531   ticl::RecoToSimCollectionWithSimClusters returnValue(productGetter_);
0532   const auto& links = makeConnections(cCCH, sCCH);
0533 
0534   const auto& scsInLayerCluster = std::get<0>(links);
0535   for (size_t lcId = 0; lcId < scsInLayerCluster.size(); ++lcId) {
0536     for (auto& scPair : scsInLayerCluster[lcId]) {
0537       LogDebug("LCToSCAssociatorByEnergyScoreImpl")
0538           << "layerCluster Id: \t" << lcId << "\t SC id: \t" << scPair.first << "\t score \t" << scPair.second << "\n";
0539       // Fill AssociationMap
0540       returnValue.insert(edm::Ref<reco::CaloClusterCollection>(cCCH, lcId),  // Ref to LC
0541                          std::make_pair(edm::Ref<SimClusterCollection>(sCCH, scPair.first),
0542                                         scPair.second)  // Pair <Ref to SC, score>
0543       );
0544     }
0545   }
0546   return returnValue;
0547 }
0548 
0549 template <typename HIT>
0550 ticl::SimToRecoCollectionWithSimClusters LCToSCAssociatorByEnergyScoreImpl<HIT>::associateSimToReco(
0551     const edm::Handle<reco::CaloClusterCollection>& cCCH, const edm::Handle<SimClusterCollection>& sCCH) const {
0552   ticl::SimToRecoCollectionWithSimClusters returnValue(productGetter_);
0553   const auto& links = makeConnections(cCCH, sCCH);
0554   const auto& lcsInSimCluster = std::get<1>(links);
0555   for (size_t scId = 0; scId < lcsInSimCluster.size(); ++scId) {
0556     for (size_t layerId = 0; layerId < lcsInSimCluster[scId].size(); ++layerId) {
0557       for (auto& lcPair : lcsInSimCluster[scId][layerId].layerClusterIdToEnergyAndScore) {
0558         returnValue.insert(
0559             edm::Ref<SimClusterCollection>(sCCH, scId),                                // Ref to SC
0560             std::make_pair(edm::Ref<reco::CaloClusterCollection>(cCCH, lcPair.first),  // Pair <Ref to LC,
0561                            std::make_pair(lcPair.second.first, lcPair.second.second))  // pair <energy, score> >
0562         );
0563       }
0564     }
0565   }
0566   return returnValue;
0567 }
0568 
0569 template class LCToSCAssociatorByEnergyScoreImpl<HGCRecHit>;
0570 template class LCToSCAssociatorByEnergyScoreImpl<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!