From f3afa7fb8dddfef5f730377ad3c9d074cb3bb44f Mon Sep 17 00:00:00 2001 From: Andreas Psallidas Date: Wed, 10 Jul 2019 12:41:36 +0200 Subject: [PATCH 1/7] correcting the mess that merge produced --- .../src/HGVHistoProducerAlgo.cc | 48 ------------------- 1 file changed, 48 deletions(-) diff --git a/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc b/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc index 42e775b273b03..62ac311843687 100644 --- a/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc +++ b/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc @@ -2156,54 +2156,6 @@ void HGVHistoProducerAlgo::multiClusters_to_CaloParticles(const Histograms& hist } } //end of loop through sim hits of current calo particle - 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) << "mclWithMaxEnergyInCP\t" << std::setw(15) - << "maxEnergyMCLinCP\t" << std::setw(20) << "CPEnergyFractionInMCL" - << "\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) << mclWithMaxEnergyInCP << "\t" - << std::setw(15) << maxEnergyMCLperlayerinCP << "\t" << std::setw(20) - << CPEnergyFractionInMCLperlayer << "\n"; - - 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 hitWithNoMCL = false; - if (cpFraction == 0.f) - continue; //hopefully this should never happen - auto hit_find_in_MCL = detIdToMultiClusterId_Map.find(cp_hitDetId); - if (hit_find_in_MCL == detIdToMultiClusterId_Map.end()) - hitWithNoMCL = 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 multiClusterId = lcPair.first; - float mclFraction = 0.f; - - if (!hitWithNoMCL) { - auto findHitIt = std::find(detIdToMultiClusterId_Map[cp_hitDetId].begin(), - detIdToMultiClusterId_Map[cp_hitDetId].end(), - HGVHistoProducerAlgo::detIdInfoInMultiCluster{multiClusterId, 0, 0.f}); - if (findHitIt != detIdToMultiClusterId_Map[cp_hitDetId].end()) - mclFraction = findHitIt->fraction; - } - if (mclFraction == 0.) { - mclFraction = -1.; - } - //Observe here that we do not divide as before by the layer cluster energy weight. We should sum first - //over all layers and divide with the total CP energy over all layers. - lcPair.second.second += (mclFraction - cpFraction) * (mclFraction - cpFraction) * hitEnergyWeight; - LogDebug("HGCalValidator") << "multiClusterId:\t" << multiClusterId << "\t" - << "mclfraction,cpfraction:\t" << mclFraction << ", " << cpFraction << "\t" - << "hitEnergyWeight:\t" << hitEnergyWeight << "\t" - << "currect score numerator:\t" << lcPair.second.second << "\n"; - } - } //end of loop through sim hits of current calo particle - if (cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore.empty()) LogDebug("HGCalValidator") << "CP Id: \t" << cpId << "\t MCL id:\t-1 " << "\t layer \t " << layerId << " Sub score in \t -1" From 058b709f9f17c4001665e2601d473e0593fff55c Mon Sep 17 00:00:00 2001 From: Andreas Psallidas Date: Wed, 10 Jul 2019 12:46:46 +0200 Subject: [PATCH 2/7] CaloParticles not from PU --- Validation/HGCalValidation/plugins/HGCalValidator.cc | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/Validation/HGCalValidation/plugins/HGCalValidator.cc b/Validation/HGCalValidation/plugins/HGCalValidator.cc index 3174fd4fb9e8f..f17310f04f979 100644 --- a/Validation/HGCalValidation/plugins/HGCalValidator.cc +++ b/Validation/HGCalValidation/plugins/HGCalValidator.cc @@ -192,6 +192,18 @@ void HGCalValidator::dqmAnalyze(const edm::Event& event, histoProducerAlgo_->fill_info_histos(histograms.histoProducerAlgo, totallayers_to_monitor_); } + //Consider CaloParticles coming from the hard scatterer, excluding the PU contribution. + std::vector caloParticlesFromHardScat; + for ( auto& it_caloPart : caloParticles) { + if (it_caloPart.g4Tracks()[0].eventId().event() != 0 or it_caloPart.g4Tracks()[0].eventId().bunchCrossing() != 0) { + LogDebug("HGCalValidator") << "Excluding CaloParticles from event: " + << it_caloPart.g4Tracks()[0].eventId().event() + << " with BX: " << it_caloPart.g4Tracks()[0].eventId().bunchCrossing() << std::endl; + continue; + } + caloParticlesFromHardScat.push_back(it_caloPart); + } + // ############################################## // fill caloparticles histograms // ############################################## From 61209fe2eb20f3514790c426220713edc18c0e68 Mon Sep 17 00:00:00 2001 From: Andreas Psallidas Date: Wed, 10 Jul 2019 18:44:04 +0200 Subject: [PATCH 3/7] caloparticles from hard scatterer --- Validation/HGCalValidation/plugins/HGCalValidator.cc | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/Validation/HGCalValidation/plugins/HGCalValidator.cc b/Validation/HGCalValidation/plugins/HGCalValidator.cc index f17310f04f979..149a8f2097c0a 100644 --- a/Validation/HGCalValidation/plugins/HGCalValidator.cc +++ b/Validation/HGCalValidation/plugins/HGCalValidator.cc @@ -236,7 +236,8 @@ void HGCalValidator::dqmAnalyze(const edm::Event& event, w, clusters, densities, - caloParticles, + // caloParticles, + caloParticlesFromHardScat, hitMap, cummatbudg, totallayers_to_monitor_, @@ -250,7 +251,10 @@ void HGCalValidator::dqmAnalyze(const edm::Event& event, if (domulticlustersPlots_) { w++; histoProducerAlgo_->fill_multi_cluster_histos( - histograms.histoProducerAlgo, w, multiClusters, caloParticles, hitMap, totallayers_to_monitor_); + histograms.histoProducerAlgo, w, multiClusters, + // caloParticles, + caloParticlesFromHardScat, + hitMap, totallayers_to_monitor_); } //General Info From fbb4c3c86dda9272c844ec8b84176c57d7fb1213 Mon Sep 17 00:00:00 2001 From: Andreas Psallidas Date: Wed, 10 Jul 2019 19:05:53 +0200 Subject: [PATCH 4/7] code format --- .../HGCalValidation/plugins/HGCalValidator.cc | 22 ++++++++++--------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/Validation/HGCalValidation/plugins/HGCalValidator.cc b/Validation/HGCalValidation/plugins/HGCalValidator.cc index 149a8f2097c0a..9be942d5ae65c 100644 --- a/Validation/HGCalValidation/plugins/HGCalValidator.cc +++ b/Validation/HGCalValidation/plugins/HGCalValidator.cc @@ -193,12 +193,12 @@ void HGCalValidator::dqmAnalyze(const edm::Event& event, } //Consider CaloParticles coming from the hard scatterer, excluding the PU contribution. - std::vector caloParticlesFromHardScat; - for ( auto& it_caloPart : caloParticles) { + std::vector caloParticlesFromHardScat; + for (auto& it_caloPart : caloParticles) { if (it_caloPart.g4Tracks()[0].eventId().event() != 0 or it_caloPart.g4Tracks()[0].eventId().bunchCrossing() != 0) { LogDebug("HGCalValidator") << "Excluding CaloParticles from event: " - << it_caloPart.g4Tracks()[0].eventId().event() - << " with BX: " << it_caloPart.g4Tracks()[0].eventId().bunchCrossing() << std::endl; + << it_caloPart.g4Tracks()[0].eventId().event() + << " with BX: " << it_caloPart.g4Tracks()[0].eventId().bunchCrossing() << std::endl; continue; } caloParticlesFromHardScat.push_back(it_caloPart); @@ -237,7 +237,7 @@ void HGCalValidator::dqmAnalyze(const edm::Event& event, clusters, densities, // caloParticles, - caloParticlesFromHardScat, + caloParticlesFromHardScat, hitMap, cummatbudg, totallayers_to_monitor_, @@ -250,11 +250,13 @@ void HGCalValidator::dqmAnalyze(const edm::Event& event, if (domulticlustersPlots_) { w++; - histoProducerAlgo_->fill_multi_cluster_histos( - histograms.histoProducerAlgo, w, multiClusters, - // caloParticles, - caloParticlesFromHardScat, - hitMap, totallayers_to_monitor_); + histoProducerAlgo_->fill_multi_cluster_histos(histograms.histoProducerAlgo, + w, + multiClusters, + // caloParticles, + caloParticlesFromHardScat, + hitMap, + totallayers_to_monitor_); } //General Info From 7e6ca2d8335f89ac1adebf85a8274c4273ab1b73 Mon Sep 17 00:00:00 2001 From: Andreas Psallidas Date: Thu, 11 Jul 2019 15:19:47 +0200 Subject: [PATCH 5/7] using emplace_back --- Validation/HGCalValidation/plugins/HGCalValidator.cc | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Validation/HGCalValidation/plugins/HGCalValidator.cc b/Validation/HGCalValidation/plugins/HGCalValidator.cc index 9be942d5ae65c..6e2fafde73460 100644 --- a/Validation/HGCalValidation/plugins/HGCalValidator.cc +++ b/Validation/HGCalValidation/plugins/HGCalValidator.cc @@ -193,15 +193,16 @@ void HGCalValidator::dqmAnalyze(const edm::Event& event, } //Consider CaloParticles coming from the hard scatterer, excluding the PU contribution. + // std::vector caloParticlesFromHardScat; std::vector caloParticlesFromHardScat; - for (auto& it_caloPart : caloParticles) { + for (const auto& it_caloPart : caloParticles) { if (it_caloPart.g4Tracks()[0].eventId().event() != 0 or it_caloPart.g4Tracks()[0].eventId().bunchCrossing() != 0) { LogDebug("HGCalValidator") << "Excluding CaloParticles from event: " << it_caloPart.g4Tracks()[0].eventId().event() << " with BX: " << it_caloPart.g4Tracks()[0].eventId().bunchCrossing() << std::endl; continue; } - caloParticlesFromHardScat.push_back(it_caloPart); + caloParticlesFromHardScat.emplace_back(it_caloPart); } // ############################################## From 32bc59523adf8ef6e6048715b1a730668cad4974 Mon Sep 17 00:00:00 2001 From: Andreas Psallidas Date: Fri, 12 Jul 2019 14:42:39 +0200 Subject: [PATCH 6/7] going with the indices --- .../interface/HGVHistoProducerAlgo.h | 8 +++-- .../HGCalValidation/plugins/HGCalValidator.cc | 32 +++++++++---------- .../src/HGVHistoProducerAlgo.cc | 32 +++++++++++-------- 3 files changed, 40 insertions(+), 32 deletions(-) diff --git a/Validation/HGCalValidation/interface/HGVHistoProducerAlgo.h b/Validation/HGCalValidation/interface/HGVHistoProducerAlgo.h index b6c647dc3dfc9..2a86d884248e1 100644 --- a/Validation/HGCalValidation/interface/HGVHistoProducerAlgo.h +++ b/Validation/HGCalValidation/interface/HGVHistoProducerAlgo.h @@ -174,11 +174,13 @@ class HGVHistoProducerAlgo { void layerClusters_to_CaloParticles(const Histograms& histograms, const reco::CaloClusterCollection& clusters, std::vector const& cP, + std::vector const& cPIndices, std::map const&, unsigned layers) const; void multiClusters_to_CaloParticles(const Histograms& histograms, const std::vector& multiClusters, std::vector const& cP, + std::vector const& cPIndices, std::map const&, unsigned layers, std::vector contimulti) const; @@ -193,6 +195,7 @@ class HGVHistoProducerAlgo { const reco::CaloClusterCollection& clusters, const Density& densities, std::vector const& cP, + std::vector const& cPIndices, std::map const&, std::map cummatbudg, unsigned layers, @@ -201,6 +204,7 @@ class HGVHistoProducerAlgo { int count, const std::vector& multiClusters, std::vector const& cP, + std::vector const& cPIndices, std::map const&, unsigned layers) const; double distance2(const double x1, const double y1, const double x2, const double y2) const; @@ -212,14 +216,14 @@ class HGVHistoProducerAlgo { struct detIdInfoInCluster { bool operator==(const detIdInfoInCluster& o) const { return clusterId == o.clusterId; }; - unsigned int clusterId; + long unsigned int clusterId; float fraction; }; struct detIdInfoInMultiCluster { bool operator==(const detIdInfoInMultiCluster& o) const { return multiclusterId == o.multiclusterId; }; unsigned int multiclusterId; - unsigned int clusterId; + long unsigned int clusterId; float fraction; }; diff --git a/Validation/HGCalValidation/plugins/HGCalValidator.cc b/Validation/HGCalValidation/plugins/HGCalValidator.cc index 6e2fafde73460..eea9bd53f821c 100644 --- a/Validation/HGCalValidation/plugins/HGCalValidator.cc +++ b/Validation/HGCalValidation/plugins/HGCalValidator.cc @@ -192,17 +192,20 @@ void HGCalValidator::dqmAnalyze(const edm::Event& event, histoProducerAlgo_->fill_info_histos(histograms.histoProducerAlgo, totallayers_to_monitor_); } - //Consider CaloParticles coming from the hard scatterer, excluding the PU contribution. - // std::vector caloParticlesFromHardScat; - std::vector caloParticlesFromHardScat; - for (const auto& it_caloPart : caloParticles) { - if (it_caloPart.g4Tracks()[0].eventId().event() != 0 or it_caloPart.g4Tracks()[0].eventId().bunchCrossing() != 0) { + 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 (caloParticles[cpId].g4Tracks()[0].eventId().event() != 0 or + caloParticles[cpId].g4Tracks()[0].eventId().bunchCrossing() != 0) { LogDebug("HGCalValidator") << "Excluding CaloParticles from event: " - << it_caloPart.g4Tracks()[0].eventId().event() - << " with BX: " << it_caloPart.g4Tracks()[0].eventId().bunchCrossing() << std::endl; + << caloParticles[cpId].g4Tracks()[0].eventId().event() + << " with BX: " << caloParticles[cpId].g4Tracks()[0].eventId().bunchCrossing() + << std::endl; continue; } - caloParticlesFromHardScat.emplace_back(it_caloPart); + cPIndices.emplace_back(cpId); } // ############################################## @@ -237,8 +240,8 @@ void HGCalValidator::dqmAnalyze(const edm::Event& event, w, clusters, densities, - // caloParticles, - caloParticlesFromHardScat, + caloParticles, + cPIndices, hitMap, cummatbudg, totallayers_to_monitor_, @@ -251,13 +254,8 @@ void HGCalValidator::dqmAnalyze(const edm::Event& event, if (domulticlustersPlots_) { w++; - histoProducerAlgo_->fill_multi_cluster_histos(histograms.histoProducerAlgo, - w, - multiClusters, - // caloParticles, - caloParticlesFromHardScat, - hitMap, - totallayers_to_monitor_); + histoProducerAlgo_->fill_multi_cluster_histos( + histograms.histoProducerAlgo, w, multiClusters, caloParticles, cPIndices, hitMap, totallayers_to_monitor_); } //General Info diff --git a/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc b/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc index 62ac311843687..fab86cdbf628c 100644 --- a/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc +++ b/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc @@ -868,10 +868,12 @@ void HGVHistoProducerAlgo::fill_cluster_histos(const Histograms& histograms, void HGVHistoProducerAlgo::layerClusters_to_CaloParticles(const Histograms& histograms, const reco::CaloClusterCollection& clusters, std::vector const& cP, + std::vector const& cPIndices, std::map const& hitMap, unsigned layers) const { auto nLayerClusters = clusters.size(); - auto nCaloParticles = cP.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; @@ -891,7 +893,7 @@ void HGVHistoProducerAlgo::layerClusters_to_CaloParticles(const Histograms& hist } } - for (unsigned int cpId = 0; cpId < nCaloParticles; ++cpId) { + for (const auto& cpId : cPIndices) { const SimClusterRefVector& simClusterRefVector = cP[cpId].simClusters(); for (const auto& it_sc : simClusterRefVector) { const SimCluster& simCluster = (*(it_sc)); @@ -1220,7 +1222,7 @@ void HGVHistoProducerAlgo::layerClusters_to_CaloParticles(const Histograms& hist histograms.h_denom_layercl_phi_perlayer.at(lcLayerId).fill(clusters[lcId].phi()); } // End of loop over LayerClusters - for (unsigned int cpId = 0; cpId < nCaloParticles; ++cpId) { + for (const auto& cpId : cPIndices) { 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; @@ -1335,6 +1337,7 @@ void HGVHistoProducerAlgo::fill_generic_cluster_histos(const Histograms& histogr const reco::CaloClusterCollection& clusters, const Density& densities, std::vector const& cP, + std::vector const& cPIndices, std::map const& hitMap, std::map cummatbudg, unsigned layers, @@ -1363,7 +1366,7 @@ 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, hitMap, layers); + layerClusters_to_CaloParticles(histograms, clusters, cP, cPIndices, hitMap, layers); //To find out the total amount of energy clustered per layer //Initialize with zeros because I see clear gives weird numbers. @@ -1374,12 +1377,12 @@ void HGVHistoProducerAlgo::fill_generic_cluster_histos(const Histograms& histogr //We need to compare with the total amount of energy coming from caloparticles double caloparteneplus = 0.; double caloparteneminus = 0.; - for (auto const caloParticle : cP) { - if (caloParticle.eta() >= 0.) { - caloparteneplus = caloparteneplus + caloParticle.energy(); + for (const auto& cpId : cPIndices) { + if (cP[cpId].eta() >= 0.) { + caloparteneplus = caloparteneplus + cP[cpId].energy(); } - if (caloParticle.eta() < 0.) { - caloparteneminus = caloparteneminus + caloParticle.energy(); + if (cP[cpId].eta() < 0.) { + caloparteneminus = caloparteneminus + cP[cpId].energy(); } } @@ -1644,11 +1647,13 @@ void HGVHistoProducerAlgo::fill_generic_cluster_histos(const Histograms& histogr void HGVHistoProducerAlgo::multiClusters_to_CaloParticles(const Histograms& histograms, const std::vector& multiClusters, std::vector const& cP, + std::vector const& cPIndices, std::map const& hitMap, unsigned layers, std::vector contimulti) const { auto nMultiClusters = multiClusters.size(); - auto nCaloParticles = cP.size(); + //Consider CaloParticles coming from the hard scatterer, excluding the PU contribution. + auto nCaloParticles = cPIndices.size(); std::unordered_map> detIdToCaloParticleId_Map; std::unordered_map> detIdToMultiClusterId_Map; @@ -1677,7 +1682,7 @@ void HGVHistoProducerAlgo::multiClusters_to_CaloParticles(const Histograms& hist } } - for (unsigned int cpId = 0; cpId < nCaloParticles; ++cpId) { + for (const auto& cpId : cPIndices) { //take sim clusters const SimClusterRefVector& simClusterRefVector = cP[cpId].simClusters(); //loop through sim clusters @@ -2074,7 +2079,7 @@ void HGVHistoProducerAlgo::multiClusters_to_CaloParticles(const Histograms& hist mclsharedenergyfrac.resize(nCaloParticles); //Loop though caloparticles - for (unsigned int cpId = 0; cpId < nCaloParticles; ++cpId) { + for (const auto& cpId : cPIndices) { for (unsigned int i = 0; i < nCaloParticles; ++i) { score3d[i].resize(nMultiClusters); mclsharedenergy[i].resize(nMultiClusters); @@ -2229,6 +2234,7 @@ void HGVHistoProducerAlgo::fill_multi_cluster_histos(const Histograms& histogram int count, const std::vector& multiClusters, std::vector const& cP, + std::vector const& cPIndices, std::map const& hitMap, unsigned layers) const { //Each event to be treated as two events: @@ -2411,7 +2417,7 @@ void HGVHistoProducerAlgo::fill_multi_cluster_histos(const Histograms& histogram histograms.h_noncontmulticlusternum.fill(tnnoncontmclmz); } - multiClusters_to_CaloParticles(histograms, multiClusters, cP, hitMap, layers, contmulti); + multiClusters_to_CaloParticles(histograms, multiClusters, cP, cPIndices, hitMap, layers, contmulti); } double HGVHistoProducerAlgo::distance2(const double x1, From f6755e4d8c8d6fa4ba3eef2f1697aa5379ef1f8b Mon Sep 17 00:00:00 2001 From: Andreas Psallidas Date: Sat, 13 Jul 2019 11:40:59 +0200 Subject: [PATCH 7/7] moving initialization outside loop, restrict iteration on multiclusters --- .../src/HGVHistoProducerAlgo.cc | 40 +++++++++++-------- 1 file changed, 23 insertions(+), 17 deletions(-) diff --git a/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc b/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc index fab86cdbf628c..96043d576b260 100644 --- a/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc +++ b/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc @@ -2078,21 +2078,25 @@ void HGVHistoProducerAlgo::multiClusters_to_CaloParticles(const Histograms& hist std::vector> mclsharedenergyfrac; mclsharedenergyfrac.resize(nCaloParticles); + for (unsigned int i = 0; i < nCaloParticles; ++i) { + score3d[i].resize(nMultiClusters); + mclsharedenergy[i].resize(nMultiClusters); + mclsharedenergyfrac[i].resize(nMultiClusters); + for (unsigned int j = 0; j < nMultiClusters; ++j) { + score3d[i][j] = 0.f; + mclsharedenergy[i][j] = 0.f; + mclsharedenergyfrac[i][j] = 0.f; + } + } + //Loop though caloparticles for (const auto& cpId : cPIndices) { - for (unsigned int i = 0; i < nCaloParticles; ++i) { - score3d[i].resize(nMultiClusters); - mclsharedenergy[i].resize(nMultiClusters); - mclsharedenergyfrac[i].resize(nMultiClusters); - for (unsigned int j = 0; j < nMultiClusters; ++j) { - score3d[i][j] = 0.f; - mclsharedenergy[i][j] = 0.f; - mclsharedenergyfrac[i][j] = 0.f; - } - } + //We need to keep the multiclusters ids that are related to + //CaloParticle under study for the final filling of the score. + std::vector cpId_mclId_related; + cpId_mclId_related.clear(); float CPenergy = 0.f; - for (unsigned int layerId = 0; layerId < layers * 2; ++layerId) { unsigned int CPNumberOfHits = cPOnLayer[cpId][layerId].hits_and_fractions.size(); //Below gives the CP energy related to multicluster per layer. @@ -2139,6 +2143,7 @@ void HGVHistoProducerAlgo::multiClusters_to_CaloParticles(const Histograms& hist float hitEnergyWeight = hit->energy() * hit->energy(); for (auto& lcPair : cPOnLayer[cpId][layerId].layerClusterIdToEnergyAndScore) { unsigned int multiClusterId = lcPair.first; + cpId_mclId_related.emplace_back(multiClusterId); float mclFraction = 0.f; if (!hitWithNoMCL) { @@ -2148,9 +2153,9 @@ void HGVHistoProducerAlgo::multiClusters_to_CaloParticles(const Histograms& hist if (findHitIt != detIdToMultiClusterId_Map[cp_hitDetId].end()) mclFraction = findHitIt->fraction; } - if (mclFraction == 0.) { - mclFraction = -1.; - } + // if (mclFraction == 0.) { + // mclFraction = -1.; + // } //Observe here that we do not divide as before by the layer cluster energy weight. We should sum first //over all layers and divide with the total CP energy over all layers. lcPair.second.second += (mclFraction - cpFraction) * (mclFraction - cpFraction) * hitEnergyWeight; @@ -2175,13 +2180,14 @@ void HGVHistoProducerAlgo::multiClusters_to_CaloParticles(const Histograms& hist float invCPEnergyWeight = 1.f / (CPenergy * CPenergy); - //Loop through multiclusters here - for (unsigned int mclId = 0; mclId < nMultiClusters; ++mclId) { + //Loop through related multiclusters here + for (unsigned int i = 0; i < cpId_mclId_related.size(); ++i) { + unsigned int mclId = cpId_mclId_related[i]; //Now time for the denominator score3d[cpId][mclId] = score3d[cpId][mclId] * invCPEnergyWeight; mclsharedenergyfrac[cpId][mclId] = (mclsharedenergy[cpId][mclId] / CPenergy); - LogDebug("HGCalValidator") << "CP Id: \t" << cpId << "\t MCL id: \t" << mclId << "\t score \t" + LogDebug("HGCalValidator") << "CP Id: \t" << cpId << "\t MCL id: \t" << mclId << "\t score \t" // << score3d[cpId][mclId] << "\t" << "invCPEnergyWeight \t" << invCPEnergyWeight << "\t" << "shared energy:\t" << mclsharedenergy[cpId][mclId] << "\t"