CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PileUp.cc
Go to the documentation of this file.
14 
17 
23 
24 #include "boost/shared_ptr.hpp"
25 
26 #include <algorithm>
27 
28 namespace edm {
29  PileUp::PileUp(ParameterSet const& pset, double averageNumber, TH1F * const histo, const bool playback) :
30  type_(pset.getParameter<std::string>("type")),
31  averageNumber_(averageNumber),
32  intAverage_(static_cast<int>(averageNumber)),
33  histo_(histo),
34  histoDistribution_(type_ == "histo"),
35  probFunctionDistribution_(type_ == "probFunction"),
36  poisson_(type_ == "poisson"),
37  fixed_(type_ == "fixed"),
38  none_(type_ == "none"),
39  productRegistry_(new SignallingProductRegistry),
40  input_(VectorInputSourceFactory::get()->makeVectorInputSource(pset, InputSourceDescription(
42  *productRegistry_,
43  boost::shared_ptr<BranchIDListHelper>(new BranchIDListHelper),
44  boost::shared_ptr<ActivityRegistry>(new ActivityRegistry),
45  -1, -1,
47  )).release()),
48  processConfiguration_(new ProcessConfiguration(std::string("@MIXING"), getReleaseVersion(), getPassID())),
49  eventPrincipal_(),
50  lumiPrincipal_(),
51  runPrincipal_(),
52  provider_(),
53  poissonDistribution_(),
54  poissonDistr_OOT_(),
55  playback_(playback),
56  sequential_(pset.getUntrackedParameter<bool>("sequential", false)),
57  samelumi_(pset.getUntrackedParameter<bool>("sameLumiBlock", false)),
58  seed_(0) {
59 
60  // Use the empty parameter set for the parameter set ID of our "@MIXING" process.
62 
63  if(pset.existsAs<std::vector<ParameterSet> >("producers", true)) {
64  std::vector<ParameterSet> producers = pset.getParameter<std::vector<ParameterSet> >("producers");
66  }
67 
68  productRegistry_->setFrozen();
69 
70  // A modified HistoryAppender must be used for unscheduled processing.
71  eventPrincipal_.reset(new EventPrincipal(input_->productRegistry(),
72  input_->branchIDListHelper(),
74  nullptr));
75 
76  if (pset.exists("nbPileupEvents")) {
77  seed_=pset.getParameter<edm::ParameterSet>("nbPileupEvents").getUntrackedParameter<int>("seed",0);
78  }
79 
80  bool DB=type_=="readDB";
81 
83  if (!rng.isAvailable()) {
84  throw cms::Exception("Configuration")
85  << "PileUp requires the RandomNumberGeneratorService\n"
86  "which is not present in the configuration file. You must add the service\n"
87  "in the configuration file or remove the modules that require it.";
88  }
89 
90  CLHEP::HepRandomEngine& engine = rng->getEngine();
91  poissonDistribution_.reset(new CLHEP::RandPoissonQ(engine, averageNumber_));
92 
93  // Get seed for the case when using user histogram or probability function
95  if(seed_ !=0) {
96  gRandom->SetSeed(seed_);
97  LogInfo("MixingModule") << " Change seed for " << type_ << " mode. The seed is set to " << seed_;
98  }
99  else {
100  gRandom->SetSeed(engine.getSeed());
101  }
102  }
103 
104 
106  throw cms::Exception("Illegal parameter value","PileUp::PileUp(ParameterSet const& pset)")
107  << "'type' parameter (a string) has a value of '" << type_ << "'.\n"
108  << "Legal values are 'poisson', 'fixed', or 'none'\n";
109  }
110 
111  if (!DB){
112  manage_OOT_ = pset.getUntrackedParameter<bool>("manage_OOT", false);
113 
114  // Check for string describing special processing. Add these here for individual cases
115 
116  PU_Study_ = false;
117  Study_type_ = "";
118 
119  Study_type_ = pset.getUntrackedParameter<std::string>("Special_Pileup_Studies", "");
120 
121  if(Study_type_ == "Fixed_ITPU_Vary_OOTPU") {
122 
123  PU_Study_ = true;
124  intFixed_ITPU_ = pset.getUntrackedParameter<int>("intFixed_ITPU", 0);
125 
126  }
127 
128  if(manage_OOT_) { // figure out what the parameters are
129 
130  // if (playback_) throw cms::Exception("Illegal parameter clash","PileUp::PileUp(ParameterSet const& pset)")
131  // << " manage_OOT option not allowed with playback ";
132 
133  std::string OOT_type = pset.getUntrackedParameter<std::string>("OOT_type");
134 
135  if(OOT_type == "Poisson" || OOT_type == "poisson") {
136  poisson_OOT_ = true;
137  poissonDistr_OOT_.reset(new CLHEP::RandPoisson(engine));
138  }
139  else if(OOT_type == "Fixed" || OOT_type == "fixed") {
140  fixed_OOT_ = true;
141  // read back the fixed number requested out-of-time
142  intFixed_OOT_ = pset.getUntrackedParameter<int>("intFixed_OOT", -1);
143  if(intFixed_OOT_ < 0) {
144  throw cms::Exception("Illegal parameter value","PileUp::PileUp(ParameterSet const& pset)")
145  << " Fixed out-of-time pileup requested, but no fixed value given ";
146  }
147  }
148  else {
149  throw cms::Exception("Illegal parameter value","PileUp::PileUp(ParameterSet const& pset)")
150  << "'OOT_type' parameter (a string) has a value of '" << OOT_type << "'.\n"
151  << "Legal values are 'poisson' or 'fixed'\n";
152  }
153  edm::LogInfo("MixingModule") <<" Out-of-time pileup will be generated with a " << OOT_type << " distribution. " ;
154  }
155  }
156 
157  }
159  input_->doBeginJob();
160  if (provider_.get() != nullptr) {
161  provider_->beginJob(*productRegistry_);
162  }
163  }
164 
165  void PileUp::endJob () {
166  if (provider_.get() != nullptr) {
167  provider_->endJob();
168  }
169  input_->doEndJob();
170  }
171 
173  if (provider_.get() != nullptr) {
174  boost::shared_ptr<RunAuxiliary> aux(new RunAuxiliary(run.runAuxiliary()));
175  runPrincipal_.reset(new RunPrincipal(aux, productRegistry_, *processConfiguration_, nullptr, 0));
176  provider_->beginRun(*runPrincipal_, setup, run.moduleCallingContext());
177  }
178  }
180  if (provider_.get() != nullptr) {
181  boost::shared_ptr<LuminosityBlockAuxiliary> aux(new LuminosityBlockAuxiliary(lumi.luminosityBlockAuxiliary()));
183  lumiPrincipal_->setRunPrincipal(runPrincipal_);
184  provider_->beginLuminosityBlock(*lumiPrincipal_, setup, lumi.moduleCallingContext());
185  }
186  }
187 
189  if (provider_.get() != nullptr) {
190  provider_->endRun(*runPrincipal_, setup, run.moduleCallingContext());
191  }
192  }
194  if (provider_.get() != nullptr) {
195  provider_->endLuminosityBlock(*lumiPrincipal_, setup, lumi.moduleCallingContext());
196  }
197  }
198 
200  if (provider_.get() != nullptr) {
201  // note: run and lumi numbers must be modified to match lumiPrincipal_
202  eventPrincipal_->setLuminosityBlockPrincipal(lumiPrincipal_);
203  eventPrincipal_->setRunAndLumiNumber(lumiPrincipal_->run(), lumiPrincipal_->luminosityBlock());
204  provider_->setupPileUpEvent(*eventPrincipal_, setup);
205  }
206  }
207 
209  //get the required parameters from DB.
211  setup.get<MixingRcd>().get(configM);
212 
213  const MixingInputConfig & config=configM->config(inputType_);
214 
215  //get the type
216  type_=config.type();
217  //set booleans
218  histoDistribution_=type_ == "histo";
219  probFunctionDistribution_=type_ == "probFunction";
220  poisson_=type_ == "poisson";
221  fixed_=type_ == "fixed";
222  none_=type_ == "none";
223 
224  if (histoDistribution_) edm::LogError("MisConfiguration")<<"type histo cannot be reloaded from DB, yet";
225 
226  if (fixed_){
228  }
229  else if (poisson_)
230  {
231  averageNumber_=config.averageNumber();
233  CLHEP::HepRandomEngine& engine = rng->getEngine();
234  poissonDistribution_.reset(new CLHEP::RandPoissonQ(engine, averageNumber_));
235  }
236  else if (probFunctionDistribution_)
237  {
238  //need to reload the histogram from DB
239  const std::vector<int> & dataProbFunctionVar = config.probFunctionVariable();
240  std::vector<double> dataProb = config.probValue();
241 
242  int varSize = (int) dataProbFunctionVar.size();
243  int probSize = (int) dataProb.size();
244 
245  if ((dataProbFunctionVar[0] != 0) || (dataProbFunctionVar[varSize - 1] != (varSize - 1)))
246  throw cms::Exception("BadProbFunction") << "Please, check the variables of the probability function! The first variable should be 0 and the difference between two variables should be 1." << std::endl;
247 
248  // Complete the vector containing the probability function data
249  // with the values "0"
250  if (probSize < varSize){
251  edm::LogWarning("MixingModule") << " The probability function data will be completed with " <<(varSize - probSize) <<" values 0.";
252 
253  for (int i=0; i<(varSize - probSize); i++) dataProb.push_back(0);
254 
255  probSize = dataProb.size();
256  edm::LogInfo("MixingModule") << " The number of the P(x) data set after adding the values 0 is " << probSize;
257  }
258 
259  // Create an histogram with the data from the probability function provided by the user
260  int xmin = (int) dataProbFunctionVar[0];
261  int xmax = (int) dataProbFunctionVar[varSize-1]+1; // need upper edge to be one beyond last value
262  int numBins = varSize;
263 
264  edm::LogInfo("MixingModule") << "An histogram will be created with " << numBins << " bins in the range ("<< xmin << "," << xmax << ")." << std::endl;
265 
266  if (histo_) delete histo_;
267  histo_ = new TH1F("h","Histo from the user's probability function",numBins,xmin,xmax);
268 
269  LogDebug("MixingModule") << "Filling histogram with the following data:" << std::endl;
270 
271  for (int j=0; j < numBins ; j++){
272  LogDebug("MixingModule") << " x = " << dataProbFunctionVar[j ]<< " P(x) = " << dataProb[j];
273  histo_->Fill(dataProbFunctionVar[j]+0.5,dataProb[j]); // assuming integer values for the bins, fill bin centers, not edges
274  }
275 
276  // Check if the histogram is normalized
277  if ( ((histo_->Integral() - 1) > 1.0e-02) && ((histo_->Integral() - 1) < -1.0e-02)){
278  throw cms::Exception("BadProbFunction") << "The probability function should be normalized!!! " << std::endl;
279  }
280  averageNumber_=histo_->GetMean();
281  }
282 
283  int oot=config.outOfTime();
284  manage_OOT_=false;
285  if (oot==1)
286  {
287  manage_OOT_=true;
288  poisson_OOT_ = false;
289  if (poissonDistr_OOT_){poissonDistr_OOT_.reset(); }
290  fixed_OOT_ = true;
291  intFixed_OOT_=config.fixedOutOfTime();
292  }
293  else if (oot==2)
294  {
295  manage_OOT_=true;
296  poisson_OOT_ = true;
297  fixed_OOT_ = false;
298  if (!poissonDistr_OOT_) {
299  //no need to trash the previous one if already there
301  CLHEP::HepRandomEngine& engine = rng->getEngine();
302  poissonDistr_OOT_.reset(new CLHEP::RandPoisson(engine));
303  }
304  }
305 
306 
307  }
309  }
310 
311  void PileUp::CalculatePileup(int MinBunch, int MaxBunch, std::vector<int>& PileupSelection, std::vector<float>& TrueNumInteractions) {
312 
313  // if we are managing the distribution of out-of-time pileup separately, select the distribution for bunch
314  // crossing zero first, save it for later.
315 
316  int nzero_crossing = -1;
317  double Fnzero_crossing = -1;
318 
319  if(manage_OOT_) {
320  if (none_){
321  nzero_crossing = 0;
322  }else if (poisson_){
323  nzero_crossing = poissonDistribution_->fire() ;
324  }else if (fixed_){
325  nzero_crossing = intAverage_ ;
327  double d = histo_->GetRandom();
328  //n = (int) floor(d + 0.5); // incorrect for bins with integer edges
329  Fnzero_crossing = d;
330  nzero_crossing = int(d);
331  }
332 
333  }
334 
335  for(int bx = MinBunch; bx < MaxBunch+1; ++bx) {
336 
337  if(manage_OOT_) {
338  if(bx==0 && !poisson_OOT_) {
339  PileupSelection.push_back(nzero_crossing) ;
340  TrueNumInteractions.push_back( nzero_crossing );
341  }
342  else{
343  if(poisson_OOT_) {
344  if(PU_Study_ && (Study_type_ == "Fixed_ITPU_Vary_OOTPU" ) && bx==0 ) {
345  PileupSelection.push_back(intFixed_ITPU_) ;
346  }
347  else{
348  PileupSelection.push_back(poissonDistr_OOT_->fire(Fnzero_crossing)) ;
349  }
350  TrueNumInteractions.push_back( Fnzero_crossing );
351  }
352  else {
353  PileupSelection.push_back(intFixed_OOT_) ;
354  TrueNumInteractions.push_back( intFixed_OOT_ );
355  }
356  }
357  }
358  else {
359  if (none_){
360  PileupSelection.push_back(0);
361  TrueNumInteractions.push_back( 0. );
362  }else if (poisson_){
363  PileupSelection.push_back(poissonDistribution_->fire());
364  TrueNumInteractions.push_back( averageNumber_ );
365  }else if (fixed_){
366  PileupSelection.push_back(intAverage_);
367  TrueNumInteractions.push_back( intAverage_ );
369  double d = histo_->GetRandom();
370  PileupSelection.push_back(int(d));
371  TrueNumInteractions.push_back( d );
372  }
373  }
374 
375  }
376  }
377 
378 
379 } //namespace edm
#define LogDebug(id)
void CalculatePileup(int MinBunch, int MaxBunch, std::vector< int > &PileupSelection, std::vector< float > &TrueNumInteractions)
Definition: PileUp.cc:311
boost::shared_ptr< RunPrincipal > runPrincipal_
Definition: PileUp.h:101
bool manage_OOT_
Definition: PileUp.h:86
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
std::string getPassID()
Definition: GetPassID.h:8
int i
Definition: DBlmapReader.cc:9
std::unique_ptr< SecondaryEventProvider > provider_
Definition: PileUp.h:102
bool existsAs(std::string const &parameterName, bool trackiness=true) const
checks if a parameter exists as a given type
Definition: ParameterSet.h:184
unsigned int inputType_
Definition: PileUp.h:76
std::unique_ptr< CLHEP::RandPoisson > poissonDistr_OOT_
Definition: PileUp.h:104
const int fixedOutOfTime() const
tuple lumi
Definition: fjr2json.py:35
bool exists(std::string const &parameterName) const
checks if a parameter exists
void endJob()
Definition: PileUp.cc:165
TH1F * histo_
Definition: PileUp.h:80
bool fixed_
Definition: PileUp.h:84
const std::vector< double > & probValue() const
const double averageNumber() const
ModuleCallingContext const * moduleCallingContext() const
PileUp(ParameterSet const &pset, double averageNumber, TH1F *const histo, const bool playback)
Definition: PileUp.cc:29
int const intAverage_
Definition: PileUp.h:79
bool probFunctionDistribution_
Definition: PileUp.h:82
std::string Study_type_
Definition: PileUp.h:91
boost::shared_ptr< LuminosityBlockPrincipal > lumiPrincipal_
Definition: PileUp.h:100
bool histoDistribution_
Definition: PileUp.h:81
bool poisson_
Definition: PileUp.h:83
void endLuminosityBlock(const edm::LuminosityBlock &lumi, const edm::EventSetup &setup)
Definition: PileUp.cc:193
double averageNumber() const
Definition: PileUp.h:47
const int outOfTime() const
std::unique_ptr< VectorInputSource > const input_
Definition: PileUp.h:97
bool poisson_OOT_
Definition: PileUp.h:87
void beginJob()
Definition: PileUp.cc:158
tuple averageNumber
set the number of pileup
bool isAvailable() const
Definition: Service.h:46
int j
Definition: DBlmapReader.cc:9
int intFixed_ITPU_
Definition: PileUp.h:94
virtual CLHEP::HepRandomEngine & getEngine() const =0
Use this to get the random number engine, this is the only function most users should call...
void beginRun(const edm::Run &run, const edm::EventSetup &setup)
Definition: PileUp.cc:172
bool PU_Study_
Definition: PileUp.h:90
void beginLuminosityBlock(const edm::LuminosityBlock &lumi, const edm::EventSetup &setup)
Definition: PileUp.cc:179
bool fixed_OOT_
Definition: PileUp.h:88
std::string getReleaseVersion()
void setupPileUpEvent(const edm::EventSetup &setup)
Definition: PileUp.cc:199
std::string type_
Definition: PileUp.h:77
LuminosityBlockAuxiliary const & luminosityBlockAuxiliary() const
boost::shared_ptr< ProductRegistry > productRegistry_
Definition: PileUp.h:96
const T & get() const
Definition: EventSetup.h:55
void endRun(const edm::Run &run, const edm::EventSetup &setup)
Definition: PileUp.cc:188
const std::vector< int > & probFunctionVariable() const
RunAuxiliary const & runAuxiliary() const
Definition: Run.h:55
boost::shared_ptr< ProcessConfiguration > processConfiguration_
Definition: PileUp.h:98
void reload(const edm::EventSetup &setup)
Definition: PileUp.cc:208
std::unique_ptr< CLHEP::RandPoissonQ > poissonDistribution_
Definition: PileUp.h:103
int seed_
Definition: PileUp.h:121
ModuleCallingContext const * moduleCallingContext() const
Definition: Run.h:135
bool none_
Definition: PileUp.h:85
tuple playback
Definition: Playback_cff.py:20
int intFixed_OOT_
Definition: PileUp.h:93
volatile std::atomic< bool > shutdown_flag false
std::unique_ptr< EventPrincipal > eventPrincipal_
Definition: PileUp.h:99
std::string type() const
static ParameterSetID emptyParameterSetID()
void setup(std::vector< TH2F > &depth, std::string name, std::string units="")
T get(const Candidate &c)
Definition: component.h:55
double averageNumber_
Definition: PileUp.h:78
Definition: Run.h:41