CMS 3D CMS Logo

CaloParticleDebugger.cc
Go to the documentation of this file.
1 //
2 // Original Author: Marco Rovere
3 // Created: Fri, 10 Nov 2017 14:39:18 GMT
4 //
5 //
6 //
7 // system include files
8 #include <memory>
9 #include <iostream>
10 #include <numeric>
11 #include <algorithm>
12 
13 // user include files
17 
22 
24 
33 
42 
43 //
44 // class declaration
45 //
46 
48 public:
49  explicit CaloParticleDebugger(const edm::ParameterSet&);
50  ~CaloParticleDebugger() override;
51 
52  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
53 
54 private:
55  void beginJob() override;
56  void analyze(const edm::Event&, const edm::EventSetup&) override;
57  void endJob() override;
58  void fillSimHits(std::map<int, float>&, const edm::Event&, const edm::EventSetup&);
65  std::vector<edm::InputTag> collectionTags_;
72  std::vector<edm::EDGetTokenT<std::vector<PCaloHit>>> collectionTagsToken_;
74  int geometryType_ = 0;
75  // ----------member data ---------------------------
76 };
77 
78 //
79 // constants, enums and typedefs
80 //
81 
82 //
83 // static data member definitions
84 //
85 
86 //
87 // constructors and destructor
88 //
90  : simTracks_(iConfig.getParameter<edm::InputTag>("simTracks")),
91  genParticles_(iConfig.getParameter<edm::InputTag>("genParticles")),
92  simVertices_(iConfig.getParameter<edm::InputTag>("simVertices")),
93  trackingParticles_(iConfig.getParameter<edm::InputTag>("trackingParticles")),
94  caloParticles_(iConfig.getParameter<edm::InputTag>("caloParticles")),
95  simClusters_(iConfig.getParameter<edm::InputTag>("simClusters")),
96  collectionTags_(iConfig.getParameter<std::vector<edm::InputTag>>("collectionTags")) {
98  simTracksToken_ = iC.consumes<std::vector<SimTrack>>(simTracks_);
99  genParticlesToken_ = iC.consumes<std::vector<reco::GenParticle>>(genParticles_);
100  simVerticesToken_ = iC.consumes<std::vector<SimVertex>>(simVertices_);
101  trackingParticlesToken_ = iC.consumes<std::vector<TrackingParticle>>(trackingParticles_);
102  caloParticlesToken_ = iC.consumes<std::vector<CaloParticle>>(caloParticles_);
103  simClustersToken_ = iC.consumes<std::vector<SimCluster>>(simClusters_);
104  for (auto const& collectionTag : collectionTags_) {
105  collectionTagsToken_.push_back(iC.consumes<std::vector<PCaloHit>>(collectionTag));
106  }
108 }
109 
111 
112 //
113 // member functions
114 //
115 
116 // ------------ method called for each event ------------
118  using namespace edm;
119  using std::begin;
120  using std::end;
121  using std::iota;
122  using std::sort;
123 
127  edm::Handle<std::vector<TrackingParticle>> trackingParticlesH;
130 
131  iEvent.getByToken(simTracksToken_, simTracksH);
132  auto const& tracks = *simTracksH.product();
133  std::vector<int> sorted_tracks_idx(tracks.size());
134  iota(begin(sorted_tracks_idx), end(sorted_tracks_idx), 0);
135  sort(begin(sorted_tracks_idx), end(sorted_tracks_idx), [&tracks](int i, int j) {
136  return tracks[i].momentum().eta() < tracks[j].momentum().eta();
137  });
138 
140  auto const& genParticles = *genParticlesH.product();
141  std::vector<int> sorted_genParticles_idx(genParticles.size());
142  iota(begin(sorted_genParticles_idx), end(sorted_genParticles_idx), 0);
143  sort(begin(sorted_genParticles_idx), end(sorted_genParticles_idx), [&genParticles](int i, int j) {
144  return genParticles[i].momentum().eta() < genParticles[j].momentum().eta();
145  });
146 
147  iEvent.getByToken(simVerticesToken_, simVerticesH);
148  auto const& vertices = *simVerticesH.product();
149  std::vector<int> sorted_vertices_idx(vertices.size());
150  iota(begin(sorted_vertices_idx), end(sorted_vertices_idx), 0);
151  sort(begin(sorted_vertices_idx), end(sorted_vertices_idx), [&vertices](int i, int j) {
152  return vertices[i].vertexId() < vertices[j].vertexId();
153  });
154 
155  iEvent.getByToken(trackingParticlesToken_, trackingParticlesH);
156  auto const& trackingpart = *trackingParticlesH.product();
157  std::vector<int> sorted_tp_idx(trackingpart.size());
158  iota(begin(sorted_tp_idx), end(sorted_tp_idx), 0);
159  sort(begin(sorted_tp_idx), end(sorted_tp_idx), [&trackingpart](int i, int j) {
160  return trackingpart[i].eta() < trackingpart[j].eta();
161  });
162 
163  iEvent.getByToken(caloParticlesToken_, caloParticlesH);
164  auto const& calopart = *caloParticlesH.product();
165  std::vector<int> sorted_cp_idx(calopart.size());
166  iota(begin(sorted_cp_idx), end(sorted_cp_idx), 0);
167  sort(begin(sorted_cp_idx), end(sorted_cp_idx), [&calopart](int i, int j) {
168  return calopart[i].eta() < calopart[j].eta();
169  });
170 
171  iEvent.getByToken(simClustersToken_, simClustersH);
172  auto const& simclusters = *simClustersH.product();
173  std::vector<int> sorted_simcl_idx(simclusters.size());
174  iota(begin(sorted_simcl_idx), end(sorted_simcl_idx), 0);
175  sort(begin(sorted_simcl_idx), end(sorted_simcl_idx), [&simclusters](int i, int j) {
176  return simclusters[i].eta() < simclusters[j].eta();
177  });
178 
179  // Let's first fill in hits information
180  std::map<int, float> detIdToTotalSimEnergy;
181  fillSimHits(detIdToTotalSimEnergy, iEvent, iSetup);
182 
183  int idx = 0;
184 
185  std::map<int, int> trackid_to_track_index;
186  LogVerbatim("CaloParticleDebuggerSimTracks") << "\n\n**Printing SimTracks information **";
187  LogVerbatim("CaloParticleDebuggerSimTracks") << "IDX\tTrackId\tPDGID\tMOMENTUM(x,y,z,E)\tVertexIdx\tGenPartIdx";
188  for (auto i : sorted_tracks_idx) {
189  auto const& t = tracks[i];
190  LogVerbatim("CaloParticleDebuggerSimTracks") << idx << "\t" << t.trackId() << "\t" << t;
191  LogVerbatim("CaloParticleDebuggerSimTracks")
192  << "Crossed Boundary: " << t.crossedBoundary() << " Position Boundary: " << t.getPositionAtBoundary()
193  << " Momentum Boundary: " << t.getMomentumAtBoundary() << " Vtx: " << t.vertIndex()
194  << " Momemtum Origin: " << t.momentum();
195  trackid_to_track_index[t.trackId()] = idx;
196  idx++;
197  }
198 
199  LogVerbatim("CaloParticleDebuggerGenParticles") << "\n\n**Printing GenParticles information **";
200  LogVerbatim("CaloParticleDebuggerGenParticles") << "IDX\tPDGID\tMOMENTUM(x,y,z)\tVertex(x,y,z)";
201  for (auto i : sorted_genParticles_idx) {
202  auto const& gp = genParticles[i];
203  LogVerbatim("CaloParticleDebuggerGenParticles")
204  << i << "\t" << gp.pdgId() << "\t" << gp.momentum() << "\t" << gp.vertex();
205  }
206 
207  LogVerbatim("CaloParticleDebuggerSimVertices") << "\n\n**Printing SimVertex information **";
208  LogVerbatim("CaloParticleDebuggerSimVertices") << "IDX\tPOSITION(x,y,z)\tPARENT_INDEX\tVERTEX_ID";
209  for (auto i : sorted_vertices_idx) {
210  auto const& v = vertices[i];
211  LogVerbatim("CaloParticleDebuggerSimVertices") << i << "\t" << v;
212  }
213 
214  LogVerbatim("CaloParticleDebuggerTrackingParticles") << "\n\n**Printing TrackingParticles information **";
215  for (auto i : sorted_tp_idx) {
216  auto const& tp = trackingpart[i];
217  LogVerbatim("CaloParticleDebuggerTrackingParticles") << i << "\t" << tp;
218  }
219 
220  LogVerbatim("CaloParticleDebuggerCaloParticles") << "\n\n**Printing CaloParticles information **";
221  idx = 0;
222  for (auto i : sorted_cp_idx) {
223  auto const& cp = calopart[i];
224  LogVerbatim("CaloParticleDebuggerCaloParticles")
225  << "\n\n"
226  << idx++ << " |Eta|: " << std::abs(cp.momentum().eta()) << "\tType: " << cp.pdgId()
227  << "\tEnergy: " << cp.energy() << "\tIdx: " << cp.g4Tracks()[0].trackId(); // << cp ;
228  double total_sim_energy = 0.;
229  double total_cp_energy = 0.;
230  LogVerbatim("CaloParticleDebuggerCaloParticles") << "--> Overall simclusters in CP: " << cp.simClusters().size();
231  // All the next mess just to print the simClusters ordered
232  auto const& simcs = cp.simClusters();
233  std::vector<int> sorted_sc_idx(simcs.size());
234  iota(begin(sorted_sc_idx), end(sorted_sc_idx), 0);
235  sort(begin(sorted_sc_idx), end(sorted_sc_idx), [&simcs](int i, int j) {
236  return simcs[i]->momentum().eta() < simcs[j]->momentum().eta();
237  });
238  for (auto i : sorted_sc_idx) {
239  LogVerbatim("CaloParticleDebuggerCaloParticles") << *(simcs[i]);
240  }
241 
242  for (auto const& sc : cp.simClusters()) {
243  for (auto const& cl : sc->hits_and_fractions()) {
244  total_sim_energy += detIdToTotalSimEnergy[cl.first] * cl.second;
245  total_cp_energy += cp.energy() * cl.second;
246  }
247  }
248  LogVerbatim("CaloParticleDebuggerCaloParticles")
249  << "--> Overall SC energy (sum using sim energies): " << total_sim_energy;
250  LogVerbatim("CaloParticleDebuggerCaloParticles")
251  << "--> Overall SC energy (sum using CaloP energies): " << total_cp_energy;
252  }
253 
254  idx = 0;
255  LogVerbatim("CaloParticleDebuggerSimClusters") << "\n\n**Printing SimClusters information **";
256  for (auto i : sorted_simcl_idx) {
257  auto const& simcl = simclusters[i];
258  LogVerbatim("CaloParticleDebuggerSimClusters")
259  << "\n\n"
260  << idx++ << " |Eta|: " << std::abs(simcl.momentum().eta()) << "\tType: " << simcl.pdgId()
261  << "\tEnergy: " << simcl.energy() << "\tKey: " << i; // << simcl ;
262  auto const& simtrack = simcl.g4Tracks()[0];
263  LogVerbatim("CaloParticleDebuggerSimClusters") << " Crossed Boundary: " << simtrack.crossedBoundary()
264  << " Position Boundary: " << simtrack.getPositionAtBoundary()
265  << " Momentum Boundary: " << simtrack.getMomentumAtBoundary();
266  double total_sim_energy = 0.;
267  LogVerbatim("CaloParticleDebuggerSimClusters") << "--> Overall simclusters's size: " << simcl.numberOfRecHits();
268  for (auto const& cl : simcl.hits_and_fractions()) {
269  total_sim_energy += detIdToTotalSimEnergy[cl.first] * cl.second;
270  }
271  LogVerbatim("CaloParticleDebuggerSimClusters") << simcl;
272  LogVerbatim("CaloParticleDebuggerSimClusters")
273  << "--> Overall SimCluster energy (sum using sim energies): " << total_sim_energy;
274  }
275 }
276 
277 // ------------ method called once each job just before starting event loop ------------
279 
280 // ------------ method called once each job just after ending the event loop ------------
282 
283 void CaloParticleDebugger::fillSimHits(std::map<int, float>& detIdToTotalSimEnergy,
284  const edm::Event& iEvent,
285  const edm::EventSetup& iSetup) {
286  // Taken needed quantities from the EventSetup
287  const auto& geom = iSetup.getHandle(geomToken_);
288  const HGCalGeometry *eegeom = nullptr, *fhgeom = nullptr, *bhgeomnew = nullptr;
289  const HcalGeometry* bhgeom = nullptr;
290  const HGCalDDDConstants* hgddd[3];
291  const HGCalTopology* hgtopo[3];
292  const HcalDDDRecConstants* hcddd = nullptr;
293  eegeom =
294  static_cast<const HGCalGeometry*>(geom->getSubdetectorGeometry(DetId::HGCalEE, ForwardSubdetector::ForwardEmpty));
295  //check if it's the new geometry
296  if (eegeom) {
297  geometryType_ = 1;
298  fhgeom = static_cast<const HGCalGeometry*>(
299  geom->getSubdetectorGeometry(DetId::HGCalHSi, ForwardSubdetector::ForwardEmpty));
300  bhgeomnew = static_cast<const HGCalGeometry*>(
301  geom->getSubdetectorGeometry(DetId::HGCalHSc, ForwardSubdetector::ForwardEmpty));
302  } else {
303  geometryType_ = 0;
304  eegeom = static_cast<const HGCalGeometry*>(geom->getSubdetectorGeometry(DetId::Forward, HGCEE));
305  fhgeom = static_cast<const HGCalGeometry*>(geom->getSubdetectorGeometry(DetId::Forward, HGCHEF));
306  bhgeom = static_cast<const HcalGeometry*>(geom->getSubdetectorGeometry(DetId::Hcal, HcalEndcap));
307  }
308  hgtopo[0] = &(eegeom->topology());
309  hgtopo[1] = &(fhgeom->topology());
310  if (bhgeomnew)
311  hgtopo[2] = &(bhgeomnew->topology());
312 
313  for (unsigned i = 0; i < 3; ++i) {
314  if (hgtopo[i])
315  hgddd[i] = &(hgtopo[i]->dddConstants());
316  }
317 
318  if (bhgeom)
319  hcddd = bhgeom->topology().dddConstants();
320 
321  // loop over the collections
322  int token = 0;
323  for (auto const& collectionTag : collectionTags_) {
325  const bool isHcal = (collectionTag.instance().find("HcalHits") != std::string::npos);
326  iEvent.getByToken(collectionTagsToken_[token++], hSimHits);
327  for (auto const& simHit : *hSimHits) {
328  DetId id(0);
329  const uint32_t simId = simHit.id();
330  if (geometryType_ == 1) {
331  //no test numbering in new geometry
332  id = simId;
333  } else if (isHcal) {
334  assert(hcddd);
335  HcalDetId hid = HcalHitRelabeller::relabel(simId, hcddd);
336  if (hid.subdet() == HcalEndcap)
337  id = hid;
338  } else {
339  int subdet, layer, cell, sec, subsec, zp;
340  HGCalTestNumbering::unpackHexagonIndex(simId, subdet, zp, layer, sec, subsec, cell);
341  const HGCalDDDConstants* ddd = hgddd[subdet - 3];
342  std::pair<int, int> recoLayerCell = ddd->simToReco(cell, layer, sec, hgtopo[subdet - 3]->detectorType());
343  cell = recoLayerCell.first;
344  layer = recoLayerCell.second;
345  // skip simhits with bad barcodes or non-existant layers
346  if (layer == -1 || simHit.geantTrackId() == 0)
347  continue;
348  id = HGCalDetId((ForwardSubdetector)subdet, zp, layer, subsec, sec, cell);
349  }
350 
351  if (DetId(0) == id)
352  continue;
353 
354  detIdToTotalSimEnergy[id.rawId()] += simHit.energy();
355  }
356  } // end of loop over InputTags
357 }
358 
359 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
362  desc.add<edm::InputTag>("simTracks", edm::InputTag("g4SimHits"));
363  desc.add<edm::InputTag>("genParticles", edm::InputTag("genParticles"));
364  desc.add<edm::InputTag>("simVertices", edm::InputTag("g4SimHits"));
365  desc.add<edm::InputTag>("trackingParticles", edm::InputTag("mix", "MergedTrackTruth"));
366  desc.add<edm::InputTag>("caloParticles", edm::InputTag("mix", "MergedCaloTruth"));
367  desc.add<edm::InputTag>("simClusters", edm::InputTag("mix", "MergedCaloTruth"));
368  desc.add<std::vector<edm::InputTag>>("collectionTags",
369  {edm::InputTag("g4SimHits", "HGCHitsEE"),
370  edm::InputTag("g4SimHits", "HGCHitsHEfront"),
371  edm::InputTag("g4SimHits", "HcalHits")});
372  descriptions.add("caloParticleDebugger", desc);
373 }
374 
375 // define this as a plug-in
EDGetTokenT< ProductType > consumes(edm::InputTag const &tag)
Log< level::Info, true > LogVerbatim
edm::EDGetTokenT< std::vector< TrackingParticle > > trackingParticlesToken_
T const * product() const
Definition: Handle.h:70
ForwardSubdetector
assert(be >=bs)
edm::EDGetTokenT< std::vector< SimCluster > > simClustersToken_
CaloParticleDebugger(const edm::ParameterSet &)
constexpr HcalSubdetector subdet() const
get the subdetector
Definition: HcalDetId.h:138
int iEvent
Definition: GenABIO.cc:224
std::vector< edm::EDGetTokenT< std::vector< PCaloHit > > > collectionTagsToken_
ConsumesCollector consumesCollector()
Use a ConsumesCollector to gather consumes information from helper functions.
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::pair< int, int > simToReco(int cell, int layer, int mod, bool half) const
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
edm::EDGetTokenT< std::vector< SimTrack > > simTracksToken_
DetId relabel(const uint32_t testId) const
ESHandle< T > getHandle(const ESGetToken< T, R > &iToken) const
Definition: EventSetup.h:130
edm::EDGetTokenT< std::vector< CaloParticle > > caloParticlesToken_
Definition: DetId.h:17
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
void fillSimHits(std::map< int, float > &, const edm::Event &, const edm::EventSetup &)
void add(std::string const &label, ParameterSetDescription const &psetDescription)
edm::EDGetTokenT< std::vector< reco::GenParticle > > genParticlesToken_
HLT enums.
const HcalDDDRecConstants * dddConstants() const
Definition: HcalTopology.h:164
void analyze(const edm::Event &, const edm::EventSetup &) override
edm::InputTag trackingParticles_
edm::ESGetToken< CaloGeometry, CaloGeometryRecord > geomToken_
const HGCalDDDConstants & dddConstants() const
Definition: HGCalTopology.h:98
static void unpackHexagonIndex(const uint32_t &idx, int &subdet, int &z, int &lay, int &wafer, int &celltyp, int &cell)
std::vector< edm::InputTag > collectionTags_
const HcalTopology & topology() const
Definition: HcalGeometry.h:111
edm::EDGetTokenT< std::vector< SimVertex > > simVerticesToken_