File indexing completed on 2023-03-17 11:25:08
0001
0002
0003
0004
0005
0006
0007
0008 #include <memory>
0009 #include <iostream>
0010 #include <numeric>
0011 #include <algorithm>
0012
0013
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
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
0076 };
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
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
0114
0115
0116
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
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();
0228 double total_sim_energy = 0.;
0229 double total_cp_energy = 0.;
0230 LogVerbatim("CaloParticleDebuggerCaloParticles") << "--> Overall simclusters in CP: " << cp.simClusters().size();
0231
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;
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
0278 void CaloParticleDebugger::beginJob() {}
0279
0280
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
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
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
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
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
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 }
0357 }
0358
0359
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
0376 DEFINE_FWK_MODULE(CaloParticleDebugger);