Base/interface/PileUp.h

0001 #ifndef Mixing_Base_PileUp_h
0002 #define Mixing_Base_PileUp_h
0003
0004 #include <memory>
0005 #include <string>
0006 #include <vector>
0007 #include <optional>
0008 #include "FWCore/Framework/interface/Frameworkfwd.h"
0009 #include "FWCore/Framework/interface/ConsumesCollector.h"
0010 #include "FWCore/Sources/interface/VectorInputSource.h"
0011 #include "FWCore/Utilities/interface/ESGetToken.h"
0012 #include "DataFormats/Provenance/interface/EventID.h"
0013 #include "FWCore/Framework/interface/EventPrincipal.h"
0014 #include "FWCore/ServiceRegistry/interface/ServiceToken.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 #include "FWCore/Utilities/interface/ExceptionCollector.h"
0017
0018 #include "TH1F.h"
0019
0020 namespace CLHEP {
0021   class RandPoissonQ;
0022   class RandPoisson;
0023   class HepRandomEngine;
0024 }  // namespace CLHEP
0025
0026 class MixingModuleConfig;
0027 class MixingRcd;
0028 class PileupRandomNumberGenerator;
0029
0030 namespace edm {
0031   class SecondaryEventProvider;
0032   class StreamID;
0033   class ProcessContext;
0034
0035   struct PileUpConfig {
0036     PileUpConfig(std::string sourcename, double averageNumber, std::unique_ptr<TH1F>& histo, const bool playback)
0037         : sourcename_(sourcename), averageNumber_(averageNumber), histo_(histo.release()), playback_(playback) {}
0038     std::string sourcename_;
0039     double averageNumber_;
0040     std::shared_ptr<TH1F> histo_;
0041     const bool playback_;
0042   };
0043
0044   class PileUp {
0045   public:
0046     explicit PileUp(ParameterSet const& pset,
0047                     const std::shared_ptr<PileUpConfig>& config,
0048                     edm::ConsumesCollector iC,
0049                     const bool mixingConfigFromDB);
0050     ~PileUp();
0051
0052     template <typename T>
0053     void readPileUp(edm::EventID const& signal,
0054                     std::vector<edm::SecondaryEventIDAndFileInfo>& ids,
0055                     T eventOperator,
0056                     int const NumPU,
0057                     StreamID const&);
0058
0059     template <typename T>
0060     void playPileUp(std::vector<edm::SecondaryEventIDAndFileInfo>::const_iterator begin,
0061                     std::vector<edm::SecondaryEventIDAndFileInfo>::const_iterator end,
0062                     std::vector<edm::SecondaryEventIDAndFileInfo>& ids,
0063                     T eventOperator);
0064
0065     template <typename T>
0066     void playOldFormatPileUp(std::vector<edm::EventID>::const_iterator begin,
0067                              std::vector<edm::EventID>::const_iterator end,
0068                              std::vector<edm::SecondaryEventIDAndFileInfo>& ids,
0069                              T eventOperator);
0070
0071     double averageNumber() const { return averageNumber_; }
0072     bool poisson() const { return poisson_; }
0073     bool doPileUp(int BX) {
0074       if (Source_type_ != "cosmics") {
0075         return none_ ? false : averageNumber_ > 0.;
0076       } else {
0077         return (BX >= minBunch_cosmics_ && BX <= maxBunch_cosmics_);
0078       }
0079     }
0080     void dropUnwantedBranches(std::vector<std::string> const& wantedBranches) {
0081       input_->dropUnwantedBranches(wantedBranches);
0082     }
0083     void beginJob(eventsetup::ESRecordsToProductResolverIndices const&);
0084     void beginStream(edm::StreamID);
0085     void endStream();
0086     void endStream(ExceptionCollector&);
0087
0088     void beginRun(const edm::Run& run, const edm::EventSetup& setup);
0089     void beginLuminosityBlock(const edm::LuminosityBlock& lumi, const edm::EventSetup& setup);
0090
0091     void endRun(const edm::Run& run, const edm::EventSetup& setup);
0092     void endLuminosityBlock(const edm::LuminosityBlock& lumi, const edm::EventSetup& setup);
0093
0094     void setupPileUpEvent(const edm::EventSetup& setup);
0095
0096     void reload(const edm::EventSetup& setup);
0097
0098     void CalculatePileup(int MinBunch,
0099                          int MaxBunch,
0100                          std::vector<int>& PileupSelection,
0101                          std::vector<float>& TrueNumInteractions,
0102                          StreamID const&);
0103
0104     //template<typename T>
0105     // void recordEventForPlayback(EventPrincipal const& eventPrincipal,
0106     //            std::vector<edm::SecondaryEventIDAndFileInfo> &ids, T& eventOperator);
0107
0108     const unsigned int& input() const { return inputType_; }
0109     void input(unsigned int s) { inputType_ = s; }
0110
0111   private:
0112     std::unique_ptr<CLHEP::RandPoissonQ> const& poissonDistribution(StreamID const& streamID);
0113     std::unique_ptr<CLHEP::RandPoisson> const& poissonDistr_OOT(StreamID const& streamID);
0114     CLHEP::HepRandomEngine* randomEngine(StreamID const& streamID);
0115     void setRandomEngine(StreamID);
0116     void setRandomEngine(LuminosityBlock const&);
0117
0118     unsigned int inputType_;
0119     std::string type_;
0120     std::string Source_type_;
0121     double averageNumber_;
0122     int const intAverage_;
0123     std::shared_ptr<TH1F> histo_;
0124     bool histoDistribution_;
0125     bool probFunctionDistribution_;
0126     bool poisson_;
0127     bool fixed_;
0128     bool none_;
0129     bool manage_OOT_;
0130     bool poisson_OOT_;
0131     bool fixed_OOT_;
0132
0133     bool PU_Study_;
0134     std::string Study_type_;
0135
0136     int intFixed_OOT_;
0137     int intFixed_ITPU_;
0138
0139     int minBunch_cosmics_;
0140     int maxBunch_cosmics_;
0141
0142     edm::ESGetToken<MixingModuleConfig, MixingRcd> configToken_;
0143     size_t fileNameHash_;
0144     std::shared_ptr<ProductRegistry> productRegistry_;
0145     std::unique_ptr<VectorInputSource> const input_;
0146     std::shared_ptr<ProcessConfiguration> processConfiguration_;
0147     std::shared_ptr<ProcessContext> processContext_;
0148     std::shared_ptr<StreamContext> streamContext_;
0149     std::unique_ptr<EventPrincipal> eventPrincipal_;
0150     std::shared_ptr<LuminosityBlockPrincipal> lumiPrincipal_;
0151     std::shared_ptr<RunPrincipal> runPrincipal_;
0152     PileupRandomNumberGenerator* randomGenerator_ = nullptr;
0153     std::optional<ServiceToken> serviceToken_;
0154     std::unique_ptr<SecondaryEventProvider> provider_;
0155     std::unique_ptr<CLHEP::RandPoissonQ> PoissonDistribution_;
0156     std::unique_ptr<CLHEP::RandPoisson> PoissonDistr_OOT_;
0157     CLHEP::HepRandomEngine* randomEngine_;
0158
0159     //TH1F *h1f;
0160     //TH1F *hprobFunction;
0161     //TFile *probFileHisto;
0162
0163     //playback info
0164     bool playback_;
0165
0166     // sequential reading
0167     bool sequential_;
0168   };
0169
0170   template <typename T>
0171   class RecordEventID {
0172   private:
0173     std::vector<edm::SecondaryEventIDAndFileInfo>& ids_;
0174     T& eventOperator_;
0175     int eventCount;
0176
0177   public:
0178     RecordEventID(std::vector<edm::SecondaryEventIDAndFileInfo>& ids, T& eventOperator)
0179         : ids_(ids), eventOperator_(eventOperator), eventCount(1) {}
0180     bool operator()(EventPrincipal const& eventPrincipal, size_t fileNameHash) {
0181       bool used = eventOperator_(eventPrincipal, eventCount);
0182       if (used) {
0183         ++eventCount;
0184         ids_.emplace_back(eventPrincipal.id(), fileNameHash);
0185       }
0186       return used;
0187     }
0188   };
0189
0190   /*! Generates events from a VectorInputSource.
0191    *  This function decides which method of VectorInputSource
0192    *  to call: sequential, random, or pre-specified.
0193    *  The ids are either ids to read or ids to store while reading.
0194    *  eventOperator has a type that matches the eventOperator in
0195    *  VectorInputSource::loopRandom.
0196    *
0197    *  The "signal" event is optionally used to restrict
0198    *  the secondary events used for pileup and mixing.
0199    */
0200   template <typename T>
0201   void PileUp::readPileUp(edm::EventID const& signal,
0202                           std::vector<edm::SecondaryEventIDAndFileInfo>& ids,
0203                           T eventOperator,
0204                           int const pileEventCnt,
0205                           StreamID const& streamID) {
0206     // One reason PileUp is responsible for recording event IDs is
0207     // that it is the one that knows how many events will be read.
0208     ids.reserve(pileEventCnt);
0209     RecordEventID<T> recorder(ids, eventOperator);
0210     int read = 0;
0211     CLHEP::HepRandomEngine* engine = (sequential_ ? nullptr : randomEngine(streamID));
0212     read = input_->loopOverEvents(*eventPrincipal_, fileNameHash_, pileEventCnt, recorder, engine, &signal);
0213     if (read != pileEventCnt)
0214       edm::LogWarning("PileUp") << "Could not read enough pileup events: only " << read << " out of " << pileEventCnt
0215                                 << " requested.";
0216   }
0217
0218   template <typename T>
0219   void PileUp::playPileUp(std::vector<edm::SecondaryEventIDAndFileInfo>::const_iterator begin,
0220                           std::vector<edm::SecondaryEventIDAndFileInfo>::const_iterator end,
0221                           std::vector<edm::SecondaryEventIDAndFileInfo>& ids,
0222                           T eventOperator) {
0223     //TrueNumInteractions.push_back( end - begin ) ;
0224     RecordEventID<T> recorder(ids, eventOperator);
0225     input_->loopSpecified(*eventPrincipal_, fileNameHash_, begin, end, recorder);
0226   }
0227
0228   template <typename T>
0229   void PileUp::playOldFormatPileUp(std::vector<edm::EventID>::const_iterator begin,
0230                                    std::vector<edm::EventID>::const_iterator end,
0231                                    std::vector<edm::SecondaryEventIDAndFileInfo>& ids,
0232                                    T eventOperator) {
0233     //TrueNumInteractions.push_back( end - begin ) ;
0234     RecordEventID<T> recorder(ids, eventOperator);
0235     input_->loopSpecified(*eventPrincipal_, fileNameHash_, begin, end, recorder);
0236   }
0237
0238 }  // namespace edm
0239
0240 #endif