Base/src/PileUp.cc

0001 #include "Mixing/Base/interface/PileUp.h"
0002 #include "DataFormats/Provenance/interface/BranchIDListHelper.h"
0003 #include "DataFormats/Provenance/interface/ModuleDescription.h"
0004 #include "DataFormats/Provenance/interface/ThinnedAssociationsHelper.h"
0005 #include "FWCore/Framework/interface/EventPrincipal.h"
0006 #include "FWCore/Framework/interface/LuminosityBlock.h"
0007 #include "FWCore/Framework/interface/Run.h"
0008 #include "FWCore/Framework/interface/SignallingProductRegistryFiller.h"
0009 #include "FWCore/Framework/interface/ESRecordsToProductResolverIndices.h"
0010 #include "FWCore/Framework/interface/ProductResolversFactory.h"
0011 #include "FWCore/ServiceRegistry/interface/ActivityRegistry.h"
0012 #include "FWCore/ServiceRegistry/interface/GlobalContext.h"
0013 #include "FWCore/ServiceRegistry/interface/ProcessContext.h"
0014 #include "FWCore/Sources/interface/VectorInputSourceDescription.h"
0015 #include "FWCore/Sources/interface/VectorInputSourceFactory.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 #include "FWCore/Utilities/interface/Exception.h"
0019 #include "FWCore/Utilities/interface/StreamID.h"
0020 #include "FWCore/Version/interface/GetReleaseVersion.h"
0021
0022 #include "FWCore/AbstractServices/interface/RandomNumberGenerator.h"
0023 #include "FWCore/ServiceRegistry/interface/Service.h"
0024
0025 #include "FWCore/Framework/interface/ESHandle.h"
0026 #include "FWCore/Framework/interface/EventSetup.h"
0027 #include "Mixing/Base/src/SecondaryEventProvider.h"
0028 #include "CondFormats/DataRecord/interface/MixingRcd.h"
0029 #include "CondFormats/RunInfo/interface/MixingModuleConfig.h"
0030
0031 #include "CLHEP/Random/RandPoissonQ.h"
0032 #include "CLHEP/Random/RandPoisson.h"
0033
0034 #include "PileupRandomNumberGenerator.h"
0035
0036 #include <algorithm>
0037 #include <memory>
0038 #include "TMath.h"
0039
0040 ////////////////////////////////////////////////////////////////////////////////
0041 /// return a random number distributed according the histogram bin contents.
0042 ///
0043 /// This routine is derived from root/hist/hist/src/TH1.cxx
0044 /*************************************************************************
0045  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
0046  * All rights reserved.                                                  *
0047  *                                                                       *
0048  * For the licensing terms see $ROOTSYS/LICENSE.                         *
0049  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
0050  *************************************************************************/
0051
0052 static Double_t GetRandom(TH1* th1, CLHEP::HepRandomEngine* rng) {
0053   Int_t nbinsx = th1->GetNbinsX();
0054   Double_t* fIntegral = th1->GetIntegral();
0055   Double_t integral = fIntegral[nbinsx];
0056
0057   if (integral == 0)
0058     return 0;
0059
0060   Double_t r1 = rng->flat();
0061   Int_t ibin = TMath::BinarySearch(nbinsx, fIntegral, r1);
0062   Double_t x = th1->GetBinLowEdge(ibin + 1);
0063   if (r1 > fIntegral[ibin])
0064     x += th1->GetBinWidth(ibin + 1) * (r1 - fIntegral[ibin]) / (fIntegral[ibin + 1] - fIntegral[ibin]);
0065   return x;
0066 }
0067 ////////////////////////////////////////////////////////////////////////////////
0068
0069 namespace edm {
0070   PileUp::PileUp(ParameterSet const& pset,
0071                  const std::shared_ptr<PileUpConfig>& config,
0072                  edm::ConsumesCollector iC,
0073                  const bool mixingConfigFromDB)
0074       : type_(pset.getParameter<std::string>("type")),
0075         Source_type_(config->sourcename_),
0076         averageNumber_(config->averageNumber_),
0077         intAverage_(static_cast<int>(averageNumber_)),
0078         histo_(std::make_shared<TH1F>(*config->histo_)),
0079         histoDistribution_(type_ == "histo"),
0080         probFunctionDistribution_(type_ == "probFunction"),
0081         poisson_(type_ == "poisson"),
0082         fixed_(type_ == "fixed"),
0083         none_(type_ == "none"),
0084         fileNameHash_(0U),
0085         productRegistry_(),
0086         input_(VectorInputSourceFactory::get()->makeVectorInputSource(
0087             pset,
0088             VectorInputSourceDescription(std::make_shared<edm::ProductRegistry>(), edm::PreallocationConfiguration()))),
0089         // hardware information is not needed for the "overlay"
0090         processConfiguration_(std::make_shared<ProcessConfiguration>(
0091             "@MIXING", getReleaseVersion(), edm::HardwareResourcesDescription())),
0092         processContext_(new ProcessContext()),
0093         eventPrincipal_(),
0094         lumiPrincipal_(),
0095         runPrincipal_(),
0096         provider_(),
0097         PoissonDistribution_(),
0098         PoissonDistr_OOT_(),
0099         randomEngine_(),
0100         playback_(config->playback_),
0101         sequential_(pset.getUntrackedParameter<bool>("sequential", false)) {
0102     if (mixingConfigFromDB) {
0103       configToken_ = iC.esConsumes<edm::Transition::BeginLuminosityBlock>();
0104     }
0105
0106     // Use the empty parameter set for the parameter set ID of our "@MIXING" process.
0107     processConfiguration_->setParameterSetID(ParameterSet::emptyParameterSetID());
0108     processContext_->setProcessConfiguration(processConfiguration_.get());
0109
0110     SignallingProductRegistryFiller filler;
0111     if (pset.existsAs<std::vector<ParameterSet> >("producers", true)) {
0112       std::vector<ParameterSet> producers = pset.getParameter<std::vector<ParameterSet> >("producers");
0113
0114       std::vector<std::string> names;
0115       names.reserve(producers.size());
0116       std::transform(producers.begin(), producers.end(), std::back_inserter(names), [](edm::ParameterSet const& iPSet) {
0117         return iPSet.getParameter<std::string>("@module_label");
0118       });
0119       auto randomGenerator = std::make_unique<PileupRandomNumberGenerator>(names);
0120       randomGenerator_ = randomGenerator.get();
0121       std::unique_ptr<edm::RandomNumberGenerator> baseGen = std::move(randomGenerator);
0122       serviceToken_ = edm::ServiceRegistry::createContaining(
0123           std::move(baseGen), edm::ServiceRegistry::instance().presentToken(), true);
0124
0125       provider_ = std::make_unique<SecondaryEventProvider>(producers, filler, processConfiguration_);
0126     }
0127     filler.addFromInput(*input_->productRegistry());
0128     filler.setFrozen();
0129     productRegistry_ = std::make_shared<ProductRegistry>(filler.moveTo());
0130
0131     // A modified HistoryAppender must be used for unscheduled processing.
0132     eventPrincipal_ = std::make_unique<EventPrincipal>(productRegistry_,
0133                                                        edm::productResolversFactory::makePrimary,
0134                                                        std::make_shared<BranchIDListHelper>(),
0135                                                        std::make_shared<ThinnedAssociationsHelper>(),
0136                                                        *processConfiguration_,
0137                                                        nullptr);
0138
0139     bool DB = type_ == "readDB";
0140
0141     if (pset.exists("nbPileupEvents")) {
0142       if (0 != pset.getParameter<edm::ParameterSet>("nbPileupEvents").getUntrackedParameter<int>("seed", 0)) {
0143         edm::LogWarning("MixingModule") << "Parameter nbPileupEvents.seed is not supported";
0144       }
0145     }
0146
0147     edm::Service<edm::RandomNumberGenerator> rng;
0148     if (!rng.isAvailable()) {
0149       throw cms::Exception("Configuration")
0150           << "PileUp requires the RandomNumberGeneratorService\n"
0151              "which is not present in the configuration file.  You must add the service\n"
0152              "in the configuration file or remove the modules that require it.";
0153     }
0154
0155     if (!(histoDistribution_ || probFunctionDistribution_ || poisson_ || fixed_ || none_) && !DB) {
0156       throw cms::Exception("Illegal parameter value", "PileUp::PileUp(ParameterSet const& pset)")
0157           << "'type' parameter (a string) has a value of '" << type_ << "'.\n"
0158           << "Legal values are 'poisson', 'fixed', or 'none'\n";
0159     }
0160
0161     if (!DB) {
0162       manage_OOT_ = pset.getUntrackedParameter<bool>("manage_OOT", false);
0163
0164       // Check for string describing special processing.  Add these here for individual cases
0165       PU_Study_ = false;
0166       Study_type_ = pset.getUntrackedParameter<std::string>("Special_Pileup_Studies", "");
0167
0168       if (Study_type_ == "Fixed_ITPU_Vary_OOTPU") {
0169         PU_Study_ = true;
0170         intFixed_ITPU_ = pset.getUntrackedParameter<int>("intFixed_ITPU", 0);
0171       }
0172
0173       if (manage_OOT_) {  // figure out what the parameters are
0174
0175         //      if (playback_) throw cms::Exception("Illegal parameter clash","PileUp::PileUp(ParameterSet const& pset)")
0176         // << " manage_OOT option not allowed with playback ";
0177
0178         std::string OOT_type = pset.getUntrackedParameter<std::string>("OOT_type");
0179
0180         if (OOT_type == "Poisson" || OOT_type == "poisson") {
0181           poisson_OOT_ = true;
0182         } else if (OOT_type == "Fixed" || OOT_type == "fixed") {
0183           fixed_OOT_ = true;
0184           // read back the fixed number requested out-of-time
0185           intFixed_OOT_ = pset.getUntrackedParameter<int>("intFixed_OOT", -1);
0186           if (intFixed_OOT_ < 0) {
0187             throw cms::Exception("Illegal parameter value", "PileUp::PileUp(ParameterSet const& pset)")
0188                 << " Fixed out-of-time pileup requested, but no fixed value given ";
0189           }
0190         } else {
0191           throw cms::Exception("Illegal parameter value", "PileUp::PileUp(ParameterSet const& pset)")
0192               << "'OOT_type' parameter (a string) has a value of '" << OOT_type << "'.\n"
0193               << "Legal values are 'poisson' or 'fixed'\n";
0194         }
0195         edm::LogInfo("MixingModule") << " Out-of-time pileup will be generated with a " << OOT_type
0196                                      << " distribution. ";
0197       }
0198     }
0199
0200     if (Source_type_ == "cosmics") {  // allow for some extra flexibility for mixing
0201       minBunch_cosmics_ = pset.getUntrackedParameter<int>("minBunch_cosmics", -1000);
0202       maxBunch_cosmics_ = pset.getUntrackedParameter<int>("maxBunch_cosmics", 1000);
0203     }
0204   }  // end of constructor
0205
0206   void PileUp::beginJob(eventsetup::ESRecordsToProductResolverIndices const& iES) {
0207     input_->doBeginJob();
0208     if (provider_.get() != nullptr) {
0209       edm::ServiceRegistry::Operate guard(*serviceToken_);
0210       GlobalContext globalContext(GlobalContext::Transition::kBeginJob, processContext_.get());
0211       provider_->beginJob(*productRegistry_, iES, globalContext);
0212     }
0213   }
0214
0215   void PileUp::beginStream(edm::StreamID) {
0216     auto iID = eventPrincipal_->streamID();  // each producer has its own workermanager, so use default streamid
0217     streamContext_ = std::make_shared<StreamContext>(iID, processContext_.get());
0218     streamContext_->setTransition(StreamContext::Transition::kBeginStream);
0219     if (provider_.get() != nullptr) {
0220       edm::ServiceRegistry::Operate guard(*serviceToken_);
0221       provider_->beginStream(iID, *streamContext_);
0222     }
0223   }
0224
0225   void PileUp::endStream() {
0226     ExceptionCollector exceptionCollector(
0227         "Multiple exceptions were thrown while executing PileUp::endStream. An exception message follows for "
0228         "each.\n");
0229     endStream(exceptionCollector);
0230
0231     if (exceptionCollector.hasThrown()) {
0232       exceptionCollector.rethrow();
0233     }
0234   }
0235
0236   void PileUp::endStream(ExceptionCollector& exceptionCollector) {
0237     if (provider_.get() != nullptr) {
0238       edm::ServiceRegistry::Operate guard(*serviceToken_);
0239       streamContext_->setTransition(StreamContext::Transition::kEndStream);
0240       provider_->endStream(streamContext_->streamID(), *streamContext_, exceptionCollector);
0241       // This is kind of strange, end of job running as part of endStream multiple times...
0242       // For the moment, I'm leaving this as is but maybe we should think about this...
0243       // I think nothing uses this code anymore anyway...
0244       GlobalContext globalContext(GlobalContext::Transition::kEndJob, processContext_.get());
0245       provider_->endJob(exceptionCollector, globalContext);
0246     }
0247     input_->doEndJob();
0248   }
0249
0250   void PileUp::beginRun(const edm::Run& run, const edm::EventSetup& setup) {
0251     if (provider_.get() != nullptr) {
0252       runPrincipal_ = std::make_shared<RunPrincipal>(
0253           productRegistry_, edm::productResolversFactory::makePrimary, *processConfiguration_, nullptr, 0);
0254       runPrincipal_->setAux(run.runAuxiliary());
0255       edm::ServiceRegistry::Operate guard(*serviceToken_);
0256       streamContext_->setTransition(StreamContext::Transition::kBeginRun);
0257       provider_->beginRun(*runPrincipal_, setup.impl(), run.moduleCallingContext(), *streamContext_);
0258     }
0259   }
0260   void PileUp::beginLuminosityBlock(const edm::LuminosityBlock& lumi, const edm::EventSetup& setup) {
0261     if (provider_.get() != nullptr) {
0262       lumiPrincipal_ = std::make_shared<LuminosityBlockPrincipal>(
0263           productRegistry_, edm::productResolversFactory::makePrimary, *processConfiguration_, nullptr, 0);
0264       lumiPrincipal_->setAux(lumi.luminosityBlockAuxiliary());
0265       lumiPrincipal_->setRunPrincipal(runPrincipal_);
0266       setRandomEngine(lumi);
0267       edm::ServiceRegistry::Operate guard(*serviceToken_);
0268       streamContext_->setTransition(StreamContext::Transition::kBeginLuminosityBlock);
0269       provider_->beginLuminosityBlock(*lumiPrincipal_, setup.impl(), lumi.moduleCallingContext(), *streamContext_);
0270     }
0271   }
0272
0273   void PileUp::endRun(const edm::Run& run, const edm::EventSetup& setup) {
0274     if (provider_.get() != nullptr) {
0275       edm::ServiceRegistry::Operate guard(*serviceToken_);
0276       streamContext_->setTransition(StreamContext::Transition::kEndRun);
0277       provider_->endRun(*runPrincipal_, setup.impl(), run.moduleCallingContext(), *streamContext_);
0278     }
0279   }
0280   void PileUp::endLuminosityBlock(const edm::LuminosityBlock& lumi, const edm::EventSetup& setup) {
0281     if (provider_.get() != nullptr) {
0282       edm::ServiceRegistry::Operate guard(*serviceToken_);
0283       streamContext_->setTransition(StreamContext::Transition::kEndLuminosityBlock);
0284       provider_->endLuminosityBlock(*lumiPrincipal_, setup.impl(), lumi.moduleCallingContext(), *streamContext_);
0285     }
0286   }
0287
0288   void PileUp::setupPileUpEvent(const edm::EventSetup& setup) {
0289     if (provider_.get() != nullptr) {
0290       // note:  run and lumi numbers must be modified to match lumiPrincipal_
0291       eventPrincipal_->setLuminosityBlockPrincipal(lumiPrincipal_.get());
0292       eventPrincipal_->setRunAndLumiNumber(lumiPrincipal_->run(), lumiPrincipal_->luminosityBlock());
0293       edm::ServiceRegistry::Operate guard(*serviceToken_);
0294       streamContext_->setTransition(StreamContext::Transition::kEvent);
0295       provider_->setupPileUpEvent(*eventPrincipal_, setup.impl(), *streamContext_);
0296     }
0297   }
0298
0299   void PileUp::reload(const edm::EventSetup& setup) {
0300     //get the required parameters from DB.
0301     const MixingInputConfig& config = setup.getData(configToken_).config(inputType_);
0302
0303     //get the type
0304     type_ = config.type();
0305     //set booleans
0306     histoDistribution_ = type_ == "histo";
0307     probFunctionDistribution_ = type_ == "probFunction";
0308     poisson_ = type_ == "poisson";
0309     fixed_ = type_ == "fixed";
0310     none_ = type_ == "none";
0311
0312     if (histoDistribution_)
0313       edm::LogError("MisConfiguration") << "type histo cannot be reloaded from DB, yet";
0314
0315     if (fixed_) {
0316       averageNumber_ = averageNumber();
0317     } else if (poisson_) {
0318       averageNumber_ = config.averageNumber();
0319       if (PoissonDistribution_) {
0320         PoissonDistribution_ = std::make_unique<CLHEP::RandPoissonQ>(PoissonDistribution_->engine(), averageNumber_);
0321       }
0322     } else if (probFunctionDistribution_) {
0323       //need to reload the histogram from DB
0324       const std::vector<int>& dataProbFunctionVar = config.probFunctionVariable();
0325       std::vector<double> dataProb = config.probValue();
0326
0327       int varSize = (int)dataProbFunctionVar.size();
0328       int probSize = (int)dataProb.size();
0329
0330       if ((dataProbFunctionVar[0] != 0) || (dataProbFunctionVar[varSize - 1] != (varSize - 1)))
0331         throw cms::Exception("BadProbFunction")
0332             << "Please, check the variables of the probability function! The first variable should be 0 and the "
0333                "difference between two variables should be 1."
0334             << std::endl;
0335
0336       // Complete the vector containing the probability  function data
0337       // with the values "0"
0338       if (probSize < varSize) {
0339         edm::LogWarning("MixingModule") << " The probability function data will be completed with "
0340                                         << (varSize - probSize) << " values 0.";
0341
0342         for (int i = 0; i < (varSize - probSize); i++)
0343           dataProb.push_back(0);
0344
0345         probSize = dataProb.size();
0346         edm::LogInfo("MixingModule") << " The number of the P(x) data set after adding the values 0 is " << probSize;
0347       }
0348
0349       // Create an histogram with the data from the probability function provided by the user
0350       int xmin = (int)dataProbFunctionVar[0];
0351       int xmax = (int)dataProbFunctionVar[varSize - 1] + 1;  // need upper edge to be one beyond last value
0352       int numBins = varSize;
0353
0354       edm::LogInfo("MixingModule") << "An histogram will be created with " << numBins << " bins in the range (" << xmin
0355                                    << "," << xmax << ")." << std::endl;
0356
0357       histo_ = std::make_shared<TH1F>("h", "Histo from the user's probability function", numBins, xmin, xmax);
0358
0359       LogDebug("MixingModule") << "Filling histogram with the following data:" << std::endl;
0360
0361       for (int j = 0; j < numBins; j++) {
0362         LogDebug("MixingModule") << " x = " << dataProbFunctionVar[j] << " P(x) = " << dataProb[j];
0363         histo_->Fill(dataProbFunctionVar[j] + 0.5,
0364                      dataProb[j]);  // assuming integer values for the bins, fill bin centers, not edges
0365       }
0366
0367       // Check if the histogram is normalized
0368       if (std::abs(histo_->Integral() - 1) > 1.0e-02) {
0369         throw cms::Exception("BadProbFunction") << "The probability function should be normalized!!! " << std::endl;
0370       }
0371       averageNumber_ = histo_->GetMean();
0372     }
0373
0374     int oot = config.outOfTime();
0375     manage_OOT_ = false;
0376     if (oot == 1) {
0377       manage_OOT_ = true;
0378       poisson_OOT_ = false;
0379       fixed_OOT_ = true;
0380       intFixed_OOT_ = config.fixedOutOfTime();
0381     } else if (oot == 2) {
0382       manage_OOT_ = true;
0383       poisson_OOT_ = true;
0384       fixed_OOT_ = false;
0385     }
0386   }
0387   PileUp::~PileUp() {}
0388
0389   std::unique_ptr<CLHEP::RandPoissonQ> const& PileUp::poissonDistribution(StreamID const& streamID) {
0390     if (!PoissonDistribution_) {
0391       CLHEP::HepRandomEngine& engine = *randomEngine(streamID);
0392       PoissonDistribution_ = std::make_unique<CLHEP::RandPoissonQ>(engine, averageNumber_);
0393     }
0394     return PoissonDistribution_;
0395   }
0396
0397   std::unique_ptr<CLHEP::RandPoisson> const& PileUp::poissonDistr_OOT(StreamID const& streamID) {
0398     if (!PoissonDistr_OOT_) {
0399       CLHEP::HepRandomEngine& engine = *randomEngine(streamID);
0400       PoissonDistr_OOT_ = std::make_unique<CLHEP::RandPoisson>(engine);
0401     }
0402     return PoissonDistr_OOT_;
0403   }
0404
0405   CLHEP::HepRandomEngine* PileUp::randomEngine(StreamID const& streamID) {
0406     if (!randomEngine_) {
0407       Service<RandomNumberGenerator> rng;
0408       randomEngine_ = &rng->getEngine(streamID);
0409     }
0410     return randomEngine_;
0411   }
0412
0413   void PileUp::setRandomEngine(StreamID streamID) { randomGenerator_->setEngine(*randomEngine(streamID)); }
0414   void PileUp::setRandomEngine(LuminosityBlock const& iLumi) {
0415     Service<RandomNumberGenerator> rng;
0416     randomGenerator_->setSeed(rng->mySeed());
0417     randomGenerator_->setEngine(rng->getEngine(iLumi.index()));
0418   }
0419
0420   void PileUp::CalculatePileup(int MinBunch,
0421                                int MaxBunch,
0422                                std::vector<int>& PileupSelection,
0423                                std::vector<float>& TrueNumInteractions,
0424                                StreamID const& streamID) {
0425     // if we are managing the distribution of out-of-time pileup separately, select the distribution for bunch
0426     // crossing zero first, save it for later.
0427
0428     int nzero_crossing = -1;
0429     double Fnzero_crossing = -1;
0430
0431     if (provider_) {
0432       setRandomEngine(streamID);
0433     }
0434
0435     if (manage_OOT_) {
0436       if (none_) {
0437         nzero_crossing = 0;
0438       } else if (poisson_) {
0439         nzero_crossing = poissonDistribution(streamID)->fire();
0440       } else if (fixed_) {
0441         nzero_crossing = intAverage_;
0442       } else if (histoDistribution_ || probFunctionDistribution_) {
0443         // RANDOM_NUMBER_ERROR
0444         // Random number should be generated by the engines from the
0445         // RandomNumberGeneratorService. This appears to use the global
0446         // engine in ROOT. This is not thread safe unless the module using
0447         // it is a one module and declares a shared resource and all
0448         // other modules using it also declare the same shared resource.
0449         // This also breaks replay.
0450         double d = GetRandom(histo_.get(), randomEngine(streamID));
0451         //n = (int) floor(d + 0.5);  // incorrect for bins with integer edges
0452         Fnzero_crossing = d;
0453         nzero_crossing = int(d);
0454       }
0455     }
0456
0457     for (int bx = MinBunch; bx < MaxBunch + 1; ++bx) {
0458       if (manage_OOT_) {
0459         if (bx == 0 && !poisson_OOT_) {
0460           PileupSelection.push_back(nzero_crossing);
0461           TrueNumInteractions.push_back(nzero_crossing);
0462         } else {
0463           if (poisson_OOT_) {
0464             if (PU_Study_ && (Study_type_ == "Fixed_ITPU_Vary_OOTPU") && bx == 0) {
0465               PileupSelection.push_back(intFixed_ITPU_);
0466             } else {
0467               PileupSelection.push_back(poissonDistr_OOT(streamID)->fire(Fnzero_crossing));
0468             }
0469             TrueNumInteractions.push_back(Fnzero_crossing);
0470           } else {
0471             PileupSelection.push_back(intFixed_OOT_);
0472             TrueNumInteractions.push_back(intFixed_OOT_);
0473           }
0474         }
0475       } else {
0476         if (none_) {
0477           PileupSelection.push_back(0);
0478           TrueNumInteractions.push_back(0.);
0479         } else if (poisson_) {
0480           PileupSelection.push_back(poissonDistribution(streamID)->fire());
0481           TrueNumInteractions.push_back(averageNumber_);
0482         } else if (fixed_) {
0483           PileupSelection.push_back(intAverage_);
0484           TrueNumInteractions.push_back(intAverage_);
0485         } else if (histoDistribution_ || probFunctionDistribution_) {
0486           // RANDOM_NUMBER_ERROR
0487           // Random number should be generated by the engines from the
0488           // RandomNumberGeneratorService. This appears to use the global
0489           // engine in ROOT. This is not thread safe unless the module using
0490           // it is a one module and declares a shared resource and all
0491           // other modules using it also declare the same shared resource.
0492           // This also breaks replay.
0493           double d = GetRandom(histo_.get(), randomEngine(streamID));
0494           PileupSelection.push_back(int(d));
0495           TrueNumInteractions.push_back(d);
0496         }
0497       }
0498     }
0499   }
0500
0501 }  //namespace edm