File indexing completed on 2024-07-28 22:48:42
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/ESRecordsToProductResolverIndices.h"
0010 #include "FWCore/ServiceRegistry/interface/ActivityRegistry.h"
0011 #include "FWCore/ServiceRegistry/interface/GlobalContext.h"
0012 #include "FWCore/ServiceRegistry/interface/ProcessContext.h"
0013 #include "FWCore/Sources/interface/VectorInputSourceDescription.h"
0014 #include "FWCore/Sources/interface/VectorInputSourceFactory.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017 #include "FWCore/Utilities/interface/Exception.h"
0018 #include "FWCore/Utilities/interface/GetPassID.h"
0019 #include "FWCore/Utilities/interface/StreamID.h"
0020 #include "FWCore/Version/interface/GetReleaseVersion.h"
0021
0022 #include "FWCore/ServiceRegistry/interface/Service.h"
0023 #include "FWCore/Utilities/interface/RandomNumberGenerator.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
0042
0043
0044
0045
0046
0047
0048
0049
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_(new SignallingProductRegistry),
0086 input_(VectorInputSourceFactory::get()
0087 ->makeVectorInputSource(
0088 pset, VectorInputSourceDescription(productRegistry_, edm::PreallocationConfiguration()))
0089 .release()),
0090 processConfiguration_(new ProcessConfiguration(std::string("@MIXING"), getReleaseVersion(), getPassID())),
0091 processContext_(new ProcessContext()),
0092 eventPrincipal_(),
0093 lumiPrincipal_(),
0094 runPrincipal_(),
0095 provider_(),
0096 PoissonDistribution_(),
0097 PoissonDistr_OOT_(),
0098 randomEngine_(),
0099 playback_(config->playback_),
0100 sequential_(pset.getUntrackedParameter<bool>("sequential", false)) {
0101 if (mixingConfigFromDB) {
0102 configToken_ = iC.esConsumes<edm::Transition::BeginLuminosityBlock>();
0103 }
0104
0105
0106 processConfiguration_->setParameterSetID(ParameterSet::emptyParameterSetID());
0107 processContext_->setProcessConfiguration(processConfiguration_.get());
0108
0109 if (pset.existsAs<std::vector<ParameterSet> >("producers", true)) {
0110 std::vector<ParameterSet> producers = pset.getParameter<std::vector<ParameterSet> >("producers");
0111
0112 std::vector<std::string> names;
0113 names.reserve(producers.size());
0114 std::transform(producers.begin(), producers.end(), std::back_inserter(names), [](edm::ParameterSet const& iPSet) {
0115 return iPSet.getParameter<std::string>("@module_label");
0116 });
0117 auto randomGenerator = std::make_unique<PileupRandomNumberGenerator>(names);
0118 randomGenerator_ = randomGenerator.get();
0119 std::unique_ptr<edm::RandomNumberGenerator> baseGen = std::move(randomGenerator);
0120 serviceToken_ = edm::ServiceRegistry::createContaining(
0121 std::move(baseGen), edm::ServiceRegistry::instance().presentToken(), true);
0122
0123 provider_ = std::make_unique<SecondaryEventProvider>(producers, *productRegistry_, processConfiguration_);
0124 }
0125
0126 productRegistry_->setFrozen();
0127
0128
0129 eventPrincipal_ = std::make_unique<EventPrincipal>(input_->productRegistry(),
0130 std::make_shared<BranchIDListHelper>(),
0131 std::make_shared<ThinnedAssociationsHelper>(),
0132 *processConfiguration_,
0133 nullptr);
0134
0135 bool DB = type_ == "readDB";
0136
0137 if (pset.exists("nbPileupEvents")) {
0138 if (0 != pset.getParameter<edm::ParameterSet>("nbPileupEvents").getUntrackedParameter<int>("seed", 0)) {
0139 edm::LogWarning("MixingModule") << "Parameter nbPileupEvents.seed is not supported";
0140 }
0141 }
0142
0143 edm::Service<edm::RandomNumberGenerator> rng;
0144 if (!rng.isAvailable()) {
0145 throw cms::Exception("Configuration")
0146 << "PileUp requires the RandomNumberGeneratorService\n"
0147 "which is not present in the configuration file. You must add the service\n"
0148 "in the configuration file or remove the modules that require it.";
0149 }
0150
0151 if (!(histoDistribution_ || probFunctionDistribution_ || poisson_ || fixed_ || none_) && !DB) {
0152 throw cms::Exception("Illegal parameter value", "PileUp::PileUp(ParameterSet const& pset)")
0153 << "'type' parameter (a string) has a value of '" << type_ << "'.\n"
0154 << "Legal values are 'poisson', 'fixed', or 'none'\n";
0155 }
0156
0157 if (!DB) {
0158 manage_OOT_ = pset.getUntrackedParameter<bool>("manage_OOT", false);
0159
0160
0161 PU_Study_ = false;
0162 Study_type_ = pset.getUntrackedParameter<std::string>("Special_Pileup_Studies", "");
0163
0164 if (Study_type_ == "Fixed_ITPU_Vary_OOTPU") {
0165 PU_Study_ = true;
0166 intFixed_ITPU_ = pset.getUntrackedParameter<int>("intFixed_ITPU", 0);
0167 }
0168
0169 if (manage_OOT_) {
0170
0171
0172
0173
0174 std::string OOT_type = pset.getUntrackedParameter<std::string>("OOT_type");
0175
0176 if (OOT_type == "Poisson" || OOT_type == "poisson") {
0177 poisson_OOT_ = true;
0178 } else if (OOT_type == "Fixed" || OOT_type == "fixed") {
0179 fixed_OOT_ = true;
0180
0181 intFixed_OOT_ = pset.getUntrackedParameter<int>("intFixed_OOT", -1);
0182 if (intFixed_OOT_ < 0) {
0183 throw cms::Exception("Illegal parameter value", "PileUp::PileUp(ParameterSet const& pset)")
0184 << " Fixed out-of-time pileup requested, but no fixed value given ";
0185 }
0186 } else {
0187 throw cms::Exception("Illegal parameter value", "PileUp::PileUp(ParameterSet const& pset)")
0188 << "'OOT_type' parameter (a string) has a value of '" << OOT_type << "'.\n"
0189 << "Legal values are 'poisson' or 'fixed'\n";
0190 }
0191 edm::LogInfo("MixingModule") << " Out-of-time pileup will be generated with a " << OOT_type
0192 << " distribution. ";
0193 }
0194 }
0195
0196 if (Source_type_ == "cosmics") {
0197 minBunch_cosmics_ = pset.getUntrackedParameter<int>("minBunch_cosmics", -1000);
0198 maxBunch_cosmics_ = pset.getUntrackedParameter<int>("maxBunch_cosmics", 1000);
0199 }
0200 }
0201
0202 void PileUp::beginJob(eventsetup::ESRecordsToProductResolverIndices const& iES) {
0203 input_->doBeginJob();
0204 if (provider_.get() != nullptr) {
0205 edm::ServiceRegistry::Operate guard(*serviceToken_);
0206 GlobalContext globalContext(GlobalContext::Transition::kBeginJob, processContext_.get());
0207 provider_->beginJob(*productRegistry_, iES, globalContext);
0208 }
0209 }
0210
0211 void PileUp::beginStream(edm::StreamID) {
0212 auto iID = eventPrincipal_->streamID();
0213 streamContext_.reset(new StreamContext(iID, processContext_.get()));
0214 streamContext_->setTransition(StreamContext::Transition::kBeginStream);
0215 if (provider_.get() != nullptr) {
0216 edm::ServiceRegistry::Operate guard(*serviceToken_);
0217 provider_->beginStream(iID, *streamContext_);
0218 }
0219 }
0220
0221 void PileUp::endStream() {
0222 ExceptionCollector exceptionCollector(
0223 "Multiple exceptions were thrown while executing PileUp::endStream. An exception message follows for "
0224 "each.\n");
0225 endStream(exceptionCollector);
0226
0227 if (exceptionCollector.hasThrown()) {
0228 exceptionCollector.rethrow();
0229 }
0230 }
0231
0232 void PileUp::endStream(ExceptionCollector& exceptionCollector) {
0233 if (provider_.get() != nullptr) {
0234 edm::ServiceRegistry::Operate guard(*serviceToken_);
0235 streamContext_->setTransition(StreamContext::Transition::kEndStream);
0236 provider_->endStream(streamContext_->streamID(), *streamContext_, exceptionCollector);
0237
0238
0239
0240 GlobalContext globalContext(GlobalContext::Transition::kEndJob, processContext_.get());
0241 provider_->endJob(exceptionCollector, globalContext);
0242 }
0243 input_->doEndJob();
0244 }
0245
0246 void PileUp::beginRun(const edm::Run& run, const edm::EventSetup& setup) {
0247 if (provider_.get() != nullptr) {
0248 runPrincipal_.reset(new RunPrincipal(productRegistry_, *processConfiguration_, nullptr, 0));
0249 runPrincipal_->setAux(run.runAuxiliary());
0250 edm::ServiceRegistry::Operate guard(*serviceToken_);
0251 streamContext_->setTransition(StreamContext::Transition::kBeginRun);
0252 provider_->beginRun(*runPrincipal_, setup.impl(), run.moduleCallingContext(), *streamContext_);
0253 }
0254 }
0255 void PileUp::beginLuminosityBlock(const edm::LuminosityBlock& lumi, const edm::EventSetup& setup) {
0256 if (provider_.get() != nullptr) {
0257 lumiPrincipal_.reset(new LuminosityBlockPrincipal(productRegistry_, *processConfiguration_, nullptr, 0));
0258 lumiPrincipal_->setAux(lumi.luminosityBlockAuxiliary());
0259 lumiPrincipal_->setRunPrincipal(runPrincipal_);
0260 setRandomEngine(lumi);
0261 edm::ServiceRegistry::Operate guard(*serviceToken_);
0262 streamContext_->setTransition(StreamContext::Transition::kBeginLuminosityBlock);
0263 provider_->beginLuminosityBlock(*lumiPrincipal_, setup.impl(), lumi.moduleCallingContext(), *streamContext_);
0264 }
0265 }
0266
0267 void PileUp::endRun(const edm::Run& run, const edm::EventSetup& setup) {
0268 if (provider_.get() != nullptr) {
0269 edm::ServiceRegistry::Operate guard(*serviceToken_);
0270 streamContext_->setTransition(StreamContext::Transition::kEndRun);
0271 provider_->endRun(*runPrincipal_, setup.impl(), run.moduleCallingContext(), *streamContext_);
0272 }
0273 }
0274 void PileUp::endLuminosityBlock(const edm::LuminosityBlock& lumi, const edm::EventSetup& setup) {
0275 if (provider_.get() != nullptr) {
0276 edm::ServiceRegistry::Operate guard(*serviceToken_);
0277 streamContext_->setTransition(StreamContext::Transition::kEndLuminosityBlock);
0278 provider_->endLuminosityBlock(*lumiPrincipal_, setup.impl(), lumi.moduleCallingContext(), *streamContext_);
0279 }
0280 }
0281
0282 void PileUp::setupPileUpEvent(const edm::EventSetup& setup) {
0283 if (provider_.get() != nullptr) {
0284
0285 eventPrincipal_->setLuminosityBlockPrincipal(lumiPrincipal_.get());
0286 eventPrincipal_->setRunAndLumiNumber(lumiPrincipal_->run(), lumiPrincipal_->luminosityBlock());
0287 edm::ServiceRegistry::Operate guard(*serviceToken_);
0288 streamContext_->setTransition(StreamContext::Transition::kEvent);
0289 provider_->setupPileUpEvent(*eventPrincipal_, setup.impl(), *streamContext_);
0290 }
0291 }
0292
0293 void PileUp::reload(const edm::EventSetup& setup) {
0294
0295 const MixingInputConfig& config = setup.getData(configToken_).config(inputType_);
0296
0297
0298 type_ = config.type();
0299
0300 histoDistribution_ = type_ == "histo";
0301 probFunctionDistribution_ = type_ == "probFunction";
0302 poisson_ = type_ == "poisson";
0303 fixed_ = type_ == "fixed";
0304 none_ = type_ == "none";
0305
0306 if (histoDistribution_)
0307 edm::LogError("MisConfiguration") << "type histo cannot be reloaded from DB, yet";
0308
0309 if (fixed_) {
0310 averageNumber_ = averageNumber();
0311 } else if (poisson_) {
0312 averageNumber_ = config.averageNumber();
0313 if (PoissonDistribution_) {
0314 PoissonDistribution_ = std::make_unique<CLHEP::RandPoissonQ>(PoissonDistribution_->engine(), averageNumber_);
0315 }
0316 } else if (probFunctionDistribution_) {
0317
0318 const std::vector<int>& dataProbFunctionVar = config.probFunctionVariable();
0319 std::vector<double> dataProb = config.probValue();
0320
0321 int varSize = (int)dataProbFunctionVar.size();
0322 int probSize = (int)dataProb.size();
0323
0324 if ((dataProbFunctionVar[0] != 0) || (dataProbFunctionVar[varSize - 1] != (varSize - 1)))
0325 throw cms::Exception("BadProbFunction")
0326 << "Please, check the variables of the probability function! The first variable should be 0 and the "
0327 "difference between two variables should be 1."
0328 << std::endl;
0329
0330
0331
0332 if (probSize < varSize) {
0333 edm::LogWarning("MixingModule") << " The probability function data will be completed with "
0334 << (varSize - probSize) << " values 0.";
0335
0336 for (int i = 0; i < (varSize - probSize); i++)
0337 dataProb.push_back(0);
0338
0339 probSize = dataProb.size();
0340 edm::LogInfo("MixingModule") << " The number of the P(x) data set after adding the values 0 is " << probSize;
0341 }
0342
0343
0344 int xmin = (int)dataProbFunctionVar[0];
0345 int xmax = (int)dataProbFunctionVar[varSize - 1] + 1;
0346 int numBins = varSize;
0347
0348 edm::LogInfo("MixingModule") << "An histogram will be created with " << numBins << " bins in the range (" << xmin
0349 << "," << xmax << ")." << std::endl;
0350
0351 histo_.reset(new TH1F("h", "Histo from the user's probability function", numBins, xmin, xmax));
0352
0353 LogDebug("MixingModule") << "Filling histogram with the following data:" << std::endl;
0354
0355 for (int j = 0; j < numBins; j++) {
0356 LogDebug("MixingModule") << " x = " << dataProbFunctionVar[j] << " P(x) = " << dataProb[j];
0357 histo_->Fill(dataProbFunctionVar[j] + 0.5,
0358 dataProb[j]);
0359 }
0360
0361
0362 if (std::abs(histo_->Integral() - 1) > 1.0e-02) {
0363 throw cms::Exception("BadProbFunction") << "The probability function should be normalized!!! " << std::endl;
0364 }
0365 averageNumber_ = histo_->GetMean();
0366 }
0367
0368 int oot = config.outOfTime();
0369 manage_OOT_ = false;
0370 if (oot == 1) {
0371 manage_OOT_ = true;
0372 poisson_OOT_ = false;
0373 fixed_OOT_ = true;
0374 intFixed_OOT_ = config.fixedOutOfTime();
0375 } else if (oot == 2) {
0376 manage_OOT_ = true;
0377 poisson_OOT_ = true;
0378 fixed_OOT_ = false;
0379 }
0380 }
0381 PileUp::~PileUp() {}
0382
0383 std::unique_ptr<CLHEP::RandPoissonQ> const& PileUp::poissonDistribution(StreamID const& streamID) {
0384 if (!PoissonDistribution_) {
0385 CLHEP::HepRandomEngine& engine = *randomEngine(streamID);
0386 PoissonDistribution_ = std::make_unique<CLHEP::RandPoissonQ>(engine, averageNumber_);
0387 }
0388 return PoissonDistribution_;
0389 }
0390
0391 std::unique_ptr<CLHEP::RandPoisson> const& PileUp::poissonDistr_OOT(StreamID const& streamID) {
0392 if (!PoissonDistr_OOT_) {
0393 CLHEP::HepRandomEngine& engine = *randomEngine(streamID);
0394 PoissonDistr_OOT_ = std::make_unique<CLHEP::RandPoisson>(engine);
0395 }
0396 return PoissonDistr_OOT_;
0397 }
0398
0399 CLHEP::HepRandomEngine* PileUp::randomEngine(StreamID const& streamID) {
0400 if (!randomEngine_) {
0401 Service<RandomNumberGenerator> rng;
0402 randomEngine_ = &rng->getEngine(streamID);
0403 }
0404 return randomEngine_;
0405 }
0406
0407 void PileUp::setRandomEngine(StreamID streamID) { randomGenerator_->setEngine(*randomEngine(streamID)); }
0408 void PileUp::setRandomEngine(LuminosityBlock const& iLumi) {
0409 Service<RandomNumberGenerator> rng;
0410 randomGenerator_->setSeed(rng->mySeed());
0411 randomGenerator_->setEngine(rng->getEngine(iLumi.index()));
0412 }
0413
0414 void PileUp::CalculatePileup(int MinBunch,
0415 int MaxBunch,
0416 std::vector<int>& PileupSelection,
0417 std::vector<float>& TrueNumInteractions,
0418 StreamID const& streamID) {
0419
0420
0421
0422 int nzero_crossing = -1;
0423 double Fnzero_crossing = -1;
0424
0425 if (provider_) {
0426 setRandomEngine(streamID);
0427 }
0428
0429 if (manage_OOT_) {
0430 if (none_) {
0431 nzero_crossing = 0;
0432 } else if (poisson_) {
0433 nzero_crossing = poissonDistribution(streamID)->fire();
0434 } else if (fixed_) {
0435 nzero_crossing = intAverage_;
0436 } else if (histoDistribution_ || probFunctionDistribution_) {
0437
0438
0439
0440
0441
0442
0443
0444 double d = GetRandom(histo_.get(), randomEngine(streamID));
0445
0446 Fnzero_crossing = d;
0447 nzero_crossing = int(d);
0448 }
0449 }
0450
0451 for (int bx = MinBunch; bx < MaxBunch + 1; ++bx) {
0452 if (manage_OOT_) {
0453 if (bx == 0 && !poisson_OOT_) {
0454 PileupSelection.push_back(nzero_crossing);
0455 TrueNumInteractions.push_back(nzero_crossing);
0456 } else {
0457 if (poisson_OOT_) {
0458 if (PU_Study_ && (Study_type_ == "Fixed_ITPU_Vary_OOTPU") && bx == 0) {
0459 PileupSelection.push_back(intFixed_ITPU_);
0460 } else {
0461 PileupSelection.push_back(poissonDistr_OOT(streamID)->fire(Fnzero_crossing));
0462 }
0463 TrueNumInteractions.push_back(Fnzero_crossing);
0464 } else {
0465 PileupSelection.push_back(intFixed_OOT_);
0466 TrueNumInteractions.push_back(intFixed_OOT_);
0467 }
0468 }
0469 } else {
0470 if (none_) {
0471 PileupSelection.push_back(0);
0472 TrueNumInteractions.push_back(0.);
0473 } else if (poisson_) {
0474 PileupSelection.push_back(poissonDistribution(streamID)->fire());
0475 TrueNumInteractions.push_back(averageNumber_);
0476 } else if (fixed_) {
0477 PileupSelection.push_back(intAverage_);
0478 TrueNumInteractions.push_back(intAverage_);
0479 } else if (histoDistribution_ || probFunctionDistribution_) {
0480
0481
0482
0483
0484
0485
0486
0487 double d = GetRandom(histo_.get(), randomEngine(streamID));
0488 PileupSelection.push_back(int(d));
0489 TrueNumInteractions.push_back(d);
0490 }
0491 }
0492 }
0493 }
0494
0495 }