CMS 3D CMS Logo
Lumi3DReWeighting.cc
Go to the documentation of this file.
1 #ifndef PhysicsTools_Utilities_interface_Lumi3DReWeighting_cc
2 #define PhysicsTools_Utilities_interface_Lumi3DReWeighting_cc
3 
19 #include "TFile.h"
20 #include "TH1.h"
21 #include "TH3.h"
22 #include "TRandom1.h"
23 #include "TRandom2.h"
24 #include "TRandom3.h"
25 #include "TStopwatch.h"
26 #include <algorithm>
27 #include <iostream>
28 #include <memory>
29 #include <string>
30 
31 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
32 #include "PhysicsTools/Utilities/interface/Lumi3DReWeighting.h"
33 #include "DataFormats/Common/interface/Handle.h"
34 #include "DataFormats/Provenance/interface/ProcessHistory.h"
35 #include "DataFormats/Provenance/interface/Provenance.h"
36 #include "FWCore/Common/interface/EventBase.h"
37 
38 using namespace edm;
39 
40 Lumi3DReWeighting::Lumi3DReWeighting(std::string generatedFile,
41  std::string dataFile,
42  std::string GenHistName = "pileup",
43  std::string DataHistName = "pileup",
44  std::string WeightOutputFile = "")
45  : generatedFileName_(generatedFile),
46  dataFileName_(dataFile),
47  GenHistName_(GenHistName),
48  DataHistName_(DataHistName),
49  weightFileName_(WeightOutputFile) {
50  generatedFile_ = std::make_shared<TFile>(generatedFileName_.c_str()); //MC distribution
51  dataFile_ = std::make_shared<TFile>(dataFileName_.c_str()); //Data distribution
52 
53  std::shared_ptr<TH1> Data_temp =
54  std::shared_ptr<TH1>((static_cast<TH1 *>(dataFile_->Get(DataHistName_.c_str())->Clone())));
55 
56  std::shared_ptr<TH1> MC_temp =
57  std::shared_ptr<TH1>((static_cast<TH1 *>(generatedFile_->Get(GenHistName_.c_str())->Clone())));
58 
59  MC_distr_ = std::shared_ptr<TH1>((static_cast<TH1 *>(generatedFile_->Get(GenHistName_.c_str())->Clone())));
60  Data_distr_ = std::shared_ptr<TH1>((static_cast<TH1 *>(dataFile_->Get(DataHistName_.c_str())->Clone())));
61 
62  // MC * data/MC = data, so the weights are data/MC:
63 
64  // normalize both histograms first
65 
66  Data_distr_->Scale(1.0 / Data_distr_->Integral());
67  MC_distr_->Scale(1.0 / MC_distr_->Integral());
68 }
69 
70 Lumi3DReWeighting::Lumi3DReWeighting(const std::vector<float> &MC_distr,
71  const std::vector<float> &Lumi_distr,
72  std::string WeightOutputFile = "") {
73  weightFileName_ = WeightOutputFile;
74 
75  // no histograms for input: use vectors
76 
77  // now, make histograms out of them:
78 
79  Int_t NMCBins = MC_distr.size();
80 
81  MC_distr_ = std::shared_ptr<TH1>(new TH1F("MC_distr", "MC dist", NMCBins, 0., float(NMCBins)));
82 
83  Int_t NDBins = Lumi_distr.size();
84 
85  Data_distr_ = std::shared_ptr<TH1>(new TH1F("Data_distr", "Data dist", NDBins, 0., float(NDBins)));
86 
87  for (int ibin = 1; ibin < NMCBins + 1; ++ibin) {
88  MC_distr_->SetBinContent(ibin, MC_distr[ibin - 1]);
89  }
90 
91  for (int ibin = 1; ibin < NDBins + 1; ++ibin) {
92  Data_distr_->SetBinContent(ibin, Lumi_distr[ibin - 1]);
93  }
94 
95  // check integrals, make sure things are normalized
96 
97  float deltaH = Data_distr_->Integral();
98  if (fabs(1.0 - deltaH) > 0.001) { //*OOPS*...
99  Data_distr_->Scale(1.0 / Data_distr_->Integral());
100  }
101  float deltaMC = MC_distr_->Integral();
102  if (fabs(1.0 - deltaMC) > 0.001) {
103  MC_distr_->Scale(1.0 / MC_distr_->Integral());
104  }
105 }
106 
107 double Lumi3DReWeighting::weight3D(int pv1, int pv2, int pv3) {
108  using std::min;
109 
110  int npm1 = min(pv1, 49);
111  int np0 = min(pv2, 49);
112  int npp1 = min(pv3, 49);
113 
114  return Weight3D_[npm1][np0][npp1];
115 }
116 
117 // This version does all of the work for you, just taking an event pointer as input
118 
119 double Lumi3DReWeighting::weight3D(const edm::EventBase &e) {
120  using std::min;
121 
122  // get pileup summary information
123 
124  Handle<std::vector<PileupSummaryInfo> > PupInfo;
125  e.getByLabel(edm::InputTag("addPileupInfo"), PupInfo);
126 
127  std::vector<PileupSummaryInfo>::const_iterator PVI;
128 
129  int npm1 = -1;
130  int np0 = -1;
131  int npp1 = -1;
132 
133  for (PVI = PupInfo->begin(); PVI != PupInfo->end(); ++PVI) {
134  int BX = PVI->getBunchCrossing();
135 
136  if (BX == -1) {
137  npm1 = PVI->getPU_NumInteractions();
138  }
139  if (BX == 0) {
140  np0 = PVI->getPU_NumInteractions();
141  }
142  if (BX == 1) {
143  npp1 = PVI->getPU_NumInteractions();
144  }
145  }
146 
147  npm1 = min(npm1, 49);
148  np0 = min(np0, 49);
149  npp1 = min(npp1, 49);
150 
151  // std::cout << " vertices " << npm1 << " " << np0 << " " << npp1 << std::endl;
152 
153  return Weight3D_[npm1][np0][npp1];
154 }
155 
156 void Lumi3DReWeighting::weight3D_set(std::string generatedFile,
157  std::string dataFile,
158  std::string GenHistName,
159  std::string DataHistName,
160  std::string WeightOutputFile) {
161  generatedFileName_ = generatedFile;
162  dataFileName_ = dataFile;
163  GenHistName_ = GenHistName;
164  DataHistName_ = DataHistName;
165  weightFileName_ = WeightOutputFile;
166 
167  std::cout << " seting values: " << generatedFileName_ << " " << dataFileName_ << " " << GenHistName_ << " "
168  << DataHistName_ << std::endl;
169 
170  generatedFile_ = std::make_shared<TFile>(generatedFileName_.c_str()); //MC distribution
171  dataFile_ = std::make_shared<TFile>(dataFileName_.c_str()); //Data distribution
172 
173  std::shared_ptr<TH1> Data_temp =
174  std::shared_ptr<TH1>((static_cast<TH1 *>(dataFile_->Get(DataHistName_.c_str())->Clone())));
175 
176  std::shared_ptr<TH1> MC_temp =
177  std::shared_ptr<TH1>((static_cast<TH1 *>(generatedFile_->Get(GenHistName_.c_str())->Clone())));
178 
179  MC_distr_ = std::shared_ptr<TH1>((static_cast<TH1 *>(generatedFile_->Get(GenHistName_.c_str())->Clone())));
180  Data_distr_ = std::shared_ptr<TH1>((static_cast<TH1 *>(dataFile_->Get(DataHistName_.c_str())->Clone())));
181 
182  // MC * data/MC = data, so the weights are data/MC:
183 
184  // normalize both histograms first
185 
186  Data_distr_->Scale(1.0 / Data_distr_->Integral());
187  MC_distr_->Scale(1.0 / MC_distr_->Integral());
188 }
189 
190 void Lumi3DReWeighting::weight3D_init(float ScaleFactor) {
191  // Scale factor is used to shift target distribution (i.e. luminosity scale) 1. = no shift
192 
193  //create histogram to write output weights, save pain of generating them again...
194 
195  TH3D *WHist = new TH3D("WHist", "3D weights", 50, 0., 50., 50, 0., 50., 50, 0., 50.);
196  TH3D *DHist = new TH3D("DHist", "3D weights", 50, 0., 50., 50, 0., 50., 50, 0., 50.);
197  TH3D *MHist = new TH3D("MHist", "3D weights", 50, 0., 50., 50, 0., 50., 50, 0., 50.);
198 
199  using std::min;
200 
201  if (MC_distr_->GetEntries() == 0) {
202  std::cout << " MC and Data distributions are not initialized! You must call the Lumi3DReWeighting constructor. "
203  << std::endl;
204  }
205 
206  // arrays for storing number of interactions
207 
208  double MC_ints[50][50][50];
209  double Data_ints[50][50][50];
210 
211  for (int i = 0; i < 50; i++) {
212  for (int j = 0; j < 50; j++) {
213  for (int k = 0; k < 50; k++) {
214  MC_ints[i][j][k] = 0.;
215  Data_ints[i][j][k] = 0.;
216  }
217  }
218  }
219 
220  double factorial[50];
221  double PowerSer[50];
222  double base = 1.;
223 
224  factorial[0] = 1.;
225  PowerSer[0] = 1.;
226 
227  for (int i = 1; i < 50; ++i) {
228  base = base * float(i);
229  factorial[i] = base;
230  }
231 
232  double x;
233  double xweight;
234  double probi, probj, probk;
235  double Expval, mean;
236  int xi;
237 
238  // Get entries for Data, MC, fill arrays:
239 
240  int NMCbin = MC_distr_->GetNbinsX();
241 
242  for (int jbin = 1; jbin < NMCbin + 1; jbin++) {
243  x = MC_distr_->GetBinCenter(jbin);
244  xweight = MC_distr_->GetBinContent(jbin); //use as weight for matrix
245 
246  //for Summer 11, we have this int feature:
247  xi = int(x);
248 
249  // Generate Poisson distribution for each value of the mean
250 
251  mean = double(xi);
252 
253  if (mean < 0.) {
254  throw cms::Exception("BadInputValue") << " Your histogram generates MC luminosity values less than zero!"
255  << " Please Check. Terminating." << std::endl;
256  }
257 
258  if (mean == 0.) {
259  Expval = 1.;
260  } else {
261  Expval = exp(-1. * mean);
262  }
263 
264  base = 1.;
265 
266  for (int i = 1; i < 50; ++i) {
267  base = base * mean;
268  PowerSer[i] = base; // PowerSer is mean^i
269  }
270 
271  // compute poisson probability for each Nvtx in weight matrix
272 
273  for (int i = 0; i < 50; i++) {
274  probi = PowerSer[i] / factorial[i] * Expval;
275  for (int j = 0; j < 50; j++) {
276  probj = PowerSer[j] / factorial[j] * Expval;
277  for (int k = 0; k < 50; k++) {
278  probk = PowerSer[k] / factorial[k] * Expval;
279  // joint probability is product of event weights multiplied by weight of input distribution bin
280  MC_ints[i][j][k] = MC_ints[i][j][k] + probi * probj * probk * xweight;
281  }
282  }
283  }
284  }
285 
286  int NDatabin = Data_distr_->GetNbinsX();
287 
288  for (int jbin = 1; jbin < NDatabin + 1; jbin++) {
289  mean = (Data_distr_->GetBinCenter(jbin)) * ScaleFactor;
290  xweight = Data_distr_->GetBinContent(jbin);
291 
292  // Generate poisson distribution for each value of the mean
293 
294  if (mean < 0.) {
295  throw cms::Exception("BadInputValue") << " Your histogram generates Data luminosity values less than zero!"
296  << " Please Check. Terminating." << std::endl;
297  }
298 
299  if (mean == 0.) {
300  Expval = 1.;
301  } else {
302  Expval = exp(-1. * mean);
303  }
304 
305  base = 1.;
306 
307  for (int i = 1; i < 50; ++i) {
308  base = base * mean;
309  PowerSer[i] = base;
310  }
311 
312  // compute poisson probability for each Nvtx in weight matrix
313 
314  for (int i = 0; i < 50; i++) {
315  probi = PowerSer[i] / factorial[i] * Expval;
316  for (int j = 0; j < 50; j++) {
317  probj = PowerSer[j] / factorial[j] * Expval;
318  for (int k = 0; k < 50; k++) {
319  probk = PowerSer[k] / factorial[k] * Expval;
320  // joint probability is product of event weights multiplied by weight of input distribution bin
321  Data_ints[i][j][k] = Data_ints[i][j][k] + probi * probj * probk * xweight;
322  }
323  }
324  }
325  }
326 
327  for (int i = 0; i < 50; i++) {
328  //if(i<5) std::cout << "i = " << i << std::endl;
329  for (int j = 0; j < 50; j++) {
330  for (int k = 0; k < 50; k++) {
331  if ((MC_ints[i][j][k]) > 0.) {
332  Weight3D_[i][j][k] = Data_ints[i][j][k] / MC_ints[i][j][k];
333  } else {
334  Weight3D_[i][j][k] = 0.;
335  }
336  WHist->SetBinContent(i + 1, j + 1, k + 1, Weight3D_[i][j][k]);
337  DHist->SetBinContent(i + 1, j + 1, k + 1, Data_ints[i][j][k]);
338  MHist->SetBinContent(i + 1, j + 1, k + 1, MC_ints[i][j][k]);
339  // if(i<5 && j<5 && k<5) std::cout << Weight3D_[i][j][k] << " " ;
340  }
341  // if(i<5 && j<5) std::cout << std::endl;
342  }
343  }
344 
345  if (!weightFileName_.empty()) {
346  std::cout << " 3D Weight Matrix initialized! " << std::endl;
347  std::cout << " Writing weights to file " << weightFileName_ << " for re-use... " << std::endl;
348 
349  TFile *outfile = new TFile(weightFileName_.c_str(), "RECREATE");
350  WHist->Write();
351  MHist->Write();
352  DHist->Write();
353  outfile->Write();
354  outfile->Close();
355  outfile->Delete();
356  }
357 
358  return;
359 }
360 
361 void Lumi3DReWeighting::weight3D_init(std::string WeightFileName) {
362  TFile *infile = new TFile(WeightFileName.c_str());
363  TH1F *WHist = (TH1F *)infile->Get("WHist");
364 
365  // Check if the histogram exists
366  if (!WHist) {
367  throw cms::Exception("HistogramNotFound") << " Could not find the histogram WHist in the file "
368  << "in the file " << WeightFileName << "." << std::endl;
369  }
370 
371  for (int i = 0; i < 50; i++) {
372  // if(i<5) std::cout << "i = " << i << std::endl;
373  for (int j = 0; j < 50; j++) {
374  for (int k = 0; k < 50; k++) {
375  Weight3D_[i][j][k] = WHist->GetBinContent(i + 1, j + 1, k + 1);
376  // if(i<5 && j<5 && k<5) std::cout << Weight3D_[i][j][k] << " ";
377  }
378  // if(i<5 && j<5) std::cout << std::endl;
379  }
380  }
381 
382  std::cout << " 3D Weight Matrix initialized! " << std::endl;
383 
384  return;
385 }
386 
387 void Lumi3DReWeighting::weight3D_init(std::string MCWeightFileName, std::string DataWeightFileName) {
388  TFile *infileMC = new TFile(MCWeightFileName.c_str());
389  TH3D *MHist = (TH3D *)infileMC->Get("MHist");
390 
391  // Check if the histogram exists
392  if (!MHist) {
393  throw cms::Exception("HistogramNotFound") << " Could not find the histogram MHist in the file "
394  << "in the file " << MCWeightFileName << "." << std::endl;
395  }
396 
397  TFile *infileD = new TFile(DataWeightFileName.c_str());
398  TH3D *DHist = (TH3D *)infileD->Get("DHist");
399 
400  // Check if the histogram exists
401  if (!DHist) {
402  throw cms::Exception("HistogramNotFound") << " Could not find the histogram DHist in the file "
403  << "in the file " << DataWeightFileName << "." << std::endl;
404  }
405 
406  for (int i = 0; i < 50; i++) {
407  for (int j = 0; j < 50; j++) {
408  for (int k = 0; k < 50; k++) {
409  Weight3D_[i][j][k] = DHist->GetBinContent(i + 1, j + 1, k + 1) / MHist->GetBinContent(i + 1, j + 1, k + 1);
410  }
411  }
412  }
413 
414  std::cout << " 3D Weight Matrix initialized! " << std::endl;
415 
416  return;
417 }
418 
419 #endif
edm::Lumi3DReWeighting::GenHistName_
std::string GenHistName_
Definition: Lumi3DReWeighting.h:59
edm::Lumi3DReWeighting::dataFile_
std::shared_ptr< TFile > dataFile_
Definition: Lumi3DReWeighting.h:63
Handle.h
PileupSummaryInfo.h
mps_fire.i
i
Definition: mps_fire.py:428
Lumi3DReWeighting.h
dqmMemoryStats.float
float
Definition: dqmMemoryStats.py:127
SiStripPI::mean
Definition: SiStripPayloadInspectorHelper.h:169
edm::Lumi3DReWeighting::Lumi3DReWeighting
Lumi3DReWeighting()
Definition: Lumi3DReWeighting.h:38
min
T min(T a, T b)
Definition: MathUtil.h:58
edm
HLT enums.
Definition: AlignableModifier.h:19
edm::min
EventID const & min(EventID const &lh, EventID const &rh)
Definition: EventID.h:116
gather_cfg.cout
cout
Definition: gather_cfg.py:144
factorial
int factorial(int n)
factorial function
Definition: L1GlobalTriggerFunctions.cc:13
edm::Lumi3DReWeighting::weight3D_set
void weight3D_set(std::string generatedFile, std::string dataFile, std::string GenHistName, std::string DataHistName, std::string WeightOutputFile)
Definition: Lumi3DReWeighting.cc:156
edm::Lumi3DReWeighting::Weight3D_
double Weight3D_[50][50][50]
Definition: Lumi3DReWeighting.h:71
EventBase.h
edm::Handle
Definition: AssociativeIterator.h:50
edm::Lumi3DReWeighting::DataHistName_
std::string DataHistName_
Definition: Lumi3DReWeighting.h:60
edm::Lumi3DReWeighting::generatedFileName_
std::string generatedFileName_
Definition: Lumi3DReWeighting.h:57
edm::Lumi3DReWeighting::Data_distr_
std::shared_ptr< TH1 > Data_distr_
Definition: Lumi3DReWeighting.h:69
ProcessHistory.h
Provenance.h
edm::Lumi3DReWeighting::generatedFile_
std::shared_ptr< TFile > generatedFile_
Definition: Lumi3DReWeighting.h:62
dqmdumpme.k
k
Definition: dqmdumpme.py:60
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
edm::Lumi3DReWeighting::weight3D_init
void weight3D_init(float Scale)
Definition: Lumi3DReWeighting.cc:190
createfilelist.int
int
Definition: createfilelist.py:10
edm::Lumi3DReWeighting::dataFileName_
std::string dataFileName_
Definition: Lumi3DReWeighting.h:58
edm::Lumi3DReWeighting::MC_distr_
std::shared_ptr< TH1 > MC_distr_
Definition: Lumi3DReWeighting.h:68
protons_cff.xi
xi
Definition: protons_cff.py:35
edm::EventBase
Definition: EventBase.h:46
Exception
Definition: hltDiff.cc:245
L1TStage2uGTEmulatorClient_cff.BX
BX
Definition: L1TStage2uGTEmulatorClient_cff.py:9
timingPdfMaker.infile
infile
Definition: timingPdfMaker.py:350
timingPdfMaker.outfile
outfile
Definition: timingPdfMaker.py:351
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
newFWLiteAna.base
base
Definition: newFWLiteAna.py:92
edm::Lumi3DReWeighting::weightFileName_
std::string weightFileName_
Definition: Lumi3DReWeighting.h:61
edm::InputTag
Definition: InputTag.h:15
edm::Lumi3DReWeighting::weight3D
double weight3D(const edm::EventBase &e)
Definition: Lumi3DReWeighting.cc:119
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
Generated for CMSSW Reference Manual by doxygen 1.8.16