PixelSeeding/interface/HitTripletEDProducerT.h

0001 #ifndef RecoTracker_PixelSeeding_HitTripletEDProducerT_h
0002 #define RecoTracker_PixelSeeding_HitTripletEDProducerT_h
0003
0004 #include "FWCore/Framework/interface/stream/EDProducer.h"
0005 #include "FWCore/Framework/interface/ConsumesCollector.h"
0006 #include "FWCore/Framework/interface/ProducesCollector.h"
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0010 #include "FWCore/Utilities/interface/EDGetToken.h"
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "DataFormats/Common/interface/Handle.h"
0013 #include "FWCore/Utilities/interface/RunningAverage.h"
0014
0015 #include "RecoTracker/TkHitPairs/interface/IntermediateHitDoublets.h"
0016 #include "RecoTracker/TkHitPairs/interface/RegionsSeedingHitSets.h"
0017 #include "RecoTracker/PixelSeeding/interface/OrderedHitTriplets.h"
0018 #include "RecoTracker/PixelSeeding/interface/IntermediateHitTriplets.h"
0019 #include "RecoTracker/PixelSeeding/interface/LayerTriplets.h"
0020
0021 #include <numeric>
0022
0023 namespace hitTripletEDProducerT {
0024   class ImplBase;
0025 }
0026
0027 template <typename T_HitTripletGenerator>
0028 class HitTripletEDProducerT : public edm::stream::EDProducer<> {
0029 public:
0030   HitTripletEDProducerT(const edm::ParameterSet& iConfig);
0031   ~HitTripletEDProducerT() override = default;
0032
0033   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0034
0035   void produce(edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0036
0037 private:
0038   edm::EDGetTokenT<IntermediateHitDoublets> doubletToken_;
0039
0040   std::unique_ptr<hitTripletEDProducerT::ImplBase> impl_;
0041 };
0042
0043 namespace hitTripletEDProducerT {
0044   class ImplBase {
0045   public:
0046     ImplBase() = default;
0047     virtual ~ImplBase() = default;
0048
0049     virtual void produces(edm::ProducesCollector) const = 0;
0050     virtual void produce(const IntermediateHitDoublets& regionDoublets,
0051                          edm::Event& iEvent,
0052                          const edm::EventSetup& iSetup) = 0;
0053
0054   protected:
0055     edm::RunningAverage localRA_;
0056   };
0057
0058   /////
0059   template <typename T_HitTripletGenerator>
0060   class ImplGeneratorBase : public ImplBase {
0061   public:
0062     ImplGeneratorBase(const edm::ParameterSet& iConfig, edm::ConsumesCollector& iC) : generator_(iConfig, iC) {}
0063     ~ImplGeneratorBase() override = default;
0064
0065   protected:
0066     T_HitTripletGenerator generator_;
0067   };
0068
0069   /////
0070   template <typename T_HitTripletGenerator, typename T_SeedingHitSets, typename T_IntermediateHitTriplets>
0071   class Impl : public ImplGeneratorBase<T_HitTripletGenerator> {
0072   public:
0073     Impl(const edm::ParameterSet& iConfig, edm::ConsumesCollector& iC)
0074         : ImplGeneratorBase<T_HitTripletGenerator>(iConfig, iC) {}
0075     ~Impl() override = default;
0076
0077     void produces(edm::ProducesCollector producesCollector) const override {
0078       T_SeedingHitSets::produces(producesCollector);
0079       T_IntermediateHitTriplets::produces(producesCollector);
0080     };
0081
0082     void produce(const IntermediateHitDoublets& regionDoublets,
0083                  edm::Event& iEvent,
0084                  const edm::EventSetup& iSetup) override {
0085       const SeedingLayerSetsHits& seedingLayerHits = regionDoublets.seedingLayerHits();
0086
0087       auto seedingHitSetsProducer = T_SeedingHitSets();
0088       auto intermediateHitTripletsProducer = T_IntermediateHitTriplets(&seedingLayerHits);
0089
0090       if (regionDoublets.empty()) {
0091         seedingHitSetsProducer.putEmpty(iEvent);
0092         intermediateHitTripletsProducer.putEmpty(iEvent);
0093         return;
0094       }
0095
0096       seedingHitSetsProducer.reserve(regionDoublets.regionSize(), this->localRA_.upper());
0097       intermediateHitTripletsProducer.reserve(regionDoublets.regionSize(), this->localRA_.upper());
0098
0099       // match-making of pair and triplet layers
0100       std::vector<LayerTriplets::LayerSetAndLayers> trilayers = LayerTriplets::layers(seedingLayerHits);
0101
0102       OrderedHitTriplets triplets;
0103       triplets.reserve(this->localRA_.upper());
0104       size_t triplets_total = 0;
0105
0106       LogDebug("HitTripletEDProducer") << "Creating triplets for " << regionDoublets.regionSize() << " regions, and "
0107                                        << trilayers.size() << " pair+3rd layers from "
0108                                        << regionDoublets.layerPairsSize() << " layer pairs";
0109
0110       for (const auto& regionLayerPairs : regionDoublets) {
0111         const TrackingRegion& region = regionLayerPairs.region();
0112
0113         auto hitCachePtr_filler_shs = seedingHitSetsProducer.beginRegion(&region, nullptr);
0114         auto hitCachePtr_filler_iht =
0115             intermediateHitTripletsProducer.beginRegion(&region, std::get<0>(hitCachePtr_filler_shs));
0116         auto hitCachePtr = std::get<0>(hitCachePtr_filler_iht);
0117
0118         LayerHitMapCache& hitCache = *hitCachePtr;
0119         hitCache.extend(regionLayerPairs.layerHitMapCache());
0120
0121         LogTrace("HitTripletEDProducer") << " starting region";
0122
0123         for (const auto& layerPair : regionLayerPairs) {
0124           LogTrace("HitTripletEDProducer")
0125               << "  starting layer pair " << layerPair.innerLayerIndex() << "," << layerPair.outerLayerIndex();
0126
0127           auto found = std::find_if(trilayers.begin(), trilayers.end(), [&](const LayerTriplets::LayerSetAndLayers& a) {
0128             return a.first[0].index() == layerPair.innerLayerIndex() &&
0129                    a.first[1].index() == layerPair.outerLayerIndex();
0130           });
0131           if (found == trilayers.end()) {
0132             auto exp = cms::Exception("LogicError") << "Did not find the layer pair from vector<pair+third layers>. "
0133                                                        "This is a sign of some internal inconsistency\n";
0134             exp << "I was looking for layer pair " << layerPair.innerLayerIndex() << "," << layerPair.outerLayerIndex()
0135                 << ". Triplets have the following pairs:\n";
0136             for (const auto& a : trilayers) {
0137               exp << " " << a.first[0].index() << "," << a.first[1].index() << ": 3rd layers";
0138               for (const auto& b : a.second) {
0139                 exp << " " << b.index();
0140               }
0141               exp << "\n";
0142             }
0143             throw exp;
0144           }
0145           const auto& thirdLayers = found->second;
0146
0147           this->generator_.hitTriplets(region,
0148                                        triplets,
0149                                        iEvent,
0150                                        iSetup,
0151                                        layerPair.doublets(),
0152                                        thirdLayers,
0153                                        intermediateHitTripletsProducer.tripletLastLayerIndexVector(),
0154                                        hitCache);
0155
0156 #ifdef EDM_ML_DEBUG
0157           LogTrace("HitTripletEDProducer")
0158               << "  created " << triplets.size() << " triplets for layer pair " << layerPair.innerLayerIndex() << ","
0159               << layerPair.outerLayerIndex() << " and 3rd layers";
0160           for (const auto& l : thirdLayers) {
0161             LogTrace("HitTripletEDProducer") << "   " << l.index();
0162           }
0163 #endif
0164
0165           triplets_total += triplets.size();
0166           seedingHitSetsProducer.fill(std::get<1>(hitCachePtr_filler_shs), triplets);
0167           intermediateHitTripletsProducer.fill(
0168               std::get<1>(hitCachePtr_filler_iht), layerPair.layerPair(), thirdLayers, triplets);
0169
0170           triplets.clear();
0171         }
0172       }
0173       this->localRA_.update(triplets_total);
0174
0175       seedingHitSetsProducer.put(iEvent);
0176       intermediateHitTripletsProducer.put(iEvent);
0177     }
0178   };
0179
0180   /////
0181   class DoNothing {
0182   public:
0183     DoNothing() {}
0184     explicit DoNothing(const SeedingLayerSetsHits* layers) {}
0185
0186     static void produces(edm::ProducesCollector) {}
0187
0188     void reserve(size_t, size_t) {}
0189
0190     auto beginRegion(const TrackingRegion*, LayerHitMapCache* ptr) { return std::make_tuple(ptr, 0); }
0191
0192     std::vector<int>* tripletLastLayerIndexVector() { return nullptr; }
0193
0194     void fill(int, const OrderedHitTriplets&) {}
0195     void fill(int,
0196               const IntermediateHitTriplets::LayerPair&,
0197               const std::vector<SeedingLayerSetsHits::SeedingLayer>&,
0198               const OrderedHitTriplets&) {}
0199
0200     void put(edm::Event& iEvent) {}
0201     void putEmpty(edm::Event& iEvent) {}
0202   };
0203
0204   /////
0205   class ImplSeedingHitSets {
0206   public:
0207     ImplSeedingHitSets() : seedingHitSets_(std::make_unique<RegionsSeedingHitSets>()) {}
0208
0209     static void produces(edm::ProducesCollector producesCollector) {
0210       producesCollector.produces<RegionsSeedingHitSets>();
0211     }
0212
0213     void reserve(size_t regionsSize, size_t localRAupper) { seedingHitSets_->reserve(regionsSize, localRAupper); }
0214
0215     auto beginRegion(const TrackingRegion* region, LayerHitMapCache* ptr) {
0216       hitCacheTmp_.clear();
0217       return std::make_tuple(&hitCacheTmp_, seedingHitSets_->beginRegion(region));
0218     }
0219
0220     void fill(RegionsSeedingHitSets::RegionFiller& filler, const OrderedHitTriplets& triplets) {
0221       for (const auto& trpl : triplets) {
0222         filler.emplace_back(trpl.inner(), trpl.middle(), trpl.outer());
0223       }
0224     }
0225
0226     void put(edm::Event& iEvent) {
0227       seedingHitSets_->shrink_to_fit();
0228       putEmpty(iEvent);
0229     }
0230     void putEmpty(edm::Event& iEvent) { iEvent.put(std::move(seedingHitSets_)); }
0231
0232   private:
0233     std::unique_ptr<RegionsSeedingHitSets> seedingHitSets_;
0234     LayerHitMapCache hitCacheTmp_;  // used if !produceIntermediateHitDoublets
0235   };
0236
0237   /////
0238   class ImplIntermediateHitTriplets {
0239   public:
0240     explicit ImplIntermediateHitTriplets(const SeedingLayerSetsHits* layers)
0241         : intermediateHitTriplets_(std::make_unique<IntermediateHitTriplets>(layers)), layers_(layers) {}
0242
0243     static void produces(edm::ProducesCollector producesCollector) {
0244       producesCollector.produces<IntermediateHitTriplets>();
0245     }
0246
0247     void reserve(size_t regionsSize, size_t localRAupper) {
0248       intermediateHitTriplets_->reserve(regionsSize, layers_->size(), localRAupper);
0249       tripletLastLayerIndex_.reserve(localRAupper);
0250     }
0251
0252     auto beginRegion(const TrackingRegion* region, LayerHitMapCache*) {
0253       auto filler = intermediateHitTriplets_->beginRegion(region);
0254       return std::make_tuple(&(filler.layerHitMapCache()), filler);
0255     }
0256
0257     std::vector<int>* tripletLastLayerIndexVector() { return &tripletLastLayerIndex_; }
0258
0259     void fill(IntermediateHitTriplets::RegionFiller& filler,
0260               const IntermediateHitTriplets::LayerPair& layerPair,
0261               const std::vector<SeedingLayerSetsHits::SeedingLayer>& thirdLayers,
0262               const OrderedHitTriplets& triplets) {
0263       if (tripletLastLayerIndex_.size() != triplets.size()) {
0264         throw cms::Exception("LogicError") << "tripletLastLayerIndex_.size() " << tripletLastLayerIndex_.size()
0265                                            << " triplets.size() " << triplets.size();
0266       }
0267       tripletPermutation_.resize(tripletLastLayerIndex_.size());
0268       std::iota(tripletPermutation_.begin(), tripletPermutation_.end(), 0);  // assign 0,1,2,...,N
0269       std::stable_sort(tripletPermutation_.begin(), tripletPermutation_.end(), [&](size_t i, size_t j) {
0270         return tripletLastLayerIndex_[i] < tripletLastLayerIndex_[j];
0271       });
0272
0273       // empty triplets need to propagate here
0274       filler.addTriplets(layerPair, thirdLayers, triplets, tripletLastLayerIndex_, tripletPermutation_);
0275       tripletLastLayerIndex_.clear();
0276     }
0277
0278     void put(edm::Event& iEvent) {
0279       intermediateHitTriplets_->shrink_to_fit();
0280       putEmpty(iEvent);
0281     }
0282     void putEmpty(edm::Event& iEvent) { iEvent.put(std::move(intermediateHitTriplets_)); }
0283
0284   private:
0285     std::unique_ptr<IntermediateHitTriplets> intermediateHitTriplets_;
0286     const SeedingLayerSetsHits* layers_;
0287     std::vector<size_t> tripletPermutation_;
0288     std::vector<int> tripletLastLayerIndex_;
0289   };
0290 }  // namespace hitTripletEDProducerT
0291
0292 template <typename T_HitTripletGenerator>
0293 HitTripletEDProducerT<T_HitTripletGenerator>::HitTripletEDProducerT(const edm::ParameterSet& iConfig)
0294     : doubletToken_(consumes<IntermediateHitDoublets>(iConfig.getParameter<edm::InputTag>("doublets"))) {
0295   const bool produceSeedingHitSets = iConfig.getParameter<bool>("produceSeedingHitSets");
0296   const bool produceIntermediateHitTriplets = iConfig.getParameter<bool>("produceIntermediateHitTriplets");
0297
0298   auto iC = consumesCollector();
0299
0300   using namespace hitTripletEDProducerT;
0301
0302   if (produceSeedingHitSets && produceIntermediateHitTriplets)
0303     impl_ = std::make_unique<Impl<T_HitTripletGenerator, ImplSeedingHitSets, ImplIntermediateHitTriplets>>(iConfig, iC);
0304   else if (produceSeedingHitSets)
0305     impl_ = std::make_unique<Impl<T_HitTripletGenerator, ImplSeedingHitSets, DoNothing>>(iConfig, iC);
0306   else if (produceIntermediateHitTriplets)
0307     impl_ = std::make_unique<Impl<T_HitTripletGenerator, DoNothing, ImplIntermediateHitTriplets>>(iConfig, iC);
0308   else
0309     throw cms::Exception("Configuration")
0310         << "HitTripletEDProducerT requires either produceIntermediateHitTriplets or produceSeedingHitSets to be True. "
0311            "If neither are needed, just remove this module from your sequence/path as it doesn't do anything useful";
0312
0313   impl_->produces(producesCollector());
0314 }
0315
0316 template <typename T_HitTripletGenerator>
0317 void HitTripletEDProducerT<T_HitTripletGenerator>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0318   edm::ParameterSetDescription desc;
0319
0320   desc.add<edm::InputTag>("doublets", edm::InputTag("hitPairEDProducer"));
0321   desc.add<bool>("produceSeedingHitSets", false);
0322   desc.add<bool>("produceIntermediateHitTriplets", false);
0323
0324   T_HitTripletGenerator::fillDescriptions(desc);
0325
0326   auto label = T_HitTripletGenerator::fillDescriptionsLabel() + std::string("EDProducerDefault");
0327   descriptions.add(label, desc);
0328 }
0329
0330 template <typename T_HitTripletGenerator>
0331 void HitTripletEDProducerT<T_HitTripletGenerator>::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0332   edm::Handle<IntermediateHitDoublets> hdoublets;
0333   iEvent.getByToken(doubletToken_, hdoublets);
0334   const auto& regionDoublets = *hdoublets;
0335
0336   const SeedingLayerSetsHits& seedingLayerHits = regionDoublets.seedingLayerHits();
0337   if (seedingLayerHits.numberOfLayersInSet() < 3) {
0338     throw cms::Exception("LogicError")
0339         << "HitTripletEDProducerT expects SeedingLayerSetsHits::numberOfLayersInSet() to be >= 3, got "
0340         << seedingLayerHits.numberOfLayersInSet()
0341         << ". This is likely caused by a configuration error of this module, HitPairEDProducer, or "
0342            "SeedingLayersEDProducer.";
0343   }
0344
0345   impl_->produce(regionDoublets, iEvent, iSetup);
0346 }
0347
0348 #endif