File indexing completed on 2025-04-22 06:27:49
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 genBarcodes_;
0062 edm::InputTag simVertices_;
0063 edm::InputTag trackingParticles_;
0064 edm::InputTag caloParticles_;
0065 edm::InputTag simClusters_;
0066 std::vector<edm::InputTag> collectionTags_;
0067 edm::EDGetTokenT<std::vector<SimTrack>> simTracksToken_;
0068 edm::EDGetTokenT<std::vector<reco::GenParticle>> genParticlesToken_;
0069 edm::EDGetTokenT<std::vector<int>> genBarcodesToken_;
0070 edm::EDGetTokenT<std::vector<SimVertex>> simVerticesToken_;
0071 edm::EDGetTokenT<std::vector<TrackingParticle>> trackingParticlesToken_;
0072 edm::EDGetTokenT<std::vector<CaloParticle>> caloParticlesToken_;
0073 edm::EDGetTokenT<std::vector<SimCluster>> simClustersToken_;
0074 std::vector<edm::EDGetTokenT<std::vector<PCaloHit>>> collectionTagsToken_;
0075 edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geomToken_;
0076 int geometryType_ = 0;
0077
0078 };
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091 CaloParticleDebugger::CaloParticleDebugger(const edm::ParameterSet& iConfig)
0092 : simTracks_(iConfig.getParameter<edm::InputTag>("simTracks")),
0093 genParticles_(iConfig.getParameter<edm::InputTag>("genParticles")),
0094 genBarcodes_(iConfig.getParameter<edm::InputTag>("genBarcodes")),
0095 simVertices_(iConfig.getParameter<edm::InputTag>("simVertices")),
0096 trackingParticles_(iConfig.getParameter<edm::InputTag>("trackingParticles")),
0097 caloParticles_(iConfig.getParameter<edm::InputTag>("caloParticles")),
0098 simClusters_(iConfig.getParameter<edm::InputTag>("simClusters")),
0099 collectionTags_(iConfig.getParameter<std::vector<edm::InputTag>>("collectionTags")) {
0100 edm::ConsumesCollector&& iC = consumesCollector();
0101 simTracksToken_ = iC.consumes<std::vector<SimTrack>>(simTracks_);
0102 genParticlesToken_ = iC.consumes<std::vector<reco::GenParticle>>(genParticles_);
0103 genBarcodesToken_ = iC.consumes<std::vector<int>>(genBarcodes_);
0104 simVerticesToken_ = iC.consumes<std::vector<SimVertex>>(simVertices_);
0105 trackingParticlesToken_ = iC.consumes<std::vector<TrackingParticle>>(trackingParticles_);
0106 caloParticlesToken_ = iC.consumes<std::vector<CaloParticle>>(caloParticles_);
0107 simClustersToken_ = iC.consumes<std::vector<SimCluster>>(simClusters_);
0108 for (auto const& collectionTag : collectionTags_) {
0109 collectionTagsToken_.push_back(iC.consumes<std::vector<PCaloHit>>(collectionTag));
0110 }
0111 geomToken_ = iC.esConsumes<CaloGeometry, CaloGeometryRecord>();
0112 }
0113
0114 CaloParticleDebugger::~CaloParticleDebugger() {}
0115
0116
0117
0118
0119
0120
0121 void CaloParticleDebugger::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0122 using namespace edm;
0123 using std::begin;
0124 using std::end;
0125
0126 edm::Handle<std::vector<SimTrack>> simTracksH;
0127 edm::Handle<std::vector<reco::GenParticle>> genParticlesH;
0128 edm::Handle<std::vector<int>> genBarcodesH;
0129 edm::Handle<std::vector<SimVertex>> simVerticesH;
0130 edm::Handle<std::vector<TrackingParticle>> trackingParticlesH;
0131 edm::Handle<std::vector<CaloParticle>> caloParticlesH;
0132 edm::Handle<std::vector<SimCluster>> simClustersH;
0133
0134 iEvent.getByToken(simTracksToken_, simTracksH);
0135 auto const& tracks = *simTracksH.product();
0136
0137 iEvent.getByToken(genParticlesToken_, genParticlesH);
0138 auto const& genParticles = *genParticlesH.product();
0139
0140 iEvent.getByToken(genBarcodesToken_, genBarcodesH);
0141 auto const& genBarcodes = *genBarcodesH.product();
0142
0143 iEvent.getByToken(simVerticesToken_, simVerticesH);
0144 auto const& vertices = *simVerticesH.product();
0145
0146 iEvent.getByToken(trackingParticlesToken_, trackingParticlesH);
0147 auto const& trackingpart = *trackingParticlesH.product();
0148
0149 iEvent.getByToken(caloParticlesToken_, caloParticlesH);
0150 auto const& calopart = *caloParticlesH.product();
0151
0152 iEvent.getByToken(simClustersToken_, simClustersH);
0153 auto const& simclusters = *simClustersH.product();
0154
0155
0156 std::map<int, float> detIdToTotalSimEnergy;
0157 fillSimHits(detIdToTotalSimEnergy, iEvent, iSetup);
0158
0159 std::map<int, int> trackid_to_track_index;
0160 LogVerbatim("CaloParticleDebuggerSimTracks") << "\n\n**Printing SimTracks information **";
0161 LogVerbatim("CaloParticleDebuggerSimTracks") << "IDX\tTrackId\tPDGID\tMOMENTUM(x,y,z,E)\tVertexIdx\tGenPartIdx";
0162 for (size_t i = 0; i < tracks.size(); ++i) {
0163 auto const& t = tracks[i];
0164 LogVerbatim("CaloParticleDebuggerSimTracks")
0165 << i << "\t" << t.trackId() << "\t" << t << " Crossed Boundary: " << t.crossedBoundary()
0166 << " Vtx: " << t.vertIndex() << " isFromBackScattering: " << t.isFromBackScattering()
0167 << " isPrimary: " << t.isPrimary() << " ParentID: " << t.getPrimaryID()
0168 << " Position Boundary: " << t.getPositionAtBoundary() << " Momentum Boundary: " << t.getMomentumAtBoundary()
0169 << " Momemtum Origin: " << t.momentum();
0170 trackid_to_track_index[t.trackId()] = i;
0171 }
0172
0173 LogVerbatim("CaloParticleDebuggerGenParticles") << "\n\n**Printing GenParticles information **";
0174 LogVerbatim("CaloParticleDebuggerGenParticles") << "BARCODE\tID\tPDGID\tMOMENTUM(x,y,z)\tVertex(x,y,z)";
0175 for (size_t i = 0; i < genParticles.size(); ++i) {
0176 auto const& gp = genParticles[i];
0177 LogVerbatim("CaloParticleDebuggerGenParticles")
0178 << genBarcodes[i] << "\t" << i << "\t" << gp.pdgId() << "\t" << gp.momentum() << "\t" << gp.vertex();
0179 }
0180
0181 LogVerbatim("CaloParticleDebuggerSimVertices") << "\n\n**Printing SimVertex information **";
0182 LogVerbatim("CaloParticleDebuggerSimVertices") << "IDX\tPOSITION(x,y,z)\tPARENT_INDEX\tVERTEX_ID";
0183 for (size_t i = 0; i < vertices.size(); ++i) {
0184 auto const& v = vertices[i];
0185 LogVerbatim("CaloParticleDebuggerSimVertices") << i << "\t" << v;
0186 }
0187
0188 LogVerbatim("CaloParticleDebuggerTrackingParticles") << "\n\n**Printing TrackingParticles information **";
0189 for (size_t i = 0; i < trackingpart.size(); ++i) {
0190 auto const& tp = trackingpart[i];
0191 LogVerbatim("CaloParticleDebuggerTrackingParticles") << i << "\t" << tp;
0192 }
0193
0194 LogVerbatim("CaloParticleDebuggerCaloParticles") << "\n\n**Printing CaloParticles information **";
0195 int totalSimCl_in_event = 0;
0196 std::set<int> simClusters_in_CaloParticles;
0197 for (size_t i = 0; i < calopart.size(); ++i) {
0198 auto const& cp = calopart[i];
0199 LogVerbatim("CaloParticleDebuggerCaloParticles") << "\n\n"
0200 << i << "\tType: " << cp.pdgId() << "\tEnergy: " << cp.energy()
0201 << "\tBarcode: " << cp.g4Tracks()[0].genpartIndex()
0202 << " G4_trackID: " << cp.g4Tracks()[0].trackId();
0203 double total_sim_energy = 0.;
0204 LogVerbatim("CaloParticleDebuggerCaloParticles") << "--> Overall simclusters in CP: " << cp.simClusters().size();
0205 totalSimCl_in_event += cp.simClusters().size();
0206
0207 auto const& simcs = cp.simClusters();
0208 for (size_t j = 0; j < simcs.size(); ++j) {
0209 LogVerbatim("CaloParticleDebuggerCaloParticles") << *(simcs[j]);
0210 simClusters_in_CaloParticles.insert(simcs[j]->g4Tracks()[0].trackId());
0211 }
0212
0213 for (auto const& sc : cp.simClusters()) {
0214 for (auto const& cl : sc->hits_and_fractions()) {
0215 total_sim_energy += detIdToTotalSimEnergy[cl.first] * cl.second;
0216 }
0217 }
0218 LogVerbatim("CaloParticleDebuggerCaloParticles")
0219 << "--> Overall SC energy (sum using sim energies): " << total_sim_energy;
0220 }
0221 LogVerbatim("CaloParticleDebuggerCaloParticles")
0222 << "--> Overall SCs attached to CPs: " << totalSimCl_in_event;
0223 LogVerbatim("CaloParticleDebuggerCaloParticles")
0224 << "--> Overall SCs in the event: " << simclusters.size();
0225 LogVerbatim("CaloParticleDebuggerCaloParticles")
0226 << "--> Orphan SCs in the event: " << simclusters.size() - totalSimCl_in_event;
0227
0228 LogVerbatim("CaloParticleDebuggerSimClusters") << "\n\n**Printing SimClusters information **";
0229 for (size_t i = 0; i < simclusters.size(); ++i) {
0230 auto const& simcl = simclusters[i];
0231 LogVerbatim("CaloParticleDebuggerSimClusters")
0232 << "\n\n"
0233 << i << "\tType: " << simcl.pdgId() << "\tEnergy: " << simcl.energy() << "\tKey: " << i;
0234 auto const& simtrack = simcl.g4Tracks()[0];
0235 LogVerbatim("CaloParticleDebuggerSimClusters") << "\n Primary GenParticleID: " << simtrack.getPrimaryID()
0236 << "\n Crossed Boundary: " << simtrack.crossedBoundary()
0237 << "\n Position Boundary: " << simtrack.getPositionAtBoundary()
0238 << "\n Momentum Boundary: " << simtrack.getMomentumAtBoundary();
0239 if (simClusters_in_CaloParticles.find(simcl.g4Tracks()[0].trackId()) == simClusters_in_CaloParticles.end()) {
0240 LogVerbatim("CaloParticleDebuggerSimClusters") << " Orphan SimCluster: " << simtrack.trackId();
0241 }
0242 double total_sim_energy = 0.;
0243 LogVerbatim("CaloParticleDebuggerSimClusters") << "--> Overall simclusters's size: " << simcl.numberOfRecHits();
0244 for (auto const& cl : simcl.hits_and_fractions()) {
0245 total_sim_energy += detIdToTotalSimEnergy[cl.first] * cl.second;
0246 }
0247 LogVerbatim("CaloParticleDebuggerSimClusters") << simcl;
0248 LogVerbatim("CaloParticleDebuggerSimClusters")
0249 << "--> Overall SimCluster energy (sum using sim energies): " << total_sim_energy;
0250 }
0251 }
0252
0253
0254 void CaloParticleDebugger::beginJob() {}
0255
0256
0257 void CaloParticleDebugger::endJob() {}
0258
0259 void CaloParticleDebugger::fillSimHits(std::map<int, float>& detIdToTotalSimEnergy,
0260 const edm::Event& iEvent,
0261 const edm::EventSetup& iSetup) {
0262
0263 const auto& geom = iSetup.getHandle(geomToken_);
0264 const HGCalGeometry *eegeom = nullptr, *fhgeom = nullptr, *bhgeomnew = nullptr;
0265 const HcalGeometry* bhgeom = nullptr;
0266 const HGCalDDDConstants* hgddd[3];
0267 const HGCalTopology* hgtopo[3];
0268 const HcalDDDRecConstants* hcddd = nullptr;
0269 eegeom =
0270 static_cast<const HGCalGeometry*>(geom->getSubdetectorGeometry(DetId::HGCalEE, ForwardSubdetector::ForwardEmpty));
0271
0272 if (eegeom) {
0273 geometryType_ = 1;
0274 fhgeom = static_cast<const HGCalGeometry*>(
0275 geom->getSubdetectorGeometry(DetId::HGCalHSi, ForwardSubdetector::ForwardEmpty));
0276 bhgeomnew = static_cast<const HGCalGeometry*>(
0277 geom->getSubdetectorGeometry(DetId::HGCalHSc, ForwardSubdetector::ForwardEmpty));
0278 } else {
0279 geometryType_ = 0;
0280 eegeom = static_cast<const HGCalGeometry*>(geom->getSubdetectorGeometry(DetId::Forward, HGCEE));
0281 fhgeom = static_cast<const HGCalGeometry*>(geom->getSubdetectorGeometry(DetId::Forward, HGCHEF));
0282 bhgeom = static_cast<const HcalGeometry*>(geom->getSubdetectorGeometry(DetId::Hcal, HcalEndcap));
0283 }
0284 hgtopo[0] = &(eegeom->topology());
0285 hgtopo[1] = &(fhgeom->topology());
0286 if (bhgeomnew)
0287 hgtopo[2] = &(bhgeomnew->topology());
0288 else
0289 hgtopo[2] = nullptr;
0290
0291 for (unsigned i = 0; i < 3; ++i) {
0292 if (hgtopo[i])
0293 hgddd[i] = &(hgtopo[i]->dddConstants());
0294 }
0295
0296 if (bhgeom)
0297 hcddd = bhgeom->topology().dddConstants();
0298
0299
0300 int token = 0;
0301 for (auto const& collectionTag : collectionTags_) {
0302 edm::Handle<std::vector<PCaloHit>> hSimHits;
0303 const bool isHcal = (collectionTag.instance().find("HcalHits") != std::string::npos);
0304 iEvent.getByToken(collectionTagsToken_[token++], hSimHits);
0305 for (auto const& simHit : *hSimHits) {
0306 DetId id(0);
0307 const uint32_t simId = simHit.id();
0308 if (geometryType_ == 1) {
0309
0310 id = simId;
0311 } else if (isHcal) {
0312 assert(hcddd);
0313 HcalDetId hid = HcalHitRelabeller::relabel(simId, hcddd);
0314 if (hid.subdet() == HcalEndcap)
0315 id = hid;
0316 } else {
0317 int subdet, layer, cell, sec, subsec, zp;
0318 HGCalTestNumbering::unpackHexagonIndex(simId, subdet, zp, layer, sec, subsec, cell);
0319 const HGCalDDDConstants* ddd = hgddd[subdet - 3];
0320 std::pair<int, int> recoLayerCell = ddd->simToReco(cell, layer, sec, hgtopo[subdet - 3]->detectorType());
0321 cell = recoLayerCell.first;
0322 layer = recoLayerCell.second;
0323
0324 if (layer == -1 || simHit.geantTrackId() == 0)
0325 continue;
0326 id = HGCalDetId((ForwardSubdetector)subdet, zp, layer, subsec, sec, cell);
0327 }
0328
0329 if (DetId(0) == id)
0330 continue;
0331
0332 detIdToTotalSimEnergy[id.rawId()] += simHit.energy();
0333 }
0334 }
0335 }
0336
0337
0338 void CaloParticleDebugger::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0339 edm::ParameterSetDescription desc;
0340 desc.add<edm::InputTag>("simTracks", edm::InputTag("g4SimHits"));
0341 desc.add<edm::InputTag>("genParticles", edm::InputTag("genParticles"));
0342 desc.add<edm::InputTag>("genBarcodes", edm::InputTag("genParticles"));
0343 desc.add<edm::InputTag>("simVertices", edm::InputTag("g4SimHits"));
0344 desc.add<edm::InputTag>("trackingParticles", edm::InputTag("mix", "MergedTrackTruth"));
0345 desc.add<edm::InputTag>("caloParticles", edm::InputTag("mix", "MergedCaloTruth"));
0346 desc.add<edm::InputTag>("simClusters", edm::InputTag("mix", "MergedCaloTruth"));
0347 desc.add<std::vector<edm::InputTag>>("collectionTags",
0348 {edm::InputTag("g4SimHits", "HGCHitsEE"),
0349 edm::InputTag("g4SimHits", "HGCHitsHEfront"),
0350 edm::InputTag("g4SimHits", "HcalHits")});
0351 descriptions.add("caloParticleDebugger", desc);
0352 }
0353
0354
0355 DEFINE_FWK_MODULE(CaloParticleDebugger);