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