LumiReWeighting.cc

Go to the documentation of this file.

#ifndef PhysicsTools_Utilities_interface_LumiReWeighting_cc
#define PhysicsTools_Utilities_interface_LumiReWeighting_cc

#include "TFile.h"
#include "TH1.h"
#include "TRandom1.h"
#include "TRandom2.h"
#include "TRandom3.h"
#include "TStopwatch.h"
#include <algorithm>
#include <iostream>
#include <memory>
#include <string>

#include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
#include "PhysicsTools/Utilities/interface/LumiReWeighting.h"
#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/Provenance/interface/ProcessHistory.h"
#include "DataFormats/Provenance/interface/Provenance.h"
#include "FWCore/Common/interface/EventBase.h"

using namespace edm;

LumiReWeighting::LumiReWeighting(std::string generatedFile,
                                 std::string dataFile,
                                 std::string GenHistName,
                                 std::string DataHistName,
                                 const edm::InputTag& PileupSumInfoInputTag)
    : generatedFileName_(generatedFile),
      dataFileName_(dataFile),
      GenHistName_(GenHistName),
      DataHistName_(DataHistName),
      pileupSumInfoTag_(PileupSumInfoInputTag) {
  generatedFile_ = std::make_shared<TFile>(generatedFileName_.c_str());  //MC distribution
  dataFile_ = std::make_shared<TFile>(dataFileName_.c_str());            //Data distribution

  Data_distr_ = std::shared_ptr<TH1>((static_cast<TH1*>(dataFile_->Get(DataHistName_.c_str())->Clone())));
  MC_distr_ = std::shared_ptr<TH1>((static_cast<TH1*>(generatedFile_->Get(GenHistName_.c_str())->Clone())));

  // MC * data/MC = data, so the weights are data/MC:

  // normalize both histograms first

  Data_distr_->Scale(1.0 / Data_distr_->Integral());
  MC_distr_->Scale(1.0 / MC_distr_->Integral());

  weights_ = std::shared_ptr<TH1>(static_cast<TH1*>(Data_distr_->Clone()));

  weights_->SetName("lumiWeights");

  TH1* den = dynamic_cast<TH1*>(MC_distr_->Clone());

  //den->Scale(1.0/ den->Integral());

  weights_->Divide(den);  // so now the average weight should be 1.0

  std::cout << " Lumi/Pileup Reweighting: Computed Weights per In-Time Nint " << std::endl;

  int NBins = weights_->GetNbinsX();

  for (int ibin = 1; ibin < NBins + 1; ++ibin) {
    std::cout << "   " << ibin - 1 << " " << weights_->GetBinContent(ibin) << std::endl;
  }

  FirstWarning_ = true;
  OldLumiSection_ = -1;
}

LumiReWeighting::LumiReWeighting(const std::vector<float>& MC_distr,
                                 const std::vector<float>& Lumi_distr,
                                 const edm::InputTag& PileupSumInfoInputTag)
    : pileupSumInfoTag_(PileupSumInfoInputTag) {
  // no histograms for input: use vectors

  // now, make histograms out of them:

  // first, check they are the same size...

  if (MC_distr.size() != Lumi_distr.size()) {
    std::cerr << "ERROR: LumiReWeighting: input vectors have different sizes. Quitting... \n";
    return;
  }

  Int_t NBins = MC_distr.size();

  MC_distr_ = std::shared_ptr<TH1>(new TH1F("MC_distr", "MC dist", NBins, -0.5, float(NBins) - 0.5));
  Data_distr_ = std::shared_ptr<TH1>(new TH1F("Data_distr", "Data dist", NBins, -0.5, float(NBins) - 0.5));

  weights_ = std::shared_ptr<TH1>(new TH1F("luminumer", "luminumer", NBins, -0.5, float(NBins) - 0.5));
  TH1* den = new TH1F("lumidenom", "lumidenom", NBins, -0.5, float(NBins) - 0.5);

  for (int ibin = 1; ibin < NBins + 1; ++ibin) {
    weights_->SetBinContent(ibin, Lumi_distr[ibin - 1]);
    Data_distr_->SetBinContent(ibin, Lumi_distr[ibin - 1]);
    den->SetBinContent(ibin, MC_distr[ibin - 1]);
    MC_distr_->SetBinContent(ibin, MC_distr[ibin - 1]);
  }

  // check integrals, make sure things are normalized

  float deltaH = weights_->Integral();
  if (fabs(1.0 - deltaH) > 0.02) {  //*OOPS*...
    weights_->Scale(1.0 / weights_->Integral());
    Data_distr_->Scale(1.0 / Data_distr_->Integral());
  }
  float deltaMC = den->Integral();
  if (fabs(1.0 - deltaMC) > 0.02) {
    den->Scale(1.0 / den->Integral());
    MC_distr_->Scale(1.0 / MC_distr_->Integral());
  }

  weights_->Divide(den);  // so now the average weight should be 1.0

  std::cout << " Lumi/Pileup Reweighting: Computed Weights per In-Time Nint " << std::endl;

  for (int ibin = 1; ibin < NBins + 1; ++ibin) {
    std::cout << "   " << ibin - 1 << " " << weights_->GetBinContent(ibin) << std::endl;
  }

  FirstWarning_ = true;
  OldLumiSection_ = -1;
}

