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 "FWCore/Framework/interface/Frameworkfwd.h"
0008 #include "FWCore/Framework/interface/ConsumesCollector.h"
0009 #include "FWCore/Sources/interface/VectorInputSource.h"
0010 #include "FWCore/Utilities/interface/ESGetToken.h"
0011 #include "DataFormats/Provenance/interface/EventID.h"
0012 #include "FWCore/Framework/interface/EventPrincipal.h"
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014
0015 #include "TRandom.h"
0016 #include "TFile.h"
0017 #include "TH1F.h"
0018
0019 class TFile;
0020 class TH1F;
0021
0022 namespace CLHEP {
0023   class RandPoissonQ;
0024   class RandPoisson;
0025   class HepRandomEngine;
0026 }  // namespace CLHEP
0027
0028 class MixingModuleConfig;
0029 class MixingRcd;
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 beginStream(edm::StreamID);
0085     void endStream();
0086
0087     void beginRun(const edm::Run& run, const edm::EventSetup& setup);
0088     void beginLuminosityBlock(const edm::LuminosityBlock& lumi, const edm::EventSetup& setup);
0089
0090     void endRun(const edm::Run& run, const edm::EventSetup& setup);
0091     void endLuminosityBlock(const edm::LuminosityBlock& lumi, const edm::EventSetup& setup);
0092
0093     void setupPileUpEvent(const edm::EventSetup& setup);
0094
0095     void reload(const edm::EventSetup& setup);
0096
0097     void CalculatePileup(int MinBunch,
0098                          int MaxBunch,
0099                          std::vector<int>& PileupSelection,
0100                          std::vector<float>& TrueNumInteractions,
0101                          StreamID const&);
0102
0103     //template<typename T>
0104     // void recordEventForPlayback(EventPrincipal const& eventPrincipal,
0105     //            std::vector<edm::SecondaryEventIDAndFileInfo> &ids, T& eventOperator);
0106
0107     const unsigned int& input() const { return inputType_; }
0108     void input(unsigned int s) { inputType_ = s; }
0109
0110   private:
0111     std::unique_ptr<CLHEP::RandPoissonQ> const& poissonDistribution(StreamID const& streamID);
0112     std::unique_ptr<CLHEP::RandPoisson> const& poissonDistr_OOT(StreamID const& streamID);
0113     CLHEP::HepRandomEngine* randomEngine(StreamID const& streamID);
0114
0115     unsigned int inputType_;
0116     std::string type_;
0117     std::string Source_type_;
0118     double averageNumber_;
0119     int const intAverage_;
0120     std::shared_ptr<TH1F> histo_;
0121     bool histoDistribution_;
0122     bool probFunctionDistribution_;
0123     bool poisson_;
0124     bool fixed_;
0125     bool none_;
0126     bool manage_OOT_;
0127     bool poisson_OOT_;
0128     bool fixed_OOT_;
0129
0130     bool PU_Study_;
0131     std::string Study_type_;
0132
0133     int intFixed_OOT_;
0134     int intFixed_ITPU_;
0135
0136     int minBunch_cosmics_;
0137     int maxBunch_cosmics_;
0138
0139     edm::ESGetToken<MixingModuleConfig, MixingRcd> configToken_;
0140     size_t fileNameHash_;
0141     std::shared_ptr<ProductRegistry> productRegistry_;
0142     std::unique_ptr<VectorInputSource> const input_;
0143     std::shared_ptr<ProcessConfiguration> processConfiguration_;
0144     std::shared_ptr<ProcessContext> processContext_;
0145     std::shared_ptr<StreamContext> streamContext_;
0146     std::unique_ptr<EventPrincipal> eventPrincipal_;
0147     std::shared_ptr<LuminosityBlockPrincipal> lumiPrincipal_;
0148     std::shared_ptr<RunPrincipal> runPrincipal_;
0149     std::unique_ptr<SecondaryEventProvider> provider_;
0150     std::unique_ptr<CLHEP::RandPoissonQ> PoissonDistribution_;
0151     std::unique_ptr<CLHEP::RandPoisson> PoissonDistr_OOT_;
0152     CLHEP::HepRandomEngine* randomEngine_;
0153
0154     //TH1F *h1f;
0155     //TH1F *hprobFunction;
0156     //TFile *probFileHisto;
0157
0158     //playback info
0159     bool playback_;
0160
0161     // sequential reading
0162     bool sequential_;
0163   };
0164
0165   template <typename T>
0166   class RecordEventID {
0167   private:
0168     std::vector<edm::SecondaryEventIDAndFileInfo>& ids_;
0169     T& eventOperator_;
0170     int eventCount;
0171
0172   public:
0173     RecordEventID(std::vector<edm::SecondaryEventIDAndFileInfo>& ids, T& eventOperator)
0174         : ids_(ids), eventOperator_(eventOperator), eventCount(1) {}
0175     bool operator()(EventPrincipal const& eventPrincipal, size_t fileNameHash) {
0176       bool used = eventOperator_(eventPrincipal, eventCount);
0177       if (used) {
0178         ++eventCount;
0179         ids_.emplace_back(eventPrincipal.id(), fileNameHash);
0180       }
0181       return used;
0182     }
0183   };
0184
0185   /*! Generates events from a VectorInputSource.
0186    *  This function decides which method of VectorInputSource
0187    *  to call: sequential, random, or pre-specified.
0188    *  The ids are either ids to read or ids to store while reading.
0189    *  eventOperator has a type that matches the eventOperator in
0190    *  VectorInputSource::loopRandom.
0191    *
0192    *  The "signal" event is optionally used to restrict
0193    *  the secondary events used for pileup and mixing.
0194    */
0195   template <typename T>
0196   void PileUp::readPileUp(edm::EventID const& signal,
0197                           std::vector<edm::SecondaryEventIDAndFileInfo>& ids,
0198                           T eventOperator,
0199                           int const pileEventCnt,
0200                           StreamID const& streamID) {
0201     // One reason PileUp is responsible for recording event IDs is
0202     // that it is the one that knows how many events will be read.
0203     ids.reserve(pileEventCnt);
0204     RecordEventID<T> recorder(ids, eventOperator);
0205     int read = 0;
0206     CLHEP::HepRandomEngine* engine = (sequential_ ? nullptr : randomEngine(streamID));
0207     read = input_->loopOverEvents(*eventPrincipal_, fileNameHash_, pileEventCnt, recorder, engine, &signal);
0208     if (read != pileEventCnt)
0209       edm::LogWarning("PileUp") << "Could not read enough pileup events: only " << read << " out of " << pileEventCnt
0210                                 << " requested.";
0211   }
0212
0213   template <typename T>
0214   void PileUp::playPileUp(std::vector<edm::SecondaryEventIDAndFileInfo>::const_iterator begin,
0215                           std::vector<edm::SecondaryEventIDAndFileInfo>::const_iterator end,
0216                           std::vector<edm::SecondaryEventIDAndFileInfo>& ids,
0217                           T eventOperator) {
0218     //TrueNumInteractions.push_back( end - begin ) ;
0219     RecordEventID<T> recorder(ids, eventOperator);
0220     input_->loopSpecified(*eventPrincipal_, fileNameHash_, begin, end, recorder);
0221   }
0222
0223   template <typename T>
0224   void PileUp::playOldFormatPileUp(std::vector<edm::EventID>::const_iterator begin,
0225                                    std::vector<edm::EventID>::const_iterator end,
0226                                    std::vector<edm::SecondaryEventIDAndFileInfo>& ids,
0227                                    T eventOperator) {
0228     //TrueNumInteractions.push_back( end - begin ) ;
0229     RecordEventID<T> recorder(ids, eventOperator);
0230     input_->loopSpecified(*eventPrincipal_, fileNameHash_, begin, end, recorder);
0231   }
0232
0233 }  // namespace edm
0234
0235 #endif