CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
BMixingModule.cc
Go to the documentation of this file.
1 // File: BMixingModule.cc
2 // Description: see BMixingModule.h
3 // Author: Ursula Berthon, LLR Palaiseau, Bill Tanenbaum
4 //
5 //--------------------------------------------
6 
16 
17 #include "TFile.h"
18 #include "TH1F.h"
19 #include <iostream>
20 const unsigned int edm::BMixingModule::maxNbSources_ =4;
21 
22 namespace
23 {
24  std::shared_ptr<edm::PileUp>
25  maybeMakePileUp(edm::ParameterSet const& ps,std::string sourceName, const int minb, const int maxb, const bool playback)
26  {
27  std::shared_ptr<edm::PileUp> pileup; // value to be returned
28  // Make sure we have a parameter named 'sourceName'
29  if (ps.exists(sourceName))
30  {
31  // We have the parameter
32  // and if we have either averageNumber or cfg by luminosity... make the PileUp
33  double averageNumber;
35  std::string histoName =" ";
36  TH1F * h = new TH1F("h","h",10,0,10);
37  std::vector<int> dataProbFunctionVar;
38  std::vector<double> dataProb;
39 
40 
41  const edm::ParameterSet & psin=ps.getParameter<edm::ParameterSet>(sourceName);
42  std::string type_=psin.getParameter<std::string>("type");
43  if (ps.exists("readDB") && ps.getParameter<bool>("readDB")){
44  //in case of DB access, do not try to load anything from the PSet, but wait for beginRun.
45  edm::LogError("BMixingModule")<<"Will read from DB: reset to a dummy PileUp object.";
46  pileup.reset(new edm::PileUp(psin,sourceName,0,0,playback));
47  return pileup;
48  }
49  if (type_!="none"){
50  if (psin.exists("nbPileupEvents")) {
51  edm::ParameterSet psin_average=psin.getParameter<edm::ParameterSet>("nbPileupEvents");
52  if (psin_average.exists("averageNumber"))
53  {
54  averageNumber=psin_average.getParameter<double>("averageNumber");
55  pileup.reset(new edm::PileUp(psin,sourceName,averageNumber,h,playback));
56  edm::LogInfo("MixingModule") <<" Created source "<<sourceName<<" with averageNumber "<<averageNumber;
57  }
58  else if (psin_average.exists("fileName") && psin_average.exists("histoName")) {
59  std::string histoFileName = psin_average.getUntrackedParameter<std::string>("fileName");
60  std::string histoName = psin_average.getUntrackedParameter<std::string>("histoName");
61 
62  TFile *infile = new TFile(histoFileName.c_str());
63  TH1F *h = (TH1F*)infile->Get(histoName.c_str());
64 
65  // Check if the histogram exists
66  if (!h) {
67  throw cms::Exception("HistogramNotFound") << " Could not find the histogram " << histoName
68  << "in the file " << histoFileName << "." << std::endl;
69  }else{
70  edm::LogInfo("MixingModule") << "Open a root file " << histoFileName << " containing the probability distribution histogram " << histoName << std::endl;
71  edm::LogInfo("MixingModule") << "The PileUp number to be added will be chosen randomly from this histogram" << std::endl;
72  }
73 
74  // Check if the histogram is normalized
75  if (std::abs(h->Integral() - 1) > 1.0e-02) throw cms::Exception("BadHistoDistribution") << "The histogram should be normalized!" << std::endl;
76 
77  // Get the averageNumber from the histo
78  averageNumber = h->GetMean();
79 
80  pileup.reset(new edm::PileUp(psin,sourceName,averageNumber,h,playback));
81  edm::LogInfo("MixingModule") <<" Created source "<<sourceName<<" with averageNumber "<<averageNumber;
82 
83  }
84  else if (psin_average.exists("probFunctionVariable") && psin_average.exists("probValue") && psin_average.exists("histoFileName"))
85  {
86  if (type_!="probFunction") {
87  edm::LogError("MisConfiguration")<<"type is set to: "<<type_<<" while parameters implies probFunction; changing.";
88  type_="probFunction";
89  }
90 
91  dataProbFunctionVar = psin_average.getParameter<std::vector<int> >("probFunctionVariable");
92  dataProb = psin_average.getParameter<std::vector<double> >("probValue");
93  histoFileName = psin_average.getUntrackedParameter<std::string>("histoFileName");
94 
95  int varSize = (int) dataProbFunctionVar.size();
96  int probSize = (int) dataProb.size();
97 
98  if ((dataProbFunctionVar[0] != 0) || (dataProbFunctionVar[varSize - 1] != (varSize - 1)))
99  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;
100 
101  // Complete the vector containing the probability function data
102  // with the values "0"
103  if (probSize < varSize){
104  edm::LogInfo("MixingModule") << " The probability function data will be completed with " <<(varSize - probSize) <<" values 0.";
105 
106  for (int i=0; i<(varSize - probSize); i++) dataProb.push_back(0);
107 
108  probSize = dataProb.size();
109  edm::LogInfo("MixingModule") << " The number of the P(x) data set after adding the values 0 is " << probSize;
110  }
111 
112  // Create an histogram with the data from the probability function provided by the user
113  int xmin = (int) dataProbFunctionVar[0];
114  int xmax = (int) dataProbFunctionVar[varSize-1]+1; // need upper edge to be one beyond last value
115  int numBins = varSize;
116 
117  edm::LogInfo("MixingModule") << "An histogram will be created with " << numBins << " bins in the range ("<< xmin << "," << xmax << ")." << std::endl;
118 
119  TH1F *hprob = new TH1F("h","Histo from the user's probability function",numBins,xmin,xmax);
120 
121  LogDebug("MixingModule") << "Filling histogram with the following data:" << std::endl;
122 
123  for (int j=0; j < numBins ; j++){
124  LogDebug("MixingModule") << " x = " << dataProbFunctionVar[j ]<< " P(x) = " << dataProb[j];
125  hprob->Fill(dataProbFunctionVar[j]+0.5,dataProb[j]); // assuming integer values for the bins, fill bin centers, not edges
126  }
127 
128  // Check if the histogram is normalized
129  if (std::abs(hprob->Integral() - 1) > 1.0e-02){
130  throw cms::Exception("BadProbFunction") << "The probability function should be normalized!!! " << std::endl;
131  }
132 
133  averageNumber = hprob->GetMean();
134 
135  // Write the created histogram into a root file
136  edm::LogInfo("MixingModule") << " The histogram created from the x, P(x) values will be written into the root file " << histoFileName;
137 
138  TFile * outfile = new TFile(histoFileName.c_str(),"RECREATE");
139  hprob->Write();
140  outfile->Write();
141  outfile->Close();
142  outfile->Delete();
143 
144  pileup.reset(new edm::PileUp(psin,sourceName,averageNumber,hprob,playback));
145  edm::LogInfo("MixingModule") <<" Created source "<<sourceName<<" with averageNumber "<<averageNumber;
146 
147  }
148  //special for pileup input
149  else if (sourceName=="input" && psin_average.exists("Lumi") && psin_average.exists("sigmaInel")) {
150  averageNumber=psin_average.getParameter<double>("Lumi")*psin_average.getParameter<double>("sigmaInel")*ps.getParameter<int>("bunchspace")/1000*3564./2808.; //FIXME
151  pileup.reset(new
152  edm::PileUp(psin,sourceName,averageNumber,h,playback));
153  edm::LogInfo("MixingModule") <<" Created source "<<sourceName<<" with minBunch,maxBunch "<<minb<<" "<<maxb;
154  edm::LogInfo("MixingModule")<<" Luminosity configuration, average number used is "<<averageNumber;
155  }
156  }
157  }
158  }
159  return pileup;
160  }
161 }
162 
163 namespace edm {
164 
165  // Constructor
167  bunchSpace_(pset.getParameter<int>("bunchspace")),
168  vertexOffset_(0),
169  minBunch_((pset.getParameter<int>("minBunch")*25)/pset.getParameter<int>("bunchspace")),
170  maxBunch_((pset.getParameter<int>("maxBunch")*25)/pset.getParameter<int>("bunchspace")),
171  mixProdStep1_(pset.getParameter<bool>("mixProdStep1")),
172  mixProdStep2_(pset.getParameter<bool>("mixProdStep2")),
173  readDB_(false)
174  {
175  if (pset.exists("readDB")) readDB_=pset.getParameter<bool>("readDB");
176 
177  playback_=pset.getUntrackedParameter<bool>("playback",false);
178 
179  if (playback_) {
180  //this could be explicitely checked
181  LogInfo("MixingModule") <<" ATTENTION:Mixing will be done in playback mode! \n"
182  <<" ATTENTION:Mixing Configuration must be the same as for the original mixing!";
183  }
184 
185  // Just for debugging print out.
186  sourceNames_.push_back("input");
187  sourceNames_.push_back("cosmics");
188  sourceNames_.push_back("beamhalo_plus");
189  sourceNames_.push_back("beamhalo_minus");
190 
191  for (size_t makeIdx = 0; makeIdx < maxNbSources_; makeIdx++ ) {
192  inputSources_.push_back(maybeMakePileUp(pset,sourceNames_[makeIdx],
194  if (inputSources_.back()) inputSources_.back()->input(makeIdx);
195  }
196  }
197 
198  // Virtual destructor needed.
200 
201  // update method call at begin run/lumi to reload the mixing configuration
203  update(setup);
204  for (size_t endIdx=0; endIdx<maxNbSources_; ++endIdx) {
205  if(inputSources_[endIdx]) inputSources_[endIdx]->beginLuminosityBlock(lumi, setup);
206  }
207  }
208 
210  update(setup);
211  for (size_t endIdx=0; endIdx<maxNbSources_; ++endIdx) {
212  if(inputSources_[endIdx]) inputSources_[endIdx]->beginRun(run, setup);
213  }
214  }
215 
217  for (size_t endIdx=0; endIdx<maxNbSources_; ++endIdx) {
218  if(inputSources_[endIdx]) inputSources_[endIdx]->endLuminosityBlock(lumi, setup);
219  }
220  }
221 
223  for (size_t endIdx=0; endIdx<maxNbSources_; ++endIdx) {
224  if(inputSources_[endIdx]) inputSources_[endIdx]->endRun(run, setup);
225  }
226  }
227 
229  if (readDB_ && parameterWatcher_.check(setup)){
230  for (size_t makeIdx = 0; makeIdx < maxNbSources_; makeIdx++ ) {
231  if (inputSources_[makeIdx]) inputSources_[makeIdx]->reload(setup);
232  }
233  reload(setup);
234  }
235  }
236 
237  // Functions that get called by framework every event
239  // Check if the signal is present in the root file
240  // for all the objects we want to mix
241  checkSignal(e);
242 
243  // Create EDProduct
245 
246  initializeEvent(e, setup);
247 
248  // Add signals
249  if (!mixProdStep1_){
250  addSignals(e,setup);
251  }
252 
253  doPileUp(e, setup);
254 
255  // Includes putting digi products into the edm::Event.
256  finalizeEvent(e, setup);
257 
258  // Put output into event (here only playback info)
259  put(e,setup);
260  }
261 
263  for (size_t dropIdx=0; dropIdx<maxNbSources_; ++dropIdx) {
264  if(inputSources_[dropIdx]) inputSources_[dropIdx]->setupPileUpEvent(setup);
265  }
266  }
267 
268  void BMixingModule::dropUnwantedBranches(std::vector<std::string> const& wantedBranches) {
269  for (size_t dropIdx=0; dropIdx<maxNbSources_; ++dropIdx) {
270  if(inputSources_[dropIdx]) inputSources_[dropIdx]->dropUnwantedBranches(wantedBranches);
271  }
272  }
273 
275  for (size_t endIdx=0; endIdx<maxNbSources_; ++endIdx) {
276  if(inputSources_[endIdx]) inputSources_[endIdx]->beginJob();
277  }
278  }
279 
281  for (size_t endIdx=0; endIdx<maxNbSources_; ++endIdx) {
282  if(inputSources_[endIdx]) inputSources_[endIdx]->endJob();
283  }
284  }
285 
286  void BMixingModule::createnewEDProduct() {std::cout << "BMixingModule::createnewEDProduct must be overwritten!" << std::endl;}
287 
288  void BMixingModule::checkSignal(const edm::Event &e) {std::cout << "BMixingModule::checkSignal must be overwritten!" << std::endl;}
289 
290  void BMixingModule::setBcrOffset () {std::cout << "BMixingModule::setBcrOffset must be overwritten!" << std::endl;} //FIXME: LogWarning
291 
292  void BMixingModule::setSourceOffset (const unsigned int s) {std::cout << "BMixingModule::setSourceOffset must be overwritten!" << std::endl;}
293 
294  void BMixingModule::doPileUp(edm::Event &e, const edm::EventSetup& c) {std::cout << "BMixingModule::doPileUp must be overwritten!" << std::endl;}
295 
296 } //edm
#define LogDebug(id)
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
edm::ESWatcher< MixingRcd > parameterWatcher_
Definition: BMixingModule.h:97
virtual void beginRun(const edm::Run &r, const edm::EventSetup &setup) override
virtual void addSignals(const edm::Event &e, const edm::EventSetup &c)
Definition: BMixingModule.h:63
virtual void produce(edm::Event &e1, const edm::EventSetup &c) override
tuple lumi
Definition: fjr2json.py:35
virtual void finalizeEvent(edm::Event &event, const edm::EventSetup &setup)
Definition: BMixingModule.h:43
BMixingModule(const edm::ParameterSet &ps)
tuple histoFileName
Definition: diJetCalib.py:108
virtual void checkSignal(const edm::Event &e)
bool exists(std::string const &parameterName) const
checks if a parameter exists
virtual void put(edm::Event &e, const edm::EventSetup &c)
Definition: BMixingModule.h:67
virtual void endRun(const edm::Run &r, const edm::EventSetup &setup) override
virtual void setBcrOffset()
virtual void createnewEDProduct()
virtual void setSourceOffset(const unsigned int s)
virtual void endJob() override
list outfile
Definition: EdgesToViz.py:91
virtual void doPileUp(edm::Event &e, const edm::EventSetup &c)
virtual void initializeEvent(const edm::Event &event, const edm::EventSetup &setup)
Definition: BMixingModule.h:40
tuple averageNumber
set the number of pileup
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int j
Definition: DBlmapReader.cc:9
static const unsigned int maxNbSources_
Definition: BMixingModule.h:86
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
virtual void beginJob() override
virtual void reload(const edm::EventSetup &setup)
Definition: BMixingModule.h:52
bool check(const edm::EventSetup &iSetup)
Definition: ESWatcher.h:57
bool const mixProdStep1_
Definition: BMixingModule.h:81
virtual ~BMixingModule()
std::vector< std::shared_ptr< PileUp > > inputSources_
Definition: BMixingModule.h:94
list infile
Definition: EdgesToViz.py:90
tuple cout
Definition: gather_cfg.py:121
void dropUnwantedBranches(std::vector< std::string > const &wantedBranches)
volatile std::atomic< bool > shutdown_flag false
std::vector< std::string > sourceNames_
Definition: BMixingModule.h:87
virtual void beginLuminosityBlock(const edm::LuminosityBlock &l, const edm::EventSetup &setup) override
void update(edm::EventSetup const &)
void setupPileUpEvent(const edm::EventSetup &setup)
void setup(std::vector< TH2F > &depth, std::string name, std::string units="")
virtual void endLuminosityBlock(const edm::LuminosityBlock &l, const edm::EventSetup &setup) override
Definition: Run.h:43