Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-04-22 06:27:49 

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 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   // ----------member data ---------------------------
0078 };
0079 
0080 //
0081 // constants, enums and typedefs
0082 //
0083 
0084 //
0085 // static data member definitions
0086 //
0087 
0088 //
0089 // constructors and destructor
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 // member functions
0118 //
0119 
0120 // ------------ method called for each event  ------------
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   // Let's first fill in hits information
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();  // << cp ;
0203     double total_sim_energy = 0.;
0204     LogVerbatim("CaloParticleDebuggerCaloParticles") << "--> Overall simclusters in CP: " << cp.simClusters().size();
0205     totalSimCl_in_event += cp.simClusters().size();
0206     // All the next mess just to print the simClusters ordered
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;  // << simcl ;
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 // ------------ method called once each job just before starting event loop  ------------
0254 void CaloParticleDebugger::beginJob() {}
0255 
0256 // ------------ method called once each job just after ending the event loop  ------------
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   // Taken needed quantities from the EventSetup
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   //check if it's the new geometry
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   // loop over the collections
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         //no test numbering in new geometry
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         // skip simhits with bad barcodes or non-existant layers
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   }  // end of loop over InputTags
0335 }
0336 
0337 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
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 // define this as a plug-in
0355 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!