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/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 /// return a random number distributed according the histogram bin contents.
0039 ///
0040 /// This routine is derived from root/hist/hist/src/TH1.cxx
0041 /*************************************************************************
0042  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
0043  * All rights reserved.                                                  *
0044  *                                                                       *
0045  * For the licensing terms see $ROOTSYS/LICENSE.                         *
0046  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
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     // Use the empty parameter set for the parameter set ID of our "@MIXING" process.
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     // A modified HistoryAppender must be used for unscheduled processing.
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       // Check for string describing special processing.  Add these here for individual cases
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_) {  // figure out what the parameters are
0155
0156         //      if (playback_) throw cms::Exception("Illegal parameter clash","PileUp::PileUp(ParameterSet const& pset)")
0157         // << " manage_OOT option not allowed with playback ";
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           // read back the fixed number requested out-of-time
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") {  // allow for some extra flexibility for mixing
0182       minBunch_cosmics_ = pset.getUntrackedParameter<int>("minBunch_cosmics", -1000);
0183       maxBunch_cosmics_ = pset.getUntrackedParameter<int>("maxBunch_cosmics", 1000);
0184     }
0185   }  // end of constructor
0186
0187   void PileUp::beginStream(edm::StreamID) {
0188     auto iID = eventPrincipal_->streamID();  // each producer has its own workermanager, so use default streamid
0189     streamContext_.reset(new StreamContext(iID, processContext_.get()));
0190     input_->doBeginJob();
0191     if (provider_.get() != nullptr) {
0192       //TODO for now, we do not support consumes from EventSetup
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       // note:  run and lumi numbers must be modified to match lumiPrincipal_
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     //get the required parameters from DB.
0244     const MixingInputConfig& config = setup.getData(configToken_).config(inputType_);
0245
0246     //get the type
0247     type_ = config.type();
0248     //set booleans
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       //need to reload the histogram from DB
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       // Complete the vector containing the probability  function data
0280       // with the values "0"
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       // Create an histogram with the data from the probability function provided by the user
0293       int xmin = (int)dataProbFunctionVar[0];
0294       int xmax = (int)dataProbFunctionVar[varSize - 1] + 1;  // need upper edge to be one beyond last value
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]);  // assuming integer values for the bins, fill bin centers, not edges
0308       }
0309
0310       // Check if the histogram is normalized
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     // if we are managing the distribution of out-of-time pileup separately, select the distribution for bunch
0362     // crossing zero first, save it for later.
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         // RANDOM_NUMBER_ERROR
0376         // Random number should be generated by the engines from the
0377         // RandomNumberGeneratorService. This appears to use the global
0378         // engine in ROOT. This is not thread safe unless the module using
0379         // it is a one module and declares a shared resource and all
0380         // other modules using it also declare the same shared resource.
0381         // This also breaks replay.
0382         double d = GetRandom(histo_.get(), randomEngine(streamID));
0383         //n = (int) floor(d + 0.5);  // incorrect for bins with integer edges
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           // RANDOM_NUMBER_ERROR
0419           // Random number should be generated by the engines from the
0420           // RandomNumberGeneratorService. This appears to use the global
0421           // engine in ROOT. This is not thread safe unless the module using
0422           // it is a one module and declares a shared resource and all
0423           // other modules using it also declare the same shared resource.
0424           // This also breaks replay.
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 }  //namespace edm