File indexing completed on 2021-09-02 03:45:54
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/SignallingProductRegistry.h"
0009 #include "FWCore/Framework/interface/ESRecordsToProxyIndices.h"
0010 #include "FWCore/ServiceRegistry/interface/ActivityRegistry.h"
0011 #include "FWCore/ServiceRegistry/interface/ProcessContext.h"
0012 #include "FWCore/Sources/interface/VectorInputSourceDescription.h"
0013 #include "FWCore/Sources/interface/VectorInputSourceFactory.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 #include "FWCore/Utilities/interface/Exception.h"
0017 #include "FWCore/Utilities/interface/GetPassID.h"
0018 #include "FWCore/Utilities/interface/StreamID.h"
0019 #include "FWCore/Version/interface/GetReleaseVersion.h"
0020
0021 #include "FWCore/ServiceRegistry/interface/Service.h"
0022 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0023
0024 #include "FWCore/Framework/interface/ESHandle.h"
0025 #include "FWCore/Framework/interface/EventSetup.h"
0026 #include "Mixing/Base/src/SecondaryEventProvider.h"
0027 #include "CondFormats/DataRecord/interface/MixingRcd.h"
0028 #include "CondFormats/RunInfo/interface/MixingModuleConfig.h"
0029
0030 #include "CLHEP/Random/RandPoissonQ.h"
0031 #include "CLHEP/Random/RandPoisson.h"
0032
0033 #include <algorithm>
0034 #include <memory>
0035 #include "TMath.h"
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049 static Double_t GetRandom(TH1* th1, CLHEP::HepRandomEngine* rng) {
0050 Int_t nbinsx = th1->GetNbinsX();
0051 Double_t* fIntegral = th1->GetIntegral();
0052 Double_t integral = fIntegral[nbinsx];
0053
0054 if (integral == 0)
0055 return 0;
0056
0057 Double_t r1 = rng->flat();
0058 Int_t ibin = TMath::BinarySearch(nbinsx, fIntegral, r1);
0059 Double_t x = th1->GetBinLowEdge(ibin + 1);
0060 if (r1 > fIntegral[ibin])
0061 x += th1->GetBinWidth(ibin + 1) * (r1 - fIntegral[ibin]) / (fIntegral[ibin + 1] - fIntegral[ibin]);
0062 return x;
0063 }
0064
0065
0066 namespace edm {
0067 PileUp::PileUp(ParameterSet const& pset,
0068 const std::shared_ptr<PileUpConfig>& config,
0069 edm::ConsumesCollector iC,
0070 const bool mixingConfigFromDB)
0071 : type_(pset.getParameter<std::string>("type")),
0072 Source_type_(config->sourcename_),
0073 averageNumber_(config->averageNumber_),
0074 intAverage_(static_cast<int>(averageNumber_)),
0075 histo_(std::make_shared<TH1F>(*config->histo_)),
0076 histoDistribution_(type_ == "histo"),
0077 probFunctionDistribution_(type_ == "probFunction"),
0078 poisson_(type_ == "poisson"),
0079 fixed_(type_ == "fixed"),
0080 none_(type_ == "none"),
0081 fileNameHash_(0U),
0082 productRegistry_(new SignallingProductRegistry),
0083 input_(VectorInputSourceFactory::get()
0084 ->makeVectorInputSource(
0085 pset, VectorInputSourceDescription(productRegistry_, edm::PreallocationConfiguration()))
0086 .release()),
0087 processConfiguration_(new ProcessConfiguration(std::string("@MIXING"), getReleaseVersion(), getPassID())),
0088 processContext_(new ProcessContext()),
0089 eventPrincipal_(),
0090 lumiPrincipal_(),
0091 runPrincipal_(),
0092 provider_(),
0093 PoissonDistribution_(),
0094 PoissonDistr_OOT_(),
0095 randomEngine_(),
0096 playback_(config->playback_),
0097 sequential_(pset.getUntrackedParameter<bool>("sequential", false)) {
0098 if (mixingConfigFromDB) {
0099 configToken_ = iC.esConsumes<edm::Transition::BeginLuminosityBlock>();
0100 }
0101
0102
0103 processConfiguration_->setParameterSetID(ParameterSet::emptyParameterSetID());
0104 processContext_->setProcessConfiguration(processConfiguration_.get());
0105
0106 if (pset.existsAs<std::vector<ParameterSet> >("producers", true)) {
0107 std::vector<ParameterSet> producers = pset.getParameter<std::vector<ParameterSet> >("producers");
0108 provider_ = std::make_unique<SecondaryEventProvider>(producers, *productRegistry_, processConfiguration_);
0109 }
0110
0111 productRegistry_->setFrozen();
0112
0113
0114 eventPrincipal_ = std::make_unique<EventPrincipal>(input_->productRegistry(),
0115 std::make_shared<BranchIDListHelper>(),
0116 std::make_shared<ThinnedAssociationsHelper>(),
0117 *processConfiguration_,
0118 nullptr);
0119
0120 bool DB = type_ == "readDB";
0121
0122 if (pset.exists("nbPileupEvents")) {
0123 if (0 != pset.getParameter<edm::ParameterSet>("nbPileupEvents").getUntrackedParameter<int>("seed", 0)) {
0124 edm::LogWarning("MixingModule") << "Parameter nbPileupEvents.seed is not supported";
0125 }
0126 }
0127
0128 edm::Service<edm::RandomNumberGenerator> rng;
0129 if (!rng.isAvailable()) {
0130 throw cms::Exception("Configuration")
0131 << "PileUp requires the RandomNumberGeneratorService\n"
0132 "which is not present in the configuration file. You must add the service\n"
0133 "in the configuration file or remove the modules that require it.";
0134 }
0135
0136 if (!(histoDistribution_ || probFunctionDistribution_ || poisson_ || fixed_ || none_) && !DB) {
0137 throw cms::Exception("Illegal parameter value", "PileUp::PileUp(ParameterSet const& pset)")
0138 << "'type' parameter (a string) has a value of '" << type_ << "'.\n"
0139 << "Legal values are 'poisson', 'fixed', or 'none'\n";
0140 }
0141
0142 if (!DB) {
0143 manage_OOT_ = pset.getUntrackedParameter<bool>("manage_OOT", false);
0144
0145
0146 PU_Study_ = false;
0147 Study_type_ = pset.getUntrackedParameter<std::string>("Special_Pileup_Studies", "");
0148
0149 if (Study_type_ == "Fixed_ITPU_Vary_OOTPU") {
0150 PU_Study_ = true;
0151 intFixed_ITPU_ = pset.getUntrackedParameter<int>("intFixed_ITPU", 0);
0152 }
0153
0154 if (manage_OOT_) {
0155
0156
0157
0158
0159 std::string OOT_type = pset.getUntrackedParameter<std::string>("OOT_type");
0160
0161 if (OOT_type == "Poisson" || OOT_type == "poisson") {
0162 poisson_OOT_ = true;
0163 } else if (OOT_type == "Fixed" || OOT_type == "fixed") {
0164 fixed_OOT_ = true;
0165
0166 intFixed_OOT_ = pset.getUntrackedParameter<int>("intFixed_OOT", -1);
0167 if (intFixed_OOT_ < 0) {
0168 throw cms::Exception("Illegal parameter value", "PileUp::PileUp(ParameterSet const& pset)")
0169 << " Fixed out-of-time pileup requested, but no fixed value given ";
0170 }
0171 } else {
0172 throw cms::Exception("Illegal parameter value", "PileUp::PileUp(ParameterSet const& pset)")
0173 << "'OOT_type' parameter (a string) has a value of '" << OOT_type << "'.\n"
0174 << "Legal values are 'poisson' or 'fixed'\n";
0175 }
0176 edm::LogInfo("MixingModule") << " Out-of-time pileup will be generated with a " << OOT_type
0177 << " distribution. ";
0178 }
0179 }
0180
0181 if (Source_type_ == "cosmics") {
0182 minBunch_cosmics_ = pset.getUntrackedParameter<int>("minBunch_cosmics", -1000);
0183 maxBunch_cosmics_ = pset.getUntrackedParameter<int>("maxBunch_cosmics", 1000);
0184 }
0185 }
0186
0187 void PileUp::beginStream(edm::StreamID) {
0188 auto iID = eventPrincipal_->streamID();
0189 streamContext_.reset(new StreamContext(iID, processContext_.get()));
0190 input_->doBeginJob();
0191 if (provider_.get() != nullptr) {
0192
0193 provider_->beginJob(*productRegistry_, eventsetup::ESRecordsToProxyIndices{{}});
0194 provider_->beginStream(iID, *streamContext_);
0195 }
0196 }
0197
0198 void PileUp::endStream() {
0199 if (provider_.get() != nullptr) {
0200 provider_->endStream(streamContext_->streamID(), *streamContext_);
0201 provider_->endJob();
0202 }
0203 input_->doEndJob();
0204 }
0205
0206 void PileUp::beginRun(const edm::Run& run, const edm::EventSetup& setup) {
0207 if (provider_.get() != nullptr) {
0208 auto aux = std::make_shared<RunAuxiliary>(run.runAuxiliary());
0209 runPrincipal_.reset(new RunPrincipal(aux, productRegistry_, *processConfiguration_, nullptr, 0));
0210 provider_->beginRun(*runPrincipal_, setup.impl(), run.moduleCallingContext(), *streamContext_);
0211 }
0212 }
0213 void PileUp::beginLuminosityBlock(const edm::LuminosityBlock& lumi, const edm::EventSetup& setup) {
0214 if (provider_.get() != nullptr) {
0215 lumiPrincipal_.reset(new LuminosityBlockPrincipal(productRegistry_, *processConfiguration_, nullptr, 0));
0216 lumiPrincipal_->setAux(lumi.luminosityBlockAuxiliary());
0217 lumiPrincipal_->setRunPrincipal(runPrincipal_);
0218 provider_->beginLuminosityBlock(*lumiPrincipal_, setup.impl(), lumi.moduleCallingContext(), *streamContext_);
0219 }
0220 }
0221
0222 void PileUp::endRun(const edm::Run& run, const edm::EventSetup& setup) {
0223 if (provider_.get() != nullptr) {
0224 provider_->endRun(*runPrincipal_, setup.impl(), run.moduleCallingContext(), *streamContext_);
0225 }
0226 }
0227 void PileUp::endLuminosityBlock(const edm::LuminosityBlock& lumi, const edm::EventSetup& setup) {
0228 if (provider_.get() != nullptr) {
0229 provider_->endLuminosityBlock(*lumiPrincipal_, setup.impl(), lumi.moduleCallingContext(), *streamContext_);
0230 }
0231 }
0232
0233 void PileUp::setupPileUpEvent(const edm::EventSetup& setup) {
0234 if (provider_.get() != nullptr) {
0235
0236 eventPrincipal_->setLuminosityBlockPrincipal(lumiPrincipal_.get());
0237 eventPrincipal_->setRunAndLumiNumber(lumiPrincipal_->run(), lumiPrincipal_->luminosityBlock());
0238 provider_->setupPileUpEvent(*eventPrincipal_, setup.impl(), *streamContext_);
0239 }
0240 }
0241
0242 void PileUp::reload(const edm::EventSetup& setup) {
0243
0244 const MixingInputConfig& config = setup.getData(configToken_).config(inputType_);
0245
0246
0247 type_ = config.type();
0248
0249 histoDistribution_ = type_ == "histo";
0250 probFunctionDistribution_ = type_ == "probFunction";
0251 poisson_ = type_ == "poisson";
0252 fixed_ = type_ == "fixed";
0253 none_ = type_ == "none";
0254
0255 if (histoDistribution_)
0256 edm::LogError("MisConfiguration") << "type histo cannot be reloaded from DB, yet";
0257
0258 if (fixed_) {
0259 averageNumber_ = averageNumber();
0260 } else if (poisson_) {
0261 averageNumber_ = config.averageNumber();
0262 if (PoissonDistribution_) {
0263 PoissonDistribution_ = std::make_unique<CLHEP::RandPoissonQ>(PoissonDistribution_->engine(), averageNumber_);
0264 }
0265 } else if (probFunctionDistribution_) {
0266
0267 const std::vector<int>& dataProbFunctionVar = config.probFunctionVariable();
0268 std::vector<double> dataProb = config.probValue();
0269
0270 int varSize = (int)dataProbFunctionVar.size();
0271 int probSize = (int)dataProb.size();
0272
0273 if ((dataProbFunctionVar[0] != 0) || (dataProbFunctionVar[varSize - 1] != (varSize - 1)))
0274 throw cms::Exception("BadProbFunction")
0275 << "Please, check the variables of the probability function! The first variable should be 0 and the "
0276 "difference between two variables should be 1."
0277 << std::endl;
0278
0279
0280
0281 if (probSize < varSize) {
0282 edm::LogWarning("MixingModule") << " The probability function data will be completed with "
0283 << (varSize - probSize) << " values 0.";
0284
0285 for (int i = 0; i < (varSize - probSize); i++)
0286 dataProb.push_back(0);
0287
0288 probSize = dataProb.size();
0289 edm::LogInfo("MixingModule") << " The number of the P(x) data set after adding the values 0 is " << probSize;
0290 }
0291
0292
0293 int xmin = (int)dataProbFunctionVar[0];
0294 int xmax = (int)dataProbFunctionVar[varSize - 1] + 1;
0295 int numBins = varSize;
0296
0297 edm::LogInfo("MixingModule") << "An histogram will be created with " << numBins << " bins in the range (" << xmin
0298 << "," << xmax << ")." << std::endl;
0299
0300 histo_.reset(new TH1F("h", "Histo from the user's probability function", numBins, xmin, xmax));
0301
0302 LogDebug("MixingModule") << "Filling histogram with the following data:" << std::endl;
0303
0304 for (int j = 0; j < numBins; j++) {
0305 LogDebug("MixingModule") << " x = " << dataProbFunctionVar[j] << " P(x) = " << dataProb[j];
0306 histo_->Fill(dataProbFunctionVar[j] + 0.5,
0307 dataProb[j]);
0308 }
0309
0310
0311 if (std::abs(histo_->Integral() - 1) > 1.0e-02) {
0312 throw cms::Exception("BadProbFunction") << "The probability function should be normalized!!! " << std::endl;
0313 }
0314 averageNumber_ = histo_->GetMean();
0315 }
0316
0317 int oot = config.outOfTime();
0318 manage_OOT_ = false;
0319 if (oot == 1) {
0320 manage_OOT_ = true;
0321 poisson_OOT_ = false;
0322 fixed_OOT_ = true;
0323 intFixed_OOT_ = config.fixedOutOfTime();
0324 } else if (oot == 2) {
0325 manage_OOT_ = true;
0326 poisson_OOT_ = true;
0327 fixed_OOT_ = false;
0328 }
0329 }
0330 PileUp::~PileUp() {}
0331
0332 std::unique_ptr<CLHEP::RandPoissonQ> const& PileUp::poissonDistribution(StreamID const& streamID) {
0333 if (!PoissonDistribution_) {
0334 CLHEP::HepRandomEngine& engine = *randomEngine(streamID);
0335 PoissonDistribution_ = std::make_unique<CLHEP::RandPoissonQ>(engine, averageNumber_);
0336 }
0337 return PoissonDistribution_;
0338 }
0339
0340 std::unique_ptr<CLHEP::RandPoisson> const& PileUp::poissonDistr_OOT(StreamID const& streamID) {
0341 if (!PoissonDistr_OOT_) {
0342 CLHEP::HepRandomEngine& engine = *randomEngine(streamID);
0343 PoissonDistr_OOT_ = std::make_unique<CLHEP::RandPoisson>(engine);
0344 }
0345 return PoissonDistr_OOT_;
0346 }
0347
0348 CLHEP::HepRandomEngine* PileUp::randomEngine(StreamID const& streamID) {
0349 if (!randomEngine_) {
0350 Service<RandomNumberGenerator> rng;
0351 randomEngine_ = &rng->getEngine(streamID);
0352 }
0353 return randomEngine_;
0354 }
0355
0356 void PileUp::CalculatePileup(int MinBunch,
0357 int MaxBunch,
0358 std::vector<int>& PileupSelection,
0359 std::vector<float>& TrueNumInteractions,
0360 StreamID const& streamID) {
0361
0362
0363
0364 int nzero_crossing = -1;
0365 double Fnzero_crossing = -1;
0366
0367 if (manage_OOT_) {
0368 if (none_) {
0369 nzero_crossing = 0;
0370 } else if (poisson_) {
0371 nzero_crossing = poissonDistribution(streamID)->fire();
0372 } else if (fixed_) {
0373 nzero_crossing = intAverage_;
0374 } else if (histoDistribution_ || probFunctionDistribution_) {
0375
0376
0377
0378
0379
0380
0381
0382 double d = GetRandom(histo_.get(), randomEngine(streamID));
0383
0384 Fnzero_crossing = d;
0385 nzero_crossing = int(d);
0386 }
0387 }
0388
0389 for (int bx = MinBunch; bx < MaxBunch + 1; ++bx) {
0390 if (manage_OOT_) {
0391 if (bx == 0 && !poisson_OOT_) {
0392 PileupSelection.push_back(nzero_crossing);
0393 TrueNumInteractions.push_back(nzero_crossing);
0394 } else {
0395 if (poisson_OOT_) {
0396 if (PU_Study_ && (Study_type_ == "Fixed_ITPU_Vary_OOTPU") && bx == 0) {
0397 PileupSelection.push_back(intFixed_ITPU_);
0398 } else {
0399 PileupSelection.push_back(poissonDistr_OOT(streamID)->fire(Fnzero_crossing));
0400 }
0401 TrueNumInteractions.push_back(Fnzero_crossing);
0402 } else {
0403 PileupSelection.push_back(intFixed_OOT_);
0404 TrueNumInteractions.push_back(intFixed_OOT_);
0405 }
0406 }
0407 } else {
0408 if (none_) {
0409 PileupSelection.push_back(0);
0410 TrueNumInteractions.push_back(0.);
0411 } else if (poisson_) {
0412 PileupSelection.push_back(poissonDistribution(streamID)->fire());
0413 TrueNumInteractions.push_back(averageNumber_);
0414 } else if (fixed_) {
0415 PileupSelection.push_back(intAverage_);
0416 TrueNumInteractions.push_back(intAverage_);
0417 } else if (histoDistribution_ || probFunctionDistribution_) {
0418
0419
0420
0421
0422
0423
0424
0425 double d = GetRandom(histo_.get(), randomEngine(streamID));
0426 PileupSelection.push_back(int(d));
0427 TrueNumInteractions.push_back(d);
0428 }
0429 }
0430 }
0431 }
0432
0433 }