TICL/test/TiclDebugger.cc

0001 //
0002 // Original Author:  Marco Rovere
0003 //         Created:  Fri May  1 07:21:02 CEST 2020
0004 //
0005 //
0006 //
0007 // system include files
0008 #include <memory>
0009 #include <iostream>
0010 #include <numeric>
0011 #include <algorithm>
0012
0013 // user include files
0014 #include "FWCore/Framework/interface/Frameworkfwd.h"
0015 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0016 #include "FWCore/Framework/interface/ConsumesCollector.h"
0017
0018 #include "FWCore/Framework/interface/Event.h"
0019 #include "FWCore/Framework/interface/EventSetup.h"
0020 #include "FWCore/Framework/interface/ESHandle.h"
0021 #include "FWCore/Framework/interface/MakerMacros.h"
0022
0023 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0024
0025 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
0026 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0027 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0028 #include "DataFormats/HGCalReco/interface/Trackster.h"
0029 #include "DataFormats/TrackReco/interface/Track.h"
0030 #include "DataFormats/Math/interface/deltaR.h"
0031
0032 #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
0033 #include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
0034 //
0035 // class declaration
0036 //
0037
0038 class TiclDebugger : public edm::one::EDAnalyzer<edm::one::WatchRuns> {
0039 public:
0040   explicit TiclDebugger(const edm::ParameterSet&);
0041   ~TiclDebugger() override;
0042
0043   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0044
0045 private:
0046   void beginJob() override;
0047   void beginRun(const edm::Run&, const edm::EventSetup&) override;
0048   void analyze(const edm::Event&, const edm::EventSetup&) override;
0049   void endRun(edm::Run const& iEvent, edm::EventSetup const&) override {}
0050   void endJob() override;
0051
0052   const edm::InputTag trackstersMerge_;
0053   const edm::InputTag tracks_;
0054   const edm::InputTag caloParticles_;
0055   const edm::InputTag layerClusters_;
0056   hgcal::RecHitTools rhtools_;
0057   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometry_token_;
0058   edm::EDGetTokenT<std::vector<ticl::Trackster>> trackstersMergeToken_;
0059   edm::EDGetTokenT<std::vector<reco::Track>> tracksToken_;
0060   edm::EDGetTokenT<std::vector<CaloParticle>> caloParticlesToken_;
0061   edm::EDGetTokenT<std::vector<reco::CaloCluster>> layerClustersToken_;
0062 };
0063
0064 TiclDebugger::TiclDebugger(const edm::ParameterSet& iConfig)
0065     : trackstersMerge_(iConfig.getParameter<edm::InputTag>("trackstersMerge")),
0066       tracks_(iConfig.getParameter<edm::InputTag>("tracks")),
0067       caloParticles_(iConfig.getParameter<edm::InputTag>("caloParticles")),
0068       layerClusters_(iConfig.getParameter<edm::InputTag>("layerClusters")),
0069       caloGeometry_token_(esConsumes<CaloGeometry, CaloGeometryRecord, edm::Transition::BeginRun>()) {
0070   edm::ConsumesCollector&& iC = consumesCollector();
0071   trackstersMergeToken_ = iC.consumes<std::vector<ticl::Trackster>>(trackstersMerge_);
0072   tracksToken_ = iC.consumes<std::vector<reco::Track>>(tracks_);
0073   caloParticlesToken_ = iC.consumes<std::vector<CaloParticle>>(caloParticles_);
0074   layerClustersToken_ = iC.consumes<std::vector<reco::CaloCluster>>(layerClusters_);
0075 }
0076
0077 TiclDebugger::~TiclDebugger() {}
0078
0079 void TiclDebugger::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0080   static const char* particle_kind[] = {"gam", "e", "mu", "pi0", "h", "h0", "?", "!"};
0081   using namespace edm;
0082   using std::begin;
0083   using std::end;
0084   using std::iota;
0085   using std::sort;
0086
0087   edm::Handle<std::vector<ticl::Trackster>> trackstersMergeH;
0088
0089   iEvent.getByToken(trackstersMergeToken_, trackstersMergeH);
0090   auto const& tracksters = *trackstersMergeH.product();
0091   std::vector<int> sorted_tracksters_idx(tracksters.size());
0092   iota(begin(sorted_tracksters_idx), end(sorted_tracksters_idx), 0);
0093   sort(begin(sorted_tracksters_idx), end(sorted_tracksters_idx), [&tracksters](int i, int j) {
0094     return tracksters[i].raw_energy() > tracksters[j].raw_energy();
0095   });
0096
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<reco::Track>> tracksH;
0102   iEvent.getByToken(tracksToken_, tracksH);
0103   const auto& tracks = *tracksH.product();
0104
0105   edm::Handle<std::vector<CaloParticle>> caloParticlesH;
0106   iEvent.getByToken(caloParticlesToken_, caloParticlesH);
0107   auto const& caloParticles = *caloParticlesH.product();
0108   std::vector<std::pair<int, float>> bestCPMatches;
0109
0110   auto bestCaloParticleMatches = [&](const ticl::Trackster& t) -> void {
0111     bestCPMatches.clear();
0112     auto idx = 0;
0113     auto separation = 0.;
0114     for (auto const& cp : caloParticles) {
0115       separation = reco::deltaR2(t.barycenter(), cp.momentum());
0116       if (separation < 0.05) {
0117         bestCPMatches.push_back(std::make_pair(idx, separation));
0118       }
0119       ++idx;
0120     }
0121   };
0122
0123   std::stringstream prob_id_str;
0124
0125   for (auto const& t : sorted_tracksters_idx) {
0126     auto const& trackster = tracksters[t];
0127     auto const& probs = trackster.id_probabilities();
0128     // Sort probs in descending order
0129     std::vector<int> sorted_probs_idx(probs.size());
0130     iota(begin(sorted_probs_idx), end(sorted_probs_idx), 0);
0131     sort(begin(sorted_probs_idx), end(sorted_probs_idx), [&probs](int i, int j) { return probs[i] > probs[j]; });
0132     // Sort edges in ascending order
0133     std::vector<int> sorted_edges_idx(trackster.edges().size());
0134     iota(begin(sorted_edges_idx), end(sorted_edges_idx), 0);
0135     sort(begin(sorted_edges_idx), end(sorted_edges_idx), [&](int i, int j) {
0136       int layers = rhtools_.lastLayer();
0137       auto const& ed_i = trackster.edges()[i];
0138       auto const& ed_j = trackster.edges()[j];
0139       auto const& cl_i_in = layerClusters[ed_i[0]].hitsAndFractions()[0].first;
0140       auto const& cl_i_out = layerClusters[ed_i[1]].hitsAndFractions()[0].first;
0141       auto const& cl_j_in = layerClusters[ed_j[0]].hitsAndFractions()[0].first;
0142       auto const& cl_j_out = layerClusters[ed_j[1]].hitsAndFractions()[0].first;
0143       auto const layer_i_in = rhtools_.getLayerWithOffset(cl_i_in) + layers * ((rhtools_.zside(cl_i_in) + 1) >> 1) - 1;
0144       auto const layer_i_out =
0145           rhtools_.getLayerWithOffset(cl_i_out) + layers * ((rhtools_.zside(cl_i_out) + 1) >> 1) - 1;
0146       auto const layer_j_in = rhtools_.getLayerWithOffset(cl_j_in) + layers * ((rhtools_.zside(cl_j_in) + 1) >> 1) - 1;
0147       auto const layer_j_out =
0148           rhtools_.getLayerWithOffset(cl_j_out) + layers * ((rhtools_.zside(cl_j_out) + 1) >> 1) - 1;
0149       if (layer_i_in != layer_j_in)
0150         return layer_i_in < layer_j_in;
0151       else
0152         return layer_i_out < layer_j_out;
0153     });
0154
0155     for (auto p_idx : sorted_probs_idx) {
0156       prob_id_str << "(" << particle_kind[p_idx] << "):" << std::fixed << std::setprecision(4) << probs[p_idx] << " ";
0157     }
0158     LogVerbatim("TICLDebugger") << "\nTrksIdx: " << t << "\n bary: " << trackster.barycenter()
0159                                 << " baryEta: " << trackster.barycenter().eta()
0160                                 << " baryPhi: " << trackster.barycenter().phi()
0161                                 << "\n raw_energy: " << trackster.raw_energy()
0162                                 << " raw_em_energy: " << trackster.raw_em_energy()
0163                                 << "\n raw_pt: " << trackster.raw_pt() << " raw_em_pt: " << trackster.raw_em_pt()
0164                                 << "\n seedIdx: " << trackster.seedIndex() << "\n Probs: " << prob_id_str.str();
0165     prob_id_str.str("");
0166     prob_id_str.clear();
0167     LogVerbatim("TICLDebugger") << "\n time: " << trackster.time() << "+/-" << trackster.timeError() << std::endl
0168                                 << " vertices: " << trackster.vertices().size() << " average usage: "
0169                                 << std::accumulate(std::begin(trackster.vertex_multiplicity()),
0170                                                    std::end(trackster.vertex_multiplicity()),
0171                                                    0.) /
0172                                        trackster.vertex_multiplicity().size()
0173                                 << std::endl;
0174     LogVerbatim("TICLDebugger") << " link connections: " << trackster.edges().size() << std::endl;
0175     auto dumpLayerCluster = [&layerClusters](hgcal::RecHitTools const& rhtools, int cluster_idx) {
0176       auto const& cluster = layerClusters[cluster_idx];
0177       const auto firstHitDetId = cluster.hitsAndFractions()[0].first;
0178       int layers = rhtools.lastLayer();
0179       int lcLayerId =
0180           rhtools.getLayerWithOffset(firstHitDetId) + layers * ((rhtools.zside(firstHitDetId) + 1) >> 1) - 1;
0181
0182       LogVerbatim("TICLDebugger") << "Idx: " << cluster_idx << "(" << lcLayerId << ", "
0183                                   << cluster.hitsAndFractions().size() << ", " << cluster.position() << ") ";
0184     };
0185     for (auto link : sorted_edges_idx) {
0186       LogVerbatim("TICLDebugger") << "(" << trackster.edges()[link][0] << ", " << trackster.edges()[link][1] << ")  ";
0187       dumpLayerCluster(rhtools_, trackster.edges()[link][0]);
0188       dumpLayerCluster(rhtools_, trackster.edges()[link][1]);
0189       LogVerbatim("TICLDebugger") << std::endl;
0190     }
0191     if (trackster.seedID().id() != 0) {
0192       auto const& track = tracks[trackster.seedIndex()];
0193       LogVerbatim("TICLDebugger") << " Seeding Track:" << std::endl;
0194       LogVerbatim("TICLDebugger") << "   p: " << track.p() << " pt: " << track.pt() << " charge: " << track.charge()
0195                                   << " eta: " << track.eta() << " outerEta: " << track.outerEta()
0196                                   << " phi: " << track.phi() << " outerPhi: " << track.outerPhi() << std::endl;
0197     }
0198     bestCaloParticleMatches(trackster);
0199     if (!bestCPMatches.empty()) {
0200       LogVerbatim("TICLDebugger") << " Best CaloParticles Matches:" << std::endl;
0201       ;
0202       for (auto const& i : bestCPMatches) {
0203         auto const& cp = caloParticles[i.first];
0204         LogVerbatim("TICLDebugger") << "   " << i.first << "(" << i.second << "):" << cp.pdgId()
0205                                     << " simCl size:" << cp.simClusters().size() << " energy:" << cp.energy()
0206                                     << " pt:" << cp.pt() << " momentum:" << cp.momentum() << std::endl;
0207       }
0208       LogVerbatim("TICLDebugger") << std::endl;
0209     }
0210   }
0211 }
0212
0213 void TiclDebugger::beginRun(edm::Run const&, edm::EventSetup const& es) {
0214   const CaloGeometry& geom = es.getData(caloGeometry_token_);
0215   rhtools_.setGeometry(geom);
0216 }
0217
0218 void TiclDebugger::beginJob() {}
0219
0220 void TiclDebugger::endJob() {}
0221
0222 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0223 void TiclDebugger::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0224   edm::ParameterSetDescription desc;
0225   desc.add<edm::InputTag>("trackstersMerge", edm::InputTag("ticlTrackstersMerge"));
0226   desc.add<edm::InputTag>("tracks", edm::InputTag("generalTracks"));
0227   desc.add<edm::InputTag>("caloParticles", edm::InputTag("mix", "MergedCaloTruth"));
0228   desc.add<edm::InputTag>("layerClusters", edm::InputTag("hgcalMergeLayerClusters"));
0229   descriptions.add("ticlDebugger", desc);
0230 }
0231
0232 // define this as a plug-in
0233 DEFINE_FWK_MODULE(TiclDebugger);