Utilities/src/LumiReWeighting.cc

0001 #ifndef PhysicsTools_Utilities_interface_LumiReWeighting_cc
0002 #define PhysicsTools_Utilities_interface_LumiReWeighting_cc
0003
0004 /**
0005   \class    LumiReWeighting LumiReWeighting.h "PhysicsTools/Utilities/interface/LumiReWeighting.h"
0006   \brief    Class to provide lumi weighting for analyzers to weight "flat-to-N" MC samples to data
0007
0008   This class will trivially take two histograms:
0009   1. The generated "flat-to-N" distributions from a given processing (or any other generated input)
0010   2. A histogram generated from the "estimatePileup" macro here:
0011
0012   https://twiki.cern.ch/twiki/bin/view/CMS/LumiCalc#How_to_use_script_estimatePileup
0013
0014   and produce weights to convert the input distribution (1) to the latter (2).
0015
0016   \author Salvatore Rappoccio, modified by Mike Hildreth
0017
0018 */
0019 #include "TFile.h"
0020 #include "TH1.h"
0021 #include "TRandom1.h"
0022 #include "TRandom2.h"
0023 #include "TRandom3.h"
0024 #include "TStopwatch.h"
0025 #include <algorithm>
0026 #include <iostream>
0027 #include <memory>
0028 #include <string>
0029
0030 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
0031 #include "PhysicsTools/Utilities/interface/LumiReWeighting.h"
0032 #include "DataFormats/Common/interface/Handle.h"
0033 #include "DataFormats/Provenance/interface/ProcessHistory.h"
0034 #include "DataFormats/Provenance/interface/Provenance.h"
0035 #include "FWCore/Common/interface/EventBase.h"
0036
0037 using namespace edm;
0038
0039 LumiReWeighting::LumiReWeighting(std::string generatedFile,
0040                                  std::string dataFile,
0041                                  std::string GenHistName,
0042                                  std::string DataHistName,
0043                                  const edm::InputTag& PileupSumInfoInputTag)
0044     : generatedFileName_(generatedFile),
0045       dataFileName_(dataFile),
0046       GenHistName_(GenHistName),
0047       DataHistName_(DataHistName),
0048       pileupSumInfoTag_(PileupSumInfoInputTag) {
0049   generatedFile_ = std::make_shared<TFile>(generatedFileName_.c_str());  //MC distribution
0050   dataFile_ = std::make_shared<TFile>(dataFileName_.c_str());            //Data distribution
0051
0052   Data_distr_ = std::shared_ptr<TH1>((static_cast<TH1*>(dataFile_->Get(DataHistName_.c_str())->Clone())));
0053   MC_distr_ = std::shared_ptr<TH1>((static_cast<TH1*>(generatedFile_->Get(GenHistName_.c_str())->Clone())));
0054
0055   // MC * data/MC = data, so the weights are data/MC:
0056
0057   // normalize both histograms first
0058
0059   Data_distr_->Scale(1.0 / Data_distr_->Integral());
0060   MC_distr_->Scale(1.0 / MC_distr_->Integral());
0061
0062   weights_ = std::shared_ptr<TH1>(static_cast<TH1*>(Data_distr_->Clone()));
0063
0064   weights_->SetName("lumiWeights");
0065
0066   TH1* den = dynamic_cast<TH1*>(MC_distr_->Clone());
0067
0068   //den->Scale(1.0/ den->Integral());
0069
0070   weights_->Divide(den);  // so now the average weight should be 1.0
0071
0072   std::cout << " Lumi/Pileup Reweighting: Computed Weights per In-Time Nint " << std::endl;
0073
0074   int NBins = weights_->GetNbinsX();
0075
0076   for (int ibin = 1; ibin < NBins + 1; ++ibin) {
0077     std::cout << "   " << ibin - 1 << " " << weights_->GetBinContent(ibin) << std::endl;
0078   }
0079
0080   FirstWarning_ = true;
0081   OldLumiSection_ = -1;
0082 }
0083
0084 LumiReWeighting::LumiReWeighting(const std::vector<float>& MC_distr,
0085                                  const std::vector<float>& Lumi_distr,
0086                                  const edm::InputTag& PileupSumInfoInputTag)
0087     : pileupSumInfoTag_(PileupSumInfoInputTag) {
0088   // no histograms for input: use vectors
0089
0090   // now, make histograms out of them:
0091
0092   // first, check they are the same size...
0093
0094   if (MC_distr.size() != Lumi_distr.size()) {
0095     std::cerr << "ERROR: LumiReWeighting: input vectors have different sizes. Quitting... \n";
0096     return;
0097   }
0098
0099   Int_t NBins = MC_distr.size();
0100
0101   MC_distr_ = std::shared_ptr<TH1>(new TH1F("MC_distr", "MC dist", NBins, -0.5, float(NBins) - 0.5));
0102   Data_distr_ = std::shared_ptr<TH1>(new TH1F("Data_distr", "Data dist", NBins, -0.5, float(NBins) - 0.5));
0103
0104   weights_ = std::shared_ptr<TH1>(new TH1F("luminumer", "luminumer", NBins, -0.5, float(NBins) - 0.5));
0105   TH1* den = new TH1F("lumidenom", "lumidenom", NBins, -0.5, float(NBins) - 0.5);
0106
0107   for (int ibin = 1; ibin < NBins + 1; ++ibin) {
0108     weights_->SetBinContent(ibin, Lumi_distr[ibin - 1]);
0109     Data_distr_->SetBinContent(ibin, Lumi_distr[ibin - 1]);
0110     den->SetBinContent(ibin, MC_distr[ibin - 1]);
0111     MC_distr_->SetBinContent(ibin, MC_distr[ibin - 1]);
0112   }
0113
0114   // check integrals, make sure things are normalized
0115
0116   float deltaH = weights_->Integral();
0117   if (fabs(1.0 - deltaH) > 0.02) {  //*OOPS*...
0118     weights_->Scale(1.0 / weights_->Integral());
0119     Data_distr_->Scale(1.0 / Data_distr_->Integral());
0120   }
0121   float deltaMC = den->Integral();
0122   if (fabs(1.0 - deltaMC) > 0.02) {
0123     den->Scale(1.0 / den->Integral());
0124     MC_distr_->Scale(1.0 / MC_distr_->Integral());
0125   }
0126
0127   weights_->Divide(den);  // so now the average weight should be 1.0
0128
0129   std::cout << " Lumi/Pileup Reweighting: Computed Weights per In-Time Nint " << std::endl;
0130
0131   for (int ibin = 1; ibin < NBins + 1; ++ibin) {
0132     std::cout << "   " << ibin - 1 << " " << weights_->GetBinContent(ibin) << std::endl;
0133   }
0134
0135   FirstWarning_ = true;
0136   OldLumiSection_ = -1;
0137 }
0138
0139 double LumiReWeighting::weight(int npv) {
0140   int bin = weights_->GetXaxis()->FindBin(npv);
0141   return weights_->GetBinContent(bin);
0142 }
0143
0144 double LumiReWeighting::weight(float npv) {
0145   int bin = weights_->GetXaxis()->FindBin(npv);
0146   return weights_->GetBinContent(bin);
0147 }
0148
0149 // This version of weight does all of the work for you, assuming you want to re-weight
0150 // using the true number of interactions in the in-time beam crossing.
0151
0152 double LumiReWeighting::weight(const edm::EventBase& e) {
0153   if (FirstWarning_) {
0154     e.processHistory();
0155     //    SetFirstFalse();
0156     FirstWarning_ = false;
0157   }
0158   Handle<std::vector<PileupSummaryInfo> > PupInfo;
0159   e.getByLabel(pileupSumInfoTag_, PupInfo);
0160
0161   std::vector<PileupSummaryInfo>::const_iterator PVI;
0162
0163   int npv = -1;
0164   for (PVI = PupInfo->begin(); PVI != PupInfo->end(); ++PVI) {
0165     int BX = PVI->getBunchCrossing();
0166
0167     if (BX == 0) {
0168       npv = PVI->getPU_NumInteractions();
0169       continue;
0170     }
0171   }
0172
0173   if (npv < 0)
0174     std::cerr << " no in-time beam crossing found\n! ";
0175
0176   return weight(npv);
0177 }
0178
0179 // Use this routine to re-weight out-of-time pileup to match the in-time distribution
0180 // As of May 2011, CMS is only sensitive to a bunch that is 50ns "late", which corresponds to
0181 // BunchCrossing +1.  So, we use that here for re-weighting.
0182
0183 double LumiReWeighting::weightOOT(const edm::EventBase& e) {
0184   //int Run = e.run();
0185   int LumiSection = e.luminosityBlock();
0186   // do some caching here, attempt to catch file boundaries
0187
0188   if (LumiSection != OldLumiSection_) {
0189     e.processHistory();  // keep the function call
0190     OldLumiSection_ = LumiSection;
0191   }
0192   // find the pileup summary information
0193
0194   Handle<std::vector<PileupSummaryInfo> > PupInfo;
0195   e.getByLabel(pileupSumInfoTag_, PupInfo);
0196
0197   std::vector<PileupSummaryInfo>::const_iterator PVI;
0198
0199   int npv = -1;
0200   int npv50ns = -1;
0201
0202   for (PVI = PupInfo->begin(); PVI != PupInfo->end(); ++PVI) {
0203     int BX = PVI->getBunchCrossing();
0204
0205     if (BX == 0) {
0206       npv = PVI->getPU_NumInteractions();
0207     }
0208
0209     if (BX == 1) {
0210       npv50ns = PVI->getPU_NumInteractions();
0211     }
0212   }
0213
0214   // Note: for the "uncorrelated" out-of-time pileup, reweighting is only done on the 50ns
0215   // "late" bunch (BX=+1), since that is basically the only one that matters in terms of
0216   // energy deposition.
0217
0218   if (npv < 0) {
0219     std::cerr << " no in-time beam crossing found\n! ";
0220     std::cerr << " Returning event weight=0\n! ";
0221     return 0.;
0222   }
0223   if (npv50ns < 0) {
0224     std::cerr << " no out-of-time beam crossing found\n! ";
0225     std::cerr << " Returning event weight=0\n! ";
0226     return 0.;
0227   }
0228
0229   int bin = weights_->GetXaxis()->FindBin(npv);
0230
0231   double inTimeWeight = weights_->GetBinContent(bin);
0232
0233   double TotalWeight = 1.0;
0234
0235   TotalWeight = inTimeWeight;
0236
0237   return TotalWeight;
0238 }
0239
0240 #endif