double LumiReWeighting::weight(int npv) {
  int bin = weights_->GetXaxis()->FindBin(npv);
  return weights_->GetBinContent(bin);
}

double LumiReWeighting::weight(float npv) {
  int bin = weights_->GetXaxis()->FindBin(npv);
  return weights_->GetBinContent(bin);
}

// This version of weight does all of the work for you, assuming you want to re-weight
// using the true number of interactions in the in-time beam crossing.

double LumiReWeighting::weight(const edm::EventBase& e) {
  if (FirstWarning_) {
    e.processHistory();
    //    SetFirstFalse();
    FirstWarning_ = false;
  }
  Handle<std::vector<PileupSummaryInfo> > PupInfo;
  e.getByLabel(pileupSumInfoTag_, PupInfo);

  std::vector<PileupSummaryInfo>::const_iterator PVI;

  int npv = -1;
  for (PVI = PupInfo->begin(); PVI != PupInfo->end(); ++PVI) {
    int BX = PVI->getBunchCrossing();

    if (BX == 0) {
      npv = PVI->getPU_NumInteractions();
      continue;
    }
  }

  if (npv < 0)
    std::cerr << " no in-time beam crossing found\n! ";

  return weight(npv);
}

// Use this routine to re-weight out-of-time pileup to match the in-time distribution
// As of May 2011, CMS is only sensitive to a bunch that is 50ns "late", which corresponds to
// BunchCrossing +1.  So, we use that here for re-weighting.

double LumiReWeighting::weightOOT(const edm::EventBase& e) {
  //int Run = e.run();
  int LumiSection = e.luminosityBlock();
  // do some caching here, attempt to catch file boundaries

  if (LumiSection != OldLumiSection_) {
    e.processHistory();  // keep the function call
    OldLumiSection_ = LumiSection;
  }
  // find the pileup summary information

  Handle<std::vector<PileupSummaryInfo> > PupInfo;
  e.getByLabel(pileupSumInfoTag_, PupInfo);

  std::vector<PileupSummaryInfo>::const_iterator PVI;

  int npv = -1;
  int npv50ns = -1;

  for (PVI = PupInfo->begin(); PVI != PupInfo->end(); ++PVI) {
    int BX = PVI->getBunchCrossing();

    if (BX == 0) {
      npv = PVI->getPU_NumInteractions();
    }

    if (BX == 1) {
      npv50ns = PVI->getPU_NumInteractions();
    }
  }

  // Note: for the "uncorrelated" out-of-time pileup, reweighting is only done on the 50ns
  // "late" bunch (BX=+1), since that is basically the only one that matters in terms of
  // energy deposition.

  if (npv < 0) {
    std::cerr << " no in-time beam crossing found\n! ";
    std::cerr << " Returning event weight=0\n! ";
    return 0.;
  }
  if (npv50ns < 0) {
    std::cerr << " no out-of-time beam crossing found\n! ";
    std::cerr << " Returning event weight=0\n! ";
    return 0.;
  }

  int bin = weights_->GetXaxis()->FindBin(npv);

  double inTimeWeight = weights_->GetBinContent(bin);

  double TotalWeight = 1.0;

  TotalWeight = inTimeWeight;

  return TotalWeight;
}

#endif