Back to home page

Project CMSSW displayed by LXR
 
 

    

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

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