CMS 3D CMS Logo
List of all members | Public Member Functions | Protected Attributes
edm::Lumi3DReWeighting Class Reference

#include <Lumi3DReWeighting.h>

Public Member Functions

 Lumi3DReWeighting ()
 
 Lumi3DReWeighting (const std::vector< float > &MC_distr, const std::vector< float > &Lumi_distr, std::string WeightOutputFile)
 
 Lumi3DReWeighting (std::string generatedFile, std::string dataFile, std::string GenHistName, std::string DataHistName, std::string WeightOutputFile)
 
double weight3D (const edm::EventBase &e)
 
double weight3D (int, int, int)
 
void weight3D_init (float Scale)
 
void weight3D_init (std::string MCFileName, std::string DataFileName)
 
void weight3D_init (std::string WeightFileName)
 
void weight3D_set (std::string generatedFile, std::string dataFile, std::string GenHistName, std::string DataHistName, std::string WeightOutputFile)
 

Protected Attributes

std::shared_ptr< TH1 > Data_distr_
 
std::shared_ptr< TFile > dataFile_
 
std::string dataFileName_
 
std::string DataHistName_
 
std::shared_ptr< TFile > generatedFile_
 
std::string generatedFileName_
 
std::string GenHistName_
 
std::shared_ptr< TH1 > MC_distr_
 
double Weight3D_ [50][50][50]
 
std::string weightFileName_
 
std::shared_ptr< TH1 > weights_
 

Detailed Description

Definition at line 26 of file Lumi3DReWeighting.h.

Constructor & Destructor Documentation

◆ Lumi3DReWeighting() [1/3]

Lumi3DReWeighting::Lumi3DReWeighting ( std::string  generatedFile,
std::string  dataFile,
std::string  GenHistName = "pileup",
std::string  DataHistName = "pileup",
std::string  WeightOutputFile = "" 
)

Definition at line 40 of file Lumi3DReWeighting.cc.

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 }

References Data_distr_, dataFile_, dataFileName_, DataHistName_, generatedFile_, generatedFileName_, GenHistName_, and MC_distr_.

◆ Lumi3DReWeighting() [2/3]

Lumi3DReWeighting::Lumi3DReWeighting ( const std::vector< float > &  MC_distr,
const std::vector< float > &  Lumi_distr,
std::string  WeightOutputFile = "" 
)

Definition at line 70 of file Lumi3DReWeighting.cc.

72  {
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 }

References Data_distr_, MC_distr_, and weightFileName_.

◆ Lumi3DReWeighting() [3/3]

edm::Lumi3DReWeighting::Lumi3DReWeighting ( )
inline

Definition at line 38 of file Lumi3DReWeighting.h.

38 {};

Member Function Documentation

◆ weight3D() [1/2]

double Lumi3DReWeighting::weight3D ( const edm::EventBase &  e)

Definition at line 119 of file Lumi3DReWeighting.cc.

119  {
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 }

References L1TStage2uGTEmulatorClient_cff::BX, MillePedeFileConverter_cfg::e, min(), edm::min(), and Weight3D_.

◆ weight3D() [2/2]

double Lumi3DReWeighting::weight3D ( int  pv1,
int  pv2,
int  pv3 
)

Definition at line 107 of file Lumi3DReWeighting.cc.

107  {
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 }

References min(), edm::min(), and Weight3D_.

◆ weight3D_init() [1/3]

void Lumi3DReWeighting::weight3D_init ( float  Scale)

Definition at line 190 of file Lumi3DReWeighting.cc.

190  {
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 }

References newFWLiteAna::base, gather_cfg::cout, Data_distr_, Exception, JetChargeProducer_cfi::exp, factorial(), dqmMemoryStats::float, mps_fire::i, createfilelist::int, dqmiolumiharvest::j, dqmdumpme::k, MC_distr_, SiStripPI::mean, min(), timingPdfMaker::outfile, Weight3D_, weightFileName_, and protons_cff::xi.

◆ weight3D_init() [2/3]

void Lumi3DReWeighting::weight3D_init ( std::string  MCFileName,
std::string  DataFileName 
)

Definition at line 387 of file Lumi3DReWeighting.cc.

387  {
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 }

References gather_cfg::cout, Exception, mps_fire::i, dqmiolumiharvest::j, dqmdumpme::k, and Weight3D_.

◆ weight3D_init() [3/3]

void Lumi3DReWeighting::weight3D_init ( std::string  WeightFileName)

