PreMixingModule/plugins/PreMixingModule.cc

0001 /** \class PreMixingModule
0002  *
0003  * PreMixingModule is the EDProducer subclass that overlays premixed
0004  * MC events on top of MC. It is similar to DataMixingModule, but
0005  * tailored for premixing use case.
0006  *
0007  ************************************************************/
0008 #include "Mixing/Base/interface/BMixingModule.h"
0009
0010 #include "FWCore/Framework/interface/ConsumesCollector.h"
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/EventPrincipal.h"
0013 #include "FWCore/Framework/interface/ModuleContextSentry.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 #include "FWCore/Utilities/interface/EDMException.h"
0017 #include "FWCore/Utilities/interface/transform.h"
0018 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0019 #include "FWCore/ServiceRegistry/interface/InternalContext.h"
0020 #include "FWCore/ServiceRegistry/interface/ModuleCallingContext.h"
0021 #include "FWCore/ServiceRegistry/interface/ParentContext.h"
0022 #include "FWCore/ServiceRegistry/interface/Service.h"
0023
0024 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0025 #include "SimDataFormats/CrossingFrame/interface/CrossingFramePlaybackInfoNew.h"
0026 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
0027 #include "SimGeneral/MixingModule/interface/PileUpEventPrincipal.h"
0028
0029 #include "SimGeneral/PreMixingModule/interface/PreMixingWorker.h"
0030 #include "SimGeneral/PreMixingModule/interface/PreMixingWorkerFactory.h"
0031 #include "PreMixingPileupCopy.h"
0032
0033 #include <CLHEP/Random/RandomEngine.h>
0034
0035 #include <functional>
0036 #include <vector>
0037
0038 namespace edm {
0039   class PreMixingModule : public BMixingModule {
0040   public:
0041     PreMixingModule(const edm::ParameterSet& ps, MixingCache::Config const* globalConf);
0042
0043     ~PreMixingModule() override = default;
0044
0045     void checkSignal(const edm::Event& e) override{};
0046     void createnewEDProduct() override {}
0047     void addSignals(const edm::Event& e, const edm::EventSetup& ES) override;
0048     void doPileUp(edm::Event& e, const edm::EventSetup& ES) override;
0049     void put(edm::Event& e, const edm::EventSetup& ES) override;
0050
0051     void initializeEvent(edm::Event const& e, edm::EventSetup const& eventSetup) override;
0052     void beginRun(edm::Run const& run, edm::EventSetup const& eventSetup) override;
0053     void beginLuminosityBlock(LuminosityBlock const& l1, EventSetup const& c) override;
0054     void endLuminosityBlock(LuminosityBlock const& l1, EventSetup const& c) override;
0055     void endRun(const edm::Run& r, const edm::EventSetup& setup) override;
0056
0057   private:
0058     class AdjustPileupDistribution {
0059     public:
0060       AdjustPileupDistribution(const edm::ParameterSet& ps)
0061           : firstRun_(ps.getParameter<unsigned int>("firstRun")),
0062             firstBinPileup_(ps.getParameter<unsigned int>("firstBinPileup")),
0063             pileupProbabilities_(ps.getParameter<std::vector<double>>("pileupProbabilities")) {
0064         for (double p : pileupProbabilities_) {
0065           if (p < 0. or p > 1.) {
0066             throw cms::Exception("Configuration") << "Invalid probability value " << p << " for firstRun " << firstRun_
0067                                                   << ". The probability must be >= 0. and <= 1.";
0068           }
0069         }
0070       }
0071
0072       edm::RunNumber_t firstRun() const { return firstRun_; }
0073       double probability(float pileup) const {
0074         unsigned int bin = static_cast<unsigned int>(pileup);
0075         if (bin < firstBinPileup_ or bin >= firstBinPileup_ + pileupProbabilities_.size()) {
0076           edm::LogWarning("PreMixingModule")
0077               << "Got pileup event with true pileup " << pileup
0078               << " that is outside of the configured pileup adjustment bounds [" << firstBinPileup_ << ", "
0079               << firstBinPileup_ + pileupProbabilities_.size() - 1 << "]. Using probability 0.";
0080           return 0.;
0081         }
0082         return pileupProbabilities_[bin - firstBinPileup_];
0083       }
0084
0085     private:
0086       edm::RunNumber_t firstRun_;
0087       unsigned int firstBinPileup_;
0088       std::vector<double> pileupProbabilities_;
0089     };
0090
0091     bool pileWorker(const edm::EventPrincipal&,
0092                     int bcr,
0093                     int EventId,
0094                     const edm::EventSetup& ES,
0095                     ModuleCallingContext const*,
0096                     AdjustPileupDistribution const* pileupAdjuster);
0097
0098     PreMixingPileupCopy puWorker_;
0099     bool addedPileup_ = false;
0100
0101     std::vector<AdjustPileupDistribution> pileupAdjusters_;
0102     std::vector<std::unique_ptr<PreMixingWorker>> workers_;
0103   };
0104
0105   PreMixingModule::PreMixingModule(const edm::ParameterSet& ps, MixingCache::Config const* globalConf)
0106       : BMixingModule(ps, globalConf),
0107         puWorker_(ps.getParameter<edm::ParameterSet>("workers").getParameter<edm::ParameterSet>("pileup"),
0108                   producesCollector(),
0109                   consumesCollector()),
0110         pileupAdjusters_(
0111             edm::vector_transform(ps.getParameter<std::vector<edm::ParameterSet>>("adjustPileupDistribution"),
0112                                   [](const auto& ps) { return AdjustPileupDistribution(ps); })) {
0113     std::sort(pileupAdjusters_.begin(), pileupAdjusters_.end(), [](const auto& a, const auto& b) {
0114       return a.firstRun() < b.firstRun();
0115     });
0116
0117     const auto& workers = ps.getParameter<edm::ParameterSet>("workers");
0118     std::vector<std::string> names = workers.getParameterNames();
0119
0120     // Hack to keep the random number sequence unchanged for migration
0121     // from DataMixingModule to PreMixingModule. To be removed in a
0122     // subsequent PR doing only that.
0123     {
0124       std::vector<std::string> tmp;
0125       auto hack = [&](const std::string& name) {
0126         auto i = std::find(names.begin(), names.end(), name);
0127         if (i != names.end()) {
0128           tmp.push_back(*i);
0129           names.erase(i);
0130         }
0131       };
0132       hack("ecal");
0133       hack("hcal");
0134       hack("strip");
0135       hack("pixel");
0136       std::copy(names.begin(), names.end(), std::back_inserter(tmp));
0137       names = std::move(tmp);
0138     }
0139
0140     for (const auto& name : names) {
0141       if (name == "pileup") {
0142         continue;
0143       }
0144       const auto& pset = workers.getParameter<edm::ParameterSet>(name);
0145       std::string type = pset.getParameter<std::string>("workerType");
0146       workers_.emplace_back(
0147           PreMixingWorkerFactory::get()->create(type, pset, producesCollector(), consumesCollector()));
0148     }
0149   }
0150
0151   void PreMixingModule::initializeEvent(const edm::Event& e, const edm::EventSetup& ES) {
0152     for (auto& w : workers_) {
0153       w->initializeEvent(e, ES);
0154     }
0155   }
0156
0157   void PreMixingModule::beginRun(edm::Run const& run, const edm::EventSetup& ES) {
0158     BMixingModule::beginRun(run, ES);
0159     for (auto& w : workers_) {
0160       w->beginRun(run, ES);
0161     }
0162   }
0163
0164   void PreMixingModule::endRun(edm::Run const& run, const edm::EventSetup& ES) {
0165     for (auto& w : workers_) {
0166       w->endRun();
0167     }
0168     BMixingModule::endRun(run, ES);
0169   }
0170
0171   void PreMixingModule::addSignals(const edm::Event& e, const edm::EventSetup& ES) {
0172     // fill in maps of hits
0173
0174     LogDebug("PreMixingModule") << "===============> adding MC signals for " << e.id();
0175
0176     for (auto& w : workers_) {
0177       w->addSignals(e, ES);
0178     }
0179
0180     addedPileup_ = false;
0181   }
0182
0183   bool PreMixingModule::pileWorker(const EventPrincipal& ep,
0184                                    int bcr,
0185                                    int eventNr,
0186                                    const edm::EventSetup& ES,
0187                                    edm::ModuleCallingContext const* mcc,
0188                                    AdjustPileupDistribution const* pileupAdjuster) {
0189     InternalContext internalContext(ep.id(), mcc);
0190     ParentContext parentContext(&internalContext);
0191     ModuleCallingContext moduleCallingContext(&moduleDescription());
0192     ModuleContextSentry moduleContextSentry(&moduleCallingContext, parentContext);
0193
0194     PileUpEventPrincipal pep(ep, &moduleCallingContext, bcr);
0195
0196     if (pileupAdjuster) {
0197       float trueNumInteractions = puWorker_.getTrueNumInteractions(pep);
0198       double prob = pileupAdjuster->probability(static_cast<unsigned int>(trueNumInteractions));
0199       edm::Service<edm::RandomNumberGenerator> rng;
0200       CLHEP::HepRandomEngine& engine = rng->getEngine(ep.streamID());
0201       if (engine.flat() > prob) {
0202         // engine.flat() should give a double in ]0,1[ range
0203         // the choice above means that "prob = 1-ulp" is treatead as 1
0204         return false;
0205       }
0206     }
0207
0208     LogDebug("PreMixingModule") << "\n===============> adding pileups from event  " << ep.id() << " for bunchcrossing "
0209                                 << bcr;
0210
0211     // Note:  setupPileUpEvent may modify the run and lumi numbers of the EventPrincipal to match that of the primary event.
0212     setupPileUpEvent(ES);
0213
0214     // check and see if we need to copy the pileup information from
0215     // secondary stream to the output stream
0216     // We only have the pileup event here, so pick the first time and store the info
0217     if (!addedPileup_) {
0218       puWorker_.addPileupInfo(pep);
0219       addedPileup_ = true;
0220     }
0221
0222     // fill in maps of hits; same code as addSignals, except now applied to the pileup events
0223
0224     for (auto& w : workers_) {
0225       w->addPileups(pep, ES);
0226     }
0227
0228     return true;
0229   }
0230
0231   void PreMixingModule::doPileUp(edm::Event& e, const edm::EventSetup& ES) {
0232     using namespace std::placeholders;
0233
0234     std::vector<edm::SecondaryEventIDAndFileInfo> recordEventID;
0235     std::vector<int> PileupList;
0236     TrueNumInteractions_.clear();
0237
0238     ModuleCallingContext const* mcc = e.moduleCallingContext();
0239
0240     AdjustPileupDistribution const* pileupAdjuster = nullptr;
0241     if (not pileupAdjusters_.empty()) {
0242       // Find the adjustment settings for the run of the signal event
0243       // the container should be small-enough to not really gain
0244       // anything with binary search
0245       auto it = std::find_if(pileupAdjusters_.rbegin(),
0246                              pileupAdjusters_.rend(),
0247                              [iRun = e.id().run()](const auto& elem) { return elem.firstRun() <= iRun; });
0248       if (it == pileupAdjusters_.rend()) {
0249         throw cms::Exception("LogicError") << "Encountered run " << e.id().run()
0250                                            << ", but the first run available in the pileup adjustment configuration is "
0251                                            << pileupAdjusters_.front().firstRun() << ". Please fix the configuration.";
0252       }
0253       pileupAdjuster = &*it;
0254     }
0255
0256     for (int bunchCrossing = minBunch_; bunchCrossing <= maxBunch_; ++bunchCrossing) {
0257       for (unsigned int isource = 0; isource < maxNbSources_; ++isource) {
0258         std::shared_ptr<PileUp> source = inputSources_[isource];
0259         if (!source || !(source->doPileUp(bunchCrossing)))
0260           continue;
0261
0262         if (isource == 0)
0263           source->CalculatePileup(minBunch_, maxBunch_, PileupList, TrueNumInteractions_, e.streamID());
0264
0265         int NumPU_Events = 0;
0266         if (isource == 0) {
0267           NumPU_Events = PileupList[bunchCrossing - minBunch_];
0268         } else {
0269           // non-minbias pileup only gets one event for now. Fix later if desired.
0270           NumPU_Events = 1;
0271         }
0272
0273         for (auto& w : workers_) {
0274           w->initializeBunchCrossing(e, ES, bunchCrossing);
0275         }
0276
0277         source->readPileUp(
0278             e.id(),
0279             recordEventID,
0280             std::bind(
0281                 &PreMixingModule::pileWorker, std::ref(*this), _1, bunchCrossing, _2, std::cref(ES), mcc, pileupAdjuster),
0282             NumPU_Events,
0283             e.streamID());
0284
0285         for (auto& w : workers_) {
0286           w->finalizeBunchCrossing(e, ES, bunchCrossing);
0287         }
0288       }
0289     }
0290   }
0291
0292   void PreMixingModule::put(edm::Event& e, const edm::EventSetup& ES) {
0293     // individual workers...
0294     // move pileup first so we have access to the information for the put step
0295     const auto& ps = puWorker_.getPileupSummaryInfo();
0296     int bunchSpacing = puWorker_.getBunchSpacing();
0297
0298     for (auto& w : workers_) {
0299       w->put(e, ES, ps, bunchSpacing);
0300     }
0301
0302     puWorker_.putPileupInfo(e);
0303   }
0304
0305   void PreMixingModule::beginLuminosityBlock(LuminosityBlock const& l1, EventSetup const& c) {
0306     BMixingModule::beginLuminosityBlock(l1, c);
0307     for (auto& w : workers_) {
0308       w->beginLuminosityBlock(l1, c);
0309     }
0310   }
0311
0312   void PreMixingModule::endLuminosityBlock(LuminosityBlock const& l1, EventSetup const& c) {
0313     BMixingModule::endLuminosityBlock(l1, c);
0314   }
0315 }  // namespace edm
0316
0317 #include "FWCore/PluginManager/interface/PluginManager.h"
0318 #include "FWCore/Framework/interface/MakerMacros.h"
0319 using edm::PreMixingModule;
0320 DEFINE_FWK_MODULE(PreMixingModule);