edm::Lumi3DReWeighting Class Reference

#include <Lumi3DReWeighting.h>

Public Member Functions
	Lumi3DReWeighting (std::string generatedFile, std::string dataFile, std::string GenHistName, std::string DataHistName, std::string WeightOutputFile)
	Lumi3DReWeighting (const std::vector< float > &MC_distr, const std::vector< float > &Lumi_distr, std::string WeightOutputFile)
	Lumi3DReWeighting ()
double	weight3D (const edm::EventBase &e)
double	weight3D (int, int, int)
void	weight3D_init (float Scale)
void	weight3D_init (std::string WeightFileName)
void	weight3D_init (std::string MCFileName, std::string DataFileName)
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.

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

    : generatedFileName_(generatedFile),
      dataFileName_(dataFile),
      GenHistName_(GenHistName),
      DataHistName_(DataHistName),
      weightFileName_(WeightOutputFile) {
  generatedFile_ = std::make_shared<TFile>(generatedFileName_.c_str());  //MC distribution
  dataFile_ = std::make_shared<TFile>(dataFileName_.c_str());            //Data distribution

  std::shared_ptr<TH1> Data_temp =
      std::shared_ptr<TH1>((static_cast<TH1 *>(dataFile_->Get(DataHistName_.c_str())->Clone())));

  std::shared_ptr<TH1> MC_temp =
      std::shared_ptr<TH1>((static_cast<TH1 *>(generatedFile_->Get(GenHistName_.c_str())->Clone())));

  MC_distr_ = std::shared_ptr<TH1>((static_cast<TH1 *>(generatedFile_->Get(GenHistName_.c_str())->Clone())));
  Data_distr_ = std::shared_ptr<TH1>((static_cast<TH1 *>(dataFile_->Get(DataHistName_.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());
}

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.

References Data_distr_, MC_distr_, and weightFileName_.

                                                                      {
  weightFileName_ = WeightOutputFile;

  // no histograms for input: use vectors

  // now, make histograms out of them:

  Int_t NMCBins = MC_distr.size();

  MC_distr_ = std::shared_ptr<TH1>(new TH1F("MC_distr", "MC dist", NMCBins, 0., float(NMCBins)));

  Int_t NDBins = Lumi_distr.size();

  Data_distr_ = std::shared_ptr<TH1>(new TH1F("Data_distr", "Data dist", NDBins, 0., float(NDBins)));

  for (int ibin = 1; ibin < NMCBins + 1; ++ibin) {
    MC_distr_->SetBinContent(ibin, MC_distr[ibin - 1]);
  }

  for (int ibin = 1; ibin < NDBins + 1; ++ibin) {
    Data_distr_->SetBinContent(ibin, Lumi_distr[ibin - 1]);
  }

  // check integrals, make sure things are normalized

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

Lumi3DReWeighting() [3/3]

edm::Lumi3DReWeighting::Lumi3DReWeighting ( )

inline

Definition at line 38 of file Lumi3DReWeighting.h.

{};

Member Function Documentation

weight3D() [1/2]

double Lumi3DReWeighting::weight3D

(

const edm::EventBase &

e

)

Definition at line 119 of file Lumi3DReWeighting.cc.

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

                                                        {
  using std::min;

  // get pileup summary information

  Handle<std::vector<PileupSummaryInfo> > PupInfo;
  e.getByLabel(edm::InputTag("addPileupInfo"), PupInfo);

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

  int npm1 = -1;
  int np0 = -1;
  int npp1 = -1;

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

    if (BX == -1) {
      npm1 = PVI->getPU_NumInteractions();
    }
    if (BX == 0) {
      np0 = PVI->getPU_NumInteractions();
    }
    if (BX == 1) {
      npp1 = PVI->getPU_NumInteractions();
    }
  }

  npm1 = min(npm1, 49);
  np0 = min(np0, 49);
  npp1 = min(npp1, 49);

  //  std::cout << " vertices " << npm1 << " " << np0 << " " << npp1 << std::endl;

  return Weight3D_[npm1][np0][npp1];
}

weight3D() [2/2]

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

Definition at line 107 of file Lumi3DReWeighting.cc.

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

                                                            {
  using std::min;

  int npm1 = min(pv1, 49);
  int np0 = min(pv2, 49);
  int npp1 = min(pv3, 49);

  return Weight3D_[npm1][np0][npp1];
}

weight3D_init() [1/3]

void Lumi3DReWeighting::weight3D_init

(

float

Scale

)

Definition at line 190 of file Lumi3DReWeighting.cc.

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

                                                       {
  // Scale factor is used to shift target distribution (i.e. luminosity scale)  1. = no shift

  //create histogram to write output weights, save pain of generating them again...

  TH3D *WHist = new TH3D("WHist", "3D weights", 50, 0., 50., 50, 0., 50., 50, 0., 50.);
  TH3D *DHist = new TH3D("DHist", "3D weights", 50, 0., 50., 50, 0., 50., 50, 0., 50.);
  TH3D *MHist = new TH3D("MHist", "3D weights", 50, 0., 50., 50, 0., 50., 50, 0., 50.);

  using std::min;

  if (MC_distr_->GetEntries() == 0) {
    std::cout << " MC and Data distributions are not initialized! You must call the Lumi3DReWeighting constructor. "
              << std::endl;
  }

  // arrays for storing number of interactions

  double MC_ints[50][50][50];
  double Data_ints[50][50][50];

  for (int i = 0; i < 50; i++) {
    for (int j = 0; j < 50; j++) {
      for (int k = 0; k < 50; k++) {
        MC_ints[i][j][k] = 0.;
        Data_ints[i][j][k] = 0.;
      }
    }
  }

  double factorial[50];
  double PowerSer[50];
  double base = 1.;

  factorial[0] = 1.;
  PowerSer[0] = 1.;

  for (int i = 1; i < 50; ++i) {
    base = base * float(i);
    factorial[i] = base;
  }

  double x;
  double xweight;
  double probi, probj, probk;
  double Expval, mean;
  int xi;

  // Get entries for Data, MC, fill arrays:

  int NMCbin = MC_distr_->GetNbinsX();

  for (int jbin = 1; jbin < NMCbin + 1; jbin++) {
    x = MC_distr_->GetBinCenter(jbin);
    xweight = MC_distr_->GetBinContent(jbin);  //use as weight for matrix

    //for Summer 11, we have this int feature:
    xi = int(x);

    // Generate Poisson distribution for each value of the mean

    mean = double(xi);

    if (mean < 0.) {
      throw cms::Exception("BadInputValue") << " Your histogram generates MC luminosity values less than zero!"
                                            << " Please Check.  Terminating." << std::endl;
    }

    if (mean == 0.) {
      Expval = 1.;
    } else {
      Expval = exp(-1. * mean);
    }

    base = 1.;

    for (int i = 1; i < 50; ++i) {
      base = base * mean;
      PowerSer[i] = base;  // PowerSer is mean^i
    }

    // compute poisson probability for each Nvtx in weight matrix

    for (int i = 0; i < 50; i++) {
      probi = PowerSer[i] / factorial[i] * Expval;
      for (int j = 0; j < 50; j++) {
        probj = PowerSer[j] / factorial[j] * Expval;
        for (int k = 0; k < 50; k++) {
          probk = PowerSer[k] / factorial[k] * Expval;
          // joint probability is product of event weights multiplied by weight of input distribution bin
          MC_ints[i][j][k] = MC_ints[i][j][k] + probi * probj * probk * xweight;
        }
      }
    }
  }

  int NDatabin = Data_distr_->GetNbinsX();

  for (int jbin = 1; jbin < NDatabin + 1; jbin++) {
    mean = (Data_distr_->GetBinCenter(jbin)) * ScaleFactor;
    xweight = Data_distr_->GetBinContent(jbin);

    // Generate poisson distribution for each value of the mean

    if (mean < 0.) {
      throw cms::Exception("BadInputValue") << " Your histogram generates Data luminosity values less than zero!"
                                            << " Please Check.  Terminating." << std::endl;
    }

    if (mean == 0.) {
      Expval = 1.;
    } else {
      Expval = exp(-1. * mean);
    }

    base = 1.;

    for (int i = 1; i < 50; ++i) {
      base = base * mean;
      PowerSer[i] = base;
    }

    // compute poisson probability for each Nvtx in weight matrix

    for (int i = 0; i < 50; i++) {
      probi = PowerSer[i] / factorial[i] * Expval;
      for (int j = 0; j < 50; j++) {
        probj = PowerSer[j] / factorial[j] * Expval;
        for (int k = 0; k < 50; k++) {
          probk = PowerSer[k] / factorial[k] * Expval;
          // joint probability is product of event weights multiplied by weight of input distribution bin
          Data_ints[i][j][k] = Data_ints[i][j][k] + probi * probj * probk * xweight;
        }
      }
    }
  }

  for (int i = 0; i < 50; i++) {
    //if(i<5) std::cout << "i = " << i << std::endl;
    for (int j = 0; j < 50; j++) {
      for (int k = 0; k < 50; k++) {
        if ((MC_ints[i][j][k]) > 0.) {
          Weight3D_[i][j][k] = Data_ints[i][j][k] / MC_ints[i][j][k];
        } else {
          Weight3D_[i][j][k] = 0.;
        }
        WHist->SetBinContent(i + 1, j + 1, k + 1, Weight3D_[i][j][k]);
        DHist->SetBinContent(i + 1, j + 1, k + 1, Data_ints[i][j][k]);
        MHist->SetBinContent(i + 1, j + 1, k + 1, MC_ints[i][j][k]);
        //  if(i<5 && j<5 && k<5) std::cout << Weight3D_[i][j][k] << " " ;
      }
      //      if(i<5 && j<5) std::cout << std::endl;
    }
  }

  if (!weightFileName_.empty()) {
    std::cout << " 3D Weight Matrix initialized! " << std::endl;
    std::cout << " Writing weights to file " << weightFileName_ << " for re-use...  " << std::endl;

    TFile *outfile = new TFile(weightFileName_.c_str(), "RECREATE");
    WHist->Write();
    MHist->Write();
    DHist->Write();
    outfile->Write();
    outfile->Close();
    outfile->Delete();
  }

  return;
}

weight3D_init() [2/3]

void Lumi3DReWeighting::weight3D_init

(

std::string

WeightFileName

)

Definition at line 361 of file Lumi3DReWeighting.cc.

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

                                                              {
  TFile *infile = new TFile(WeightFileName.c_str());
  TH1F *WHist = (TH1F *)infile->Get("WHist");

  // Check if the histogram exists
  if (!WHist) {
    throw cms::Exception("HistogramNotFound") << " Could not find the histogram WHist in the file "
                                              << "in the file " << WeightFileName << "." << std::endl;
  }

  for (int i = 0; i < 50; i++) {
    //    if(i<5) std::cout << "i = " << i << std::endl;
    for (int j = 0; j < 50; j++) {
      for (int k = 0; k < 50; k++) {
        Weight3D_[i][j][k] = WHist->GetBinContent(i + 1, j + 1, k + 1);
        //  if(i<5 && j<5 && k<5) std::cout << Weight3D_[i][j][k] << " ";
      }
      //      if(i<5 && j<5) std::cout << std::endl;
    }
  }

  std::cout << " 3D Weight Matrix initialized! " << std::endl;

  return;
}

weight3D_init() [3/3]

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

Definition at line 387 of file Lumi3DReWeighting.cc.

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

                                                                                              {
  TFile *infileMC = new TFile(MCWeightFileName.c_str());
  TH3D *MHist = (TH3D *)infileMC->Get("MHist");

  // Check if the histogram exists
  if (!MHist) {
    throw cms::Exception("HistogramNotFound") << " Could not find the histogram MHist in the file "
                                              << "in the file " << MCWeightFileName << "." << std::endl;
  }

  TFile *infileD = new TFile(DataWeightFileName.c_str());
  TH3D *DHist = (TH3D *)infileD->Get("DHist");

  // Check if the histogram exists
  if (!DHist) {
    throw cms::Exception("HistogramNotFound") << " Could not find the histogram DHist in the file "
                                              << "in the file " << DataWeightFileName << "." << std::endl;
  }

  for (int i = 0; i < 50; i++) {
    for (int j = 0; j < 50; j++) {
      for (int k = 0; k < 50; k++) {
        Weight3D_[i][j][k] = DHist->GetBinContent(i + 1, j + 1, k + 1) / MHist->GetBinContent(i + 1, j + 1, k + 1);
      }
    }
  }

  std::cout << " 3D Weight Matrix initialized! " << std::endl;

  return;
}

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.

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

                                                                 {
  generatedFileName_ = generatedFile;
  dataFileName_ = dataFile;
  GenHistName_ = GenHistName;
  DataHistName_ = DataHistName;
  weightFileName_ = WeightOutputFile;

  std::cout << " seting values: " << generatedFileName_ << " " << dataFileName_ << " " << GenHistName_ << " "
            << DataHistName_ << std::endl;

  generatedFile_ = std::make_shared<TFile>(generatedFileName_.c_str());  //MC distribution
  dataFile_ = std::make_shared<TFile>(dataFileName_.c_str());            //Data distribution

  std::shared_ptr<TH1> Data_temp =
      std::shared_ptr<TH1>((static_cast<TH1 *>(dataFile_->Get(DataHistName_.c_str())->Clone())));

  std::shared_ptr<TH1> MC_temp =
      std::shared_ptr<TH1>((static_cast<TH1 *>(generatedFile_->Get(GenHistName_.c_str())->Clone())));

  MC_distr_ = std::shared_ptr<TH1>((static_cast<TH1 *>(generatedFile_->Get(GenHistName_.c_str())->Clone())));
  Data_distr_ = std::shared_ptr<TH1>((static_cast<TH1 *>(dataFile_->Get(DataHistName_.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());
}

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.