File indexing completed on 2024-04-06 12:24:22
0001 #ifndef PhysicsTools_Utilities_interface_Lumi3DReWeighting_cc
0002 #define PhysicsTools_Utilities_interface_Lumi3DReWeighting_cc
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019 #include "TFile.h"
0020 #include "TH1.h"
0021 #include "TH3.h"
0022 #include "TRandom1.h"
0023 #include "TRandom2.h"
0024 #include "TRandom3.h"
0025 #include "TStopwatch.h"
0026 #include <algorithm>
0027 #include <iostream>
0028 #include <memory>
0029 #include <string>
0030
0031 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
0032 #include "PhysicsTools/Utilities/interface/Lumi3DReWeighting.h"
0033 #include "DataFormats/Common/interface/Handle.h"
0034 #include "DataFormats/Provenance/interface/ProcessHistory.h"
0035 #include "DataFormats/Provenance/interface/Provenance.h"
0036 #include "FWCore/Common/interface/EventBase.h"
0037
0038 using namespace edm;
0039
0040 Lumi3DReWeighting::Lumi3DReWeighting(std::string generatedFile,
0041 std::string dataFile,
0042 std::string GenHistName = "pileup",
0043 std::string DataHistName = "pileup",
0044 std::string WeightOutputFile = "")
0045 : generatedFileName_(generatedFile),
0046 dataFileName_(dataFile),
0047 GenHistName_(GenHistName),
0048 DataHistName_(DataHistName),
0049 weightFileName_(WeightOutputFile) {
0050 generatedFile_ = std::make_shared<TFile>(generatedFileName_.c_str());
0051 dataFile_ = std::make_shared<TFile>(dataFileName_.c_str());
0052
0053 std::shared_ptr<TH1> Data_temp =
0054 std::shared_ptr<TH1>((static_cast<TH1 *>(dataFile_->Get(DataHistName_.c_str())->Clone())));
0055
0056 std::shared_ptr<TH1> MC_temp =
0057 std::shared_ptr<TH1>((static_cast<TH1 *>(generatedFile_->Get(GenHistName_.c_str())->Clone())));
0058
0059 MC_distr_ = std::shared_ptr<TH1>((static_cast<TH1 *>(generatedFile_->Get(GenHistName_.c_str())->Clone())));
0060 Data_distr_ = std::shared_ptr<TH1>((static_cast<TH1 *>(dataFile_->Get(DataHistName_.c_str())->Clone())));
0061
0062
0063
0064
0065
0066 Data_distr_->Scale(1.0 / Data_distr_->Integral());
0067 MC_distr_->Scale(1.0 / MC_distr_->Integral());
0068 }
0069
0070 Lumi3DReWeighting::Lumi3DReWeighting(const std::vector<float> &MC_distr,
0071 const std::vector<float> &Lumi_distr,
0072 std::string WeightOutputFile = "") {
0073 weightFileName_ = WeightOutputFile;
0074
0075
0076
0077
0078
0079 Int_t NMCBins = MC_distr.size();
0080
0081 MC_distr_ = std::shared_ptr<TH1>(new TH1F("MC_distr", "MC dist", NMCBins, 0., float(NMCBins)));
0082
0083 Int_t NDBins = Lumi_distr.size();
0084
0085 Data_distr_ = std::shared_ptr<TH1>(new TH1F("Data_distr", "Data dist", NDBins, 0., float(NDBins)));
0086
0087 for (int ibin = 1; ibin < NMCBins + 1; ++ibin) {
0088 MC_distr_->SetBinContent(ibin, MC_distr[ibin - 1]);
0089 }
0090
0091 for (int ibin = 1; ibin < NDBins + 1; ++ibin) {
0092 Data_distr_->SetBinContent(ibin, Lumi_distr[ibin - 1]);
0093 }
0094
0095
0096
0097 float deltaH = Data_distr_->Integral();
0098 if (fabs(1.0 - deltaH) > 0.001) {
0099 Data_distr_->Scale(1.0 / Data_distr_->Integral());
0100 }
0101 float deltaMC = MC_distr_->Integral();
0102 if (fabs(1.0 - deltaMC) > 0.001) {
0103 MC_distr_->Scale(1.0 / MC_distr_->Integral());
0104 }
0105 }
0106
0107 double Lumi3DReWeighting::weight3D(int pv1, int pv2, int pv3) {
0108 using std::min;
0109
0110 int npm1 = min(pv1, 49);
0111 int np0 = min(pv2, 49);
0112 int npp1 = min(pv3, 49);
0113
0114 return Weight3D_[npm1][np0][npp1];
0115 }
0116
0117
0118
0119 double Lumi3DReWeighting::weight3D(const edm::EventBase &e) {
0120 using std::min;
0121
0122
0123
0124 Handle<std::vector<PileupSummaryInfo> > PupInfo;
0125 e.getByLabel(edm::InputTag("addPileupInfo"), PupInfo);
0126
0127 std::vector<PileupSummaryInfo>::const_iterator PVI;
0128
0129 int npm1 = -1;
0130 int np0 = -1;
0131 int npp1 = -1;
0132
0133 for (PVI = PupInfo->begin(); PVI != PupInfo->end(); ++PVI) {
0134 int BX = PVI->getBunchCrossing();
0135
0136 if (BX == -1) {
0137 npm1 = PVI->getPU_NumInteractions();
0138 }
0139 if (BX == 0) {
0140 np0 = PVI->getPU_NumInteractions();
0141 }
0142 if (BX == 1) {
0143 npp1 = PVI->getPU_NumInteractions();
0144 }
0145 }
0146
0147 npm1 = min(npm1, 49);
0148 np0 = min(np0, 49);
0149 npp1 = min(npp1, 49);
0150
0151
0152
0153 return Weight3D_[npm1][np0][npp1];
0154 }
0155
0156 void Lumi3DReWeighting::weight3D_set(std::string generatedFile,
0157 std::string dataFile,
0158 std::string GenHistName,
0159 std::string DataHistName,
0160 std::string WeightOutputFile) {
0161 generatedFileName_ = generatedFile;
0162 dataFileName_ = dataFile;
0163 GenHistName_ = GenHistName;
0164 DataHistName_ = DataHistName;
0165 weightFileName_ = WeightOutputFile;
0166
0167 std::cout << " seting values: " << generatedFileName_ << " " << dataFileName_ << " " << GenHistName_ << " "
0168 << DataHistName_ << std::endl;
0169
0170 generatedFile_ = std::make_shared<TFile>(generatedFileName_.c_str());
0171 dataFile_ = std::make_shared<TFile>(dataFileName_.c_str());
0172
0173 std::shared_ptr<TH1> Data_temp =
0174 std::shared_ptr<TH1>((static_cast<TH1 *>(dataFile_->Get(DataHistName_.c_str())->Clone())));
0175
0176 std::shared_ptr<TH1> MC_temp =
0177 std::shared_ptr<TH1>((static_cast<TH1 *>(generatedFile_->Get(GenHistName_.c_str())->Clone())));
0178
0179 MC_distr_ = std::shared_ptr<TH1>((static_cast<TH1 *>(generatedFile_->Get(GenHistName_.c_str())->Clone())));
0180 Data_distr_ = std::shared_ptr<TH1>((static_cast<TH1 *>(dataFile_->Get(DataHistName_.c_str())->Clone())));
0181
0182
0183
0184
0185
0186 Data_distr_->Scale(1.0 / Data_distr_->Integral());
0187 MC_distr_->Scale(1.0 / MC_distr_->Integral());
0188 }
0189
0190 void Lumi3DReWeighting::weight3D_init(float ScaleFactor) {
0191
0192
0193
0194
0195 TH3D *WHist = new TH3D("WHist", "3D weights", 50, 0., 50., 50, 0., 50., 50, 0., 50.);
0196 TH3D *DHist = new TH3D("DHist", "3D weights", 50, 0., 50., 50, 0., 50., 50, 0., 50.);
0197 TH3D *MHist = new TH3D("MHist", "3D weights", 50, 0., 50., 50, 0., 50., 50, 0., 50.);
0198
0199 using std::min;
0200
0201 if (MC_distr_->GetEntries() == 0) {
0202 std::cout << " MC and Data distributions are not initialized! You must call the Lumi3DReWeighting constructor. "
0203 << std::endl;
0204 }
0205
0206
0207
0208 double MC_ints[50][50][50];
0209 double Data_ints[50][50][50];
0210
0211 for (int i = 0; i < 50; i++) {
0212 for (int j = 0; j < 50; j++) {
0213 for (int k = 0; k < 50; k++) {
0214 MC_ints[i][j][k] = 0.;
0215 Data_ints[i][j][k] = 0.;
0216 }
0217 }
0218 }
0219
0220 double factorial[50];
0221 double PowerSer[50];
0222 double base = 1.;
0223
0224 factorial[0] = 1.;
0225 PowerSer[0] = 1.;
0226
0227 for (int i = 1; i < 50; ++i) {
0228 base = base * float(i);
0229 factorial[i] = base;
0230 }
0231
0232 double x;
0233 double xweight;
0234 double probi, probj, probk;
0235 double Expval, mean;
0236 int xi;
0237
0238
0239
0240 int NMCbin = MC_distr_->GetNbinsX();
0241
0242 for (int jbin = 1; jbin < NMCbin + 1; jbin++) {
0243 x = MC_distr_->GetBinCenter(jbin);
0244 xweight = MC_distr_->GetBinContent(jbin);
0245
0246
0247 xi = int(x);
0248
0249
0250
0251 mean = double(xi);
0252
0253 if (mean < 0.) {
0254 throw cms::Exception("BadInputValue") << " Your histogram generates MC luminosity values less than zero!"
0255 << " Please Check. Terminating." << std::endl;
0256 }
0257
0258 if (mean == 0.) {
0259 Expval = 1.;
0260 } else {
0261 Expval = exp(-1. * mean);
0262 }
0263
0264 base = 1.;
0265
0266 for (int i = 1; i < 50; ++i) {
0267 base = base * mean;
0268 PowerSer[i] = base;
0269 }
0270
0271
0272
0273 for (int i = 0; i < 50; i++) {
0274 probi = PowerSer[i] / factorial[i] * Expval;
0275 for (int j = 0; j < 50; j++) {
0276 probj = PowerSer[j] / factorial[j] * Expval;
0277 for (int k = 0; k < 50; k++) {
0278 probk = PowerSer[k] / factorial[k] * Expval;
0279
0280 MC_ints[i][j][k] = MC_ints[i][j][k] + probi * probj * probk * xweight;
0281 }
0282 }
0283 }
0284 }
0285
0286 int NDatabin = Data_distr_->GetNbinsX();
0287
0288 for (int jbin = 1; jbin < NDatabin + 1; jbin++) {
0289 mean = (Data_distr_->GetBinCenter(jbin)) * ScaleFactor;
0290 xweight = Data_distr_->GetBinContent(jbin);
0291
0292
0293
0294 if (mean < 0.) {
0295 throw cms::Exception("BadInputValue") << " Your histogram generates Data luminosity values less than zero!"
0296 << " Please Check. Terminating." << std::endl;
0297 }
0298
0299 if (mean == 0.) {
0300 Expval = 1.;
0301 } else {
0302 Expval = exp(-1. * mean);
0303 }
0304
0305 base = 1.;
0306
0307 for (int i = 1; i < 50; ++i) {
0308 base = base * mean;
0309 PowerSer[i] = base;
0310 }
0311
0312
0313
0314 for (int i = 0; i < 50; i++) {
0315 probi = PowerSer[i] / factorial[i] * Expval;
0316 for (int j = 0; j < 50; j++) {
0317 probj = PowerSer[j] / factorial[j] * Expval;
0318 for (int k = 0; k < 50; k++) {
0319 probk = PowerSer[k] / factorial[k] * Expval;
0320
0321 Data_ints[i][j][k] = Data_ints[i][j][k] + probi * probj * probk * xweight;
0322 }
0323 }
0324 }
0325 }
0326
0327 for (int i = 0; i < 50; i++) {
0328
0329 for (int j = 0; j < 50; j++) {
0330 for (int k = 0; k < 50; k++) {
0331 if ((MC_ints[i][j][k]) > 0.) {
0332 Weight3D_[i][j][k] = Data_ints[i][j][k] / MC_ints[i][j][k];
0333 } else {
0334 Weight3D_[i][j][k] = 0.;
0335 }
0336 WHist->SetBinContent(i + 1, j + 1, k + 1, Weight3D_[i][j][k]);
0337 DHist->SetBinContent(i + 1, j + 1, k + 1, Data_ints[i][j][k]);
0338 MHist->SetBinContent(i + 1, j + 1, k + 1, MC_ints[i][j][k]);
0339
0340 }
0341
0342 }
0343 }
0344
0345 if (!weightFileName_.empty()) {
0346 std::cout << " 3D Weight Matrix initialized! " << std::endl;
0347 std::cout << " Writing weights to file " << weightFileName_ << " for re-use... " << std::endl;
0348
0349 TFile *outfile = new TFile(weightFileName_.c_str(), "RECREATE");
0350 WHist->Write();
0351 MHist->Write();
0352 DHist->Write();
0353 outfile->Write();
0354 outfile->Close();
0355 outfile->Delete();
0356 }
0357
0358 return;
0359 }
0360
0361 void Lumi3DReWeighting::weight3D_init(std::string WeightFileName) {
0362 TFile *infile = new TFile(WeightFileName.c_str());
0363 TH1F *WHist = (TH1F *)infile->Get("WHist");
0364
0365
0366 if (!WHist) {
0367 throw cms::Exception("HistogramNotFound") << " Could not find the histogram WHist in the file "
0368 << "in the file " << WeightFileName << "." << std::endl;
0369 }
0370
0371 for (int i = 0; i < 50; i++) {
0372
0373 for (int j = 0; j < 50; j++) {
0374 for (int k = 0; k < 50; k++) {
0375 Weight3D_[i][j][k] = WHist->GetBinContent(i + 1, j + 1, k + 1);
0376
0377 }
0378
0379 }
0380 }
0381
0382 std::cout << " 3D Weight Matrix initialized! " << std::endl;
0383
0384 return;
0385 }
0386
0387 void Lumi3DReWeighting::weight3D_init(std::string MCWeightFileName, std::string DataWeightFileName) {
0388 TFile *infileMC = new TFile(MCWeightFileName.c_str());
0389 TH3D *MHist = (TH3D *)infileMC->Get("MHist");
0390
0391
0392 if (!MHist) {
0393 throw cms::Exception("HistogramNotFound") << " Could not find the histogram MHist in the file "
0394 << "in the file " << MCWeightFileName << "." << std::endl;
0395 }
0396
0397 TFile *infileD = new TFile(DataWeightFileName.c_str());
0398 TH3D *DHist = (TH3D *)infileD->Get("DHist");
0399
0400
0401 if (!DHist) {
0402 throw cms::Exception("HistogramNotFound") << " Could not find the histogram DHist in the file "
0403 << "in the file " << DataWeightFileName << "." << std::endl;
0404 }
0405
0406 for (int i = 0; i < 50; i++) {
0407 for (int j = 0; j < 50; j++) {
0408 for (int k = 0; k < 50; k++) {
0409 Weight3D_[i][j][k] = DHist->GetBinContent(i + 1, j + 1, k + 1) / MHist->GetBinContent(i + 1, j + 1, k + 1);
0410 }
0411 }
0412 }
0413
0414 std::cout << " 3D Weight Matrix initialized! " << std::endl;
0415
0416 return;
0417 }
0418
0419 #endif