diff --git a/RecoHGCal/TICL/plugins/HGCDoublet.cc b/RecoHGCal/TICL/plugins/HGCDoublet.cc index cac2c4d1a52df..cec1229a6b8f0 100644 --- a/RecoHGCal/TICL/plugins/HGCDoublet.cc +++ b/RecoHGCal/TICL/plugins/HGCDoublet.cc @@ -110,9 +110,10 @@ int HGCDoublet::areAligned(double xi, auto minCosPointing_sq = minCosPointing * minCosPointing; bool isWithinLimitsPointing = (dot_pointing_sq > minCosPointing_sq * mag_pointing_sq * mag2sq); if (debug) { - LogDebug("HGCDoublet") << "-- Are Aligned -- dot_pointing_sq: " << dot_pointing * dot_pointing + LogDebug("HGCDoublet") << "Pointing direction: " << pointingDir << std::endl; + LogDebug("HGCDoublet") << "-- Are Aligned -- dot_pointing_sq: " << dot_pointing_sq << " mag_pointing_sq: " << mag_pointing_sq << " mag2sq: " << mag2sq - << " isWithinLimits: " << isWithinLimitsPointing << std::endl; + << " isWithinLimitsPointing: " << isWithinLimitsPointing << std::endl; } // by squaring cosTheta and multiplying by the squares of the magnitudes // an equivalent comparison is made without the division and square root which are costly FP ops. diff --git a/Validation/Configuration/python/hgcalSimValid_cff.py b/Validation/Configuration/python/hgcalSimValid_cff.py index 2fc76100faa36..b0c5a299ccf7c 100644 --- a/Validation/Configuration/python/hgcalSimValid_cff.py +++ b/Validation/Configuration/python/hgcalSimValid_cff.py @@ -17,6 +17,9 @@ from Validation.HGCalValidation.ticlPFValidation_cfi import ticlPFValidation hgcalTiclPFValidation = cms.Sequence(ticlPFValidation) +from Validation.HGCalValidation.ticlTrackstersEdgesValidation_cfi import ticlTrackstersEdgesValidation +hgcalTiclTrackstersEdgesValidationSequence = cms.Sequence(ticlTrackstersEdgesValidation) + hgcalValidatorSequence = cms.Sequence(hgcalValidator) hgcalPFJetValidation = _hgcalPFJetValidation.clone(BenchmarkLabel = 'PFJetValidation/HGCAlCompWithGenJet', VariablePtBins=[10., 30., 80., 120., 250., 600.], @@ -38,6 +41,8 @@ + hgcalHitValidationSequence + hgcalValidatorSequence + hgcalTiclPFValidation + #Currently commented out until trackster edges are saved +# + hgcalTiclTrackstersEdgesValidationSequence + hgcalPFJetValidation) _hfnose_hgcalAssociatorsTask = hgcalAssociators.copy() diff --git a/Validation/HGCalValidation/interface/HGCalValidator.h b/Validation/HGCalValidation/interface/HGCalValidator.h index 2e603edcc6a98..f617830125c1d 100644 --- a/Validation/HGCalValidation/interface/HGCalValidator.h +++ b/Validation/HGCalValidation/interface/HGCalValidator.h @@ -60,6 +60,7 @@ class HGCalValidator : public DQMGlobalEDAnalyzer { std::unordered_map const&) const; protected: + edm::ESGetToken caloGeomToken_; edm::InputTag label_lcl; std::vector label_mcl; edm::InputTag associator_; @@ -67,9 +68,9 @@ class HGCalValidator : public DQMGlobalEDAnalyzer { const bool SaveGeneralInfo_; const bool doCaloParticlePlots_; const bool doCaloParticleSelection_; - const bool dosimclustersPlots_; - const bool dolayerclustersPlots_; - const bool domulticlustersPlots_; + const bool doSimClustersPlots_; + const bool doLayerClustersPlots_; + const bool doMultiClustersPlots_; std::vector label_clustersmask; const edm::FileInPath cummatbudinxo_; diff --git a/Validation/HGCalValidation/interface/HGVHistoProducerAlgo.h b/Validation/HGCalValidation/interface/HGVHistoProducerAlgo.h index a25a1f0d4eb64..f6359fd9294b1 100644 --- a/Validation/HGCalValidation/interface/HGVHistoProducerAlgo.h +++ b/Validation/HGCalValidation/interface/HGVHistoProducerAlgo.h @@ -149,6 +149,7 @@ struct HGVHistoProducerAlgoHistograms { std::vector h_numMerge_multicl_phi; std::vector h_sharedenergy_multicl2caloparticle; std::vector h_sharedenergy_caloparticle2multicl; + std::vector h_sharedenergy_caloparticle2multicl_assoc; std::vector h_sharedenergy_multicl2caloparticle_vs_eta; std::vector h_sharedenergy_multicl2caloparticle_vs_phi; std::vector h_sharedenergy_caloparticle2multicl_vs_eta; diff --git a/Validation/HGCalValidation/plugins/HGCalValidator.cc b/Validation/HGCalValidation/plugins/HGCalValidator.cc index 3bc004f7df673..e662bd1ff7b32 100644 --- a/Validation/HGCalValidation/plugins/HGCalValidator.cc +++ b/Validation/HGCalValidation/plugins/HGCalValidator.cc @@ -9,16 +9,17 @@ using namespace std; using namespace edm; HGCalValidator::HGCalValidator(const edm::ParameterSet& pset) - : label_lcl(pset.getParameter("label_lcl")), + : caloGeomToken_(esConsumes()), + label_lcl(pset.getParameter("label_lcl")), label_mcl(pset.getParameter>("label_mcl")), associator_(pset.getUntrackedParameter("associator")), associatorSim_(pset.getUntrackedParameter("associatorSim")), SaveGeneralInfo_(pset.getUntrackedParameter("SaveGeneralInfo")), doCaloParticlePlots_(pset.getUntrackedParameter("doCaloParticlePlots")), doCaloParticleSelection_(pset.getUntrackedParameter("doCaloParticleSelection")), - dosimclustersPlots_(pset.getUntrackedParameter("dosimclustersPlots")), - dolayerclustersPlots_(pset.getUntrackedParameter("dolayerclustersPlots")), - domulticlustersPlots_(pset.getUntrackedParameter("domulticlustersPlots")), + doSimClustersPlots_(pset.getUntrackedParameter("doSimClustersPlots")), + doLayerClustersPlots_(pset.getUntrackedParameter("doLayerClustersPlots")), + doMultiClustersPlots_(pset.getUntrackedParameter("doMultiClustersPlots")), label_clustersmask(pset.getParameter>("LayerClustersInputMask")), cummatbudinxo_(pset.getParameter("cummatbudinxo")) { //In this way we can easily generalize to associations between other objects also. @@ -115,7 +116,7 @@ void HGCalValidator::bookHistograms(DQMStore::IBooker& ibook, } //Booking histograms concerning with simclusters - if (dosimclustersPlots_) { + if (doSimClustersPlots_) { ibook.cd(); ibook.setCurrentFolder(dirName_ + "simClusters/ClusterLevel"); histoProducerAlgo_->bookSimClusterHistos( @@ -149,7 +150,7 @@ void HGCalValidator::bookHistograms(DQMStore::IBooker& ibook, } //if for simcluster plots //Booking histograms concerning with hgcal layer clusters - if (dolayerclustersPlots_) { + if (doLayerClustersPlots_) { ibook.cd(); ibook.setCurrentFolder(dirName_ + "hgcalLayerClusters/ClusterLevel"); histoProducerAlgo_->bookClusterHistos_ClusterLevel(ibook, @@ -192,7 +193,7 @@ void HGCalValidator::bookHistograms(DQMStore::IBooker& ibook, ibook.setCurrentFolder(dirName); //Booking histograms concerning for hgcal multi clusters - if (domulticlustersPlots_) { + if (doMultiClustersPlots_) { histoProducerAlgo_->bookMultiClusterHistos(ibook, histograms.histoProducerAlgo, totallayers_to_monitor_); } } //end of booking multiclusters loop @@ -241,8 +242,7 @@ void HGCalValidator::dqmAnalyze(const edm::Event& event, event.getByToken(label_cp_effic, caloParticleHandle); std::vector const& caloParticles = *caloParticleHandle; - edm::ESHandle geom; - setup.get().get(geom); + edm::ESHandle geom = setup.getHandle(caloGeomToken_); tools_->setGeometry(*geom); histoProducerAlgo_->setRecHitTools(tools_); @@ -310,7 +310,7 @@ void HGCalValidator::dqmAnalyze(const edm::Event& event, // ############################################## // fill simcluster histograms // ############################################## - if (dosimclustersPlots_) { + if (doSimClustersPlots_) { histoProducerAlgo_->fill_simcluster_histos( histograms.histoProducerAlgo, simclusters, totallayers_to_monitor_, thicknesses_to_monitor_); @@ -347,7 +347,7 @@ void HGCalValidator::dqmAnalyze(const edm::Event& event, // fill layercluster histograms // ############################################## int w = 0; //counter counting the number of sets of histograms - if (dolayerclustersPlots_) { + if (doLayerClustersPlots_) { histoProducerAlgo_->fill_generic_cluster_histos(histograms.histoProducerAlgo, w, clusterHandle, @@ -377,7 +377,7 @@ void HGCalValidator::dqmAnalyze(const edm::Event& event, // fill multicluster histograms // ############################################## for (unsigned int wml = 0; wml < label_mclTokens.size(); wml++) { - if (domulticlustersPlots_) { + if (doMultiClustersPlots_) { edm::Handle> multiClusterHandle; event.getByToken(label_mclTokens[wml], multiClusterHandle); const std::vector& multiClusters = *multiClusterHandle; diff --git a/Validation/HGCalValidation/plugins/TICLTrackstersEdgesValidation.cc b/Validation/HGCalValidation/plugins/TICLTrackstersEdgesValidation.cc new file mode 100644 index 0000000000000..792808e343161 --- /dev/null +++ b/Validation/HGCalValidation/plugins/TICLTrackstersEdgesValidation.cc @@ -0,0 +1,329 @@ +#include +#include +#include + +// user include files +#include "DataFormats/Math/interface/Point3D.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "DQMServices/Core/interface/DQMGlobalEDAnalyzer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "FWCore/Utilities/interface/transform.h" +#include "DataFormats/HGCalReco/interface/Trackster.h" +#include "DataFormats/HGCalReco/interface/TICLSeedingRegion.h" + +#include "DataFormats/CaloRecHit/interface/CaloCluster.h" +#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h" + +using namespace ticl; + +struct Histogram_TICLTrackstersEdgesValidation { + dqm::reco::MonitorElement* number_; + dqm::reco::MonitorElement *raw_energy_, *raw_energy_1plusLC_; + dqm::reco::MonitorElement *regr_energy_, *regr_energy_1plusLC_; + dqm::reco::MonitorElement *raw_energy_vs_regr_energy_, *raw_energy_vs_regr_energy_1plusLC_; + dqm::reco::MonitorElement *id_prob_, *id_prob_1plusLC_; + dqm::reco::MonitorElement* delta_energy_; + dqm::reco::MonitorElement* delta_energy_relative_; + dqm::reco::MonitorElement* delta_energy_vs_energy_; + dqm::reco::MonitorElement* delta_energy_vs_layer_; + dqm::reco::MonitorElement* delta_energy_relative_vs_layer_; + dqm::reco::MonitorElement* delta_layer_; + dqm::reco::MonitorElement* ingoing_links_vs_layer_; + dqm::reco::MonitorElement* outgoing_links_vs_layer_; + // For the definition of the angles, read http://hgcal.web.cern.ch/hgcal/Reconstruction/Tutorial/ + dqm::reco::MonitorElement* angle_alpha_; + dqm::reco::MonitorElement* angle_alpha_alternative_; + dqm::reco::MonitorElement* angle_beta_; + std::vector angle_beta_byLayer_; + std::vector angle_beta_w_byLayer_; +}; + +using Histograms_TICLTrackstersEdgesValidation = + std::unordered_map; + +class TICLTrackstersEdgesValidation : public DQMGlobalEDAnalyzer { +public: + explicit TICLTrackstersEdgesValidation(const edm::ParameterSet&); + ~TICLTrackstersEdgesValidation() override; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + +private: + void bookHistograms(DQMStore::IBooker&, + edm::Run const&, + edm::EventSetup const&, + Histograms_TICLTrackstersEdgesValidation&) const override; + + void dqmAnalyze(edm::Event const&, + edm::EventSetup const&, + Histograms_TICLTrackstersEdgesValidation const&) const override; + void dqmBeginRun(edm::Run const&, edm::EventSetup const&, Histograms_TICLTrackstersEdgesValidation&) const override; + + edm::ESGetToken caloGeomToken_; + std::string folder_; + std::vector trackstersCollectionsNames_; + std::vector>> tracksterTokens_; + edm::EDGetTokenT> layerClustersToken_; + edm::EDGetTokenT> ticlSeedingGlobalToken_; + edm::EDGetTokenT> ticlSeedingTrkToken_; + mutable hgcal::RecHitTools rhtools_; +}; + +TICLTrackstersEdgesValidation::TICLTrackstersEdgesValidation(const edm::ParameterSet& iConfig) + : caloGeomToken_(esConsumes()), + folder_(iConfig.getParameter("folder")) { + tracksterTokens_ = edm::vector_transform(iConfig.getParameter>("tracksterCollections"), + [this](edm::InputTag const& tag) { + trackstersCollectionsNames_.emplace_back(tag.label()); + return consumes>(tag); + }); + layerClustersToken_ = consumes>(iConfig.getParameter("layerClusters")); + ticlSeedingGlobalToken_ = + consumes>(iConfig.getParameter("ticlSeedingGlobal")); + ticlSeedingTrkToken_ = + consumes>(iConfig.getParameter("ticlSeedingTrk")); +} + +TICLTrackstersEdgesValidation::~TICLTrackstersEdgesValidation() {} + +void TICLTrackstersEdgesValidation::dqmAnalyze(edm::Event const& iEvent, + edm::EventSetup const& iSetup, + Histograms_TICLTrackstersEdgesValidation const& histos) const { + edm::Handle> layerClustersH; + iEvent.getByToken(layerClustersToken_, layerClustersH); + auto const& layerClusters = *layerClustersH.product(); + + edm::Handle> ticlSeedingGlobalH; + iEvent.getByToken(ticlSeedingGlobalToken_, ticlSeedingGlobalH); + auto const& ticlSeedingGlobal = *ticlSeedingGlobalH.product(); + + edm::Handle> ticlSeedingTrkH; + iEvent.getByToken(ticlSeedingTrkToken_, ticlSeedingTrkH); + auto const& ticlSeedingTrk = *ticlSeedingTrkH.product(); + + for (const auto& trackster_token : tracksterTokens_) { + edm::Handle> trackster_h; + iEvent.getByToken(trackster_token, trackster_h); + auto numberOfTracksters = trackster_h->size(); + //using .at() as [] is not const + const auto& histo = histos.at(trackster_token.index()); + histo.number_->Fill(numberOfTracksters); + for (unsigned int i = 0; i < numberOfTracksters; ++i) { + const auto& thisTrackster = trackster_h->at(i); + + // The following plots should be moved to HGVHistoProducerAlgo + // when we get rid of the MultiClusters and use only Tracksters + histo.raw_energy_->Fill(thisTrackster.raw_energy()); + histo.regr_energy_->Fill(thisTrackster.regressed_energy()); + histo.raw_energy_vs_regr_energy_->Fill(thisTrackster.regressed_energy(), thisTrackster.raw_energy()); + const auto& probs = thisTrackster.id_probabilities(); + std::vector sorted_probs_idx(probs.size()); + std::iota(begin(sorted_probs_idx), end(sorted_probs_idx), 0); + std::sort(begin(sorted_probs_idx), end(sorted_probs_idx), [&probs](int i, int j) { return probs[i] > probs[j]; }); + histo.id_prob_->Fill(sorted_probs_idx[0]); + if (!thisTrackster.vertices().empty()) { + histo.raw_energy_1plusLC_->Fill(thisTrackster.raw_energy()); + histo.regr_energy_1plusLC_->Fill(thisTrackster.regressed_energy()); + histo.raw_energy_vs_regr_energy_1plusLC_->Fill(thisTrackster.regressed_energy(), thisTrackster.raw_energy()); + histo.id_prob_1plusLC_->Fill(sorted_probs_idx[0]); + } + + // Plots on edges + for (const auto& edge : thisTrackster.edges()) { + auto& ic = layerClusters[edge[0]]; + auto& oc = layerClusters[edge[1]]; + auto const& cl_in = ic.hitsAndFractions()[0].first; + auto const& cl_out = oc.hitsAndFractions()[0].first; + auto const layer_in = rhtools_.getLayerWithOffset(cl_in); + auto const layer_out = rhtools_.getLayerWithOffset(cl_out); + histo.delta_energy_->Fill(oc.energy() - ic.energy()); + histo.delta_energy_relative_->Fill((oc.energy() - ic.energy()) / ic.energy()); + histo.delta_energy_vs_energy_->Fill(oc.energy() - ic.energy(), ic.energy()); + histo.delta_energy_vs_layer_->Fill(layer_in, oc.energy() - ic.energy()); + histo.delta_energy_relative_vs_layer_->Fill(layer_in, (oc.energy() - ic.energy()) / ic.energy()); + histo.delta_layer_->Fill(layer_out - layer_in); + + // Alpha angles + const auto& outer_outer_pos = oc.position(); + const auto& outer_inner_pos = ic.position(); + const auto& seed = thisTrackster.seedIndex(); + auto seedGlobalPos = math::XYZPoint( + ticlSeedingGlobal[0].origin.x(), ticlSeedingGlobal[0].origin.y(), ticlSeedingGlobal[0].origin.z()); + auto seedDirectionPos = outer_inner_pos; + if (thisTrackster.seedID().id() != 0) { + // Seed to trackster association is, at present, rather convoluted. + for (auto const& s : ticlSeedingTrk) { + if (s.index == seed) { + seedGlobalPos = math::XYZPoint(s.origin.x(), s.origin.y(), s.origin.z()); + seedDirectionPos = + math::XYZPoint(s.directionAtOrigin.x(), s.directionAtOrigin.y(), s.directionAtOrigin.z()); + break; + } + } + } + + auto alpha = (outer_inner_pos - seedGlobalPos).Dot(outer_outer_pos - outer_inner_pos) / + sqrt((outer_inner_pos - seedGlobalPos).Mag2() * (outer_outer_pos - outer_inner_pos).Mag2()); + auto alpha_alternative = (outer_outer_pos - seedGlobalPos).Dot(seedDirectionPos) / + sqrt((outer_outer_pos - seedGlobalPos).Mag2() * seedDirectionPos.Mag2()); + histo.angle_alpha_->Fill(alpha); + histo.angle_alpha_alternative_->Fill(alpha_alternative); + + // Beta angle is usually computed using 2 edges. Another inner loop + // is therefore needed. + std::vector> innerDoublets; + std::vector> outerDoublets; + for (const auto& otherEdge : thisTrackster.edges()) { + if (otherEdge[1] == edge[0]) { + innerDoublets.push_back(otherEdge); + } + if (edge[1] == otherEdge[0]) { + outerDoublets.push_back(otherEdge); + } + } + + histo.ingoing_links_vs_layer_->Fill(layer_in, innerDoublets.size()); + histo.outgoing_links_vs_layer_->Fill(layer_out, outerDoublets.size()); + for (const auto& inner : innerDoublets) { + const auto& inner_ic = layerClusters[inner[0]]; + const auto& inner_inner_pos = inner_ic.position(); + auto beta = (outer_inner_pos - inner_inner_pos).Dot(outer_outer_pos - inner_inner_pos) / + sqrt((outer_inner_pos - inner_inner_pos).Mag2() * (outer_outer_pos - inner_inner_pos).Mag2()); + histo.angle_beta_->Fill(beta); + histo.angle_beta_byLayer_[layer_in]->Fill(beta); + histo.angle_beta_w_byLayer_[layer_in]->Fill(beta, ic.energy()); + } + } + } + } +} + +void TICLTrackstersEdgesValidation::bookHistograms(DQMStore::IBooker& ibook, + edm::Run const& run, + edm::EventSetup const& iSetup, + Histograms_TICLTrackstersEdgesValidation& histos) const { + float eMin = 0., eThresh = 70., eMax = 500; + float eWidth[] = {0.5, 2.}; + std::vector eBins; + float eVal = eMin; + while (eVal <= eThresh) { + eBins.push_back(eVal); + eVal += eWidth[0]; + } + while (eVal < eMax) { + eVal += eWidth[1]; + eBins.push_back(eVal); + } + int eNBins = eBins.size() - 1; + + TString onePlusLC[] = {"1plus LC", "for tracksters with at least one LC"}; + TString trkers = "Tracksters"; + static const char* particle_kind[] = { + "photon", "electron", "muon", "neutral_pion", "charged_hadron", "neutral_hadron", "ambiguous", "unknown"}; + auto nCategory = sizeof(particle_kind) / sizeof(*particle_kind); + int labelIndex = 0; + for (const auto& trackster_token : tracksterTokens_) { + auto& histo = histos[trackster_token.index()]; + ibook.setCurrentFolder(folder_ + "HGCalValidator/" + trackstersCollectionsNames_[labelIndex]); + histo.number_ = ibook.book1D( + "Number of Tracksters per Event", "Number of Tracksters per Event;# Tracksters;Events", 250, 0., 250.); + // The following plots should be moved to HGVHistoProducerAlgo + // when we get rid of the MultiClusters and use only Tracksters + histo.raw_energy_ = ibook.book1D("Raw Energy", "Raw Energy;Raw Energy [GeV];" + trkers, eNBins, &eBins[0]); + histo.regr_energy_ = + ibook.book1D("Regressed Energy", "Regressed Energy;Regressed Energy [GeV];" + trkers, eNBins, &eBins[0]); + histo.raw_energy_vs_regr_energy_ = ibook.book2D("Raw Energy vs Regressed Energy", + "Raw vs Regressed Energy;Regressed Energy [GeV];Raw Energy [GeV]", + eNBins, + &eBins[0], + eNBins, + &eBins[0]); + histo.id_prob_ = + ibook.book1D("ID probability", "ID probability;category;Max ID probability", nCategory, 0, nCategory); + histo.raw_energy_1plusLC_ = ibook.book1D( + "Raw Energy " + onePlusLC[0], "Raw Energy " + onePlusLC[1] + ";Raw Energy [GeV];" + trkers, eNBins, &eBins[0]); + histo.regr_energy_1plusLC_ = ibook.book1D("Regressed Energy " + onePlusLC[0], + "Regressed Energy " + onePlusLC[1] + ";Regressed Energy [GeV];" + trkers, + eNBins, + &eBins[0]); + histo.raw_energy_vs_regr_energy_1plusLC_ = + ibook.book2D("Raw Energy vs Regressed Energy " + onePlusLC[0], + "Raw vs Regressed Energy " + onePlusLC[1] + ";Regressed Energy [GeV];Raw Energy [GeV]", + eNBins, + &eBins[0], + eNBins, + &eBins[0]); + histo.id_prob_1plusLC_ = ibook.book1D("ID probability " + onePlusLC[0], + "ID probability " + onePlusLC[1] + ";category;Max ID probability", + nCategory, + 0, + nCategory); + for (int iBin = 0; iBin < histo.id_prob_->getNbinsX(); iBin++) { + histo.id_prob_->setBinLabel(iBin + 1, particle_kind[iBin]); + histo.id_prob_1plusLC_->setBinLabel(iBin + 1, particle_kind[iBin]); + } + // Plots on edges + histo.delta_energy_ = ibook.book1D("Delta energy", "Delta Energy (O-I)", 800, -20., 20.); + histo.delta_energy_relative_ = + ibook.book1D("Relative Delta energy", "Relative Delta Energy (O-I)/I", 200, -10., 10.); + histo.delta_energy_vs_energy_ = + ibook.book2D("Energy vs Delta Energy", "Energy (I) vs Delta Energy (O-I)", 800, -20., 20., 200, 0., 20.); + histo.delta_energy_vs_layer_ = ibook.book2D( + "Delta Energy (O-I) vs Layer Number (I)", "Delta Energy (O-I) vs Layer Number (I)", 50, 0., 50., 800, -20., 20.); + histo.delta_energy_relative_vs_layer_ = ibook.book2D("Relative Delta Energy (O-I)_I vs Layer Number (I)", + "Relative Delta Energy (O-I)_I vs Layer Number (I)", + 50, + 0., + 50., + 200, + -10., + 10.); + histo.ingoing_links_vs_layer_ = + ibook.book2D("Ingoing links Layer Number", "Ingoing links vs Layer Number", 50, 0., 50., 40, 0., 40.); + histo.outgoing_links_vs_layer_ = + ibook.book2D("Outgoing links vs Layer Number", "Outgoing links vs Layer Number", 50, 0., 50., 40, 0., 40.); + histo.delta_layer_ = ibook.book1D("Delta Layer", "Delta Layer", 10, 0., 10.); + histo.angle_alpha_ = ibook.book1D("cosAngle Alpha", "cosAngle Alpha", 200, -1., 1.); + histo.angle_beta_ = ibook.book1D("cosAngle Beta", "cosAngle Beta", 200, -1., 1.); + histo.angle_alpha_alternative_ = ibook.book1D("cosAngle Alpha Alternative", "Angle Alpha Alternative", 200, 0., 1.); + for (int layer = 0; layer < 50; layer++) { + auto layerstr = std::to_string(layer + 1); + if (layerstr.length() < 2) + layerstr.insert(0, 2 - layerstr.length(), '0'); + histo.angle_beta_byLayer_.push_back( + ibook.book1D("cosAngle Beta on Layer " + layerstr, "cosAngle Beta on Layer " + layerstr, 200, -1., 1.)); + histo.angle_beta_w_byLayer_.push_back(ibook.book1D( + "cosAngle Beta Weighted on Layer " + layerstr, "cosAngle Beta Weighted on Layer " + layerstr, 200, -1., 1.)); + } + labelIndex++; + } +} + +void TICLTrackstersEdgesValidation::dqmBeginRun(edm::Run const& run, + edm::EventSetup const& iSetup, + Histograms_TICLTrackstersEdgesValidation& histograms) const { + edm::ESHandle geom = iSetup.getHandle(caloGeomToken_); + rhtools_.setGeometry(*geom); +} + +void TICLTrackstersEdgesValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + std::vector source_vector{edm::InputTag("ticlTrackstersTrk"), + edm::InputTag("ticlTrackstersTrkEM"), + edm::InputTag("ticlTrackstersEM"), + edm::InputTag("ticlTrackstersHAD"), + edm::InputTag("ticlTrackstersMerge")}; + desc.add>("tracksterCollections", source_vector); + desc.add("layerClusters", edm::InputTag("hgcalLayerClusters")); + desc.add("ticlSeedingGlobal", edm::InputTag("ticlSeedingGlobal")); + desc.add("ticlSeedingTrk", edm::InputTag("ticlSeedingTrk")); + desc.add("folder", "HGCAL/"); + descriptions.add("ticlTrackstersEdgesValidationDefault", desc); +} + +DEFINE_FWK_MODULE(TICLTrackstersEdgesValidation); diff --git a/Validation/HGCalValidation/python/CaloParticleSelectionForEfficiency_cfi.py b/Validation/HGCalValidation/python/CaloParticleSelectionForEfficiency_cfi.py index f3e9667d8c40f..5e84ad75752b6 100644 --- a/Validation/HGCalValidation/python/CaloParticleSelectionForEfficiency_cfi.py +++ b/Validation/HGCalValidation/python/CaloParticleSelectionForEfficiency_cfi.py @@ -2,7 +2,7 @@ CaloParticleSelectionForEfficiency = cms.PSet( ptMinCP = cms.double(0.5), - ptMaxCP = cms.double(250.), + ptMaxCP = cms.double(300.), minRapidityCP = cms.double(-3.1), maxRapidityCP = cms.double(3.1), #--z position of the origin vertex less than lipCP diff --git a/Validation/HGCalValidation/python/HGCalValidator_cfi.py b/Validation/HGCalValidation/python/HGCalValidator_cfi.py index 8a6925bcb1f6d..1e75ef6be8d07 100644 --- a/Validation/HGCalValidation/python/HGCalValidator_cfi.py +++ b/Validation/HGCalValidation/python/HGCalValidator_cfi.py @@ -35,11 +35,11 @@ #Select caloParticles for efficiency or pass through doCaloParticleSelection = cms.untracked.bool(True), #SimCluster related plots - dosimclustersPlots = cms.untracked.bool(True), + doSimClustersPlots = cms.untracked.bool(True), #Layer Cluster related plots - dolayerclustersPlots = cms.untracked.bool(True), + doLayerClustersPlots = cms.untracked.bool(True), #Multi Cluster related plots - domulticlustersPlots = cms.untracked.bool(True), + doMultiClustersPlots = cms.untracked.bool(True), #The cumulative material budget in front of each layer. To be more specific, it #is the material budget just in front of the active material (not including it). diff --git a/Validation/HGCalValidation/python/hgcalPlots.py b/Validation/HGCalValidation/python/hgcalPlots.py index 627126d4719af..009f55a6c05a3 100644 --- a/Validation/HGCalValidation/python/hgcalPlots.py +++ b/Validation/HGCalValidation/python/hgcalPlots.py @@ -2493,8 +2493,11 @@ def append_hgcalCaloParticlesPlots(files, collection = '-211', name_collection = "Rec-matched Hits Sum Energy vs layer"] dqmfolder = "DQMData/Run 1/HGCAL/Run summary/HGCalValidator/SelectedCaloParticles/" + collection - print(dqmfolder) templateFile = ROOT.TFile.Open(files[0]) # assuming all files have same structure + if not gDirectory.GetDirectory(dqmfolder): + print("Error: GeneralInfo directory %s not found in DQM file, exit"%dqmfolder) + return hgcalTrackstersPlotter + keys = gDirectory.GetDirectory(dqmfolder,True).GetListOfKeys() key = keys[0] while key: @@ -2511,7 +2514,7 @@ def append_hgcalCaloParticlesPlots(files, collection = '-211', name_collection = ncols=1) if name in list_2D_histos : - pg= PlotOnSideGroup(fileName.Data(), + pg= PlotOnSideGroup(plotName.Data(), Plot(name, xtitle=obj.GetXaxis().GetTitle(), ytitle=obj.GetYaxis().GetTitle(), drawCommand = "COLZ", @@ -2533,6 +2536,123 @@ def append_hgcalCaloParticlesPlots(files, collection = '-211', name_collection = return hgcalCaloParticlesPlotter +#================================================================================================= +def create_hgcalTrackstersPlotter(files, collection = 'ticlTrackstersMerge', name_collection = "MultiClustersMerge"): + grouped = {"cosAngle Beta": PlotGroup("cosAngle_Beta_per_layer",[],ncols=10), "cosAngle Beta Weighted": PlotGroup("cosAngle_Beta_Weighted_per_layer",[],ncols=10)} + groupingFlag = " on Layer " + + hgcalTrackstersPlotter = Plotter() + dqmfolder = "DQMData/Run 1/HGCAL/Run summary/HGCalValidator/" + collection + #_multiplicity_tracksters_numberOfEventsHistogram = dqmfolder+"/Number of Trackster per Event" + + _common["ymin"] = 0.0 + _common["staty"] = 0.85 + templateFile = ROOT.TFile.Open(files[0]) # assuming all files have same structure + if not gDirectory.GetDirectory(dqmfolder): + print("Error: GeneralInfo directory %s not found in DQM file, exit"%dqmfolder) + return hgcalTrackstersPlotter + + keys = gDirectory.GetDirectory(dqmfolder,True).GetListOfKeys() + key = keys[0] + while key: + obj = key.ReadObj() + name = obj.GetName() + plotName = TString(name) + plotName.ReplaceAll(" ","_") + + if groupingFlag in name: + for group in grouped: + if group+groupingFlag in name: + grouped[group].append(Plot(name, + xtitle=obj.GetXaxis().GetTitle(), ytitle=obj.GetYaxis().GetTitle(), + **_common) + ) + else: + pg = None + if obj.InheritsFrom("TH2"): + pg = PlotOnSideGroup(plotName.Data(), + Plot(name, + xtitle=obj.GetXaxis().GetTitle(), ytitle=obj.GetYaxis().GetTitle(), + drawCommand = "COLZ", + **_common), + ncols=1) + else: + pg = PlotGroup(plotName.Data(), + [Plot(name, + xtitle=obj.GetXaxis().GetTitle(), ytitle=obj.GetYaxis().GetTitle(), + drawCommand = "COLZ", # ineffective for TH1 + **_common) + ], + ncols=1, legendDh=-0.03 * len(files)) + + hgcalTrackstersPlotter.append(name_collection+"_TICLDebugger", + [dqmfolder], PlotFolder(pg, + loopSubFolders=False, + purpose=PlotPurpose.Timing, page="MultiClusters", section=name_collection) + #numberOfEventsHistogram=_multiplicity_tracksters_numberOfEventsHistogram) + ) + + key = keys.After(key) + + for group in grouped: + hgcalTrackstersPlotter.append(name_collection+"_TICLDebugger", + [dqmfolder], PlotFolder(grouped[group], + loopSubFolders=False, + purpose=PlotPurpose.Timing, page="MultiClusters", section=name_collection) + #numberOfEventsHistogram=_multiplicity_tracksters_numberOfEventsHistogram) + ) + + templateFile.Close() + + return hgcalTrackstersPlotter + +#================================================================================================= +_common_Calo = {"stat": False, "drawStyle": "hist", "staty": 0.65, "ymin": 0.0, "ylog": False} + +hgcalCaloParticlesPlotter = Plotter() + +def append_hgcalCaloParticlesPlots(files, collection = '-211', name_collection = "pion-"): + dqmfolder = "DQMData/Run 1/HGCAL/Run summary/HGCalValidator/SelectedCaloParticles/" + collection + print(dqmfolder) +# _common["ymin"] = 0.0 + templateFile = ROOT.TFile.Open(files[0]) # assuming all files have same structure + keys = gDirectory.GetDirectory(dqmfolder,True).GetListOfKeys() + key = keys[0] + while key: + obj = key.ReadObj() + name = obj.GetName() + plotName = TString(name) + plotName.ReplaceAll(" ","_") + pg= PlotGroup(plotName.Data(),[ + Plot(name, + xtitle=obj.GetXaxis().GetTitle(), ytitle=obj.GetYaxis().GetTitle(), + drawCommand = "", # may want to customize for TH2 (colz, etc.) + normalizeToNumberOfEvents = True, **_common_Calo) + ], + ncols=1) + + if obj.InheritsFrom("TH2"): + pg= PlotOnSideGroup(plotName.Data(), + Plot(name, + xtitle=obj.GetXaxis().GetTitle(), ytitle=obj.GetYaxis().GetTitle(), + drawCommand = "COLZ", + normalizeToNumberOfEvents = True, **_common_Calo), + ncols=1) + + hgcalCaloParticlesPlotter.append("CaloParticles_"+name_collection, [ + dqmfolder + ], PlotFolder( + pg, + loopSubFolders=False, + purpose=PlotPurpose.Timing, page="CaloParticles", section=name_collection) + ) + + key = keys.After(key) + + templateFile.Close() + + return hgcalCaloParticlesPlotter + #================================================================================================= # hitValidation def _hgcalHitFolders(dirName="HGCalSimHitsV/HGCalEESensitive"): diff --git a/Validation/HGCalValidation/python/ticlTrackstersEdgesValidation_cfi.py b/Validation/HGCalValidation/python/ticlTrackstersEdgesValidation_cfi.py new file mode 100644 index 0000000000000..3bbf5c8607edc --- /dev/null +++ b/Validation/HGCalValidation/python/ticlTrackstersEdgesValidation_cfi.py @@ -0,0 +1,4 @@ +import FWCore.ParameterSet.Config as cms + +from Validation.HGCalValidation.ticlTrackstersEdgesValidationDefault_cfi import ticlTrackstersEdgesValidationDefault as _ticlTrackstersEdgesValidationDefault +ticlTrackstersEdgesValidation = _ticlTrackstersEdgesValidationDefault.clone() diff --git a/Validation/HGCalValidation/scripts/makeHGCalValidationPlots.py b/Validation/HGCalValidation/scripts/makeHGCalValidationPlots.py index 3f7eaa6bef211..5d77728a0dc8d 100755 --- a/Validation/HGCalValidation/scripts/makeHGCalValidationPlots.py +++ b/Validation/HGCalValidation/scripts/makeHGCalValidationPlots.py @@ -13,8 +13,7 @@ trackstersIters = ["ticlMultiClustersFromTrackstersMerge", "ticlMultiClustersFromTrackstersMIP", "ticlMultiClustersFromTrackstersTrk","ticlMultiClustersFromTrackstersTrkEM", - "ticlMultiClustersFromTrackstersEM", "ticlMultiClustersFromTrackstersHAD", - "ticlMultiClustersFromTrackstersDummy"] + "ticlMultiClustersFromTrackstersEM", "ticlMultiClustersFromTrackstersHAD"] simClustersGeneralLabel = 'simClusters' layerClustersGeneralLabel = 'hgcalLayerClusters' @@ -74,6 +73,11 @@ def main(opts): tracksterCollection = i_iter.replace("ticlMultiClustersFromTracksters","ticlTracksters") hgcalPlots.append_hgcalMultiClustersPlots(i_iter, tracksterCollection) val.doPlots(hgcmulticlus, plotterDrawArgs=drawArgs) + # TICLTrackstersEdges plots + for i_iter in trackstersIters : + tracksterCollection = i_iter.replace("ticlMultiClustersFromTracksters","ticlTracksters") + hgctracksters = [hgcalPlots.create_hgcalTrackstersPlotter(sample.files(), tracksterCollection, tracksterCollection)] + val.doPlots(hgctracksters, plotterDrawArgs=drawArgs) elif (opts.collection == caloParticlesLabel): particletypes = {"pion-":"-211", "pion+":"211", "pion0": "111", "muon-": "-13", "muon+":"13", @@ -135,6 +139,11 @@ def main(opts): tracksterCollection = i_iter.replace("ticlMultiClustersFromTracksters","ticlTracksters") hgcalPlots.append_hgcalMultiClustersPlots(i_iter, tracksterCollection) val.doPlots(hgcmulticlus, plotterDrawArgs=drawArgs) + #TICLTrackstersEdges plots + for i_iter in trackstersIters : + tracksterCollection = i_iter.replace("ticlMultiClustersFromTracksters","ticlTracksters") + hgctracksters = [hgcalPlots.create_hgcalTrackstersPlotter(sample.files(), tracksterCollection, tracksterCollection)] + val.doPlots(hgctracksters, plotterDrawArgs=drawArgs) if opts.no_html: diff --git a/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc b/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc index 35500a2cb63ee..759950175b85d 100644 --- a/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc +++ b/Validation/HGCalValidation/src/HGVHistoProducerAlgo.cc @@ -1030,6 +1030,12 @@ void HGVHistoProducerAlgo::bookMultiClusterHistos(DQMStore::IBooker& ibook, nintSharedEneFrac_, minMCLSharedEneFrac_, maxMCLSharedEneFrac_)); + histograms.h_sharedenergy_caloparticle2multicl_assoc.push_back( + ibook.book1D("SharedEnergy_caloparticle2multicl_assoc", + "Shared Energy of Associated CaloParticle per Multi Cluster", + nintSharedEneFrac_, + minMCLSharedEneFrac_, + maxMCLSharedEneFrac_)); histograms.h_sharedenergy_caloparticle2multicl_vs_eta.push_back( ibook.bookProfile("SharedEnergy_caloparticle2multicl_vs_eta", "Shared Energy of CaloParticle vs #eta per best Multi Cluster", @@ -1250,13 +1256,13 @@ void HGVHistoProducerAlgo::fill_caloparticle_histos(const Histograms& histograms std::map totenergy_layer; for (auto const& sc : caloparticle.simClusters()) { + LogDebug("HGCalValidator") << " This sim cluster has " << sc->hits_and_fractions().size() << " simHits and " + << sc->energy() << " energy. " << std::endl; simHits += sc->hits_and_fractions().size(); - for (auto const& h_and_f : sc->hits_and_fractions()) { const auto hitDetId = h_and_f.first; int layerId = recHitTools_->getLayerWithOffset(hitDetId) + layers * ((recHitTools_->zside(hitDetId) + 1) >> 1) - 1; - // set to 0 if matched RecHit not found int layerId_matched_min = 999; int layerId_matched_max = 0; @@ -1279,6 +1285,8 @@ void HGVHistoProducerAlgo::fill_caloparticle_histos(const Histograms& histograms if (caloparticle.simClusters().size() == 1) histograms.h_caloparticle_nHits_matched_energy_layer_1SimCl.at(pdgid)->Fill(layerId, hit->energy() * h_and_f.second); + } else { + LogDebug("HGCalValidator") << " matched to RecHit NOT found !" << std::endl; } minLayerId = std::min(minLayerId, layerId); @@ -1286,6 +1294,7 @@ void HGVHistoProducerAlgo::fill_caloparticle_histos(const Histograms& histograms minLayerId_matched = std::min(minLayerId_matched, layerId_matched_min); maxLayerId_matched = std::max(maxLayerId_matched, layerId_matched_max); } + LogDebug("HGCalValidator") << std::endl; } histograms.h_caloparticle_firstlayer.at(pdgid)->Fill(minLayerId); histograms.h_caloparticle_lastlayer.at(pdgid)->Fill(maxLayerId); @@ -2817,6 +2826,7 @@ void HGVHistoProducerAlgo::multiClusters_to_CaloParticles(const Histograms& hist multiClusters[bestmclId].energy() / CPenergy); histograms.h_sharedenergy_caloparticle2multicl_vs_phi[count]->Fill(cP[cpId].g4Tracks()[0].momentum().phi(), multiClusters[bestmclId].energy() / CPenergy); + histograms.h_sharedenergy_caloparticle2multicl_assoc[count]->Fill(mclsharedenergyfrac[cpId][bestmclId]); } if (assocDup >= 2) { auto match = std::find_if(std::begin(score3d[cpId]), std::end(score3d[cpId]), is_assoc); diff --git a/Validation/RecoTrack/python/plotting/plotting.py b/Validation/RecoTrack/python/plotting/plotting.py index db39adce17885..e09b6530e6c5d 100644 --- a/Validation/RecoTrack/python/plotting/plotting.py +++ b/Validation/RecoTrack/python/plotting/plotting.py @@ -1967,10 +1967,11 @@ def _doStats(h, col, dy): if self._fit: st.SetOptFit(0o010) st.SetOptStat(1001) + st.SetOptStat(1110) st.SetX1NDC(startingX) st.SetX2NDC(startingX+0.3) st.SetY1NDC(startingY+dy) - st.SetY2NDC(startingY+dy+0.15) + st.SetY2NDC(startingY+dy+0.12) st.SetTextColor(col) dy = 0.0 @@ -1979,7 +1980,7 @@ def _doStats(h, col, dy): dy += self._statyadjust[i] _doStats(h, _plotStylesColor[i], dy) - dy -= 0.19 + dy -= 0.16 def _normalize(self): """Normalise histograms to unit area"""