Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:23:54 

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