File indexing completed on 2021-05-01 06:57:28
0001 #define DEBUG false
0002 #if DEBUG
0003
0004
0005 #pragma GCC diagnostic push
0006 #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
0007 #endif
0008
0009
0010 #include <boost/graph/adjacency_list.hpp>
0011 #include <boost/graph/breadth_first_search.hpp>
0012 #include <boost/graph/depth_first_search.hpp>
0013 #include <boost/graph/graphviz.hpp>
0014
0015 #if DEBUG
0016 #pragma GCC diagnostic pop
0017 #endif
0018
0019 #include <iterator>
0020 #include <memory>
0021
0022 #include <numeric> // for std::accumulate
0023 #include <unordered_map>
0024
0025 #include "FWCore/Framework/interface/ConsumesCollector.h"
0026 #include "FWCore/Framework/interface/ESWatcher.h"
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/EventSetup.h"
0029 #include "FWCore/Framework/interface/ProducesCollector.h"
0030 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0032 #include "FWCore/Utilities/interface/InputTag.h"
0033
0034 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0035 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0036 #include "DataFormats/HcalDetId/interface/HcalTestNumbering.h"
0037 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0038 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0039 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0040
0041 #include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
0042 #include "SimDataFormats/CaloAnalysis/interface/CaloParticleFwd.h"
0043 #include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
0044 #include "SimDataFormats/CaloAnalysis/interface/SimClusterFwd.h"
0045 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0046 #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h"
0047 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0048 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0049
0050 #include "SimGeneral/MixingModule/interface/DigiAccumulatorMixMod.h"
0051 #include "SimGeneral/MixingModule/interface/DigiAccumulatorMixModFactory.h"
0052 #include "SimGeneral/MixingModule/interface/PileUpEventPrincipal.h"
0053 #include "SimGeneral/TrackingAnalysis/interface/EncodedTruthId.h"
0054
0055 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0056 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0057 #include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h"
0058 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0059 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0060
0061 namespace {
0062 using Index_t = unsigned;
0063 using Barcode_t = int;
0064 const std::string messageCategoryGraph_("CaloTruthAccumulatorGraphProducer");
0065 }
0066
0067 using boost::add_edge;
0068 using boost::adjacency_list;
0069 using boost::directedS;
0070 using boost::edge;
0071 using boost::edge_weight;
0072 using boost::edge_weight_t;
0073 using boost::listS;
0074 using boost::property;
0075 using boost::vecS;
0076 using boost::vertex;
0077 using boost::vertex_name;
0078 using boost::vertex_name_t;
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106 struct EdgeProperty {
0107 EdgeProperty(const SimTrack *t, int h, int c) : simTrack(t), simHits(h), cumulative_simHits(c) {}
0108 const SimTrack *simTrack;
0109 int simHits;
0110 int cumulative_simHits;
0111 };
0112
0113 struct VertexProperty {
0114 VertexProperty() : simTrack(nullptr), cumulative_simHits(0) {}
0115 VertexProperty(const SimTrack *t, int c) : simTrack(t), cumulative_simHits(c) {}
0116 VertexProperty(const VertexProperty &other)
0117 : simTrack(other.simTrack), cumulative_simHits(other.cumulative_simHits) {}
0118 const SimTrack *simTrack;
0119 int cumulative_simHits;
0120 };
0121
0122 using EdgeParticleClustersProperty = property<edge_weight_t, EdgeProperty>;
0123 using VertexMotherParticleProperty = property<vertex_name_t, VertexProperty>;
0124 using DecayChain = adjacency_list<listS, vecS, directedS, VertexMotherParticleProperty, EdgeParticleClustersProperty>;
0125
0126 class CaloTruthAccumulator : public DigiAccumulatorMixMod {
0127 public:
0128 explicit CaloTruthAccumulator(const edm::ParameterSet &config, edm::ProducesCollector, edm::ConsumesCollector &iC);
0129
0130 private:
0131 void initializeEvent(const edm::Event &event, const edm::EventSetup &setup) override;
0132 void accumulate(const edm::Event &event, const edm::EventSetup &setup) override;
0133 void accumulate(const PileUpEventPrincipal &event, const edm::EventSetup &setup, edm::StreamID const &) override;
0134 void finalizeEvent(edm::Event &event, const edm::EventSetup &setup) override;
0135
0136
0137 template <class T>
0138 void accumulateEvent(const T &event,
0139 const edm::EventSetup &setup,
0140 const edm::Handle<edm::HepMCProduct> &hepMCproduct);
0141
0142
0143
0144 template <class T>
0145 void fillSimHits(std::vector<std::pair<DetId, const PCaloHit *>> &returnValue,
0146 std::unordered_map<int, std::map<int, float>> &simTrackDetIdEnergyMap,
0147 const T &event,
0148 const edm::EventSetup &setup);
0149
0150 const std::string messageCategory_;
0151
0152 std::unordered_map<Index_t, float> m_detIdToTotalSimEnergy;
0153 std::unordered_multimap<Barcode_t, Index_t> m_simHitBarcodeToIndex;
0154
0155
0156
0157
0158
0159
0160 const unsigned int maximumPreviousBunchCrossing_;
0161
0162
0163
0164
0165
0166 const unsigned int maximumSubsequentBunchCrossing_;
0167
0168 const edm::InputTag simTrackLabel_;
0169 const edm::InputTag simVertexLabel_;
0170 edm::Handle<std::vector<SimTrack>> hSimTracks;
0171 edm::Handle<std::vector<SimVertex>> hSimVertices;
0172
0173 std::vector<edm::InputTag> collectionTags_;
0174 edm::InputTag genParticleLabel_;
0175
0176 edm::InputTag hepMCproductLabel_;
0177 const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geomToken_;
0178 edm::ESWatcher<CaloGeometryRecord> geomWatcher_;
0179
0180 const double minEnergy_, maxPseudoRapidity_;
0181 const bool premixStage1_;
0182
0183 public:
0184 struct OutputCollections {
0185 std::unique_ptr<SimClusterCollection> pSimClusters;
0186 std::unique_ptr<CaloParticleCollection> pCaloParticles;
0187 };
0188
0189 struct calo_particles {
0190 std::vector<uint32_t> sc_start_;
0191 std::vector<uint32_t> sc_stop_;
0192
0193 void swap(calo_particles &oth) {
0194 sc_start_.swap(oth.sc_start_);
0195 sc_stop_.swap(oth.sc_stop_);
0196 }
0197
0198 void clear() {
0199 sc_start_.clear();
0200 sc_stop_.clear();
0201 }
0202 };
0203
0204 private:
0205 const HGCalTopology *hgtopo_[3] = {nullptr, nullptr, nullptr};
0206 const HGCalDDDConstants *hgddd_[3] = {nullptr, nullptr, nullptr};
0207 const HcalDDDRecConstants *hcddd_ = nullptr;
0208 OutputCollections output_;
0209 calo_particles m_caloParticles;
0210
0211 int geometryType_;
0212 bool doHGCAL;
0213 };
0214
0215
0216
0217 namespace {
0218 template <typename Edge, typename Graph, typename Visitor>
0219 void accumulateSimHits_edge(Edge &e, const Graph &g, Visitor *v) {
0220 auto const edge_property = get(edge_weight, g, e);
0221 v->total_simHits += edge_property.simHits;
0222 IfLogDebug(DEBUG, messageCategoryGraph_)
0223 << " Examining edges " << e << " --> particle " << edge_property.simTrack->type() << "("
0224 << edge_property.simTrack->trackId() << ")"
0225 << " with SimClusters: " << edge_property.simHits << " Accumulated SimClusters: " << v->total_simHits
0226 << std::endl;
0227 }
0228 template <typename Vertex, typename Graph>
0229 void print_vertex(Vertex &u, const Graph &g) {
0230 auto const vertex_property = get(vertex_name, g, u);
0231 IfLogDebug(DEBUG, messageCategoryGraph_) << " At " << u;
0232
0233 if (vertex_property.simTrack)
0234 IfLogDebug(DEBUG, messageCategoryGraph_) << " [" << vertex_property.simTrack->type() << "]"
0235 << "(" << vertex_property.simTrack->trackId() << ")";
0236 IfLogDebug(DEBUG, messageCategoryGraph_) << std::endl;
0237 }
0238
0239
0240 #if DEBUG
0241 std::string graphviz_vertex(const VertexProperty &v) {
0242 std::ostringstream oss;
0243 oss << "{id: " << (v.simTrack ? v.simTrack->trackId() : 0) << ",\\ntype: " << (v.simTrack ? v.simTrack->type() : 0)
0244 << ",\\nchits: " << v.cumulative_simHits << "}";
0245 return oss.str();
0246 }
0247
0248 std::string graphviz_edge(const EdgeProperty &e) {
0249 std::ostringstream oss;
0250 oss << "[" << (e.simTrack ? e.simTrack->trackId() : 0) << "," << (e.simTrack ? e.simTrack->type() : 0) << ","
0251 << e.simHits << "," << e.cumulative_simHits << "]";
0252 return oss.str();
0253 }
0254 #endif
0255
0256 class SimHitsAccumulator_dfs_visitor : public boost::default_dfs_visitor {
0257 public:
0258 int total_simHits = 0;
0259 template <typename Edge, typename Graph>
0260 void examine_edge(Edge e, const Graph &g) {
0261 accumulateSimHits_edge(e, g, this);
0262 }
0263 template <typename Edge, typename Graph>
0264 void finish_edge(Edge e, const Graph &g) {
0265 auto const edge_property = get(edge_weight, g, e);
0266 auto src = source(e, g);
0267 auto trg = target(e, g);
0268 auto cumulative = edge_property.simHits + get(vertex_name, g, trg).cumulative_simHits +
0269 (get(vertex_name, g, src).simTrack ? get(vertex_name, g, src).cumulative_simHits
0270 : 0);
0271
0272 auto const src_vertex_property = get(vertex_name, g, src);
0273 put(get(vertex_name, const_cast<Graph &>(g)), src, VertexProperty(src_vertex_property.simTrack, cumulative));
0274 put(get(edge_weight, const_cast<Graph &>(g)),
0275 e,
0276 EdgeProperty(edge_property.simTrack, edge_property.simHits, cumulative));
0277 IfLogDebug(DEBUG, messageCategoryGraph_)
0278 << " Finished edge: " << e << " Track id: " << get(edge_weight, g, e).simTrack->trackId()
0279 << " has accumulated " << cumulative << " hits" << std::endl;
0280 IfLogDebug(DEBUG, messageCategoryGraph_) << " SrcVtx: " << src << "\t" << get(vertex_name, g, src).simTrack
0281 << "\t" << get(vertex_name, g, src).cumulative_simHits << std::endl;
0282 IfLogDebug(DEBUG, messageCategoryGraph_) << " TrgVtx: " << trg << "\t" << get(vertex_name, g, trg).simTrack
0283 << "\t" << get(vertex_name, g, trg).cumulative_simHits << std::endl;
0284 }
0285 };
0286
0287 using Selector = std::function<bool(EdgeProperty &)>;
0288
0289 class CaloParticle_dfs_visitor : public boost::default_dfs_visitor {
0290 public:
0291 CaloParticle_dfs_visitor(CaloTruthAccumulator::OutputCollections &output,
0292 CaloTruthAccumulator::calo_particles &caloParticles,
0293 std::unordered_multimap<Barcode_t, Index_t> &simHitBarcodeToIndex,
0294 std::unordered_map<int, std::map<int, float>> &simTrackDetIdEnergyMap,
0295 Selector selector)
0296 : output_(output),
0297 caloParticles_(caloParticles),
0298 simHitBarcodeToIndex_(simHitBarcodeToIndex),
0299 simTrackDetIdEnergyMap_(simTrackDetIdEnergyMap),
0300 selector_(selector) {}
0301 template <typename Vertex, typename Graph>
0302 void discover_vertex(Vertex u, const Graph &g) {
0303
0304
0305
0306 print_vertex(u, g);
0307 auto const vertex_property = get(vertex_name, g, u);
0308 if (!vertex_property.simTrack)
0309 return;
0310 auto trackIdx = vertex_property.simTrack->trackId();
0311 IfLogDebug(DEBUG, messageCategoryGraph_)
0312 << " Found " << simHitBarcodeToIndex_.count(trackIdx) << " associated simHits" << std::endl;
0313 if (simHitBarcodeToIndex_.count(trackIdx)) {
0314 output_.pSimClusters->emplace_back(*vertex_property.simTrack);
0315 auto &simcluster = output_.pSimClusters->back();
0316 std::unordered_map<uint32_t, float> acc_energy;
0317 for (auto const &hit_and_energy : simTrackDetIdEnergyMap_[trackIdx]) {
0318 acc_energy[hit_and_energy.first] += hit_and_energy.second;
0319 }
0320 for (auto const &hit_and_energy : acc_energy) {
0321 simcluster.addRecHitAndFraction(hit_and_energy.first, hit_and_energy.second);
0322 }
0323 }
0324 }
0325 template <typename Edge, typename Graph>
0326 void examine_edge(Edge e, const Graph &g) {
0327 auto src = source(e, g);
0328 auto vertex_property = get(vertex_name, g, src);
0329 if (src == 0 or (vertex_property.simTrack == nullptr)) {
0330 auto edge_property = get(edge_weight, g, e);
0331 IfLogDebug(DEBUG, messageCategoryGraph_) << "Considering CaloParticle: " << edge_property.simTrack->trackId();
0332 if (selector_(edge_property)) {
0333 IfLogDebug(DEBUG, messageCategoryGraph_) << "Adding CaloParticle: " << edge_property.simTrack->trackId();
0334 output_.pCaloParticles->emplace_back(*(edge_property.simTrack));
0335 caloParticles_.sc_start_.push_back(output_.pSimClusters->size());
0336 }
0337 }
0338 }
0339
0340 template <typename Edge, typename Graph>
0341 void finish_edge(Edge e, const Graph &g) {
0342 auto src = source(e, g);
0343 auto vertex_property = get(vertex_name, g, src);
0344 if (src == 0 or (vertex_property.simTrack == nullptr)) {
0345 auto edge_property = get(edge_weight, g, e);
0346 if (selector_(edge_property)) {
0347 caloParticles_.sc_stop_.push_back(output_.pSimClusters->size());
0348 }
0349 }
0350 }
0351
0352 private:
0353 CaloTruthAccumulator::OutputCollections &output_;
0354 CaloTruthAccumulator::calo_particles &caloParticles_;
0355 std::unordered_multimap<Barcode_t, Index_t> &simHitBarcodeToIndex_;
0356 std::unordered_map<int, std::map<int, float>> &simTrackDetIdEnergyMap_;
0357 Selector selector_;
0358 };
0359 }
0360
0361 CaloTruthAccumulator::CaloTruthAccumulator(const edm::ParameterSet &config,
0362 edm::ProducesCollector producesCollector,
0363 edm::ConsumesCollector &iC)
0364 : messageCategory_("CaloTruthAccumulator"),
0365 maximumPreviousBunchCrossing_(config.getParameter<unsigned int>("maximumPreviousBunchCrossing")),
0366 maximumSubsequentBunchCrossing_(config.getParameter<unsigned int>("maximumSubsequentBunchCrossing")),
0367 simTrackLabel_(config.getParameter<edm::InputTag>("simTrackCollection")),
0368 simVertexLabel_(config.getParameter<edm::InputTag>("simVertexCollection")),
0369 collectionTags_(),
0370 genParticleLabel_(config.getParameter<edm::InputTag>("genParticleCollection")),
0371 hepMCproductLabel_(config.getParameter<edm::InputTag>("HepMCProductLabel")),
0372 geomToken_(iC.esConsumes()),
0373 minEnergy_(config.getParameter<double>("MinEnergy")),
0374 maxPseudoRapidity_(config.getParameter<double>("MaxPseudoRapidity")),
0375 premixStage1_(config.getParameter<bool>("premixStage1")),
0376 geometryType_(-1),
0377 doHGCAL(config.getParameter<bool>("doHGCAL")) {
0378 producesCollector.produces<SimClusterCollection>("MergedCaloTruth");
0379 producesCollector.produces<CaloParticleCollection>("MergedCaloTruth");
0380 if (premixStage1_) {
0381 producesCollector.produces<std::vector<std::pair<unsigned int, float>>>("MergedCaloTruth");
0382 }
0383
0384 iC.consumes<std::vector<SimTrack>>(simTrackLabel_);
0385 iC.consumes<std::vector<SimVertex>>(simVertexLabel_);
0386 iC.consumes<std::vector<reco::GenParticle>>(genParticleLabel_);
0387 iC.consumes<std::vector<int>>(genParticleLabel_);
0388 iC.consumes<std::vector<int>>(hepMCproductLabel_);
0389
0390
0391 const edm::ParameterSet &simHitCollectionConfig = config.getParameterSet("simHitCollections");
0392 std::vector<std::string> parameterNames = simHitCollectionConfig.getParameterNames();
0393
0394 for (auto const ¶meterName : parameterNames) {
0395 std::vector<edm::InputTag> tags = simHitCollectionConfig.getParameter<std::vector<edm::InputTag>>(parameterName);
0396 collectionTags_.insert(collectionTags_.end(), tags.begin(), tags.end());
0397 }
0398
0399 for (auto const &collectionTag : collectionTags_) {
0400 iC.consumes<std::vector<PCaloHit>>(collectionTag);
0401 }
0402 }
0403
0404 void CaloTruthAccumulator::initializeEvent(edm::Event const &event, edm::EventSetup const &setup) {
0405 output_.pSimClusters = std::make_unique<SimClusterCollection>();
0406 output_.pCaloParticles = std::make_unique<CaloParticleCollection>();
0407
0408 m_detIdToTotalSimEnergy.clear();
0409
0410 if (geomWatcher_.check(setup)) {
0411 auto const &geom = setup.getData(geomToken_);
0412 const HGCalGeometry *eegeom = nullptr, *fhgeom = nullptr, *bhgeomnew = nullptr;
0413 const HcalGeometry *bhgeom = nullptr;
0414 bhgeom = static_cast<const HcalGeometry *>(geom.getSubdetectorGeometry(DetId::Hcal, HcalEndcap));
0415
0416 if (doHGCAL) {
0417 eegeom = static_cast<const HGCalGeometry *>(
0418 geom.getSubdetectorGeometry(DetId::HGCalEE, ForwardSubdetector::ForwardEmpty));
0419
0420 if (eegeom) {
0421 geometryType_ = 1;
0422 fhgeom = static_cast<const HGCalGeometry *>(
0423 geom.getSubdetectorGeometry(DetId::HGCalHSi, ForwardSubdetector::ForwardEmpty));
0424 bhgeomnew = static_cast<const HGCalGeometry *>(
0425 geom.getSubdetectorGeometry(DetId::HGCalHSc, ForwardSubdetector::ForwardEmpty));
0426 } else {
0427 geometryType_ = 0;
0428 eegeom = static_cast<const HGCalGeometry *>(geom.getSubdetectorGeometry(DetId::Forward, HGCEE));
0429 fhgeom = static_cast<const HGCalGeometry *>(geom.getSubdetectorGeometry(DetId::Forward, HGCHEF));
0430 bhgeom = static_cast<const HcalGeometry *>(geom.getSubdetectorGeometry(DetId::Hcal, HcalEndcap));
0431 }
0432 hgtopo_[0] = &(eegeom->topology());
0433 hgtopo_[1] = &(fhgeom->topology());
0434 if (bhgeomnew)
0435 hgtopo_[2] = &(bhgeomnew->topology());
0436
0437 for (unsigned i = 0; i < 3; ++i) {
0438 if (hgtopo_[i])
0439 hgddd_[i] = &(hgtopo_[i]->dddConstants());
0440 }
0441 }
0442
0443 if (bhgeom) {
0444 hcddd_ = bhgeom->topology().dddConstants();
0445 }
0446 }
0447 }
0448
0449
0450
0451
0452
0453
0454
0455
0456 void CaloTruthAccumulator::accumulate(edm::Event const &event, edm::EventSetup const &setup) {
0457 edm::Handle<edm::HepMCProduct> hepmc;
0458 event.getByLabel(hepMCproductLabel_, hepmc);
0459
0460 edm::LogInfo(messageCategory_) << " CaloTruthAccumulator::accumulate (signal)";
0461 accumulateEvent(event, setup, hepmc);
0462 }
0463
0464 void CaloTruthAccumulator::accumulate(PileUpEventPrincipal const &event,
0465 edm::EventSetup const &setup,
0466 edm::StreamID const &) {
0467 if (event.bunchCrossing() >= -static_cast<int>(maximumPreviousBunchCrossing_) &&
0468 event.bunchCrossing() <= static_cast<int>(maximumSubsequentBunchCrossing_)) {
0469
0470 edm::Handle<edm::HepMCProduct> hepmc;
0471 edm::LogInfo(messageCategory_) << " CaloTruthAccumulator::accumulate (pileup) bunchCrossing="
0472 << event.bunchCrossing();
0473 accumulateEvent(event, setup, hepmc);
0474 } else {
0475 edm::LogInfo(messageCategory_) << "Skipping pileup event for bunch crossing " << event.bunchCrossing();
0476 }
0477 }
0478
0479 void CaloTruthAccumulator::finalizeEvent(edm::Event &event, edm::EventSetup const &setup) {
0480 edm::LogInfo(messageCategory_) << "Adding " << output_.pSimClusters->size() << " SimParticles and "
0481 << output_.pCaloParticles->size() << " CaloParticles to the event.";
0482
0483
0484
0485
0486
0487
0488 if (premixStage1_) {
0489 auto totalEnergies = std::make_unique<std::vector<std::pair<unsigned int, float>>>();
0490 totalEnergies->reserve(m_detIdToTotalSimEnergy.size());
0491 std::copy(m_detIdToTotalSimEnergy.begin(), m_detIdToTotalSimEnergy.end(), std::back_inserter(*totalEnergies));
0492 std::sort(totalEnergies->begin(), totalEnergies->end());
0493 event.put(std::move(totalEnergies), "MergedCaloTruth");
0494 } else {
0495 for (auto &sc : *(output_.pSimClusters)) {
0496 auto hitsAndEnergies = sc.hits_and_fractions();
0497 sc.clearHitsAndFractions();
0498 sc.clearHitsEnergy();
0499 for (auto &hAndE : hitsAndEnergies) {
0500 const float totalenergy = m_detIdToTotalSimEnergy[hAndE.first];
0501 float fraction = 0.;
0502 if (totalenergy > 0)
0503 fraction = hAndE.second / totalenergy;
0504 else
0505 edm::LogWarning(messageCategory_)
0506 << "TotalSimEnergy for hit " << hAndE.first << " is 0! The fraction for this hit cannot be computed.";
0507 sc.addRecHitAndFraction(hAndE.first, fraction);
0508 sc.addHitEnergy(hAndE.second);
0509 }
0510 }
0511 }
0512
0513
0514 auto scHandle = event.put(std::move(output_.pSimClusters), "MergedCaloTruth");
0515
0516
0517 for (unsigned i = 0; i < output_.pCaloParticles->size(); ++i) {
0518 auto &cp = (*output_.pCaloParticles)[i];
0519 for (unsigned j = m_caloParticles.sc_start_[i]; j < m_caloParticles.sc_stop_[i]; ++j) {
0520 edm::Ref<SimClusterCollection> ref(scHandle, j);
0521 cp.addSimCluster(ref);
0522 }
0523 }
0524
0525 event.put(std::move(output_.pCaloParticles), "MergedCaloTruth");
0526
0527 calo_particles().swap(m_caloParticles);
0528
0529 std::unordered_map<Index_t, float>().swap(m_detIdToTotalSimEnergy);
0530 std::unordered_multimap<Barcode_t, Index_t>().swap(m_simHitBarcodeToIndex);
0531 }
0532
0533 template <class T>
0534 void CaloTruthAccumulator::accumulateEvent(const T &event,
0535 const edm::EventSetup &setup,
0536 const edm::Handle<edm::HepMCProduct> &hepMCproduct) {
0537 edm::Handle<std::vector<reco::GenParticle>> hGenParticles;
0538 edm::Handle<std::vector<int>> hGenParticleIndices;
0539
0540 event.getByLabel(simTrackLabel_, hSimTracks);
0541 event.getByLabel(simVertexLabel_, hSimVertices);
0542
0543 event.getByLabel(genParticleLabel_, hGenParticles);
0544 event.getByLabel(genParticleLabel_, hGenParticleIndices);
0545
0546 std::vector<std::pair<DetId, const PCaloHit *>> simHitPointers;
0547 std::unordered_map<int, std::map<int, float>> simTrackDetIdEnergyMap;
0548 fillSimHits(simHitPointers, simTrackDetIdEnergyMap, event, setup);
0549
0550
0551 m_simHitBarcodeToIndex.clear();
0552 for (unsigned int i = 0; i < simHitPointers.size(); ++i) {
0553 m_simHitBarcodeToIndex.emplace(simHitPointers[i].second->geantTrackId(), i);
0554 }
0555
0556 auto const &tracks = *hSimTracks;
0557 auto const &vertices = *hSimVertices;
0558 std::unordered_map<int, int> trackid_to_track_index;
0559 DecayChain decay;
0560 int idx = 0;
0561
0562 IfLogDebug(DEBUG, messageCategory_) << " TRACKS" << std::endl;
0563 for (auto const &t : tracks) {
0564 IfLogDebug(DEBUG, messageCategory_) << " " << idx << "\t" << t.trackId() << "\t" << t << std::endl;
0565 trackid_to_track_index[t.trackId()] = idx;
0566 idx++;
0567 }
0568
0569
0570
0571
0572
0573
0574
0575
0576
0577
0578
0579
0580
0581
0582
0583
0584
0585
0586
0587
0588
0589
0590
0591
0592
0593
0594
0595 idx = 0;
0596 std::vector<int> used_sim_tracks(tracks.size(), 0);
0597 std::vector<int> collapsed_vertices(vertices.size(), 0);
0598 IfLogDebug(DEBUG, messageCategory_) << " VERTICES" << std::endl;
0599 for (auto const &v : vertices) {
0600 IfLogDebug(DEBUG, messageCategory_) << " " << idx++ << "\t" << v << std::endl;
0601 if (v.parentIndex() != -1) {
0602 auto trk_idx = trackid_to_track_index[v.parentIndex()];
0603 auto origin_vtx = tracks[trk_idx].vertIndex();
0604 if (used_sim_tracks[trk_idx]) {
0605
0606
0607
0608 collapsed_vertices[v.vertexId()] = used_sim_tracks[trk_idx];
0609 continue;
0610 }
0611
0612 if (collapsed_vertices[origin_vtx])
0613 origin_vtx = collapsed_vertices[origin_vtx];
0614 add_edge(origin_vtx,
0615 v.vertexId(),
0616 EdgeProperty(&tracks[trk_idx], simTrackDetIdEnergyMap[v.parentIndex()].size(), 0),
0617 decay);
0618 used_sim_tracks[trk_idx] = v.vertexId();
0619 }
0620 }
0621
0622 auto const &vertexMothersProp = get(vertex_name, decay);
0623
0624
0625 int offset = vertices.size();
0626 for (size_t i = 0; i < tracks.size(); ++i) {
0627 if (!used_sim_tracks[i]) {
0628 auto origin_vtx = tracks[i].vertIndex();
0629
0630 if (collapsed_vertices[origin_vtx])
0631 origin_vtx = collapsed_vertices[origin_vtx];
0632 add_edge(
0633 origin_vtx, offset, EdgeProperty(&tracks[i], simTrackDetIdEnergyMap[tracks[i].trackId()].size(), 0), decay);
0634
0635
0636
0637 put(vertexMothersProp, offset, VertexProperty(&tracks[i], 0));
0638 offset++;
0639 }
0640 }
0641 for (auto const &v : vertices) {
0642 if (v.parentIndex() != -1) {
0643
0644 if (collapsed_vertices[v.vertexId()])
0645 continue;
0646 put(vertexMothersProp, v.vertexId(), VertexProperty(&tracks[trackid_to_track_index[v.parentIndex()]], 0));
0647 }
0648 }
0649 SimHitsAccumulator_dfs_visitor vis;
0650 depth_first_search(decay, visitor(vis));
0651 CaloParticle_dfs_visitor caloParticleCreator(
0652 output_,
0653 m_caloParticles,
0654 m_simHitBarcodeToIndex,
0655 simTrackDetIdEnergyMap,
0656 [&](EdgeProperty &edge_property) -> bool {
0657
0658
0659
0660
0661 return (edge_property.cumulative_simHits != 0 and !edge_property.simTrack->noGenpart() and
0662 edge_property.simTrack->momentum().E() > minEnergy_ and
0663 std::abs(edge_property.simTrack->momentum().Eta()) < maxPseudoRapidity_);
0664 });
0665 depth_first_search(decay, visitor(caloParticleCreator));
0666
0667 #if DEBUG
0668 boost::write_graphviz(std::cout,
0669 decay,
0670 make_label_writer(make_transform_value_property_map(&graphviz_vertex, get(vertex_name, decay))),
0671 make_label_writer(make_transform_value_property_map(&graphviz_edge, get(edge_weight, decay))));
0672 #endif
0673 }
0674
0675 template <class T>
0676 void CaloTruthAccumulator::fillSimHits(std::vector<std::pair<DetId, const PCaloHit *>> &returnValue,
0677 std::unordered_map<int, std::map<int, float>> &simTrackDetIdEnergyMap,
0678 const T &event,
0679 const edm::EventSetup &setup) {
0680 for (auto const &collectionTag : collectionTags_) {
0681 edm::Handle<std::vector<PCaloHit>> hSimHits;
0682 const bool isHcal = (collectionTag.instance().find("HcalHits") != std::string::npos);
0683 event.getByLabel(collectionTag, hSimHits);
0684
0685 for (auto const &simHit : *hSimHits) {
0686 DetId id(0);
0687
0688
0689 if (doHGCAL) {
0690 const uint32_t simId = simHit.id();
0691 if (geometryType_ == 1) {
0692
0693 id = simId;
0694 } else if (isHcal) {
0695 HcalDetId hid = HcalHitRelabeller::relabel(simId, hcddd_);
0696 if (hid.subdet() == HcalEndcap)
0697 id = hid;
0698 } else {
0699 int subdet, layer, cell, sec, subsec, zp;
0700 HGCalTestNumbering::unpackHexagonIndex(simId, subdet, zp, layer, sec, subsec, cell);
0701 const HGCalDDDConstants *ddd = hgddd_[subdet - 3];
0702 std::pair<int, int> recoLayerCell = ddd->simToReco(cell, layer, sec, hgtopo_[subdet - 3]->detectorType());
0703 cell = recoLayerCell.first;
0704 layer = recoLayerCell.second;
0705
0706 if (layer == -1 || simHit.geantTrackId() == 0)
0707 continue;
0708 id = HGCalDetId((ForwardSubdetector)subdet, zp, layer, subsec, sec, cell);
0709 }
0710 } else {
0711 id = simHit.id();
0712
0713 if (isHcal) {
0714 HcalDetId hid = HcalHitRelabeller::relabel(simHit.id(), hcddd_);
0715 id = hid;
0716 }
0717 }
0718
0719 if (id == DetId(0)) {
0720 continue;
0721 }
0722 if (simHit.geantTrackId() == 0) {
0723 continue;
0724 }
0725
0726 returnValue.emplace_back(id, &simHit);
0727 simTrackDetIdEnergyMap[simHit.geantTrackId()][id.rawId()] += simHit.energy();
0728 m_detIdToTotalSimEnergy[id.rawId()] += simHit.energy();
0729 }
0730 }
0731 }
0732
0733
0734 DEFINE_DIGI_ACCUMULATOR(CaloTruthAccumulator);