Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:32:35 

0001 #include <string>
0002 #include <unordered_map>
0003 #include <numeric>
0004 
0005 // user include files
0006 #include "DataFormats/Math/interface/Point3D.h"
0007 #include "FWCore/Framework/interface/Frameworkfwd.h"
0008 #include "DQMServices/Core/interface/DQMGlobalEDAnalyzer.h"
0009 
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/MakerMacros.h"
0012 
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 
0015 #include "FWCore/Utilities/interface/transform.h"
0016 #include "DataFormats/HGCalReco/interface/Trackster.h"
0017 #include "DataFormats/HGCalReco/interface/TICLSeedingRegion.h"
0018 
0019 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
0020 #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
0021 
0022 using namespace ticl;
0023 
0024 struct Histogram_TICLTrackstersEdgesValidation {
0025   dqm::reco::MonitorElement* number_;
0026   dqm::reco::MonitorElement *raw_energy_, *raw_energy_1plusLC_;
0027   dqm::reco::MonitorElement *regr_energy_, *regr_energy_1plusLC_;
0028   dqm::reco::MonitorElement *raw_energy_vs_regr_energy_, *raw_energy_vs_regr_energy_1plusLC_;
0029   dqm::reco::MonitorElement *id_prob_, *id_prob_1plusLC_;
0030   dqm::reco::MonitorElement* delta_energy_;
0031   dqm::reco::MonitorElement* delta_energy_relative_;
0032   dqm::reco::MonitorElement* delta_energy_vs_energy_;
0033   dqm::reco::MonitorElement* delta_energy_vs_layer_;
0034   dqm::reco::MonitorElement* delta_energy_relative_vs_layer_;
0035   dqm::reco::MonitorElement* delta_layer_;
0036   dqm::reco::MonitorElement* ingoing_links_vs_layer_;
0037   dqm::reco::MonitorElement* outgoing_links_vs_layer_;
0038   // For the definition of the angles, read http://hgcal.web.cern.ch/hgcal/Reconstruction/Tutorial/
0039   dqm::reco::MonitorElement* angle_alpha_;
0040   dqm::reco::MonitorElement* angle_alpha_alternative_;
0041   dqm::reco::MonitorElement* angle_beta_;
0042   std::vector<dqm::reco::MonitorElement*> angle_beta_byLayer_;
0043   std::vector<dqm::reco::MonitorElement*> angle_beta_w_byLayer_;
0044 };
0045 
0046 using Histograms_TICLTrackstersEdgesValidation =
0047     std::unordered_map<unsigned int, Histogram_TICLTrackstersEdgesValidation>;
0048 
0049 class TICLTrackstersEdgesValidation : public DQMGlobalEDAnalyzer<Histograms_TICLTrackstersEdgesValidation> {
0050 public:
0051   explicit TICLTrackstersEdgesValidation(const edm::ParameterSet&);
0052   ~TICLTrackstersEdgesValidation() override;
0053 
0054   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0055 
0056 private:
0057   void bookHistograms(DQMStore::IBooker&,
0058                       edm::Run const&,
0059                       edm::EventSetup const&,
0060                       Histograms_TICLTrackstersEdgesValidation&) const override;
0061 
0062   void dqmAnalyze(edm::Event const&,
0063                   edm::EventSetup const&,
0064                   Histograms_TICLTrackstersEdgesValidation const&) const override;
0065   void dqmBeginRun(edm::Run const&, edm::EventSetup const&, Histograms_TICLTrackstersEdgesValidation&) const override;
0066 
0067   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeomToken_;
0068   std::string folder_;
0069   std::vector<std::string> trackstersCollectionsNames_;
0070   std::vector<edm::EDGetTokenT<std::vector<Trackster>>> tracksterTokens_;
0071   edm::EDGetTokenT<std::vector<reco::CaloCluster>> layerClustersToken_;
0072   edm::EDGetTokenT<std::vector<TICLSeedingRegion>> ticlSeedingGlobalToken_;
0073   edm::EDGetTokenT<std::vector<TICLSeedingRegion>> ticlSeedingTrkToken_;
0074   mutable hgcal::RecHitTools rhtools_;
0075 };
0076 
0077 TICLTrackstersEdgesValidation::TICLTrackstersEdgesValidation(const edm::ParameterSet& iConfig)
0078     : caloGeomToken_(esConsumes<CaloGeometry, CaloGeometryRecord>()),
0079       folder_(iConfig.getParameter<std::string>("folder")) {
0080   tracksterTokens_ = edm::vector_transform(iConfig.getParameter<std::vector<edm::InputTag>>("tracksterCollections"),
0081                                            [this](edm::InputTag const& tag) {
0082                                              trackstersCollectionsNames_.emplace_back(tag.label());
0083                                              return consumes<std::vector<Trackster>>(tag);
0084                                            });
0085   layerClustersToken_ = consumes<std::vector<reco::CaloCluster>>(iConfig.getParameter<edm::InputTag>("layerClusters"));
0086   ticlSeedingGlobalToken_ =
0087       consumes<std::vector<TICLSeedingRegion>>(iConfig.getParameter<edm::InputTag>("ticlSeedingGlobal"));
0088   ticlSeedingTrkToken_ =
0089       consumes<std::vector<TICLSeedingRegion>>(iConfig.getParameter<edm::InputTag>("ticlSeedingTrk"));
0090 }
0091 
0092 TICLTrackstersEdgesValidation::~TICLTrackstersEdgesValidation() {}
0093 
0094 void TICLTrackstersEdgesValidation::dqmAnalyze(edm::Event const& iEvent,
0095                                                edm::EventSetup const& iSetup,
0096                                                Histograms_TICLTrackstersEdgesValidation const& histos) const {
0097   edm::Handle<std::vector<reco::CaloCluster>> layerClustersH;
0098   iEvent.getByToken(layerClustersToken_, layerClustersH);
0099   auto const& layerClusters = *layerClustersH.product();
0100 
0101   edm::Handle<std::vector<TICLSeedingRegion>> ticlSeedingGlobalH;
0102   iEvent.getByToken(ticlSeedingGlobalToken_, ticlSeedingGlobalH);
0103   auto const& ticlSeedingGlobal = *ticlSeedingGlobalH.product();
0104 
0105   edm::Handle<std::vector<TICLSeedingRegion>> ticlSeedingTrkH;
0106   iEvent.getByToken(ticlSeedingTrkToken_, ticlSeedingTrkH);
0107   auto const& ticlSeedingTrk = *ticlSeedingTrkH.product();
0108 
0109   for (const auto& trackster_token : tracksterTokens_) {
0110     edm::Handle<std::vector<Trackster>> trackster_h;
0111     iEvent.getByToken(trackster_token, trackster_h);
0112     auto numberOfTracksters = trackster_h->size();
0113     //using .at() as [] is not const
0114     const auto& histo = histos.at(trackster_token.index());
0115     histo.number_->Fill(numberOfTracksters);
0116     for (unsigned int i = 0; i < numberOfTracksters; ++i) {
0117       const auto& thisTrackster = trackster_h->at(i);
0118 
0119       // The following plots should be moved to HGVHistoProducerAlgo
0120       // when we get rid of the MultiClusters and use only Tracksters
0121       histo.raw_energy_->Fill(thisTrackster.raw_energy());
0122       histo.regr_energy_->Fill(thisTrackster.regressed_energy());
0123       histo.raw_energy_vs_regr_energy_->Fill(thisTrackster.regressed_energy(), thisTrackster.raw_energy());
0124       const auto& probs = thisTrackster.id_probabilities();
0125       std::vector<int> sorted_probs_idx(probs.size());
0126       std::iota(begin(sorted_probs_idx), end(sorted_probs_idx), 0);
0127       std::sort(begin(sorted_probs_idx), end(sorted_probs_idx), [&probs](int i, int j) { return probs[i] > probs[j]; });
0128       histo.id_prob_->Fill(sorted_probs_idx[0]);
0129       if (!thisTrackster.vertices().empty()) {
0130         histo.raw_energy_1plusLC_->Fill(thisTrackster.raw_energy());
0131         histo.regr_energy_1plusLC_->Fill(thisTrackster.regressed_energy());
0132         histo.raw_energy_vs_regr_energy_1plusLC_->Fill(thisTrackster.regressed_energy(), thisTrackster.raw_energy());
0133         histo.id_prob_1plusLC_->Fill(sorted_probs_idx[0]);
0134       }
0135 
0136       // Plots on edges
0137       for (const auto& edge : thisTrackster.edges()) {
0138         auto& ic = layerClusters[edge[0]];
0139         auto& oc = layerClusters[edge[1]];
0140         auto const& cl_in = ic.hitsAndFractions()[0].first;
0141         auto const& cl_out = oc.hitsAndFractions()[0].first;
0142         auto const layer_in = rhtools_.getLayerWithOffset(cl_in);
0143         auto const layer_out = rhtools_.getLayerWithOffset(cl_out);
0144         histo.delta_energy_->Fill(oc.energy() - ic.energy());
0145         histo.delta_energy_relative_->Fill((oc.energy() - ic.energy()) / ic.energy());
0146         histo.delta_energy_vs_energy_->Fill(oc.energy() - ic.energy(), ic.energy());
0147         histo.delta_energy_vs_layer_->Fill(layer_in, oc.energy() - ic.energy());
0148         histo.delta_energy_relative_vs_layer_->Fill(layer_in, (oc.energy() - ic.energy()) / ic.energy());
0149         histo.delta_layer_->Fill(layer_out - layer_in);
0150 
0151         // Alpha angles
0152         const auto& outer_outer_pos = oc.position();
0153         const auto& outer_inner_pos = ic.position();
0154         const auto& seed = thisTrackster.seedIndex();
0155         auto seedGlobalPos = math::XYZPoint(
0156             ticlSeedingGlobal[0].origin.x(), ticlSeedingGlobal[0].origin.y(), ticlSeedingGlobal[0].origin.z());
0157         auto seedDirectionPos = outer_inner_pos;
0158         if (thisTrackster.seedID().id() != 0) {
0159           // Seed to trackster association is, at present, rather convoluted.
0160           for (auto const& s : ticlSeedingTrk) {
0161             if (s.index == seed) {
0162               seedGlobalPos = math::XYZPoint(s.origin.x(), s.origin.y(), s.origin.z());
0163               seedDirectionPos =
0164                   math::XYZPoint(s.directionAtOrigin.x(), s.directionAtOrigin.y(), s.directionAtOrigin.z());
0165               break;
0166             }
0167           }
0168         }
0169 
0170         auto alpha = (outer_inner_pos - seedGlobalPos).Dot(outer_outer_pos - outer_inner_pos) /
0171                      sqrt((outer_inner_pos - seedGlobalPos).Mag2() * (outer_outer_pos - outer_inner_pos).Mag2());
0172         auto alpha_alternative = (outer_outer_pos - seedGlobalPos).Dot(seedDirectionPos) /
0173                                  sqrt((outer_outer_pos - seedGlobalPos).Mag2() * seedDirectionPos.Mag2());
0174         histo.angle_alpha_->Fill(alpha);
0175         histo.angle_alpha_alternative_->Fill(alpha_alternative);
0176 
0177         // Beta angle is usually computed using 2 edges. Another inner loop
0178         // is therefore needed.
0179         std::vector<std::array<unsigned int, 2>> innerDoublets;
0180         std::vector<std::array<unsigned int, 2>> outerDoublets;
0181         for (const auto& otherEdge : thisTrackster.edges()) {
0182           if (otherEdge[1] == edge[0]) {
0183             innerDoublets.push_back(otherEdge);
0184           }
0185           if (edge[1] == otherEdge[0]) {
0186             outerDoublets.push_back(otherEdge);
0187           }
0188         }
0189 
0190         histo.ingoing_links_vs_layer_->Fill(layer_in, innerDoublets.size());
0191         histo.outgoing_links_vs_layer_->Fill(layer_out, outerDoublets.size());
0192         for (const auto& inner : innerDoublets) {
0193           const auto& inner_ic = layerClusters[inner[0]];
0194           const auto& inner_inner_pos = inner_ic.position();
0195           auto beta = (outer_inner_pos - inner_inner_pos).Dot(outer_outer_pos - inner_inner_pos) /
0196                       sqrt((outer_inner_pos - inner_inner_pos).Mag2() * (outer_outer_pos - inner_inner_pos).Mag2());
0197           histo.angle_beta_->Fill(beta);
0198           histo.angle_beta_byLayer_[layer_in]->Fill(beta);
0199           histo.angle_beta_w_byLayer_[layer_in]->Fill(beta, ic.energy());
0200         }
0201       }
0202     }
0203   }
0204 }
0205 
0206 void TICLTrackstersEdgesValidation::bookHistograms(DQMStore::IBooker& ibook,
0207                                                    edm::Run const& run,
0208                                                    edm::EventSetup const& iSetup,
0209                                                    Histograms_TICLTrackstersEdgesValidation& histos) const {
0210   float eMin = 0., eThresh = 70., eMax = 500;
0211   float eWidth[] = {0.5, 2.};
0212   std::vector<float> eBins;
0213   float eVal = eMin;
0214   while (eVal <= eThresh) {
0215     eBins.push_back(eVal);
0216     eVal += eWidth[0];
0217   }
0218   while (eVal < eMax) {
0219     eVal += eWidth[1];
0220     eBins.push_back(eVal);
0221   }
0222   int eNBins = eBins.size() - 1;
0223 
0224   TString onePlusLC[] = {"1plus LC", "for tracksters with at least one LC"};
0225   TString trkers = "Tracksters";
0226   static const char* particle_kind[] = {
0227       "photon", "electron", "muon", "neutral_pion", "charged_hadron", "neutral_hadron", "ambiguous", "unknown"};
0228   auto nCategory = sizeof(particle_kind) / sizeof(*particle_kind);
0229   int labelIndex = 0;
0230   for (const auto& trackster_token : tracksterTokens_) {
0231     auto& histo = histos[trackster_token.index()];
0232     ibook.setCurrentFolder(folder_ + "HGCalValidator/" + trackstersCollectionsNames_[labelIndex]);
0233     histo.number_ = ibook.book1D(
0234         "Number of Tracksters per Event", "Number of Tracksters per Event;# Tracksters;Events", 250, 0., 250.);
0235     // The following plots should be moved to HGVHistoProducerAlgo
0236     // when we get rid of the MultiClusters and use only Tracksters
0237     histo.raw_energy_ = ibook.book1D("Raw Energy", "Raw Energy;Raw Energy [GeV];" + trkers, eNBins, &eBins[0]);
0238     histo.regr_energy_ =
0239         ibook.book1D("Regressed Energy", "Regressed Energy;Regressed Energy [GeV];" + trkers, eNBins, &eBins[0]);
0240     histo.raw_energy_vs_regr_energy_ = ibook.book2D("Raw Energy vs Regressed Energy",
0241                                                     "Raw vs Regressed Energy;Regressed Energy [GeV];Raw Energy [GeV]",
0242                                                     eNBins,
0243                                                     &eBins[0],
0244                                                     eNBins,
0245                                                     &eBins[0]);
0246     histo.id_prob_ =
0247         ibook.book1D("ID probability", "ID probability;category;Max ID probability", nCategory, 0, nCategory);
0248     histo.raw_energy_1plusLC_ = ibook.book1D(
0249         "Raw Energy " + onePlusLC[0], "Raw Energy " + onePlusLC[1] + ";Raw Energy [GeV];" + trkers, eNBins, &eBins[0]);
0250     histo.regr_energy_1plusLC_ = ibook.book1D("Regressed Energy " + onePlusLC[0],
0251                                               "Regressed Energy " + onePlusLC[1] + ";Regressed Energy [GeV];" + trkers,
0252                                               eNBins,
0253                                               &eBins[0]);
0254     histo.raw_energy_vs_regr_energy_1plusLC_ =
0255         ibook.book2D("Raw Energy vs Regressed Energy " + onePlusLC[0],
0256                      "Raw vs Regressed Energy " + onePlusLC[1] + ";Regressed Energy [GeV];Raw Energy [GeV]",
0257                      eNBins,
0258                      &eBins[0],
0259                      eNBins,
0260                      &eBins[0]);
0261     histo.id_prob_1plusLC_ = ibook.book1D("ID probability " + onePlusLC[0],
0262                                           "ID probability " + onePlusLC[1] + ";category;Max ID probability",
0263                                           nCategory,
0264                                           0,
0265                                           nCategory);
0266     for (int iBin = 0; iBin < histo.id_prob_->getNbinsX(); iBin++) {
0267       histo.id_prob_->setBinLabel(iBin + 1, particle_kind[iBin]);
0268       histo.id_prob_1plusLC_->setBinLabel(iBin + 1, particle_kind[iBin]);
0269     }
0270     // Plots on edges
0271     histo.delta_energy_ = ibook.book1D("Delta energy", "Delta Energy (O-I)", 800, -20., 20.);
0272     histo.delta_energy_relative_ =
0273         ibook.book1D("Relative Delta energy", "Relative Delta Energy (O-I)/I", 200, -10., 10.);
0274     histo.delta_energy_vs_energy_ =
0275         ibook.book2D("Energy vs Delta Energy", "Energy (I) vs Delta Energy (O-I)", 800, -20., 20., 200, 0., 20.);
0276     histo.delta_energy_vs_layer_ = ibook.book2D(
0277         "Delta Energy (O-I) vs Layer Number (I)", "Delta Energy (O-I) vs Layer Number (I)", 50, 0., 50., 800, -20., 20.);
0278     histo.delta_energy_relative_vs_layer_ = ibook.book2D("Relative Delta Energy (O-I)_I vs Layer Number (I)",
0279                                                          "Relative Delta Energy (O-I)_I vs Layer Number (I)",
0280                                                          50,
0281                                                          0.,
0282                                                          50.,
0283                                                          200,
0284                                                          -10.,
0285                                                          10.);
0286     histo.ingoing_links_vs_layer_ =
0287         ibook.book2D("Ingoing links Layer Number", "Ingoing links vs Layer Number", 50, 0., 50., 40, 0., 40.);
0288     histo.outgoing_links_vs_layer_ =
0289         ibook.book2D("Outgoing links vs Layer Number", "Outgoing links vs Layer Number", 50, 0., 50., 40, 0., 40.);
0290     histo.delta_layer_ = ibook.book1D("Delta Layer", "Delta Layer", 10, 0., 10.);
0291     histo.angle_alpha_ = ibook.book1D("cosAngle Alpha", "cosAngle Alpha", 200, -1., 1.);
0292     histo.angle_beta_ = ibook.book1D("cosAngle Beta", "cosAngle Beta", 200, -1., 1.);
0293     histo.angle_alpha_alternative_ = ibook.book1D("cosAngle Alpha Alternative", "Angle Alpha Alternative", 200, 0., 1.);
0294     for (int layer = 0; layer < 50; layer++) {
0295       auto layerstr = std::to_string(layer + 1);
0296       if (layerstr.length() < 2)
0297         layerstr.insert(0, 2 - layerstr.length(), '0');
0298       histo.angle_beta_byLayer_.push_back(
0299           ibook.book1D("cosAngle Beta on Layer " + layerstr, "cosAngle Beta on Layer " + layerstr, 200, -1., 1.));
0300       histo.angle_beta_w_byLayer_.push_back(ibook.book1D(
0301           "cosAngle Beta Weighted on Layer " + layerstr, "cosAngle Beta Weighted on Layer " + layerstr, 200, -1., 1.));
0302     }
0303     labelIndex++;
0304   }
0305 }
0306 
0307 void TICLTrackstersEdgesValidation::dqmBeginRun(edm::Run const& run,
0308                                                 edm::EventSetup const& iSetup,
0309                                                 Histograms_TICLTrackstersEdgesValidation& histograms) const {
0310   edm::ESHandle<CaloGeometry> geom = iSetup.getHandle(caloGeomToken_);
0311   rhtools_.setGeometry(*geom);
0312 }
0313 
0314 void TICLTrackstersEdgesValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0315   edm::ParameterSetDescription desc;
0316   std::vector<edm::InputTag> source_vector{edm::InputTag("ticlTrackstersTrk"),
0317                                            edm::InputTag("ticlTrackstersTrkEM"),
0318                                            edm::InputTag("ticlTrackstersEM"),
0319                                            edm::InputTag("ticlTrackstersHAD"),
0320                                            edm::InputTag("ticlTrackstersMerge")};
0321   desc.add<std::vector<edm::InputTag>>("tracksterCollections", source_vector);
0322   desc.add<edm::InputTag>("layerClusters", edm::InputTag("hgcalMergeLayerClusters"));
0323   desc.add<edm::InputTag>("ticlSeedingGlobal", edm::InputTag("ticlSeedingGlobal"));
0324   desc.add<edm::InputTag>("ticlSeedingTrk", edm::InputTag("ticlSeedingTrk"));
0325   desc.add<std::string>("folder", "HGCAL/");
0326   descriptions.add("ticlTrackstersEdgesValidationDefault", desc);
0327 }
0328 
0329 DEFINE_FWK_MODULE(TICLTrackstersEdgesValidation);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!