CMS 3D CMS Logo

PileUp.cc
Go to the documentation of this file.
20 
23 
29 
30 #include "CLHEP/Random/RandPoissonQ.h"
31 #include "CLHEP/Random/RandPoisson.h"
32 
33 #include <algorithm>
34 #include <memory>
35 #include "TMath.h"
36 
41 /*************************************************************************
42  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
43  * All rights reserved. *
44  * *
45  * For the licensing terms see $ROOTSYS/LICENSE. *
46  * For the list of contributors see $ROOTSYS/README/CREDITS. *
47  *************************************************************************/
48 
49 static Double_t GetRandom(TH1* th1, CLHEP::HepRandomEngine* rng) {
50  Int_t nbinsx = th1->GetNbinsX();
51  Double_t* fIntegral = th1->GetIntegral();
52  Double_t integral = fIntegral[nbinsx];
53 
54  if (integral == 0)
55  return 0;
56 
57  Double_t r1 = rng->flat();
58  Int_t ibin = TMath::BinarySearch(nbinsx, fIntegral, r1);
59  Double_t x = th1->GetBinLowEdge(ibin + 1);
60  if (r1 > fIntegral[ibin])
61  x += th1->GetBinWidth(ibin + 1) * (r1 - fIntegral[ibin]) / (fIntegral[ibin + 1] - fIntegral[ibin]);
62  return x;
63 }
65 
66 namespace edm {
68  const std::shared_ptr<PileUpConfig>& config,
70  const bool mixingConfigFromDB)
71  : type_(pset.getParameter<std::string>("type")),
72  Source_type_(config->sourcename_),
73  averageNumber_(config->averageNumber_),
74  intAverage_(static_cast<int>(averageNumber_)),
75  histo_(std::make_shared<TH1F>(*config->histo_)),
76  histoDistribution_(type_ == "histo"),
77  probFunctionDistribution_(type_ == "probFunction"),
78  poisson_(type_ == "poisson"),
79  fixed_(type_ == "fixed"),
80  none_(type_ == "none"),
81  fileNameHash_(0U),
82  productRegistry_(new SignallingProductRegistry),
84  ->makeVectorInputSource(
86  .release()),
87  processConfiguration_(new ProcessConfiguration(std::string("@MIXING"), getReleaseVersion(), getPassID())),
88  processContext_(new ProcessContext()),
89  eventPrincipal_(),
90  lumiPrincipal_(),
91  runPrincipal_(),
92  provider_(),
93  PoissonDistribution_(),
94  PoissonDistr_OOT_(),
95  randomEngine_(),
96  playback_(config->playback_),
97  sequential_(pset.getUntrackedParameter<bool>("sequential", false)) {
98  if (mixingConfigFromDB) {
100  }
101 
102  // Use the empty parameter set for the parameter set ID of our "@MIXING" process.
104  processContext_->setProcessConfiguration(processConfiguration_.get());
105 
106  if (pset.existsAs<std::vector<ParameterSet> >("producers", true)) {
107  std::vector<ParameterSet> producers = pset.getParameter<std::vector<ParameterSet> >("producers");
108  provider_ = std::make_unique<SecondaryEventProvider>(producers, *productRegistry_, processConfiguration_);
109  }
110 
111  productRegistry_->setFrozen();
112 
113  // A modified HistoryAppender must be used for unscheduled processing.
114  eventPrincipal_ = std::make_unique<EventPrincipal>(input_->productRegistry(),
115  std::make_shared<BranchIDListHelper>(),
116  std::make_shared<ThinnedAssociationsHelper>(),
118  nullptr);
119 
120  bool DB = type_ == "readDB";
121 
122  if (pset.exists("nbPileupEvents")) {
123  if (0 != pset.getParameter<edm::ParameterSet>("nbPileupEvents").getUntrackedParameter<int>("seed", 0)) {
124  edm::LogWarning("MixingModule") << "Parameter nbPileupEvents.seed is not supported";
125  }
126  }
127 
129  if (!rng.isAvailable()) {
130  throw cms::Exception("Configuration")
131  << "PileUp requires the RandomNumberGeneratorService\n"
132  "which is not present in the configuration file. You must add the service\n"
133  "in the configuration file or remove the modules that require it.";
134  }
135 
137  throw cms::Exception("Illegal parameter value", "PileUp::PileUp(ParameterSet const& pset)")
138  << "'type' parameter (a string) has a value of '" << type_ << "'.\n"
139  << "Legal values are 'poisson', 'fixed', or 'none'\n";
140  }
141 
142  if (!DB) {
143  manage_OOT_ = pset.getUntrackedParameter<bool>("manage_OOT", false);
144 
145  // Check for string describing special processing. Add these here for individual cases
146  PU_Study_ = false;
147  Study_type_ = pset.getUntrackedParameter<std::string>("Special_Pileup_Studies", "");
148 
149  if (Study_type_ == "Fixed_ITPU_Vary_OOTPU") {
150  PU_Study_ = true;
151  intFixed_ITPU_ = pset.getUntrackedParameter<int>("intFixed_ITPU", 0);
152  }
153 
154  if (manage_OOT_) { // figure out what the parameters are
155 
156  // if (playback_) throw cms::Exception("Illegal parameter clash","PileUp::PileUp(ParameterSet const& pset)")
157  // << " manage_OOT option not allowed with playback ";
158 
159  std::string OOT_type = pset.getUntrackedParameter<std::string>("OOT_type");
160 
161  if (OOT_type == "Poisson" || OOT_type == "poisson") {
162  poisson_OOT_ = true;
163  } else if (OOT_type == "Fixed" || OOT_type == "fixed") {
164  fixed_OOT_ = true;
165  // read back the fixed number requested out-of-time
166  intFixed_OOT_ = pset.getUntrackedParameter<int>("intFixed_OOT", -1);
167  if (intFixed_OOT_ < 0) {
168  throw cms::Exception("Illegal parameter value", "PileUp::PileUp(ParameterSet const& pset)")
169  << " Fixed out-of-time pileup requested, but no fixed value given ";
170  }
171  } else {
172  throw cms::Exception("Illegal parameter value", "PileUp::PileUp(ParameterSet const& pset)")
173  << "'OOT_type' parameter (a string) has a value of '" << OOT_type << "'.\n"
174  << "Legal values are 'poisson' or 'fixed'\n";
175  }
176  edm::LogInfo("MixingModule") << " Out-of-time pileup will be generated with a " << OOT_type
177  << " distribution. ";
178  }
179  }
180 
181  if (Source_type_ == "cosmics") { // allow for some extra flexibility for mixing
182  minBunch_cosmics_ = pset.getUntrackedParameter<int>("minBunch_cosmics", -1000);
183  maxBunch_cosmics_ = pset.getUntrackedParameter<int>("maxBunch_cosmics", 1000);
184  }
185  } // end of constructor
186 
188  auto iID = eventPrincipal_->streamID(); // each producer has its own workermanager, so use default streamid
189  streamContext_.reset(new StreamContext(iID, processContext_.get()));
190  input_->doBeginJob();
191  if (provider_.get() != nullptr) {
192  //TODO for now, we do not support consumes from EventSetup
194  provider_->beginStream(iID, *streamContext_);
195  }
196  }
197 
199  if (provider_.get() != nullptr) {
200  provider_->endStream(streamContext_->streamID(), *streamContext_);
201  provider_->endJob();
202  }
203  input_->doEndJob();
204  }
205 
207  if (provider_.get() != nullptr) {
209  runPrincipal_->setAux(run.runAuxiliary());
210  provider_->beginRun(*runPrincipal_, setup.impl(), run.moduleCallingContext(), *streamContext_);
211  }
212  }
214  if (provider_.get() != nullptr) {
216  lumiPrincipal_->setAux(lumi.luminosityBlockAuxiliary());
217  lumiPrincipal_->setRunPrincipal(runPrincipal_);
218  provider_->beginLuminosityBlock(*lumiPrincipal_, setup.impl(), lumi.moduleCallingContext(), *streamContext_);
219  }
220  }
221 
223  if (provider_.get() != nullptr) {
224  provider_->endRun(*runPrincipal_, setup.impl(), run.moduleCallingContext(), *streamContext_);
225  }
226  }
228  if (provider_.get() != nullptr) {
229  provider_->endLuminosityBlock(*lumiPrincipal_, setup.impl(), lumi.moduleCallingContext(), *streamContext_);
230  }
231  }
232 
234  if (provider_.get() != nullptr) {
235  // note: run and lumi numbers must be modified to match lumiPrincipal_
236  eventPrincipal_->setLuminosityBlockPrincipal(lumiPrincipal_.get());
237  eventPrincipal_->setRunAndLumiNumber(lumiPrincipal_->run(), lumiPrincipal_->luminosityBlock());
238  provider_->setupPileUpEvent(*eventPrincipal_, setup.impl(), *streamContext_);
239  }
240  }
241 
243  //get the required parameters from DB.
244  const MixingInputConfig& config = setup.getData(configToken_).config(inputType_);
245 
246  //get the type
247  type_ = config.type();
248  //set booleans
249  histoDistribution_ = type_ == "histo";
250  probFunctionDistribution_ = type_ == "probFunction";
251  poisson_ = type_ == "poisson";
252  fixed_ = type_ == "fixed";
253  none_ = type_ == "none";
254 
255  if (histoDistribution_)
256  edm::LogError("MisConfiguration") << "type histo cannot be reloaded from DB, yet";
257 
258  if (fixed_) {
260  } else if (poisson_) {
261  averageNumber_ = config.averageNumber();
262  if (PoissonDistribution_) {
263  PoissonDistribution_ = std::make_unique<CLHEP::RandPoissonQ>(PoissonDistribution_->engine(), averageNumber_);
264  }
265  } else if (probFunctionDistribution_) {
266  //need to reload the histogram from DB
267  const std::vector<int>& dataProbFunctionVar = config.probFunctionVariable();
268  std::vector<double> dataProb = config.probValue();
269 
270  int varSize = (int)dataProbFunctionVar.size();
271  int probSize = (int)dataProb.size();
272 
273  if ((dataProbFunctionVar[0] != 0) || (dataProbFunctionVar[varSize - 1] != (varSize - 1)))
274  throw cms::Exception("BadProbFunction")
275  << "Please, check the variables of the probability function! The first variable should be 0 and the "
276  "difference between two variables should be 1."
277  << std::endl;
278 
279  // Complete the vector containing the probability function data
280  // with the values "0"
281  if (probSize < varSize) {
282  edm::LogWarning("MixingModule") << " The probability function data will be completed with "
283  << (varSize - probSize) << " values 0.";
284 
285  for (int i = 0; i < (varSize - probSize); i++)
286  dataProb.push_back(0);
287 
288  probSize = dataProb.size();
289  edm::LogInfo("MixingModule") << " The number of the P(x) data set after adding the values 0 is " << probSize;
290  }
291 
292  // Create an histogram with the data from the probability function provided by the user
293  int xmin = (int)dataProbFunctionVar[0];
294  int xmax = (int)dataProbFunctionVar[varSize - 1] + 1; // need upper edge to be one beyond last value
295  int numBins = varSize;
296 
297  edm::LogInfo("MixingModule") << "An histogram will be created with " << numBins << " bins in the range (" << xmin
298  << "," << xmax << ")." << std::endl;
299 
300  histo_.reset(new TH1F("h", "Histo from the user's probability function", numBins, xmin, xmax));
301 
302  LogDebug("MixingModule") << "Filling histogram with the following data:" << std::endl;
303 
304  for (int j = 0; j < numBins; j++) {
305  LogDebug("MixingModule") << " x = " << dataProbFunctionVar[j] << " P(x) = " << dataProb[j];
306  histo_->Fill(dataProbFunctionVar[j] + 0.5,
307  dataProb[j]); // assuming integer values for the bins, fill bin centers, not edges
308  }
309 
310  // Check if the histogram is normalized
311  if (std::abs(histo_->Integral() - 1) > 1.0e-02) {
312  throw cms::Exception("BadProbFunction") << "The probability function should be normalized!!! " << std::endl;
313  }
314  averageNumber_ = histo_->GetMean();
315  }
316 
317  int oot = config.outOfTime();
318  manage_OOT_ = false;
319  if (oot == 1) {
320  manage_OOT_ = true;
321  poisson_OOT_ = false;
322  fixed_OOT_ = true;
323  intFixed_OOT_ = config.fixedOutOfTime();
324  } else if (oot == 2) {
325  manage_OOT_ = true;
326  poisson_OOT_ = true;
327  fixed_OOT_ = false;
328  }
329  }
331 
332  std::unique_ptr<CLHEP::RandPoissonQ> const& PileUp::poissonDistribution(StreamID const& streamID) {
333  if (!PoissonDistribution_) {
334  CLHEP::HepRandomEngine& engine = *randomEngine(streamID);
335  PoissonDistribution_ = std::make_unique<CLHEP::RandPoissonQ>(engine, averageNumber_);
336  }
337  return PoissonDistribution_;
338  }
339 
340  std::unique_ptr<CLHEP::RandPoisson> const& PileUp::poissonDistr_OOT(StreamID const& streamID) {
341  if (!PoissonDistr_OOT_) {
342  CLHEP::HepRandomEngine& engine = *randomEngine(streamID);
343  PoissonDistr_OOT_ = std::make_unique<CLHEP::RandPoisson>(engine);
344  }
345  return PoissonDistr_OOT_;
346  }
347 
348  CLHEP::HepRandomEngine* PileUp::randomEngine(StreamID const& streamID) {
349  if (!randomEngine_) {
351  randomEngine_ = &rng->getEngine(streamID);
352  }
353  return randomEngine_;
354  }
355 
356  void PileUp::CalculatePileup(int MinBunch,
357  int MaxBunch,
358  std::vector<int>& PileupSelection,
359  std::vector<float>& TrueNumInteractions,
360  StreamID const& streamID) {
361  // if we are managing the distribution of out-of-time pileup separately, select the distribution for bunch
362  // crossing zero first, save it for later.
363 
364  int nzero_crossing = -1;
365  double Fnzero_crossing = -1;
366 
367  if (manage_OOT_) {
368  if (none_) {
369  nzero_crossing = 0;
370  } else if (poisson_) {
371  nzero_crossing = poissonDistribution(streamID)->fire();
372  } else if (fixed_) {
373  nzero_crossing = intAverage_;
375  // RANDOM_NUMBER_ERROR
376  // Random number should be generated by the engines from the
377  // RandomNumberGeneratorService. This appears to use the global
378  // engine in ROOT. This is not thread safe unless the module using
379  // it is a one module and declares a shared resource and all
380  // other modules using it also declare the same shared resource.
381  // This also breaks replay.
382  double d = GetRandom(histo_.get(), randomEngine(streamID));
383  //n = (int) floor(d + 0.5); // incorrect for bins with integer edges
384  Fnzero_crossing = d;
385  nzero_crossing = int(d);
386  }
387  }
388 
389  for (int bx = MinBunch; bx < MaxBunch + 1; ++bx) {
390  if (manage_OOT_) {
391  if (bx == 0 && !poisson_OOT_) {
392  PileupSelection.push_back(nzero_crossing);
393  TrueNumInteractions.push_back(nzero_crossing);
394  } else {
395  if (poisson_OOT_) {
396  if (PU_Study_ && (Study_type_ == "Fixed_ITPU_Vary_OOTPU") && bx == 0) {
397  PileupSelection.push_back(intFixed_ITPU_);
398  } else {
399  PileupSelection.push_back(poissonDistr_OOT(streamID)->fire(Fnzero_crossing));
400  }
401  TrueNumInteractions.push_back(Fnzero_crossing);
402  } else {
403  PileupSelection.push_back(intFixed_OOT_);
404  TrueNumInteractions.push_back(intFixed_OOT_);
405  }
406  }
407  } else {
408  if (none_) {
409  PileupSelection.push_back(0);
410  TrueNumInteractions.push_back(0.);
411  } else if (poisson_) {
412  PileupSelection.push_back(poissonDistribution(streamID)->fire());
413  TrueNumInteractions.push_back(averageNumber_);
414  } else if (fixed_) {
415  PileupSelection.push_back(intAverage_);
416  TrueNumInteractions.push_back(intAverage_);
418  // RANDOM_NUMBER_ERROR
419  // Random number should be generated by the engines from the
420  // RandomNumberGeneratorService. This appears to use the global
421  // engine in ROOT. This is not thread safe unless the module using
422  // it is a one module and declares a shared resource and all
423  // other modules using it also declare the same shared resource.
424  // This also breaks replay.
425  double d = GetRandom(histo_.get(), randomEngine(streamID));
426  PileupSelection.push_back(int(d));
427  TrueNumInteractions.push_back(d);
428  }
429  }
430  }
431  }
432 
433 } //namespace edm
bool manage_OOT_
Definition: PileUp.h:126
std::string getPassID()
Definition: GetPassID.h:7
std::unique_ptr< SecondaryEventProvider > provider_
Definition: PileUp.h:149
OOT_type
manage out-of-time pileup setting this to True means that the out-of-time pileup will have a differen...
std::shared_ptr< ProcessConfiguration > processConfiguration_
Definition: PileUp.h:143
double averageNumber() const
Definition: PileUp.h:72
unsigned int inputType_
Definition: PileUp.h:115
static Double_t GetRandom(TH1 *th1, CLHEP::HepRandomEngine *rng)
Definition: PileUp.cc:49
void beginStream(edm::StreamID)
Definition: PileUp.cc:187
bool fixed_
Definition: PileUp.h:124
Definition: config.py:1
Log< level::Error, false > LogError
PileUp(ParameterSet const &pset, const std::shared_ptr< PileUpConfig > &config, edm::ConsumesCollector iC, const bool mixingConfigFromDB)
Definition: PileUp.cc:67
std::unique_ptr< CLHEP::RandPoissonQ > const & poissonDistribution(StreamID const &streamID)
Definition: PileUp.cc:332
int const intAverage_
Definition: PileUp.h:119
std::shared_ptr< TH1F > histo_
Definition: PileUp.h:120
bool probFunctionDistribution_
Definition: PileUp.h:122
int maxBunch_cosmics_
Definition: PileUp.h:137
std::string Study_type_
Definition: PileUp.h:131
T getUntrackedParameter(std::string const &, T const &) const
bool histoDistribution_
Definition: PileUp.h:121
bool poisson_
Definition: PileUp.h:123
void endLuminosityBlock(const edm::LuminosityBlock &lumi, const edm::EventSetup &setup)
Definition: PileUp.cc:227
std::unique_ptr< VectorInputSource > const input_
Definition: PileUp.h:142
bool poisson_OOT_
Definition: PileUp.h:127
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int intFixed_ITPU_
Definition: PileUp.h:134
void beginRun(const edm::Run &run, const edm::EventSetup &setup)
Definition: PileUp.cc:206
void CalculatePileup(int MinBunch, int MaxBunch, std::vector< int > &PileupSelection, std::vector< float > &TrueNumInteractions, StreamID const &)
Definition: PileUp.cc:356
bool PU_Study_
Definition: PileUp.h:130
void beginLuminosityBlock(const edm::LuminosityBlock &lumi, const edm::EventSetup &setup)
Definition: PileUp.cc:213
bool fixed_OOT_
Definition: PileUp.h:128
d
Definition: ztail.py:151
int minBunch_cosmics_
Definition: PileUp.h:136
Integral< F, X >::type integral(const F &f)
Definition: Integral.h:70
std::unique_ptr< CLHEP::RandPoisson > const & poissonDistr_OOT(StreamID const &streamID)
Definition: PileUp.cc:340
std::unique_ptr< CLHEP::RandPoisson > PoissonDistr_OOT_
Definition: PileUp.h:151
Log< level::Info, false > LogInfo
std::string getReleaseVersion()
CLHEP::HepRandomEngine * randomEngine_
Definition: PileUp.h:152
void setupPileUpEvent(const edm::EventSetup &setup)
Definition: PileUp.cc:233
std::shared_ptr< LuminosityBlockPrincipal > lumiPrincipal_
Definition: PileUp.h:147
std::string type_
Definition: PileUp.h:116
void endRun(const edm::Run &run, const edm::EventSetup &setup)
Definition: PileUp.cc:222
std::shared_ptr< ProductRegistry > productRegistry_
Definition: PileUp.h:141
std::string Source_type_
Definition: PileUp.h:117
void reload(const edm::EventSetup &setup)
Definition: PileUp.cc:242
HLT enums.
bool none_
Definition: PileUp.h:125
T const & get(Event const &event, InputTag const &tag) noexcept(false)
Definition: Event.h:680
CLHEP::HepRandomEngine * randomEngine(StreamID const &streamID)
Definition: PileUp.cc:348
float x
std::unique_ptr< CLHEP::RandPoissonQ > PoissonDistribution_
Definition: PileUp.h:150
bool isAvailable() const
Definition: Service.h:40
void endStream()
Definition: PileUp.cc:198
int intFixed_OOT_
Definition: PileUp.h:133
std::shared_ptr< ProcessContext > processContext_
Definition: PileUp.h:144
Log< level::Warning, false > LogWarning
std::unique_ptr< EventPrincipal > eventPrincipal_
Definition: PileUp.h:146
std::shared_ptr< RunPrincipal > runPrincipal_
Definition: PileUp.h:148
static ParameterSetID emptyParameterSetID()
Definition: ParameterSet.cc:94
std::shared_ptr< StreamContext > streamContext_
Definition: PileUp.h:145
double averageNumber_
Definition: PileUp.h:118
Definition: Run.h:45
edm::ESGetToken< MixingModuleConfig, MixingRcd > configToken_
Definition: PileUp.h:139
#define LogDebug(id)