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