diff --git a/DataFormats/HGCRecHit/src/classes_def.xml b/DataFormats/HGCRecHit/src/classes_def.xml index f2d3e91a3401e..13a1f5fdd1750 100644 --- a/DataFormats/HGCRecHit/src/classes_def.xml +++ b/DataFormats/HGCRecHit/src/classes_def.xml @@ -35,4 +35,8 @@ + + + + diff --git a/RecoLocalCalo/HGCalRecProducers/plugins/HGCalRecHitMapProducer.cc b/RecoLocalCalo/HGCalRecProducers/plugins/HGCalRecHitMapProducer.cc new file mode 100644 index 0000000000000..0c461379a3bef --- /dev/null +++ b/RecoLocalCalo/HGCalRecProducers/plugins/HGCalRecHitMapProducer.cc @@ -0,0 +1,64 @@ +// user include files +#include + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/global/EDProducer.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" + +#include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h" + +class HGCalRecHitMapProducer : public edm::global::EDProducer<> { +public: + HGCalRecHitMapProducer(const edm::ParameterSet&); + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + + void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override; + +private: + const edm::EDGetTokenT hits_ee_token_; + const edm::EDGetTokenT hits_fh_token_; + const edm::EDGetTokenT hits_bh_token_; +}; + +DEFINE_FWK_MODULE(HGCalRecHitMapProducer); + +HGCalRecHitMapProducer::HGCalRecHitMapProducer(const edm::ParameterSet& ps) + : hits_ee_token_(consumes(ps.getParameter("EEInput"))), + hits_fh_token_(consumes(ps.getParameter("FHInput"))), + hits_bh_token_(consumes(ps.getParameter("BHInput"))) { + produces>(); +} + +void HGCalRecHitMapProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("EEInput", {"HGCalRecHit", "HGCEERecHits"}); + desc.add("FHInput", {"HGCalRecHit", "HGCHEFRecHits"}); + desc.add("BHInput", {"HGCalRecHit", "HGCHEBRecHits"}); + descriptions.add("hgcalRecHitMapProducer", desc); +} + +void HGCalRecHitMapProducer::produce(edm::StreamID, edm::Event& evt, const edm::EventSetup& es) const { + auto hitMap = std::make_unique>(); + const auto& ee_hits = evt.get(hits_ee_token_); + const auto& fh_hits = evt.get(hits_fh_token_); + const auto& bh_hits = evt.get(hits_bh_token_); + + for (const auto& hit : ee_hits) { + hitMap->emplace(hit.detid(), &hit); + } + + for (const auto& hit : fh_hits) { + hitMap->emplace(hit.detid(), &hit); + } + + for (const auto& hit : bh_hits) { + hitMap->emplace(hit.detid(), &hit); + } + evt.put(std::move(hitMap)); +} diff --git a/SimCalorimetry/HGCalAssociatorProducers/plugins/BuildFile.xml b/SimCalorimetry/HGCalAssociatorProducers/plugins/BuildFile.xml new file mode 100644 index 0000000000000..776bbfac5b8cd --- /dev/null +++ b/SimCalorimetry/HGCalAssociatorProducers/plugins/BuildFile.xml @@ -0,0 +1,7 @@ + + + + + + + diff --git a/SimCalorimetry/HGCalAssociatorProducers/plugins/LayerClusterAssociatorByEnergyScoreImpl.cc b/SimCalorimetry/HGCalAssociatorProducers/plugins/LayerClusterAssociatorByEnergyScoreImpl.cc new file mode 100644 index 0000000000000..2a13fde87e2b7 --- /dev/null +++ b/SimCalorimetry/HGCalAssociatorProducers/plugins/LayerClusterAssociatorByEnergyScoreImpl.cc @@ -0,0 +1,526 @@ +// Original Author: Marco Rovere +// + +#include "LayerClusterAssociatorByEnergyScoreImpl.h" + +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h" +#include "SimDataFormats/CaloAnalysis/interface/SimCluster.h" + +LayerClusterAssociatorByEnergyScoreImpl::LayerClusterAssociatorByEnergyScoreImpl( + edm::EDProductGetter const& productGetter, + bool hardScatterOnly, + std::shared_ptr recHitTools, + const std::map*& hitMap) + : hardScatterOnly_(hardScatterOnly), recHitTools_(recHitTools), hitMap_(hitMap), productGetter_(&productGetter) { + layers_ = recHitTools_->lastLayerBH(); +} + +hgcal::association LayerClusterAssociatorByEnergyScoreImpl::makeConnections( + const edm::Handle& cCCH, const edm::Handle& cPCH) const { + // 1. Extract collections and filter CaloParticles, if required + const auto& clusters = *cCCH.product(); + const auto& caloParticles = *cPCH.product(); + auto nLayerClusters = clusters.size(); + //Consider CaloParticles coming from the hard scatterer, excluding the PU contribution. + auto nCaloParticles = caloParticles.size(); + std::vector cPIndices; + //Consider CaloParticles coming from the hard scatterer + //excluding the PU contribution and save the indices. + for (unsigned int cpId = 0; cpId < nCaloParticles; ++cpId) { + if (hardScatterOnly_ && (caloParticles[cpId].g4Tracks()[0].eventId().event() != 0 or + caloParticles[cpId].g4Tracks()[0].eventId().bunchCrossing() != 0)) { + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") + << "Excluding CaloParticles from event: " << caloParticles[cpId].g4Tracks()[0].eventId().event() + << " with BX: " << caloParticles[cpId].g4Tracks()[0].eventId().bunchCrossing() << std::endl; + continue; + } + cPIndices.emplace_back(cpId); + } + nCaloParticles = cPIndices.size(); + + // Initialize cPOnLayer. To be returned outside, since it contains the + // information to compute the CaloParticle-To-LayerCluster score. + hgcal::caloParticleToLayerCluster cPOnLayer; + cPOnLayer.resize(nCaloParticles); + for (unsigned int i = 0; i < nCaloParticles; ++i) { + cPOnLayer[i].resize(layers_ * 2); + for (unsigned int j = 0; j < layers_ * 2; ++j) { + cPOnLayer[i][j].caloParticleId = i; + cPOnLayer[i][j].energy = 0.f; + cPOnLayer[i][j].hits_and_fractions.clear(); + //cPOnLayer[i][j].layerClusterIdToEnergyAndScore.reserve(nLayerClusters); // Not necessary but may improve performance + } + } + + // Fill detIdToCaloParticleId_Map and update cPOnLayer + std::unordered_map> detIdToCaloParticleId_Map; + for (const auto& cpId : cPIndices) { + const SimClusterRefVector& simClusterRefVector = caloParticles[cpId].simClusters(); + for (const auto& it_sc : simClusterRefVector) { + const SimCluster& simCluster = (*(it_sc)); + const auto& hits_and_fractions = simCluster.hits_and_fractions(); + for (const auto& it_haf : hits_and_fractions) { + const auto hitid = (it_haf.first); + const auto cpLayerId = + recHitTools_->getLayerWithOffset(hitid) + layers_ * ((recHitTools_->zside(hitid) + 1) >> 1) - 1; + const auto itcheck = hitMap_->find(hitid); + if (itcheck != hitMap_->end()) { + auto hit_find_it = detIdToCaloParticleId_Map.find(hitid); + if (hit_find_it == detIdToCaloParticleId_Map.end()) { + detIdToCaloParticleId_Map[hitid] = std::vector(); + detIdToCaloParticleId_Map[hitid].emplace_back(cpId, it_haf.second); + } else { + auto findHitIt = std::find(detIdToCaloParticleId_Map[hitid].begin(), + detIdToCaloParticleId_Map[hitid].end(), + hgcal::detIdInfoInCluster{cpId, it_haf.second}); + if (findHitIt != detIdToCaloParticleId_Map[hitid].end()) { + findHitIt->fraction += it_haf.second; + } else { + detIdToCaloParticleId_Map[hitid].emplace_back(cpId, it_haf.second); + } + } + const HGCRecHit* hit = itcheck->second; + cPOnLayer[cpId][cpLayerId].energy += it_haf.second * hit->energy(); + // We need to compress the hits and fractions in order to have a + // reasonable score between CP and LC. Imagine, for example, that a + // CP has detID X used by 2 SimClusters with different fractions. If + // a single LC uses X with fraction 1 and is compared to the 2 + // contributions separately, it will be assigned a score != 0, which + // is wrong. + auto& haf = cPOnLayer[cpId][cpLayerId].hits_and_fractions; + auto found = std::find_if( + std::begin(haf), std::end(haf), [&hitid](const std::pair& v) { return v.first == hitid; }); + if (found != haf.end()) { + found->second += it_haf.second; + } else { + cPOnLayer[cpId][cpLayerId].hits_and_fractions.emplace_back(hitid, it_haf.second); + } + } + } + } + } + +#ifdef EDM_ML_DEBUG + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << "cPOnLayer INFO" << std::endl; + for (size_t cp = 0; cp < cPOnLayer.size(); ++cp) { + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << "For CaloParticle Idx: " << cp << " we have: " << std::endl; + for (size_t cpp = 0; cpp < cPOnLayer[cp].size(); ++cpp) { + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << " On Layer: " << cpp << " we have:" << std::endl; + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") + << " CaloParticleIdx: " << cPOnLayer[cp][cpp].caloParticleId << std::endl; + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") + << " Energy: " << cPOnLayer[cp][cpp].energy << std::endl; + double tot_energy = 0.; + for (auto const& haf : cPOnLayer[cp][cpp].hits_and_fractions) { + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") + << " Hits/fraction/energy: " << (uint32_t)haf.first << "/" << haf.second << "/" + << haf.second * hitMap_->at(haf.first)->energy() << std::endl; + tot_energy += haf.second * hitMap_->at(haf.first)->energy(); + } + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << " Tot Sum haf: " << tot_energy << std::endl; + for (auto const& lc : cPOnLayer[cp][cpp].layerClusterIdToEnergyAndScore) { + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << " lcIdx/energy/score: " << lc.first << "/" + << lc.second.first << "/" << lc.second.second << std::endl; + } + } + } + + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << "detIdToCaloParticleId_Map INFO" << std::endl; + for (auto const& cp : detIdToCaloParticleId_Map) { + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") + << "For detId: " << (uint32_t)cp.first + << " we have found the following connections with CaloParticles:" << std::endl; + for (auto const& cpp : cp.second) { + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") + << " CaloParticle Id: " << cpp.clusterId << " with fraction: " << cpp.fraction + << " and energy: " << cpp.fraction * hitMap_->at(cp.first)->energy() << std::endl; + } + } +#endif + + // Fill detIdToLayerClusterId_Map and cpsInLayerCluster; update cPOnLayer + std::unordered_map> detIdToLayerClusterId_Map; + // this contains the ids of the caloparticles contributing with at least one + // hit to the layer cluster and the reconstruction error. To be returned + // since this contains the information to compute the + // LayerCluster-To-CaloParticle score. + hgcal::layerClusterToCaloParticle cpsInLayerCluster; + cpsInLayerCluster.resize(nLayerClusters); + + for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) { + const std::vector>& hits_and_fractions = clusters[lcId].hitsAndFractions(); + unsigned int numberOfHitsInLC = hits_and_fractions.size(); + const auto firstHitDetId = hits_and_fractions[0].first; + int lcLayerId = + recHitTools_->getLayerWithOffset(firstHitDetId) + layers_ * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1; + + for (unsigned int hitId = 0; hitId < numberOfHitsInLC; hitId++) { + const auto rh_detid = hits_and_fractions[hitId].first; + const auto rhFraction = hits_and_fractions[hitId].second; + + auto hit_find_in_LC = detIdToLayerClusterId_Map.find(rh_detid); + if (hit_find_in_LC == detIdToLayerClusterId_Map.end()) { + detIdToLayerClusterId_Map[rh_detid] = std::vector(); + } + detIdToLayerClusterId_Map[rh_detid].emplace_back(lcId, rhFraction); + + auto hit_find_in_CP = detIdToCaloParticleId_Map.find(rh_detid); + + if (hit_find_in_CP != detIdToCaloParticleId_Map.end()) { + const auto itcheck = hitMap_->find(rh_detid); + const HGCRecHit* hit = itcheck->second; + for (auto& h : hit_find_in_CP->second) { + cPOnLayer[h.clusterId][lcLayerId].layerClusterIdToEnergyAndScore[lcId].first += h.fraction * hit->energy(); + cpsInLayerCluster[lcId].emplace_back(h.clusterId, 0.f); + } + } + } // End loop over hits on a LayerCluster + } // End of loop over LayerClusters + +#ifdef EDM_ML_DEBUG + for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) { + const auto& hits_and_fractions = clusters[lcId].hitsAndFractions(); + unsigned int numberOfHitsInLC = hits_and_fractions.size(); + const auto firstHitDetId = hits_and_fractions[0].first; + int lcLayerId = + recHitTools_->getLayerWithOffset(firstHitDetId) + layers_ * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1; + + // This vector will store, for each hit in the Layercluster, the index of + // the CaloParticle that contributed the most, in terms of energy, to it. + // Special values are: + // + // -2 --> the reconstruction fraction of the RecHit is 0 (used in the past to monitor Halo Hits) + // -3 --> same as before with the added condition that no CaloParticle has been linked to this RecHit + // -1 --> the reco fraction is >0, but no CaloParticle has been linked to it + // >=0 --> index of the linked CaloParticle + std::vector hitsToCaloParticleId(numberOfHitsInLC); + // This will store the index of the CaloParticle linked to the LayerCluster that has the most number of hits in common. + int maxCPId_byNumberOfHits = -1; + // This will store the maximum number of shared hits between a Layercluster andd a CaloParticle + unsigned int maxCPNumberOfHitsInLC = 0; + // This will store the index of the CaloParticle linked to the LayerCluster that has the most energy in common. + int maxCPId_byEnergy = -1; + // This will store the maximum number of shared energy between a Layercluster and a CaloParticle + float maxEnergySharedLCandCP = 0.f; + // This will store the fraction of the LayerCluster energy shared with the best(energy) CaloParticle: e_shared/lc_energy + float energyFractionOfLCinCP = 0.f; + // This will store the fraction of the CaloParticle energy shared with the LayerCluster: e_shared/cp_energy + float energyFractionOfCPinLC = 0.f; + std::unordered_map occurrencesCPinLC; + unsigned int numberOfNoiseHitsInLC = 0; + std::unordered_map CPEnergyInLC; + + for (unsigned int hitId = 0; hitId < numberOfHitsInLC; hitId++) { + const auto rh_detid = hits_and_fractions[hitId].first; + const auto rhFraction = hits_and_fractions[hitId].second; + + auto hit_find_in_CP = detIdToCaloParticleId_Map.find(rh_detid); + + // if the fraction is zero or the hit does not belong to any calo + // particle, set the caloparticleId for the hit to -1 this will + // contribute to the number of noise hits + + // MR Remove the case in which the fraction is 0, since this could be a + // real hit that has been marked as halo. + if (rhFraction == 0.) { + hitsToCaloParticleId[hitId] = -2; + } + if (hit_find_in_CP == detIdToCaloParticleId_Map.end()) { + hitsToCaloParticleId[hitId] -= 1; + } else { + const auto itcheck = hitMap_->find(rh_detid); + const HGCRecHit* hit = itcheck->second; + auto maxCPEnergyInLC = 0.f; + auto maxCPId = -1; + for (auto& h : hit_find_in_CP->second) { + CPEnergyInLC[h.clusterId] += h.fraction * hit->energy(); + // Keep track of which CaloParticle ccontributed the most, in terms + // of energy, to this specific LayerCluster. + if (CPEnergyInLC[h.clusterId] > maxCPEnergyInLC) { + maxCPEnergyInLC = CPEnergyInLC[h.clusterId]; + maxCPId = h.clusterId; + } + } + hitsToCaloParticleId[hitId] = maxCPId; + } + } // End loop over hits on a LayerCluster + + for (const auto& c : hitsToCaloParticleId) { + if (c < 0) { + numberOfNoiseHitsInLC++; + } else { + occurrencesCPinLC[c]++; + } + } + + for (const auto& c : occurrencesCPinLC) { + if (c.second > maxCPNumberOfHitsInLC) { + maxCPId_byNumberOfHits = c.first; + maxCPNumberOfHitsInLC = c.second; + } + } + + for (const auto& c : CPEnergyInLC) { + if (c.second > maxEnergySharedLCandCP) { + maxCPId_byEnergy = c.first; + maxEnergySharedLCandCP = c.second; + } + } + + float totalCPEnergyOnLayer = 0.f; + if (maxCPId_byEnergy >= 0) { + totalCPEnergyOnLayer = cPOnLayer[maxCPId_byEnergy][lcLayerId].energy; + energyFractionOfCPinLC = maxEnergySharedLCandCP / totalCPEnergyOnLayer; + if (clusters[lcId].energy() > 0.f) { + energyFractionOfLCinCP = maxEnergySharedLCandCP / clusters[lcId].energy(); + } + } + + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << std::setw(10) << "LayerId:" + << "\t" << std::setw(12) << "layerCluster" + << "\t" << std::setw(10) << "lc energy" + << "\t" << std::setw(5) << "nhits" + << "\t" << std::setw(12) << "noise hits" + << "\t" << std::setw(22) << "maxCPId_byNumberOfHits" + << "\t" << std::setw(8) << "nhitsCP" + << "\t" << std::setw(13) << "maxCPId_byEnergy" + << "\t" << std::setw(20) << "maxEnergySharedLCandCP" + << "\t" << std::setw(22) << "totalCPEnergyOnLayer" + << "\t" << std::setw(22) << "energyFractionOfLCinCP" + << "\t" << std::setw(25) << "energyFractionOfCPinLC" + << "\t" + << "\n"; + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") + << std::setw(10) << lcLayerId << "\t" << std::setw(12) << lcId << "\t" << std::setw(10) + << clusters[lcId].energy() << "\t" << std::setw(5) << numberOfHitsInLC << "\t" << std::setw(12) + << numberOfNoiseHitsInLC << "\t" << std::setw(22) << maxCPId_byNumberOfHits << "\t" << std::setw(8) + << maxCPNumberOfHitsInLC << "\t" << std::setw(13) << maxCPId_byEnergy << "\t" << std::setw(20) + << maxEnergySharedLCandCP << "\t" << std::setw(22) << totalCPEnergyOnLayer << "\t" << std::setw(22) + << energyFractionOfLCinCP << "\t" << std::setw(25) << energyFractionOfCPinLC << "\n"; + } // End of loop over LayerClusters + + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << "Improved cPOnLayer INFO" << std::endl; + for (size_t cp = 0; cp < cPOnLayer.size(); ++cp) { + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << "For CaloParticle Idx: " << cp << " we have: " << std::endl; + for (size_t cpp = 0; cpp < cPOnLayer[cp].size(); ++cpp) { + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << " On Layer: " << cpp << " we have:" << std::endl; + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") + << " CaloParticleIdx: " << cPOnLayer[cp][cpp].caloParticleId << std::endl; + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") + << " Energy: " << cPOnLayer[cp][cpp].energy << std::endl; + double tot_energy = 0.; + for (auto const& haf : cPOnLayer[cp][cpp].hits_and_fractions) { + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") + << " Hits/fraction/energy: " << (uint32_t)haf.first << "/" << haf.second << "/" + << haf.second * hitMap_->at(haf.first)->energy() << std::endl; + tot_energy += haf.second * hitMap_->at(haf.first)->energy(); + } + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << " Tot Sum haf: " << tot_energy << std::endl; + for (auto const& lc : cPOnLayer[cp][cpp].layerClusterIdToEnergyAndScore) { + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << " lcIdx/energy/score: " << lc.first << "/" + << lc.second.first << "/" << lc.second.second << std::endl; + } + } + } + + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << "Improved detIdToCaloParticleId_Map INFO" << std::endl; + for (auto const& cp : detIdToCaloParticleId_Map) { + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") + << "For detId: " << (uint32_t)cp.first + << " we have found the following connections with CaloParticles:" << std::endl; + for (auto const& cpp : cp.second) { + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") + << " CaloParticle Id: " << cpp.clusterId << " with fraction: " << cpp.fraction + << " and energy: " << cpp.fraction * hitMap_->at(cp.first)->energy() << std::endl; + } + } +#endif + + // Update cpsInLayerCluster; compute the score LayerCluster-to-CaloParticle, + // together with the returned AssociationMap + for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) { + // find the unique caloparticles id contributing to the layer clusters + std::sort(cpsInLayerCluster[lcId].begin(), cpsInLayerCluster[lcId].end()); + auto last = std::unique(cpsInLayerCluster[lcId].begin(), cpsInLayerCluster[lcId].end()); + cpsInLayerCluster[lcId].erase(last, cpsInLayerCluster[lcId].end()); + const auto& hits_and_fractions = clusters[lcId].hitsAndFractions(); + unsigned int numberOfHitsInLC = hits_and_fractions.size(); + // If a reconstructed LayerCluster has energy 0 but is linked to a + // CaloParticle, assigned score 1 + if (clusters[lcId].energy() == 0. && !cpsInLayerCluster[lcId].empty()) { + for (auto& cpPair : cpsInLayerCluster[lcId]) { + cpPair.second = 1.; + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << "layerClusterId : \t " << lcId << "\t CP id : \t" + << cpPair.first << "\t score \t " << cpPair.second << "\n"; + } + continue; + } + + // Compute the correct normalization + float invLayerClusterEnergyWeight = 0.f; + for (auto const& haf : clusters[lcId].hitsAndFractions()) { + invLayerClusterEnergyWeight += + (haf.second * hitMap_->at(haf.first)->energy()) * (haf.second * hitMap_->at(haf.first)->energy()); + } + invLayerClusterEnergyWeight = 1.f / invLayerClusterEnergyWeight; + for (unsigned int i = 0; i < numberOfHitsInLC; ++i) { + DetId rh_detid = hits_and_fractions[i].first; + float rhFraction = hits_and_fractions[i].second; + + bool hitWithNoCP = (detIdToCaloParticleId_Map.find(rh_detid) == detIdToCaloParticleId_Map.end()); + + auto itcheck = hitMap_->find(rh_detid); + const HGCRecHit* hit = itcheck->second; + float hitEnergyWeight = hit->energy() * hit->energy(); + + for (auto& cpPair : cpsInLayerCluster[lcId]) { + float cpFraction = 0.f; + if (!hitWithNoCP) { + auto findHitIt = std::find(detIdToCaloParticleId_Map[rh_detid].begin(), + detIdToCaloParticleId_Map[rh_detid].end(), + hgcal::detIdInfoInCluster{cpPair.first, 0.f}); + if (findHitIt != detIdToCaloParticleId_Map[rh_detid].end()) + cpFraction = findHitIt->fraction; + } + cpPair.second += + (rhFraction - cpFraction) * (rhFraction - cpFraction) * hitEnergyWeight * invLayerClusterEnergyWeight; + } + } // End of loop over Hits within a LayerCluster +#ifdef EDM_ML_DEBUG + if (cpsInLayerCluster[lcId].empty()) + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << "layerCluster Id: \t" << lcId << "\tCP id:\t-1 " + << "\t score \t-1" + << "\n"; +#endif + } // End of loop over LayerClusters + + // Compute the CaloParticle-To-LayerCluster score + for (const auto& cpId : cPIndices) { + for (unsigned int layerId = 0; layerId < layers_ * 2; ++layerId) { + unsigned int CPNumberOfHits = cPOnLayer[cpId][layerId].hits_and_fractions.size(); + if (CPNumberOfHits == 0) + continue; +#ifdef EDM_ML_DEBUG + int lcWithMaxEnergyInCP = -1; + float maxEnergyLCinCP = 0.f; + float CPenergy = cPOnLayer[cpId][layerId].energy; + float CPEnergyFractionInLC = 0.f; + for (auto& lc : cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore) { + if (lc.second.first > maxEnergyLCinCP) { + maxEnergyLCinCP = lc.second.first; + lcWithMaxEnergyInCP = lc.first; + } + } + if (CPenergy > 0.f) + CPEnergyFractionInLC = maxEnergyLCinCP / CPenergy; + + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") + << std::setw(8) << "LayerId:\t" << std::setw(12) << "caloparticle\t" << std::setw(15) << "cp total energy\t" + << std::setw(15) << "cpEnergyOnLayer\t" << std::setw(14) << "CPNhitsOnLayer\t" << std::setw(18) + << "lcWithMaxEnergyInCP\t" << std::setw(15) << "maxEnergyLCinCP\t" << std::setw(20) << "CPEnergyFractionInLC" + << "\n"; + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") + << std::setw(8) << layerId << "\t" << std::setw(12) << cpId << "\t" << std::setw(15) + << caloParticles[cpId].energy() << "\t" << std::setw(15) << CPenergy << "\t" << std::setw(14) + << CPNumberOfHits << "\t" << std::setw(18) << lcWithMaxEnergyInCP << "\t" << std::setw(15) << maxEnergyLCinCP + << "\t" << std::setw(20) << CPEnergyFractionInLC << "\n"; +#endif + // Compute the correct normalization + float invCPEnergyWeight = 0.f; + for (auto const& haf : cPOnLayer[cpId][layerId].hits_and_fractions) { + invCPEnergyWeight += std::pow(haf.second * hitMap_->at(haf.first)->energy(), 2); + } + invCPEnergyWeight = 1.f / invCPEnergyWeight; + for (unsigned int i = 0; i < CPNumberOfHits; ++i) { + auto& cp_hitDetId = cPOnLayer[cpId][layerId].hits_and_fractions[i].first; + auto& cpFraction = cPOnLayer[cpId][layerId].hits_and_fractions[i].second; + + bool hitWithNoLC = false; + if (cpFraction == 0.f) + continue; //hopefully this should never happen + auto hit_find_in_LC = detIdToLayerClusterId_Map.find(cp_hitDetId); + if (hit_find_in_LC == detIdToLayerClusterId_Map.end()) + hitWithNoLC = true; + auto itcheck = hitMap_->find(cp_hitDetId); + const HGCRecHit* hit = itcheck->second; + float hitEnergyWeight = hit->energy() * hit->energy(); + for (auto& lcPair : cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore) { + unsigned int layerClusterId = lcPair.first; + float lcFraction = 0.f; + + if (!hitWithNoLC) { + auto findHitIt = std::find(detIdToLayerClusterId_Map[cp_hitDetId].begin(), + detIdToLayerClusterId_Map[cp_hitDetId].end(), + hgcal::detIdInfoInCluster{layerClusterId, 0.f}); + if (findHitIt != detIdToLayerClusterId_Map[cp_hitDetId].end()) + lcFraction = findHitIt->fraction; + } + lcPair.second.second += + (lcFraction - cpFraction) * (lcFraction - cpFraction) * hitEnergyWeight * invCPEnergyWeight; +#ifdef EDM_ML_DEBUG + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") + << "cpDetId:\t" << (uint32_t)cp_hitDetId << "\tlayerClusterId:\t" << layerClusterId << "\t" + << "lcfraction,cpfraction:\t" << lcFraction << ", " << cpFraction << "\t" + << "hitEnergyWeight:\t" << hitEnergyWeight << "\t" + << "current score:\t" << lcPair.second.second << "\t" + << "invCPEnergyWeight:\t" << invCPEnergyWeight << "\n"; +#endif + } // End of loop over LayerClusters linked to hits of this CaloParticle + } // End of loop over hits of CaloParticle on a Layer +#ifdef EDM_ML_DEBUG + if (cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore.empty()) + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") << "CP Id: \t" << cpId << "\tLC id:\t-1 " + << "\t score \t-1" + << "\n"; + + for (const auto& lcPair : cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore) { + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") + << "CP Id: \t" << cpId << "\t LC id: \t" << lcPair.first << "\t score \t" << lcPair.second.second + << "\t shared energy:\t" << lcPair.second.first << "\t shared energy fraction:\t" + << (lcPair.second.first / CPenergy) << "\n"; + } +#endif + } + } + return {cpsInLayerCluster, cPOnLayer}; +} + +hgcal::RecoToSimCollection LayerClusterAssociatorByEnergyScoreImpl::associateRecoToSim( + const edm::Handle& cCCH, const edm::Handle& cPCH) const { + hgcal::RecoToSimCollection returnValue(productGetter_); + const auto& links = makeConnections(cCCH, cPCH); + + const auto& cpsInLayerCluster = std::get<0>(links); + for (size_t lcId = 0; lcId < cpsInLayerCluster.size(); ++lcId) { + for (auto& cpPair : cpsInLayerCluster[lcId]) { + LogDebug("LayerClusterAssociatorByEnergyScoreImpl") + << "layerCluster Id: \t" << lcId << "\t CP id: \t" << cpPair.first << "\t score \t" << cpPair.second << "\n"; + // Fill AssociationMap + returnValue.insert(edm::Ref(cCCH, lcId), // Ref to LC + std::make_pair(edm::Ref(cPCH, cpPair.first), + cpPair.second) // Pair + ); + } + } + return returnValue; +} + +hgcal::SimToRecoCollection LayerClusterAssociatorByEnergyScoreImpl::associateSimToReco( + const edm::Handle& cCCH, const edm::Handle& cPCH) const { + hgcal::SimToRecoCollection returnValue(productGetter_); + const auto& links = makeConnections(cCCH, cPCH); + const auto& cPOnLayer = std::get<1>(links); + for (size_t cpId = 0; cpId < cPOnLayer.size(); ++cpId) { + for (size_t layerId = 0; layerId < cPOnLayer[cpId].size(); ++layerId) { + for (auto& lcPair : cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore) { + returnValue.insert( + edm::Ref(cPCH, cpId), // Ref to CP + std::make_pair(edm::Ref(cCCH, lcPair.first), // Pair > + ); + } + } + } + return returnValue; +} diff --git a/SimCalorimetry/HGCalAssociatorProducers/plugins/LayerClusterAssociatorByEnergyScoreImpl.h b/SimCalorimetry/HGCalAssociatorProducers/plugins/LayerClusterAssociatorByEnergyScoreImpl.h new file mode 100644 index 0000000000000..c55af520753a5 --- /dev/null +++ b/SimCalorimetry/HGCalAssociatorProducers/plugins/LayerClusterAssociatorByEnergyScoreImpl.h @@ -0,0 +1,68 @@ +// Original Author: Marco Rovere + +#include +#include +#include +#include // shared_ptr + +#include "DataFormats/ForwardDetId/interface/HGCalDetId.h" +#include "DataFormats/HGCRecHit/interface/HGCRecHit.h" +#include "SimDataFormats/Associations/interface/LayerClusterToCaloParticleAssociator.h" +#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h" + +namespace edm { + class EDProductGetter; +} + +namespace hgcal { + struct detIdInfoInCluster { + bool operator==(const detIdInfoInCluster &o) const { return clusterId == o.clusterId; }; + long unsigned int clusterId; + float fraction; + detIdInfoInCluster(long unsigned int cId, float fr) { + clusterId = cId; + fraction = fr; + } + }; + + struct detIdInfoInMultiCluster { + bool operator==(const detIdInfoInMultiCluster &o) const { return multiclusterId == o.multiclusterId; }; + unsigned int multiclusterId; + long unsigned int clusterId; + float fraction; + }; + + struct caloParticleOnLayer { + unsigned int caloParticleId; + float energy = 0; + std::vector> hits_and_fractions; + std::unordered_map> layerClusterIdToEnergyAndScore; + }; + + typedef std::vector>> layerClusterToCaloParticle; + typedef std::vector> caloParticleToLayerCluster; + typedef std::tuple association; +} // namespace hgcal + +class LayerClusterAssociatorByEnergyScoreImpl : public hgcal::LayerClusterToCaloParticleAssociatorBaseImpl { +public: + explicit LayerClusterAssociatorByEnergyScoreImpl(edm::EDProductGetter const &, + bool, + std::shared_ptr, + const std::map *&); + + hgcal::RecoToSimCollection associateRecoToSim(const edm::Handle &cCH, + const edm::Handle &cPCH) const override; + + hgcal::SimToRecoCollection associateSimToReco(const edm::Handle &cCH, + const edm::Handle &cPCH) const override; + +private: + const bool hardScatterOnly_; + std::shared_ptr recHitTools_; + const std::map *hitMap_; + unsigned layers_; + edm::EDProductGetter const *productGetter_; + hgcal::association makeConnections(const edm::Handle &cCH, + const edm::Handle &cPCH) const; +}; diff --git a/SimCalorimetry/HGCalAssociatorProducers/plugins/LayerClusterAssociatorByEnergyScoreProducer.cc b/SimCalorimetry/HGCalAssociatorProducers/plugins/LayerClusterAssociatorByEnergyScoreProducer.cc new file mode 100644 index 0000000000000..9cff1cd72ed8f --- /dev/null +++ b/SimCalorimetry/HGCalAssociatorProducers/plugins/LayerClusterAssociatorByEnergyScoreProducer.cc @@ -0,0 +1,64 @@ +// Original author: Marco Rovere + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/global/EDProducer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/EDGetToken.h" + +#include "SimDataFormats/Associations/interface/LayerClusterToCaloParticleAssociator.h" +#include "LayerClusterAssociatorByEnergyScoreImpl.h" + +class LayerClusterAssociatorByEnergyScoreProducer : public edm::global::EDProducer<> { +public: + explicit LayerClusterAssociatorByEnergyScoreProducer(const edm::ParameterSet &); + ~LayerClusterAssociatorByEnergyScoreProducer() override; + + static void fillDescriptions(edm::ConfigurationDescriptions &descriptions); + +private: + void produce(edm::StreamID, edm::Event &, const edm::EventSetup &) const override; + edm::EDGetTokenT> hitMapToken_; + const bool hardScatterOnly_; + std::shared_ptr rhtools_; +}; + +LayerClusterAssociatorByEnergyScoreProducer::LayerClusterAssociatorByEnergyScoreProducer(const edm::ParameterSet &ps) + : hitMapToken_(consumes>(ps.getParameter("hitMapTag"))), + hardScatterOnly_(ps.getParameter("hardScatterOnly")) { + rhtools_.reset(new hgcal::RecHitTools()); + + // Register the product + produces(); +} + +LayerClusterAssociatorByEnergyScoreProducer::~LayerClusterAssociatorByEnergyScoreProducer() {} + +void LayerClusterAssociatorByEnergyScoreProducer::produce(edm::StreamID, + edm::Event &iEvent, + const edm::EventSetup &es) const { + rhtools_->getEventSetup(es); + edm::Handle> hitMapHandle; + iEvent.getByToken(hitMapToken_, hitMapHandle); + const std::map *hitMap = &*hitMapHandle; + + auto impl = std::make_unique( + iEvent.productGetter(), hardScatterOnly_, rhtools_, hitMap); + auto toPut = std::make_unique(std::move(impl)); + iEvent.put(std::move(toPut)); +} + +void LayerClusterAssociatorByEnergyScoreProducer::fillDescriptions(edm::ConfigurationDescriptions &cfg) { + edm::ParameterSetDescription desc; + desc.add("hitMapTag", edm::InputTag("hgcRecHitMapProducer")); + desc.add("hardScatterOnly", true); + + cfg.add("layerClusterAssociatorByEnergyScore", desc); +} + +//define this as a plug-in +DEFINE_FWK_MODULE(LayerClusterAssociatorByEnergyScoreProducer); diff --git a/SimCalorimetry/HGCalSimProducers/python/hgcHitAssociation_cfi.py b/SimCalorimetry/HGCalSimProducers/python/hgcHitAssociation_cfi.py new file mode 100644 index 0000000000000..fb7ecd1f31074 --- /dev/null +++ b/SimCalorimetry/HGCalSimProducers/python/hgcHitAssociation_cfi.py @@ -0,0 +1,3 @@ +from SimCalorimetry.HGCalAssociatorProducers.layerClusterAssociatorByEnergyScore_cfi import layerClusterAssociatorByEnergyScore as lcAssocByEnergyScoreProducer + +from RecoLocalCalo.HGCalRecProducers.hgcalRecHitMapProducer_cfi import hgcalRecHitMapProducer as hgcRecHitMapProducer diff --git a/SimDataFormats/Associations/BuildFile.xml b/SimDataFormats/Associations/BuildFile.xml index 4ef4e81b0a335..355c655076159 100644 --- a/SimDataFormats/Associations/BuildFile.xml +++ b/SimDataFormats/Associations/BuildFile.xml @@ -1,5 +1,7 @@ + + diff --git a/SimDataFormats/Associations/interface/LayerClusterToCaloParticleAssociator.h b/SimDataFormats/Associations/interface/LayerClusterToCaloParticleAssociator.h new file mode 100644 index 0000000000000..7a0a1ada18f9e --- /dev/null +++ b/SimDataFormats/Associations/interface/LayerClusterToCaloParticleAssociator.h @@ -0,0 +1,48 @@ +#ifndef SimDataFormats_Associations_LayerClusterToCaloParticleAssociator_h +#define SimDataFormats_Associations_LayerClusterToCaloParticleAssociator_h +// Original Author: Marco Rovere + +// system include files +#include + +// user include files + +#include "SimDataFormats/Associations/interface/LayerClusterToCaloParticleAssociatorBaseImpl.h" + +// forward declarations + +namespace hgcal { + + class LayerClusterToCaloParticleAssociator { + public: + LayerClusterToCaloParticleAssociator(std::unique_ptr); + LayerClusterToCaloParticleAssociator() = default; + LayerClusterToCaloParticleAssociator(LayerClusterToCaloParticleAssociator &&) = default; + LayerClusterToCaloParticleAssociator &operator=(LayerClusterToCaloParticleAssociator &&) = default; + ~LayerClusterToCaloParticleAssociator() = default; + + // ---------- const member functions --------------------- + /// Associate a LayerCluster to CaloParticles + hgcal::RecoToSimCollection associateRecoToSim(const edm::Handle &cCCH, + const edm::Handle &cPCH) const { + return m_impl->associateRecoToSim(cCCH, cPCH); + }; + + /// Associate a CaloParticle to LayerClusters + hgcal::SimToRecoCollection associateSimToReco(const edm::Handle &cCCH, + const edm::Handle &cPCH) const { + return m_impl->associateSimToReco(cCCH, cPCH); + } + + private: + LayerClusterToCaloParticleAssociator(const LayerClusterToCaloParticleAssociator &) = delete; // stop default + + const LayerClusterToCaloParticleAssociator &operator=(const LayerClusterToCaloParticleAssociator &) = + delete; // stop default + + // ---------- member data -------------------------------- + std::unique_ptr m_impl; + }; +} // namespace hgcal + +#endif diff --git a/SimDataFormats/Associations/interface/LayerClusterToCaloParticleAssociatorBaseImpl.h b/SimDataFormats/Associations/interface/LayerClusterToCaloParticleAssociatorBaseImpl.h new file mode 100644 index 0000000000000..8ff66876ba3e8 --- /dev/null +++ b/SimDataFormats/Associations/interface/LayerClusterToCaloParticleAssociatorBaseImpl.h @@ -0,0 +1,45 @@ +#ifndef SimDataFormats_Associations_LayerClusterToCaloParticleAssociatorBaseImpl_h +#define SimDataFormats_Associations_LayerClusterToCaloParticleAssociatorBaseImpl_h + +/** \class LayerClusterToCaloParticleAssociatorBaseImpl + * + * Base class for LayerClusterToCaloParticleAssociators. Methods take as input + * the handle of LayerClusters and the CaloParticle collections and return an + * AssociationMap (oneToManyWithQuality) + * + * \author Marco Rovere + */ + +#include "DataFormats/Common/interface/Handle.h" +#include "DataFormats/Common/interface/AssociationMap.h" +#include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h" + +#include "SimDataFormats/CaloAnalysis/interface/CaloParticleFwd.h" + +namespace hgcal { + + typedef edm::AssociationMap< + edm::OneToManyWithQualityGeneric>> + SimToRecoCollection; + typedef edm::AssociationMap< + edm::OneToManyWithQualityGeneric> + RecoToSimCollection; + + class LayerClusterToCaloParticleAssociatorBaseImpl { + public: + /// Constructor + LayerClusterToCaloParticleAssociatorBaseImpl(); + /// Destructor + virtual ~LayerClusterToCaloParticleAssociatorBaseImpl(); + + /// Associate a LayerCluster to CaloParticles + virtual hgcal::RecoToSimCollection associateRecoToSim(const edm::Handle &cCH, + const edm::Handle &cPCH) const; + + /// Associate a CaloParticle to LayerClusters + virtual hgcal::SimToRecoCollection associateSimToReco(const edm::Handle &cCH, + const edm::Handle &cPCH) const; + }; +} // namespace hgcal + +#endif diff --git a/SimDataFormats/Associations/src/LayerClusterToCaloParticleAssociator.cc b/SimDataFormats/Associations/src/LayerClusterToCaloParticleAssociator.cc new file mode 100644 index 0000000000000..38d74d2b4e12f --- /dev/null +++ b/SimDataFormats/Associations/src/LayerClusterToCaloParticleAssociator.cc @@ -0,0 +1,7 @@ +// Original Author: Marco Rovere + +#include "SimDataFormats/Associations/interface/LayerClusterToCaloParticleAssociator.h" + +hgcal::LayerClusterToCaloParticleAssociator::LayerClusterToCaloParticleAssociator( + std::unique_ptr ptr) + : m_impl(std::move(ptr)) {} diff --git a/SimDataFormats/Associations/src/LayerClusterToCaloParticleAssociatorBaseImpl.cc b/SimDataFormats/Associations/src/LayerClusterToCaloParticleAssociatorBaseImpl.cc new file mode 100644 index 0000000000000..ccc25fdd7bac4 --- /dev/null +++ b/SimDataFormats/Associations/src/LayerClusterToCaloParticleAssociatorBaseImpl.cc @@ -0,0 +1,19 @@ +// Original Author: Marco Rovere + +#include "SimDataFormats/Associations/interface/LayerClusterToCaloParticleAssociatorBaseImpl.h" + +namespace hgcal { + LayerClusterToCaloParticleAssociatorBaseImpl::LayerClusterToCaloParticleAssociatorBaseImpl(){}; + LayerClusterToCaloParticleAssociatorBaseImpl::~LayerClusterToCaloParticleAssociatorBaseImpl(){}; + + hgcal::RecoToSimCollection LayerClusterToCaloParticleAssociatorBaseImpl::associateRecoToSim( + const edm::Handle &cCCH, const edm::Handle &cPCH) const { + return hgcal::RecoToSimCollection(); + } + + hgcal::SimToRecoCollection LayerClusterToCaloParticleAssociatorBaseImpl::associateSimToReco( + const edm::Handle &cCCH, const edm::Handle &cPCH) const { + return hgcal::SimToRecoCollection(); + } + +} // namespace hgcal diff --git a/SimDataFormats/Associations/src/classes.h b/SimDataFormats/Associations/src/classes.h index 509321977d5be..20edf2497569e 100644 --- a/SimDataFormats/Associations/src/classes.h +++ b/SimDataFormats/Associations/src/classes.h @@ -6,6 +6,7 @@ #include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h" #include "SimDataFormats/Associations/interface/VertexAssociation.h" #include "SimDataFormats/Associations/interface/VertexToTrackingVertexAssociator.h" +#include "SimDataFormats/Associations/interface/LayerClusterToCaloParticleAssociator.h" namespace SimDataFormats_Associations { struct SimDataFormats_Associations { @@ -16,6 +17,8 @@ namespace SimDataFormats_Associations { edm::Wrapper dummy4; + edm::Wrapper dummy5; + reco::VertexSimToRecoCollection vstrc; reco::VertexSimToRecoCollection::const_iterator vstrci; edm::Wrapper wvstrc; diff --git a/SimDataFormats/Associations/src/classes_def.xml b/SimDataFormats/Associations/src/classes_def.xml index 19266c9ba4d5c..5080eb4b92b75 100644 --- a/SimDataFormats/Associations/src/classes_def.xml +++ b/SimDataFormats/Associations/src/classes_def.xml @@ -9,6 +9,9 @@ + + + diff --git a/Utilities/ReleaseScripts/scripts/duplicateReflexLibrarySearch.py b/Utilities/ReleaseScripts/scripts/duplicateReflexLibrarySearch.py index 68ea71032bb47..a167a3c701509 100755 --- a/Utilities/ReleaseScripts/scripts/duplicateReflexLibrarySearch.py +++ b/Utilities/ReleaseScripts/scripts/duplicateReflexLibrarySearch.py @@ -78,6 +78,7 @@ {'TauReco' : ['reco::PFJetRef']}, {'JetReco' : ['reco::.*Jet','reco::.*Jet(Collection|Ref)']}, {'HGCDigi' : ['HGCSample']}, + {'HGCRecHit' : ['constHGCRecHit','HGCRecHit']}, {'SiPixelObjects' : ['SiPixelQuality.*']}, ] diff --git a/Validation/Configuration/python/autoValidation.py b/Validation/Configuration/python/autoValidation.py index 24e74a58bd6a4..cf6b26a8f979a 100644 --- a/Validation/Configuration/python/autoValidation.py +++ b/Validation/Configuration/python/autoValidation.py @@ -14,7 +14,7 @@ 'miniAODValidation' : ['prevalidationMiniAOD','validationMiniAOD','validationHarvestingMiniAOD'], 'standardValidation' : ['prevalidation','validation','validationHarvesting'], 'standardValidationNoHLT' : ['prevalidationNoHLT','validationNoHLT','validationHarvestingNoHLT'], - 'HGCalValidation' : ['', 'globalValidationHGCal', 'hgcalValidatorPostProcessor'], + 'HGCalValidation' : ['globalPrevalidationHGCal', 'globalValidationHGCal', 'hgcalValidatorPostProcessor'], 'MTDValidation' : ['', 'globalValidationMTD', 'mtdValidationPostProcessor'], 'OuterTrackerValidation' : ['', 'globalValidationOuterTracker', 'postValidationOuterTracker'], 'ecalValidation_phase2' : ['', 'validationECALPhase2', ''], diff --git a/Validation/Configuration/python/globalValidation_cff.py b/Validation/Configuration/python/globalValidation_cff.py index d46e656c8bc9d..ee633de079a2f 100644 --- a/Validation/Configuration/python/globalValidation_cff.py +++ b/Validation/Configuration/python/globalValidation_cff.py @@ -182,6 +182,7 @@ ) globalValidationHGCal = cms.Sequence(hgcalValidation) +globalPrevalidationHGCal = cms.Sequence(hgcalAssociators) globalValidationMTD = cms.Sequence() diff --git a/Validation/Configuration/python/hgcalSimValid_cff.py b/Validation/Configuration/python/hgcalSimValid_cff.py index 432682dfd04f4..250a5e9b26ee5 100644 --- a/Validation/Configuration/python/hgcalSimValid_cff.py +++ b/Validation/Configuration/python/hgcalSimValid_cff.py @@ -1,5 +1,6 @@ import FWCore.ParameterSet.Config as cms +from SimCalorimetry.HGCalSimProducers.hgcHitAssociation_cfi import * from Validation.HGCalValidation.simhitValidation_cff import * from Validation.HGCalValidation.digiValidation_cff import * from Validation.HGCalValidation.rechitValidation_cff import * @@ -16,6 +17,9 @@ VariablePtBins=[10., 30., 80., 120., 250., 600.], DeltaPtOvPtHistoParameter = dict(EROn=True,EREtaMax=3.0, EREtaMin=1.6, slicingOn=True)) +hgcalAssociators = cms.Task(hgcRecHitMapProducer + , lcAssocByEnergyScoreProducer) + hgcalValidation = cms.Sequence(hgcalSimHitValidationEE + hgcalSimHitValidationHEF + hgcalSimHitValidationHEB diff --git a/Validation/HGCalValidation/interface/HGCalValidator.h b/Validation/HGCalValidation/interface/HGCalValidator.h index 19579f00c7073..01ec4b3ecf2de 100644 --- a/Validation/HGCalValidation/interface/HGCalValidator.h +++ b/Validation/HGCalValidation/interface/HGCalValidator.h @@ -27,6 +27,8 @@ #include "Validation/HGCalValidation/interface/CaloParticleSelector.h" #include "RecoLocalCalo/HGCalRecProducers/interface/HGCalClusteringAlgoBase.h" +#include "SimDataFormats/Associations/interface/LayerClusterToCaloParticleAssociator.h" + class PileupSummaryInfo; struct HGCalValidatorHistograms { @@ -74,6 +76,7 @@ class HGCalValidator : public DQMGlobalEDAnalyzer { edm::EDGetTokenT recHitsFH_; edm::EDGetTokenT recHitsBH_; edm::EDGetTokenT density_; + edm::EDGetTokenT LCAssocByEnergyScoreProducer_; std::unique_ptr histoProducerAlgo_; private: diff --git a/Validation/HGCalValidation/interface/HGVHistoProducerAlgo.h b/Validation/HGCalValidation/interface/HGVHistoProducerAlgo.h index a4f6fb33c5654..ab257bd1c9772 100644 --- a/Validation/HGCalValidation/interface/HGVHistoProducerAlgo.h +++ b/Validation/HGCalValidation/interface/HGVHistoProducerAlgo.h @@ -23,6 +23,7 @@ #include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h" #include "SimDataFormats/Vertex/interface/SimVertex.h" #include "RecoLocalCalo/HGCalRecProducers/interface/HGCalClusteringAlgoBase.h" +#include "SimDataFormats/Associations/interface/LayerClusterToCaloParticleAssociator.h" #include "DQMServices/Core/interface/DQMStore.h" @@ -162,13 +163,17 @@ class HGVHistoProducerAlgo { std::vector thicknesses, std::string pathtomatbudfile); void bookMultiClusterHistos(DQMStore::IBooker& ibook, Histograms& histograms, unsigned layers); - void layerClusters_to_CaloParticles(const Histograms& histograms, - const reco::CaloClusterCollection& clusters, - std::vector const& cP, - std::vector const& cPIndices, - std::vector const& cPSelectedIndices, - std::map const&, - unsigned layers) const; + void layerClusters_to_CaloParticles( + const Histograms& histograms, + edm::Handle clusterHandle, + const reco::CaloClusterCollection& clusters, + edm::Handle> caloParticleHandle, + std::vector const& cP, + std::vector const& cPIndices, + std::vector const& cPSelectedIndices, + std::map const&, + unsigned layers, + const edm::Handle& LCAssocByEnergyScoreHandle) const; void multiClusters_to_CaloParticles(const Histograms& histograms, int count, const std::vector& multiClusters, @@ -183,17 +188,21 @@ class HGVHistoProducerAlgo { const CaloParticle& caloparticle, std::vector const& simVertices) const; void fill_cluster_histos(const Histograms& histograms, int count, const reco::CaloCluster& cluster) const; - void fill_generic_cluster_histos(const Histograms& histograms, - int count, - const reco::CaloClusterCollection& clusters, - const Density& densities, - std::vector const& cP, - std::vector const& cPIndices, - std::vector const& cPSelectedIndices, - std::map const&, - std::map cummatbudg, - unsigned layers, - std::vector thicknesses) const; + void fill_generic_cluster_histos( + const Histograms& histograms, + int count, + edm::Handle clusterHandle, + const reco::CaloClusterCollection& clusters, + const Density& densities, + edm::Handle> caloParticleHandle, + std::vector const& cP, + std::vector const& cPIndices, + std::vector const& cPSelectedIndices, + std::map const&, + std::map cummatbudg, + unsigned layers, + std::vector thicknesses, + edm::Handle& LCAssocByEnergyScoreHandle) const; void fill_multi_cluster_histos(const Histograms& histograms, int count, const std::vector& multiClusters, diff --git a/Validation/HGCalValidation/plugins/HGCalValidator.cc b/Validation/HGCalValidation/plugins/HGCalValidator.cc index 8752281a0996f..824a4f5749fc5 100644 --- a/Validation/HGCalValidation/plugins/HGCalValidator.cc +++ b/Validation/HGCalValidation/plugins/HGCalValidator.cc @@ -36,6 +36,9 @@ HGCalValidator::HGCalValidator(const edm::ParameterSet& pset) label_mclTokens.push_back(consumes>(itag)); } + LCAssocByEnergyScoreProducer_ = + consumes(edm::InputTag("lcAssocByEnergyScoreProducer")); + cpSelector = CaloParticleSelector(pset.getParameter("ptMinCP"), pset.getParameter("ptMaxCP"), pset.getParameter("minRapidityCP"), @@ -187,6 +190,9 @@ void HGCalValidator::dqmAnalyze(const edm::Event& event, edm::Handle recHitHandleBH; event.getByToken(recHitsBH_, recHitHandleBH); + edm::Handle LCAssocByEnergyScoreHandle; + event.getByToken(LCAssocByEnergyScoreProducer_, LCAssocByEnergyScoreHandle); + std::map hitMap; fillHitMap(hitMap, *recHitHandleEE, *recHitHandleFH, *recHitHandleBH); @@ -236,15 +242,18 @@ void HGCalValidator::dqmAnalyze(const edm::Event& event, if (dolayerclustersPlots_) { histoProducerAlgo_->fill_generic_cluster_histos(histograms.histoProducerAlgo, w, + clusterHandle, clusters, densities, + caloParticleHandle, caloParticles, cPIndices, selected_cPeff, hitMap, cummatbudg, totallayers_to_monitor_, - thicknesses_to_monitor_); + thicknesses_to_monitor_, + LCAssocByEnergyScoreHandle); for (unsigned int layerclusterIndex = 0; layerclusterIndex < clusters.size(); layerclusterIndex++) { histoProducerAlgo_->fill_cluster_histos(histograms.histoProducerAlgo, w, clusters[layerclusterIndex]); diff --git a/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc b/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc index 86e1ae1aa318f..3359f339acd65 100644 --- a/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc +++ b/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc @@ -808,36 +808,22 @@ void HGVHistoProducerAlgo::fill_cluster_histos(const Histograms& histograms, histograms.h_cluster_eta[count]->Fill(eta); } -void HGVHistoProducerAlgo::layerClusters_to_CaloParticles(const Histograms& histograms, - const reco::CaloClusterCollection& clusters, - std::vector const& cP, - std::vector const& cPIndices, - std::vector const& cPSelectedIndices, - std::map const& hitMap, - unsigned layers) const { +void HGVHistoProducerAlgo::layerClusters_to_CaloParticles( + const Histograms& histograms, + edm::Handle clusterHandle, + const reco::CaloClusterCollection& clusters, + edm::Handle> caloParticleHandle, + std::vector const& cP, + std::vector const& cPIndices, + std::vector const& cPSelectedIndices, + std::map const& hitMap, + unsigned layers, + const edm::Handle& LCAssocByEnergyScoreHandle) const { auto nLayerClusters = clusters.size(); - //Consider CaloParticles coming from the hard scatterer, excluding the PU contribution. - auto nCaloParticles = cPIndices.size(); std::unordered_map> detIdToCaloParticleId_Map; std::unordered_map> detIdToLayerClusterId_Map; - // this contains the ids of the caloparticles contributing with at least one hit to the layer cluster and the reconstruction error - std::vector>> cpsInLayerCluster; - cpsInLayerCluster.resize(nLayerClusters); - - std::unordered_map> cPOnLayer; - // Initialization of cPOnLayer - for (unsigned int i = 0; i < nCaloParticles; ++i) { - auto cpIndex = cPIndices[i]; - cPOnLayer[cpIndex].resize(layers * 2); - for (unsigned int j = 0; j < layers * 2; ++j) { - cPOnLayer[cpIndex][j].caloParticleId = cpIndex; - cPOnLayer[cpIndex][j].energy = 0.f; - cPOnLayer[cpIndex][j].hits_and_fractions.clear(); - } - } - // The association has to be done in an all-vs-all fashion. // For this reason we use the full set of caloParticles, with the only filter on bx for (const auto& cpId : cPIndices) { @@ -847,10 +833,8 @@ void HGVHistoProducerAlgo::layerClusters_to_CaloParticles(const Histograms& hist const auto& hits_and_fractions = simCluster.hits_and_fractions(); for (const auto& it_haf : hits_and_fractions) { DetId hitid = (it_haf.first); - int cpLayerId = recHitTools_->getLayerWithOffset(hitid) + layers * ((recHitTools_->zside(hitid) + 1) >> 1) - 1; std::map::const_iterator itcheck = hitMap.find(hitid); if (itcheck != hitMap.end()) { - const HGCRecHit* hit = itcheck->second; auto hit_find_it = detIdToCaloParticleId_Map.find(hitid); if (hit_find_it == detIdToCaloParticleId_Map.end()) { detIdToCaloParticleId_Map[hitid] = std::vector(); @@ -867,57 +851,11 @@ void HGVHistoProducerAlgo::layerClusters_to_CaloParticles(const Histograms& hist HGVHistoProducerAlgo::detIdInfoInCluster{cpId, it_haf.second}); } } - cPOnLayer[cpId][cpLayerId].energy += it_haf.second * hit->energy(); - // We need to compress the hits and fractions in order to have a - // reasonable score between CP and LC. Imagine, for example, that a - // CP has detID X used by 2 SimClusters with different fractions. If - // a single LC uses X with fraction 1 and is compared to the 2 - // contributions separately, it will be assigned a score != 0, which - // is wrong. - auto& haf = cPOnLayer[cpId][cpLayerId].hits_and_fractions; - auto found = std::find_if( - std::begin(haf), std::end(haf), [&hitid](const std::pair& v) { return v.first == hitid; }); - if (found != haf.end()) { - found->second += it_haf.second; - } else { - cPOnLayer[cpId][cpLayerId].hits_and_fractions.emplace_back(hitid, it_haf.second); - } } } } } - LogDebug("HGCalValidator") << "cPOnLayer INFO" << std::endl; - for (size_t cp = 0; cp < cPOnLayer.size(); ++cp) { - LogDebug("HGCalValidator") << "For CaloParticle Idx: " << cp << " we have: " << std::endl; - for (size_t cpp = 0; cpp < cPOnLayer[cp].size(); ++cpp) { - LogDebug("HGCalValidator") << " On Layer: " << cpp << " we have:" << std::endl; - LogDebug("HGCalValidator") << " CaloParticleIdx: " << cPOnLayer[cp][cpp].caloParticleId << std::endl; - LogDebug("HGCalValidator") << " Energy: " << cPOnLayer[cp][cpp].energy << std::endl; - double tot_energy = 0.; - for (auto const& haf : cPOnLayer[cp][cpp].hits_and_fractions) { - LogDebug("HGCalValidator") << " Hits/fraction/energy: " << (uint32_t)haf.first << "/" << haf.second << "/" - << haf.second * hitMap.at(haf.first)->energy() << std::endl; - tot_energy += haf.second * hitMap.at(haf.first)->energy(); - } - LogDebug("HGCalValidator") << " Tot Sum haf: " << tot_energy << std::endl; - for (auto const& lc : cPOnLayer[cp][cpp].layerClusterIdToEnergyAndScore) { - LogDebug("HGCalValidator") << " lcIdx/energy/score: " << lc.first << "/" << lc.second.first << "/" - << lc.second.second << std::endl; - } - } - } - - LogDebug("HGCalValidator") << "detIdToCaloParticleId_Map INFO" << std::endl; - for (auto const& cp : detIdToCaloParticleId_Map) { - LogDebug("HGCalValidator") << "For detId: " << (uint32_t)cp.first - << " we have found the following connections with CaloParticles:" << std::endl; - for (auto const& cpp : cp.second) { - LogDebug("HGCalValidator") << " CaloParticle Id: " << cpp.clusterId << " with fraction: " << cpp.fraction - << " and energy: " << cpp.fraction * hitMap.at(cp.first)->energy() << std::endl; - } - } - for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) { const std::vector>& hits_and_fractions = clusters[lcId].hitsAndFractions(); unsigned int numberOfHitsInLC = hits_and_fractions.size(); @@ -935,23 +873,8 @@ void HGVHistoProducerAlgo::layerClusters_to_CaloParticles(const Histograms& hist int lcLayerId = recHitTools_->getLayerWithOffset(firstHitDetId) + layers * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1; - // This will store the index of the CaloParticle linked to the LayerCluster that has the most number of hits in common. - int maxCPId_byNumberOfHits = -1; - // This will store the maximum number of shared hits between a Layercluster andd a CaloParticle - unsigned int maxCPNumberOfHitsInLC = 0; - // This will store the index of the CaloParticle linked to the LayerCluster that has the most energy in common. - // - int maxCPId_byEnergy = -1; - // This will store the maximum number of shared energy between a Layercluster and a CaloParticle - float maxEnergySharedLCandCP = 0.f; - // This will store the fraction of the LayerCluster energy shared with the best(energy) CaloParticle: e_shared/lc_energy - float energyFractionOfLCinCP = 0.f; // This will store the fraction of the CaloParticle energy shared with the LayerCluster: e_shared/cp_energy - float energyFractionOfCPinLC = 0.f; - std::unordered_map occurrencesCPinLC; std::unordered_map CPEnergyInLC; - unsigned int numberOfNoiseHitsInLC = 0; - unsigned int numberOfHaloHitsInLC = 0; for (unsigned int hitId = 0; hitId < numberOfHitsInLC; hitId++) { DetId rh_detid = hits_and_fractions[hitId].first; @@ -976,7 +899,6 @@ void HGVHistoProducerAlgo::layerClusters_to_CaloParticles(const Histograms& hist // real hit that has been marked as halo. if (rhFraction == 0.) { hitsToCaloParticleId[hitId] = -2; - numberOfHaloHitsInLC++; } if (hit_find_in_CP == detIdToCaloParticleId_Map.end()) { hitsToCaloParticleId[hitId] -= 1; @@ -985,10 +907,7 @@ void HGVHistoProducerAlgo::layerClusters_to_CaloParticles(const Histograms& hist auto maxCPId = -1; for (auto& h : hit_find_in_CP->second) { CPEnergyInLC[h.clusterId] += h.fraction * hit->energy(); - cPOnLayer[h.clusterId][lcLayerId].layerClusterIdToEnergyAndScore[lcId].first += h.fraction * hit->energy(); - cPOnLayer[h.clusterId][lcLayerId].layerClusterIdToEnergyAndScore[lcId].second = FLT_MAX; - cpsInLayerCluster[lcId].emplace_back(std::make_pair(h.clusterId, FLT_MAX)); - // Keep track of which CaloParticle ccontributed the most, in terms + // Keep track of which CaloParticle contributed the most, in terms // of energy, to this specific LayerCluster. if (CPEnergyInLC[h.clusterId] > maxCPEnergyInLC) { maxCPEnergyInLC = CPEnergyInLC[h.clusterId]; @@ -1001,169 +920,55 @@ void HGVHistoProducerAlgo::layerClusters_to_CaloParticles(const Histograms& hist hitsToCaloParticleId[hitId] > 0. ? 0. : hitsToCaloParticleId[hitId]); } // End loop over hits on a LayerCluster - for (auto& c : hitsToCaloParticleId) { - if (c < 0) { - numberOfNoiseHitsInLC++; - } else { - occurrencesCPinLC[c]++; - } - } - - for (auto& c : occurrencesCPinLC) { - if (c.second > maxCPNumberOfHitsInLC) { - maxCPId_byNumberOfHits = c.first; - maxCPNumberOfHitsInLC = c.second; - } - } - - for (auto& c : CPEnergyInLC) { - if (c.second > maxEnergySharedLCandCP) { - maxCPId_byEnergy = c.first; - maxEnergySharedLCandCP = c.second; - } - } - float totalCPEnergyOnLayer = 0.f; - if (maxCPId_byEnergy >= 0) { - totalCPEnergyOnLayer = cPOnLayer[maxCPId_byEnergy][lcLayerId].energy; - energyFractionOfCPinLC = maxEnergySharedLCandCP / totalCPEnergyOnLayer; - if (clusters[lcId].energy() > 0.f) { - energyFractionOfLCinCP = maxEnergySharedLCandCP / clusters[lcId].energy(); - } - } - LogDebug("HGCalValidator") << std::setw(10) << "LayerId:" - << "\t" << std::setw(12) << "layerCluster" - << "\t" << std::setw(10) << "lc energy" - << "\t" << std::setw(5) << "nhits" - << "\t" << std::setw(12) << "noise hits" - << "\t" << std::setw(22) << "maxCPId_byNumberOfHits" - << "\t" << std::setw(8) << "nhitsCP" - << "\t" << std::setw(13) << "maxCPId_byEnergy" - << "\t" << std::setw(20) << "maxEnergySharedLCandCP" - << "\t" << std::setw(22) << "totalCPEnergyOnLayer" - << "\t" << std::setw(22) << "energyFractionOfLCinCP" - << "\t" << std::setw(25) << "energyFractionOfCPinLC" - << "\t" - << "\n"; - LogDebug("HGCalValidator") << std::setw(10) << lcLayerId << "\t" << std::setw(12) << lcId << "\t" << std::setw(10) - << clusters[lcId].energy() << "\t" << std::setw(5) << numberOfHitsInLC << "\t" - << std::setw(12) << numberOfNoiseHitsInLC << "\t" << std::setw(22) - << maxCPId_byNumberOfHits << "\t" << std::setw(8) << maxCPNumberOfHitsInLC << "\t" - << std::setw(13) << maxCPId_byEnergy << "\t" << std::setw(20) << maxEnergySharedLCandCP - << "\t" << std::setw(22) << totalCPEnergyOnLayer << "\t" << std::setw(22) - << energyFractionOfLCinCP << "\t" << std::setw(25) << energyFractionOfCPinLC << "\n"; } // End of loop over LayerClusters - LogDebug("HGCalValidator") << "Improved cPOnLayer INFO" << std::endl; - for (size_t cp = 0; cp < cPOnLayer.size(); ++cp) { - LogDebug("HGCalValidator") << "For CaloParticle Idx: " << cp << " we have: " << std::endl; - for (size_t cpp = 0; cpp < cPOnLayer[cp].size(); ++cpp) { - LogDebug("HGCalValidator") << " On Layer: " << cpp << " we have:" << std::endl; - LogDebug("HGCalValidator") << " CaloParticleIdx: " << cPOnLayer[cp][cpp].caloParticleId << std::endl; - LogDebug("HGCalValidator") << " Energy: " << cPOnLayer[cp][cpp].energy << std::endl; - double tot_energy = 0.; - for (auto const& haf : cPOnLayer[cp][cpp].hits_and_fractions) { - LogDebug("HGCalValidator") << " Hits/fraction/energy: " << (uint32_t)haf.first << "/" << haf.second << "/" - << haf.second * hitMap.at(haf.first)->energy() << std::endl; - tot_energy += haf.second * hitMap.at(haf.first)->energy(); - } - LogDebug("HGCalValidator") << " Tot Sum haf: " << tot_energy << std::endl; - for (auto const& lc : cPOnLayer[cp][cpp].layerClusterIdToEnergyAndScore) { - LogDebug("HGCalValidator") << " lcIdx/energy/score: " << lc.first << "/" << lc.second.first << "/" - << lc.second.second << std::endl; - } - } - } - - LogDebug("HGCalValidator") << "Improved detIdToCaloParticleId_Map INFO" << std::endl; - for (auto const& cp : detIdToCaloParticleId_Map) { - LogDebug("HGCalValidator") << "For detId: " << (uint32_t)cp.first - << " we have found the following connections with CaloParticles:" << std::endl; - for (auto const& cpp : cp.second) { - LogDebug("HGCalValidator") << " CaloParticle Id: " << cpp.clusterId << " with fraction: " << cpp.fraction - << " and energy: " << cpp.fraction * hitMap.at(cp.first)->energy() << std::endl; - } - } - + hgcal::RecoToSimCollection cpsInLayerClusterMap = + LCAssocByEnergyScoreHandle->associateRecoToSim(clusterHandle, caloParticleHandle); + hgcal::SimToRecoCollection cPOnLayerMap = + LCAssocByEnergyScoreHandle->associateSimToReco(clusterHandle, caloParticleHandle); // Here we do fill the plots to compute the different metrics linked to // reco-level, namely fake-rate an merge-rate. In this loop we should *not* // restrict only to the selected caloParaticles. for (unsigned int lcId = 0; lcId < nLayerClusters; ++lcId) { - // find the unique caloparticles id contributing to the layer clusters - std::sort(cpsInLayerCluster[lcId].begin(), cpsInLayerCluster[lcId].end()); - auto last = std::unique(cpsInLayerCluster[lcId].begin(), cpsInLayerCluster[lcId].end()); - cpsInLayerCluster[lcId].erase(last, cpsInLayerCluster[lcId].end()); const std::vector>& hits_and_fractions = clusters[lcId].hitsAndFractions(); - unsigned int numberOfHitsInLC = hits_and_fractions.size(); - auto firstHitDetId = hits_and_fractions[0].first; - int lcLayerId = + const auto firstHitDetId = hits_and_fractions[0].first; + const int lcLayerId = recHitTools_->getLayerWithOffset(firstHitDetId) + layers * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1; - // If a reconstructed LayerCluster has energy 0 but is linked to a CaloParticle, assigned score 1 - if (clusters[lcId].energy() == 0. && !cpsInLayerCluster[lcId].empty()) { - for (auto& cpPair : cpsInLayerCluster[lcId]) { - cpPair.second = 1.; - LogDebug("HGCalValidator") << "layerCluster Id: \t" << lcId << "\t CP id: \t" << cpPair.first << "\t score \t" - << cpPair.second << "\n"; + histograms.h_denom_layercl_eta_perlayer.at(lcLayerId)->Fill(clusters[lcId].eta()); + histograms.h_denom_layercl_phi_perlayer.at(lcLayerId)->Fill(clusters[lcId].phi()); + // + const edm::Ref lcRef(clusterHandle, lcId); + const auto& cpsIt = cpsInLayerClusterMap.find(lcRef); + if (cpsIt == cpsInLayerClusterMap.end()) + continue; + + const auto& cps = cpsIt->val; + if (clusters[lcId].energy() == 0. && !cps.empty()) { + for (const auto& cpPair : cps) { histograms.h_score_layercl2caloparticle_perlayer.at(lcLayerId)->Fill(cpPair.second); } continue; } - - // Compute the correct normalization - float invLayerClusterEnergyWeight = 0.f; - for (auto const& haf : clusters[lcId].hitsAndFractions()) { - invLayerClusterEnergyWeight += - (haf.second * hitMap.at(haf.first)->energy()) * (haf.second * hitMap.at(haf.first)->energy()); - } - invLayerClusterEnergyWeight = 1.f / invLayerClusterEnergyWeight; - - for (unsigned int i = 0; i < numberOfHitsInLC; ++i) { - DetId rh_detid = hits_and_fractions[i].first; - float rhFraction = hits_and_fractions[i].second; - bool hitWithNoCP = false; - - auto hit_find_in_CP = detIdToCaloParticleId_Map.find(rh_detid); - if (hit_find_in_CP == detIdToCaloParticleId_Map.end()) - hitWithNoCP = true; - auto itcheck = hitMap.find(rh_detid); - const HGCRecHit* hit = itcheck->second; - float hitEnergyWeight = hit->energy() * hit->energy(); - - for (auto& cpPair : cpsInLayerCluster[lcId]) { - float cpFraction = 0.f; - if (!hitWithNoCP) { - auto findHitIt = std::find(detIdToCaloParticleId_Map[rh_detid].begin(), - detIdToCaloParticleId_Map[rh_detid].end(), - HGVHistoProducerAlgo::detIdInfoInCluster{cpPair.first, 0.f}); - if (findHitIt != detIdToCaloParticleId_Map[rh_detid].end()) - cpFraction = findHitIt->fraction; - } - if (cpPair.second == FLT_MAX) { - cpPair.second = 0.f; - } - cpPair.second += - (rhFraction - cpFraction) * (rhFraction - cpFraction) * hitEnergyWeight * invLayerClusterEnergyWeight; - } - } // End of loop over Hits within a LayerCluster - - if (cpsInLayerCluster[lcId].empty()) - LogDebug("HGCalValidator") << "layerCluster Id: \t" << lcId << "\tCP id:\t-1 " - << "\t score \t-1" - << "\n"; - - for (auto& cpPair : cpsInLayerCluster[lcId]) { - LogDebug("HGCalValidator") << "layerCluster Id: \t" << lcId << "\t CP id: \t" << cpPair.first << "\t score \t" - << cpPair.second << "\n"; + for (const auto& cpPair : cps) { + LogDebug("HGCalValidator") << "layerCluster Id: \t" << lcId << "\t CP id: \t" << cpPair.first.index() + << "\t score \t" << cpPair.second << std::endl; histograms.h_score_layercl2caloparticle_perlayer.at(lcLayerId)->Fill(cpPair.second); - auto const& cp_linked = cPOnLayer[cpPair.first][lcLayerId].layerClusterIdToEnergyAndScore[lcId]; + auto const& cp_linked = + std::find_if(std::begin(cPOnLayerMap[cpPair.first]), + std::end(cPOnLayerMap[cpPair.first]), + [&lcRef](const std::pair, std::pair>& p) { + return p.first == lcRef; + }); + if (cp_linked == + cPOnLayerMap[cpPair.first].end()) // This should never happen by construction of the association maps + continue; histograms.h_sharedenergy_layercl2caloparticle_perlayer.at(lcLayerId)->Fill( - cp_linked.first / clusters[lcId].energy(), clusters[lcId].energy()); + cp_linked->second.first / clusters[lcId].energy(), clusters[lcId].energy()); histograms.h_energy_vs_score_layercl2caloparticle_perlayer.at(lcLayerId)->Fill( - cpPair.second, cp_linked.first / clusters[lcId].energy()); + cpPair.second, cp_linked->second.first / clusters[lcId].energy()); } - - auto assoc = std::count_if(std::begin(cpsInLayerCluster[lcId]), - std::end(cpsInLayerCluster[lcId]), - [](const auto& obj) { return obj.second < ScoreCutLCtoCP_; }); + const auto assoc = + std::count_if(std::begin(cps), std::end(cps), [](const auto& obj) { return obj.second < ScoreCutLCtoCP_; }); if (assoc) { histograms.h_num_layercl_eta_perlayer.at(lcLayerId)->Fill(clusters[lcId].eta()); histograms.h_num_layercl_phi_perlayer.at(lcLayerId)->Fill(clusters[lcId].phi()); @@ -1171,119 +976,85 @@ void HGVHistoProducerAlgo::layerClusters_to_CaloParticles(const Histograms& hist histograms.h_numMerge_layercl_eta_perlayer.at(lcLayerId)->Fill(clusters[lcId].eta()); histograms.h_numMerge_layercl_phi_perlayer.at(lcLayerId)->Fill(clusters[lcId].phi()); } - auto best = std::min_element(std::begin(cpsInLayerCluster[lcId]), - std::end(cpsInLayerCluster[lcId]), - [](const auto& obj1, const auto& obj2) { return obj1.second < obj2.second; }); - auto const& best_cp_linked = cPOnLayer[best->first][lcLayerId].layerClusterIdToEnergyAndScore[lcId]; + const auto& best = std::min_element( + std::begin(cps), std::end(cps), [](const auto& obj1, const auto& obj2) { return obj1.second < obj2.second; }); + const auto& best_cp_linked = + std::find_if(std::begin(cPOnLayerMap[best->first]), + std::end(cPOnLayerMap[best->first]), + [&lcRef](const std::pair, std::pair>& p) { + return p.first == lcRef; + }); + if (best_cp_linked == + cPOnLayerMap[best->first].end()) // This should never happen by construction of the association maps + continue; histograms.h_sharedenergy_layercl2caloparticle_vs_eta_perlayer.at(lcLayerId)->Fill( - clusters[lcId].eta(), best_cp_linked.first / clusters[lcId].energy()); + clusters[lcId].eta(), best_cp_linked->second.first / clusters[lcId].energy()); histograms.h_sharedenergy_layercl2caloparticle_vs_phi_perlayer.at(lcLayerId)->Fill( - clusters[lcId].phi(), best_cp_linked.first / clusters[lcId].energy()); + clusters[lcId].phi(), best_cp_linked->second.first / clusters[lcId].energy()); } - histograms.h_denom_layercl_eta_perlayer.at(lcLayerId)->Fill(clusters[lcId].eta()); - histograms.h_denom_layercl_phi_perlayer.at(lcLayerId)->Fill(clusters[lcId].phi()); } // End of loop over LayerClusters // Here we do fill the plots to compute the different metrics linked to - // gen-level, namely efficiency an duplicate. In this loop we should restrict + // gen-level, namely efficiency and duplicate. In this loop we should restrict // only to the selected caloParaticles. for (const auto& cpId : cPSelectedIndices) { - for (unsigned int layerId = 0; layerId < layers * 2; ++layerId) { - unsigned int CPNumberOfHits = cPOnLayer[cpId][layerId].hits_and_fractions.size(); - float CPenergy = cPOnLayer[cpId][layerId].energy; - if (CPNumberOfHits == 0) - continue; - int lcWithMaxEnergyInCP = -1; - float maxEnergyLCinCP = 0.f; - float CPEnergyFractionInLC = 0.f; - for (auto& lc : cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore) { - if (lc.second.first > maxEnergyLCinCP) { - maxEnergyLCinCP = lc.second.first; - lcWithMaxEnergyInCP = lc.first; - } - } - if (CPenergy > 0.f) - CPEnergyFractionInLC = maxEnergyLCinCP / CPenergy; + const edm::Ref cpRef(caloParticleHandle, cpId); + const auto& lcsIt = cPOnLayerMap.find(cpRef); - LogDebug("HGCalValidator") << std::setw(8) << "LayerId:\t" << std::setw(12) << "caloparticle\t" << std::setw(15) - << "cp total energy\t" << std::setw(15) << "cpEnergyOnLayer\t" << std::setw(14) - << "CPNhitsOnLayer\t" << std::setw(18) << "lcWithMaxEnergyInCP\t" << std::setw(15) - << "maxEnergyLCinCP\t" << std::setw(20) << "CPEnergyFractionInLC" - << "\n"; - LogDebug("HGCalValidator") << std::setw(8) << layerId << "\t" << std::setw(12) << cpId << "\t" << std::setw(15) - << cP[cpId].energy() << "\t" << std::setw(15) << CPenergy << "\t" << std::setw(14) - << CPNumberOfHits << "\t" << std::setw(18) << lcWithMaxEnergyInCP << "\t" - << std::setw(15) << maxEnergyLCinCP << "\t" << std::setw(20) << CPEnergyFractionInLC - << "\n"; + std::map cPEnergyOnLayer; + for (unsigned int layerId = 0; layerId < layers * 2; ++layerId) + cPEnergyOnLayer[layerId] = 0; - // Compute the correct normalization - float invCPEnergyWeight = 0.f; - for (auto const& haf : cPOnLayer[cpId][layerId].hits_and_fractions) { - invCPEnergyWeight += - (haf.second * hitMap.at(haf.first)->energy()) * (haf.second * hitMap.at(haf.first)->energy()); + const SimClusterRefVector& simClusterRefVector = cP[cpId].simClusters(); + for (const auto& it_sc : simClusterRefVector) { + const SimCluster& simCluster = (*(it_sc)); + const auto& hits_and_fractions = simCluster.hits_and_fractions(); + for (const auto& it_haf : hits_and_fractions) { + const DetId hitid = (it_haf.first); + const int cpLayerId = + recHitTools_->getLayerWithOffset(hitid) + layers * ((recHitTools_->zside(hitid) + 1) >> 1) - 1; + std::map::const_iterator itcheck = hitMap.find(hitid); + if (itcheck != hitMap.end()) { + const HGCRecHit* hit = itcheck->second; + cPEnergyOnLayer[cpLayerId] += it_haf.second * hit->energy(); + } } - invCPEnergyWeight = 1.f / invCPEnergyWeight; - - for (unsigned int i = 0; i < CPNumberOfHits; ++i) { - auto& cp_hitDetId = cPOnLayer[cpId][layerId].hits_and_fractions[i].first; - auto& cpFraction = cPOnLayer[cpId][layerId].hits_and_fractions[i].second; + } - bool hitWithNoLC = false; - if (cpFraction == 0.f) - continue; //hopefully this should never happen - auto hit_find_in_LC = detIdToLayerClusterId_Map.find(cp_hitDetId); - if (hit_find_in_LC == detIdToLayerClusterId_Map.end()) - hitWithNoLC = true; - auto itcheck = hitMap.find(cp_hitDetId); - const HGCRecHit* hit = itcheck->second; - float hitEnergyWeight = hit->energy() * hit->energy(); - for (auto& lcPair : cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore) { - unsigned int layerClusterId = lcPair.first; - float lcFraction = 0.f; - - if (!hitWithNoLC) { - auto findHitIt = std::find(detIdToLayerClusterId_Map[cp_hitDetId].begin(), - detIdToLayerClusterId_Map[cp_hitDetId].end(), - HGVHistoProducerAlgo::detIdInfoInCluster{layerClusterId, 0.f}); - if (findHitIt != detIdToLayerClusterId_Map[cp_hitDetId].end()) - lcFraction = findHitIt->fraction; - } - // if (lcFraction == 0.) { - // lcFraction = -1.; - // } - if (lcPair.second.second == FLT_MAX) { - lcPair.second.second = 0.f; - } - lcPair.second.second += - (lcFraction - cpFraction) * (lcFraction - cpFraction) * hitEnergyWeight * invCPEnergyWeight; - LogDebug("HGCalValidator") << "cpDetId:\t" << (uint32_t)cp_hitDetId << "\tlayerClusterId:\t" << layerClusterId - << "\t" - << "lcfraction,cpfraction:\t" << lcFraction << ", " << cpFraction << "\t" - << "hitEnergyWeight:\t" << hitEnergyWeight << "\t" - << "current score:\t" << lcPair.second.second << "\t" - << "invCPEnergyWeight:\t" << invCPEnergyWeight << "\n"; - } // End of loop over LayerClusters linked to hits of this CaloParticle - } // End of loop over hits of CaloParticle on a Layer + for (unsigned int layerId = 0; layerId < layers * 2; ++layerId) { + if (!cPEnergyOnLayer[layerId]) + continue; - if (cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore.empty()) - LogDebug("HGCalValidator") << "CP Id: \t" << cpId << "\tLC id:\t-1 " - << "\t score \t-1" - << "\n"; + histograms.h_denom_caloparticle_eta_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().eta()); + histograms.h_denom_caloparticle_phi_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().phi()); - for (auto& lcPair : cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore) { - LogDebug("HGCalValidator") << "CP Id: \t" << cpId << "\t LC id: \t" << lcPair.first << "\t score \t" - << lcPair.second.second << "\t" - << "shared energy:\t" << lcPair.second.first << "\t" - << "shared energy fraction:\t" << (lcPair.second.first / CPenergy) << "\n"; + if (lcsIt == cPOnLayerMap.end()) + continue; + const auto& lcs = lcsIt->val; + + auto getLCLayerId = [&](const unsigned int lcId) { + const std::vector>& hits_and_fractions = clusters[lcId].hitsAndFractions(); + const auto firstHitDetId = hits_and_fractions[0].first; + const unsigned int lcLayerId = recHitTools_->getLayerWithOffset(firstHitDetId) + + layers * ((recHitTools_->zside(firstHitDetId) + 1) >> 1) - 1; + return lcLayerId; + }; + + for (const auto& lcPair : lcs) { + if (getLCLayerId(lcPair.first.index()) != layerId) + continue; histograms.h_score_caloparticle2layercl_perlayer.at(layerId)->Fill(lcPair.second.second); - histograms.h_sharedenergy_caloparticle2layercl_perlayer.at(layerId)->Fill(lcPair.second.first / CPenergy, - CPenergy); - histograms.h_energy_vs_score_caloparticle2layercl_perlayer.at(layerId)->Fill(lcPair.second.second, - lcPair.second.first / CPenergy); + histograms.h_sharedenergy_caloparticle2layercl_perlayer.at(layerId)->Fill( + lcPair.second.first / cPEnergyOnLayer[layerId], cPEnergyOnLayer[layerId]); + histograms.h_energy_vs_score_caloparticle2layercl_perlayer.at(layerId)->Fill( + lcPair.second.second, lcPair.second.first / cPEnergyOnLayer[layerId]); } - auto assoc = std::count_if(std::begin(cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore), - std::end(cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore), - [](const auto& obj) { return obj.second.second < ScoreCutCPtoLC_; }); + const auto assoc = std::count_if(std::begin(lcs), std::end(lcs), [&](const auto& obj) { + if (getLCLayerId(obj.first.index()) != layerId) + return false; + else + return obj.second.second < ScoreCutCPtoLC_; + }); if (assoc) { histograms.h_num_caloparticle_eta_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().eta()); histograms.h_num_caloparticle_phi_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().phi()); @@ -1291,32 +1062,38 @@ void HGVHistoProducerAlgo::layerClusters_to_CaloParticles(const Histograms& hist histograms.h_numDup_caloparticle_eta_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().eta()); histograms.h_numDup_caloparticle_phi_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().phi()); } - auto best = std::min_element( - std::begin(cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore), - std::end(cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore), - [](const auto& obj1, const auto& obj2) { return obj1.second.second < obj2.second.second; }); + const auto best = std::min_element(std::begin(lcs), std::end(lcs), [&](const auto& obj1, const auto& obj2) { + if (getLCLayerId(obj1.first.index()) != layerId) + return false; + else if (getLCLayerId(obj2.first.index()) == layerId) + return obj1.second.second < obj2.second.second; + else + return true; + }); histograms.h_sharedenergy_caloparticle2layercl_vs_eta_perlayer.at(layerId)->Fill( - cP[cpId].g4Tracks()[0].momentum().eta(), best->second.first / CPenergy); + cP[cpId].g4Tracks()[0].momentum().eta(), best->second.first / cPEnergyOnLayer[layerId]); histograms.h_sharedenergy_caloparticle2layercl_vs_phi_perlayer.at(layerId)->Fill( - cP[cpId].g4Tracks()[0].momentum().phi(), best->second.first / CPenergy); + cP[cpId].g4Tracks()[0].momentum().phi(), best->second.first / cPEnergyOnLayer[layerId]); } - histograms.h_denom_caloparticle_eta_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().eta()); - histograms.h_denom_caloparticle_phi_perlayer.at(layerId)->Fill(cP[cpId].g4Tracks()[0].momentum().phi()); } } } -void HGVHistoProducerAlgo::fill_generic_cluster_histos(const Histograms& histograms, - int count, - const reco::CaloClusterCollection& clusters, - const Density& densities, - std::vector const& cP, - std::vector const& cPIndices, - std::vector const& cPSelectedIndices, - std::map const& hitMap, - std::map cummatbudg, - unsigned layers, - std::vector thicknesses) const { +void HGVHistoProducerAlgo::fill_generic_cluster_histos( + const Histograms& histograms, + int count, + edm::Handle clusterHandle, + const reco::CaloClusterCollection& clusters, + const Density& densities, + edm::Handle> caloParticleHandle, + std::vector const& cP, + std::vector const& cPIndices, + std::vector const& cPSelectedIndices, + std::map const& hitMap, + std::map cummatbudg, + unsigned layers, + std::vector thicknesses, + edm::Handle& LCAssocByEnergyScoreHandle) const { //Each event to be treated as two events: an event in +ve endcap, //plus another event in -ve endcap. In this spirit there will be //a layer variable (layerid) that maps the layers in : @@ -1341,7 +1118,16 @@ void HGVHistoProducerAlgo::fill_generic_cluster_histos(const Histograms& histogr tnlcpthplus.insert(std::pair("mixed", 0)); tnlcpthminus.insert(std::pair("mixed", 0)); - layerClusters_to_CaloParticles(histograms, clusters, cP, cPIndices, cPSelectedIndices, hitMap, layers); + layerClusters_to_CaloParticles(histograms, + clusterHandle, + clusters, + caloParticleHandle, + cP, + cPIndices, + cPSelectedIndices, + hitMap, + layers, + LCAssocByEnergyScoreHandle); //To find out the total amount of energy clustered per layer //Initialize with zeros because I see clear gives weird numbers.