Definition at line 361 of file Lumi3DReWeighting.cc.

361  {
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 }

References gather_cfg::cout, Exception, mps_fire::i, timingPdfMaker::infile, dqmiolumiharvest::j, dqmdumpme::k, and Weight3D_.

◆ weight3D_set()

void Lumi3DReWeighting::weight3D_set ( std::string  generatedFile,
std::string  dataFile,
std::string  GenHistName,
std::string  DataHistName,
std::string  WeightOutputFile 
)

Definition at line 156 of file Lumi3DReWeighting.cc.

160  {
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 }

References gather_cfg::cout, Data_distr_, dataFile_, dataFileName_, DataHistName_, generatedFile_, generatedFileName_, GenHistName_, MC_distr_, and weightFileName_.

Member Data Documentation

◆ Data_distr_

std::shared_ptr<TH1> edm::Lumi3DReWeighting::Data_distr_
protected

Definition at line 69 of file Lumi3DReWeighting.h.

Referenced by Lumi3DReWeighting(), weight3D_init(), and weight3D_set().

◆ dataFile_

std::shared_ptr<TFile> edm::Lumi3DReWeighting::dataFile_
protected

Definition at line 63 of file Lumi3DReWeighting.h.

Referenced by Lumi3DReWeighting(), and weight3D_set().

◆ dataFileName_

std::string edm::Lumi3DReWeighting::dataFileName_
protected

Definition at line 58 of file Lumi3DReWeighting.h.

Referenced by Lumi3DReWeighting(), and weight3D_set().

◆ DataHistName_

std::string edm::Lumi3DReWeighting::DataHistName_
protected

Definition at line 60 of file Lumi3DReWeighting.h.

Referenced by Lumi3DReWeighting(), and weight3D_set().

◆ generatedFile_

std::shared_ptr<TFile> edm::Lumi3DReWeighting::generatedFile_
protected

Definition at line 62 of file Lumi3DReWeighting.h.

Referenced by Lumi3DReWeighting(), and weight3D_set().

◆ generatedFileName_

std::string edm::Lumi3DReWeighting::generatedFileName_
protected

Definition at line 57 of file Lumi3DReWeighting.h.

Referenced by Lumi3DReWeighting(), and weight3D_set().

◆ GenHistName_

std::string edm::Lumi3DReWeighting::GenHistName_
protected

Definition at line 59 of file Lumi3DReWeighting.h.

Referenced by Lumi3DReWeighting(), and weight3D_set().

◆ MC_distr_

std::shared_ptr<TH1> edm::Lumi3DReWeighting::MC_distr_
protected

Definition at line 68 of file Lumi3DReWeighting.h.

Referenced by Lumi3DReWeighting(), weight3D_init(), and weight3D_set().

◆ Weight3D_

double edm::Lumi3DReWeighting::Weight3D_[50][50][50]
protected

Definition at line 71 of file Lumi3DReWeighting.h.

Referenced by weight3D(), and weight3D_init().

◆ weightFileName_

std::string edm::Lumi3DReWeighting::weightFileName_
protected

Definition at line 61 of file Lumi3DReWeighting.h.

Referenced by Lumi3DReWeighting(), weight3D_init(), and weight3D_set().

◆ weights_

std::shared_ptr<TH1> edm::Lumi3DReWeighting::weights_
protected

Definition at line 64 of file Lumi3DReWeighting.h.

edm::Lumi3DReWeighting::GenHistName_
std::string GenHistName_
Definition: Lumi3DReWeighting.h:59
edm::Lumi3DReWeighting::dataFile_
std::shared_ptr< TFile > dataFile_
Definition: Lumi3DReWeighting.h:63
mps_fire.i
i
Definition: mps_fire.py:428
dqmMemoryStats.float
float
Definition: dqmMemoryStats.py:127
SiStripPI::mean
Definition: SiStripPayloadInspectorHelper.h:169
min
T min(T a, T b)
Definition: MathUtil.h:58
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_
double Weight3D_[50][50][50]
Definition: Lumi3DReWeighting.h:71
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
edm::Lumi3DReWeighting::generatedFile_
std::shared_ptr< TFile > generatedFile_
Definition: Lumi3DReWeighting.h:62
dqmdumpme.k
k
Definition: dqmdumpme.py:60
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
Exception
Definition: hltDiff.cc:245
genVertex_cff.x
x
Definition: genVertex_cff.py:12
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
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
Generated for CMSSW Reference Manual by doxygen 1.8.16