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