Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:25:08 

0001 //
0002 // Original Author:  Marco Rovere
0003 //         Created:  Fri, 10 Nov 2017 14:39:18 GMT
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/ForwardDetId/interface/HGCalDetId.h"
0026 #include "SimDataFormats/Track/interface/SimTrack.h"
0027 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0028 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0029 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0030 #include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
0031 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0032 #include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
0033 
0034 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0035 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0036 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0037 #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h"
0038 #include "Geometry/CaloTopology/interface/HGCalTopology.h"
0039 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0040 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0041 #include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h"
0042 
0043 //
0044 // class declaration
0045 //
0046 
0047 class CaloParticleDebugger : public edm::one::EDAnalyzer<> {
0048 public:
0049   explicit CaloParticleDebugger(const edm::ParameterSet&);
0050   ~CaloParticleDebugger() override;
0051 
0052   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0053 
0054 private:
0055   void beginJob() override;
0056   void analyze(const edm::Event&, const edm::EventSetup&) override;
0057   void endJob() override;
0058   void fillSimHits(std::map<int, float>&, const edm::Event&, const edm::EventSetup&);
0059   edm::InputTag simTracks_;
0060   edm::InputTag genParticles_;
0061   edm::InputTag simVertices_;
0062   edm::InputTag trackingParticles_;
0063   edm::InputTag caloParticles_;
0064   edm::InputTag simClusters_;
0065   std::vector<edm::InputTag> collectionTags_;
0066   edm::EDGetTokenT<std::vector<SimTrack>> simTracksToken_;
0067   edm::EDGetTokenT<std::vector<reco::GenParticle>> genParticlesToken_;
0068   edm::EDGetTokenT<std::vector<SimVertex>> simVerticesToken_;
0069   edm::EDGetTokenT<std::vector<TrackingParticle>> trackingParticlesToken_;
0070   edm::EDGetTokenT<std::vector<CaloParticle>> caloParticlesToken_;
0071   edm::EDGetTokenT<std::vector<SimCluster>> simClustersToken_;
0072   std::vector<edm::EDGetTokenT<std::vector<PCaloHit>>> collectionTagsToken_;
0073   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geomToken_;
0074   int geometryType_ = 0;
0075   // ----------member data ---------------------------
0076 };
0077 
0078 //
0079 // constants, enums and typedefs
0080 //
0081 
0082 //
0083 // static data member definitions
0084 //
0085 
0086 //
0087 // constructors and destructor
0088 //
0089 CaloParticleDebugger::CaloParticleDebugger(const edm::ParameterSet& iConfig)
0090     : simTracks_(iConfig.getParameter<edm::InputTag>("simTracks")),
0091       genParticles_(iConfig.getParameter<edm::InputTag>("genParticles")),
0092       simVertices_(iConfig.getParameter<edm::InputTag>("simVertices")),
0093       trackingParticles_(iConfig.getParameter<edm::InputTag>("trackingParticles")),
0094       caloParticles_(iConfig.getParameter<edm::InputTag>("caloParticles")),
0095       simClusters_(iConfig.getParameter<edm::InputTag>("simClusters")),
0096       collectionTags_(iConfig.getParameter<std::vector<edm::InputTag>>("collectionTags")) {
0097   edm::ConsumesCollector&& iC = consumesCollector();
0098   simTracksToken_ = iC.consumes<std::vector<SimTrack>>(simTracks_);
0099   genParticlesToken_ = iC.consumes<std::vector<reco::GenParticle>>(genParticles_);
0100   simVerticesToken_ = iC.consumes<std::vector<SimVertex>>(simVertices_);
0101   trackingParticlesToken_ = iC.consumes<std::vector<TrackingParticle>>(trackingParticles_);
0102   caloParticlesToken_ = iC.consumes<std::vector<CaloParticle>>(caloParticles_);
0103   simClustersToken_ = iC.consumes<std::vector<SimCluster>>(simClusters_);
0104   for (auto const& collectionTag : collectionTags_) {
0105     collectionTagsToken_.push_back(iC.consumes<std::vector<PCaloHit>>(collectionTag));
0106   }
0107   geomToken_ = iC.esConsumes<CaloGeometry, CaloGeometryRecord>();
0108 }
0109 
0110 CaloParticleDebugger::~CaloParticleDebugger() {}
0111 
0112 //
0113 // member functions
0114 //
0115 
0116 // ------------ method called for each event  ------------
0117 void CaloParticleDebugger::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0118   using namespace edm;
0119   using std::begin;
0120   using std::end;
0121   using std::iota;
0122   using std::sort;
0123 
0124   edm::Handle<std::vector<SimTrack>> simTracksH;
0125   edm::Handle<std::vector<reco::GenParticle>> genParticlesH;
0126   edm::Handle<std::vector<SimVertex>> simVerticesH;
0127   edm::Handle<std::vector<TrackingParticle>> trackingParticlesH;
0128   edm::Handle<std::vector<CaloParticle>> caloParticlesH;
0129   edm::Handle<std::vector<SimCluster>> simClustersH;
0130 
0131   iEvent.getByToken(simTracksToken_, simTracksH);
0132   auto const& tracks = *simTracksH.product();
0133   std::vector<int> sorted_tracks_idx(tracks.size());
0134   iota(begin(sorted_tracks_idx), end(sorted_tracks_idx), 0);
0135   sort(begin(sorted_tracks_idx), end(sorted_tracks_idx), [&tracks](int i, int j) {
0136     return tracks[i].momentum().eta() < tracks[j].momentum().eta();
0137   });
0138 
0139   iEvent.getByToken(genParticlesToken_, genParticlesH);
0140   auto const& genParticles = *genParticlesH.product();
0141   std::vector<int> sorted_genParticles_idx(genParticles.size());
0142   iota(begin(sorted_genParticles_idx), end(sorted_genParticles_idx), 0);
0143   sort(begin(sorted_genParticles_idx), end(sorted_genParticles_idx), [&genParticles](int i, int j) {
0144     return genParticles[i].momentum().eta() < genParticles[j].momentum().eta();
0145   });
0146 
0147   iEvent.getByToken(simVerticesToken_, simVerticesH);
0148   auto const& vertices = *simVerticesH.product();
0149   std::vector<int> sorted_vertices_idx(vertices.size());
0150   iota(begin(sorted_vertices_idx), end(sorted_vertices_idx), 0);
0151   sort(begin(sorted_vertices_idx), end(sorted_vertices_idx), [&vertices](int i, int j) {
0152     return vertices[i].vertexId() < vertices[j].vertexId();
0153   });
0154 
0155   iEvent.getByToken(trackingParticlesToken_, trackingParticlesH);
0156   auto const& trackingpart = *trackingParticlesH.product();
0157   std::vector<int> sorted_tp_idx(trackingpart.size());
0158   iota(begin(sorted_tp_idx), end(sorted_tp_idx), 0);
0159   sort(begin(sorted_tp_idx), end(sorted_tp_idx), [&trackingpart](int i, int j) {
0160     return trackingpart[i].eta() < trackingpart[j].eta();
0161   });
0162 
0163   iEvent.getByToken(caloParticlesToken_, caloParticlesH);
0164   auto const& calopart = *caloParticlesH.product();
0165   std::vector<int> sorted_cp_idx(calopart.size());
0166   iota(begin(sorted_cp_idx), end(sorted_cp_idx), 0);
0167   sort(begin(sorted_cp_idx), end(sorted_cp_idx), [&calopart](int i, int j) {
0168     return calopart[i].eta() < calopart[j].eta();
0169   });
0170 
0171   iEvent.getByToken(simClustersToken_, simClustersH);
0172   auto const& simclusters = *simClustersH.product();
0173   std::vector<int> sorted_simcl_idx(simclusters.size());
0174   iota(begin(sorted_simcl_idx), end(sorted_simcl_idx), 0);
0175   sort(begin(sorted_simcl_idx), end(sorted_simcl_idx), [&simclusters](int i, int j) {
0176     return simclusters[i].eta() < simclusters[j].eta();
0177   });
0178 
0179   // Let's first fill in hits information
0180   std::map<int, float> detIdToTotalSimEnergy;
0181   fillSimHits(detIdToTotalSimEnergy, iEvent, iSetup);
0182 
0183   int idx = 0;
0184 
0185   std::map<int, int> trackid_to_track_index;
0186   LogVerbatim("CaloParticleDebuggerSimTracks") << "\n\n**Printing SimTracks information **";
0187   LogVerbatim("CaloParticleDebuggerSimTracks") << "IDX\tTrackId\tPDGID\tMOMENTUM(x,y,z,E)\tVertexIdx\tGenPartIdx";
0188   for (auto i : sorted_tracks_idx) {
0189     auto const& t = tracks[i];
0190     LogVerbatim("CaloParticleDebuggerSimTracks") << idx << "\t" << t.trackId() << "\t" << t;
0191     LogVerbatim("CaloParticleDebuggerSimTracks")
0192         << "Crossed  Boundary: " << t.crossedBoundary() << "  Position Boundary: " << t.getPositionAtBoundary()
0193         << "  Momentum Boundary: " << t.getMomentumAtBoundary() << "  Vtx:               " << t.vertIndex()
0194         << "  Momemtum Origin:   " << t.momentum();
0195     trackid_to_track_index[t.trackId()] = idx;
0196     idx++;
0197   }
0198 
0199   LogVerbatim("CaloParticleDebuggerGenParticles") << "\n\n**Printing GenParticles information **";
0200   LogVerbatim("CaloParticleDebuggerGenParticles") << "IDX\tPDGID\tMOMENTUM(x,y,z)\tVertex(x,y,z)";
0201   for (auto i : sorted_genParticles_idx) {
0202     auto const& gp = genParticles[i];
0203     LogVerbatim("CaloParticleDebuggerGenParticles")
0204         << i << "\t" << gp.pdgId() << "\t" << gp.momentum() << "\t" << gp.vertex();
0205   }
0206 
0207   LogVerbatim("CaloParticleDebuggerSimVertices") << "\n\n**Printing SimVertex information **";
0208   LogVerbatim("CaloParticleDebuggerSimVertices") << "IDX\tPOSITION(x,y,z)\tPARENT_INDEX\tVERTEX_ID";
0209   for (auto i : sorted_vertices_idx) {
0210     auto const& v = vertices[i];
0211     LogVerbatim("CaloParticleDebuggerSimVertices") << i << "\t" << v;
0212   }
0213 
0214   LogVerbatim("CaloParticleDebuggerTrackingParticles") << "\n\n**Printing TrackingParticles information **";
0215   for (auto i : sorted_tp_idx) {
0216     auto const& tp = trackingpart[i];
0217     LogVerbatim("CaloParticleDebuggerTrackingParticles") << i << "\t" << tp;
0218   }
0219 
0220   LogVerbatim("CaloParticleDebuggerCaloParticles") << "\n\n**Printing CaloParticles information **";
0221   idx = 0;
0222   for (auto i : sorted_cp_idx) {
0223     auto const& cp = calopart[i];
0224     LogVerbatim("CaloParticleDebuggerCaloParticles")
0225         << "\n\n"
0226         << idx++ << " |Eta|: " << std::abs(cp.momentum().eta()) << "\tType: " << cp.pdgId()
0227         << "\tEnergy: " << cp.energy() << "\tIdx: " << cp.g4Tracks()[0].trackId();  // << cp ;
0228     double total_sim_energy = 0.;
0229     double total_cp_energy = 0.;
0230     LogVerbatim("CaloParticleDebuggerCaloParticles") << "--> Overall simclusters in CP: " << cp.simClusters().size();
0231     // All the next mess just to print the simClusters ordered
0232     auto const& simcs = cp.simClusters();
0233     std::vector<int> sorted_sc_idx(simcs.size());
0234     iota(begin(sorted_sc_idx), end(sorted_sc_idx), 0);
0235     sort(begin(sorted_sc_idx), end(sorted_sc_idx), [&simcs](int i, int j) {
0236       return simcs[i]->momentum().eta() < simcs[j]->momentum().eta();
0237     });
0238     for (auto i : sorted_sc_idx) {
0239       LogVerbatim("CaloParticleDebuggerCaloParticles") << *(simcs[i]);
0240     }
0241 
0242     for (auto const& sc : cp.simClusters()) {
0243       for (auto const& cl : sc->hits_and_fractions()) {
0244         total_sim_energy += detIdToTotalSimEnergy[cl.first] * cl.second;
0245         total_cp_energy += cp.energy() * cl.second;
0246       }
0247     }
0248     LogVerbatim("CaloParticleDebuggerCaloParticles")
0249         << "--> Overall SC energy (sum using sim energies): " << total_sim_energy;
0250     LogVerbatim("CaloParticleDebuggerCaloParticles")
0251         << "--> Overall SC energy (sum using CaloP energies): " << total_cp_energy;
0252   }
0253 
0254   idx = 0;
0255   LogVerbatim("CaloParticleDebuggerSimClusters") << "\n\n**Printing SimClusters information **";
0256   for (auto i : sorted_simcl_idx) {
0257     auto const& simcl = simclusters[i];
0258     LogVerbatim("CaloParticleDebuggerSimClusters")
0259         << "\n\n"
0260         << idx++ << " |Eta|: " << std::abs(simcl.momentum().eta()) << "\tType: " << simcl.pdgId()
0261         << "\tEnergy: " << simcl.energy() << "\tKey: " << i;  // << simcl ;
0262     auto const& simtrack = simcl.g4Tracks()[0];
0263     LogVerbatim("CaloParticleDebuggerSimClusters") << "  Crossed  Boundary: " << simtrack.crossedBoundary()
0264                                                    << "  Position Boundary: " << simtrack.getPositionAtBoundary()
0265                                                    << "  Momentum Boundary: " << simtrack.getMomentumAtBoundary();
0266     double total_sim_energy = 0.;
0267     LogVerbatim("CaloParticleDebuggerSimClusters") << "--> Overall simclusters's size: " << simcl.numberOfRecHits();
0268     for (auto const& cl : simcl.hits_and_fractions()) {
0269       total_sim_energy += detIdToTotalSimEnergy[cl.first] * cl.second;
0270     }
0271     LogVerbatim("CaloParticleDebuggerSimClusters") << simcl;
0272     LogVerbatim("CaloParticleDebuggerSimClusters")
0273         << "--> Overall SimCluster energy (sum using sim energies): " << total_sim_energy;
0274   }
0275 }
0276 
0277 // ------------ method called once each job just before starting event loop  ------------
0278 void CaloParticleDebugger::beginJob() {}
0279 
0280 // ------------ method called once each job just after ending the event loop  ------------
0281 void CaloParticleDebugger::endJob() {}
0282 
0283 void CaloParticleDebugger::fillSimHits(std::map<int, float>& detIdToTotalSimEnergy,
0284                                        const edm::Event& iEvent,
0285                                        const edm::EventSetup& iSetup) {
0286   // Taken needed quantities from the EventSetup
0287   const auto& geom = iSetup.getHandle(geomToken_);
0288   const HGCalGeometry *eegeom = nullptr, *fhgeom = nullptr, *bhgeomnew = nullptr;
0289   const HcalGeometry* bhgeom = nullptr;
0290   const HGCalDDDConstants* hgddd[3];
0291   const HGCalTopology* hgtopo[3];
0292   const HcalDDDRecConstants* hcddd = nullptr;
0293   eegeom =
0294       static_cast<const HGCalGeometry*>(geom->getSubdetectorGeometry(DetId::HGCalEE, ForwardSubdetector::ForwardEmpty));
0295   //check if it's the new geometry
0296   if (eegeom) {
0297     geometryType_ = 1;
0298     fhgeom = static_cast<const HGCalGeometry*>(
0299         geom->getSubdetectorGeometry(DetId::HGCalHSi, ForwardSubdetector::ForwardEmpty));
0300     bhgeomnew = static_cast<const HGCalGeometry*>(
0301         geom->getSubdetectorGeometry(DetId::HGCalHSc, ForwardSubdetector::ForwardEmpty));
0302   } else {
0303     geometryType_ = 0;
0304     eegeom = static_cast<const HGCalGeometry*>(geom->getSubdetectorGeometry(DetId::Forward, HGCEE));
0305     fhgeom = static_cast<const HGCalGeometry*>(geom->getSubdetectorGeometry(DetId::Forward, HGCHEF));
0306     bhgeom = static_cast<const HcalGeometry*>(geom->getSubdetectorGeometry(DetId::Hcal, HcalEndcap));
0307   }
0308   hgtopo[0] = &(eegeom->topology());
0309   hgtopo[1] = &(fhgeom->topology());
0310   if (bhgeomnew)
0311     hgtopo[2] = &(bhgeomnew->topology());
0312 
0313   for (unsigned i = 0; i < 3; ++i) {
0314     if (hgtopo[i])
0315       hgddd[i] = &(hgtopo[i]->dddConstants());
0316   }
0317 
0318   if (bhgeom)
0319     hcddd = bhgeom->topology().dddConstants();
0320 
0321   // loop over the collections
0322   int token = 0;
0323   for (auto const& collectionTag : collectionTags_) {
0324     edm::Handle<std::vector<PCaloHit>> hSimHits;
0325     const bool isHcal = (collectionTag.instance().find("HcalHits") != std::string::npos);
0326     iEvent.getByToken(collectionTagsToken_[token++], hSimHits);
0327     for (auto const& simHit : *hSimHits) {
0328       DetId id(0);
0329       const uint32_t simId = simHit.id();
0330       if (geometryType_ == 1) {
0331         //no test numbering in new geometry
0332         id = simId;
0333       } else if (isHcal) {
0334         assert(hcddd);
0335         HcalDetId hid = HcalHitRelabeller::relabel(simId, hcddd);
0336         if (hid.subdet() == HcalEndcap)
0337           id = hid;
0338       } else {
0339         int subdet, layer, cell, sec, subsec, zp;
0340         HGCalTestNumbering::unpackHexagonIndex(simId, subdet, zp, layer, sec, subsec, cell);
0341         const HGCalDDDConstants* ddd = hgddd[subdet - 3];
0342         std::pair<int, int> recoLayerCell = ddd->simToReco(cell, layer, sec, hgtopo[subdet - 3]->detectorType());
0343         cell = recoLayerCell.first;
0344         layer = recoLayerCell.second;
0345         // skip simhits with bad barcodes or non-existant layers
0346         if (layer == -1 || simHit.geantTrackId() == 0)
0347           continue;
0348         id = HGCalDetId((ForwardSubdetector)subdet, zp, layer, subsec, sec, cell);
0349       }
0350 
0351       if (DetId(0) == id)
0352         continue;
0353 
0354       detIdToTotalSimEnergy[id.rawId()] += simHit.energy();
0355     }
0356   }  // end of loop over InputTags
0357 }
0358 
0359 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0360 void CaloParticleDebugger::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0361   edm::ParameterSetDescription desc;
0362   desc.add<edm::InputTag>("simTracks", edm::InputTag("g4SimHits"));
0363   desc.add<edm::InputTag>("genParticles", edm::InputTag("genParticles"));
0364   desc.add<edm::InputTag>("simVertices", edm::InputTag("g4SimHits"));
0365   desc.add<edm::InputTag>("trackingParticles", edm::InputTag("mix", "MergedTrackTruth"));
0366   desc.add<edm::InputTag>("caloParticles", edm::InputTag("mix", "MergedCaloTruth"));
0367   desc.add<edm::InputTag>("simClusters", edm::InputTag("mix", "MergedCaloTruth"));
0368   desc.add<std::vector<edm::InputTag>>("collectionTags",
0369                                        {edm::InputTag("g4SimHits", "HGCHitsEE"),
0370                                         edm::InputTag("g4SimHits", "HGCHitsHEfront"),
0371                                         edm::InputTag("g4SimHits", "HcalHits")});
0372   descriptions.add("caloParticleDebugger", desc);
0373 }
0374 
0375 // define this as a plug-in
0376 DEFINE_FWK_MODULE(CaloParticleDebugger);
[ 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!