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
boost::shared_ptr< TH1 >	Data_distr_
boost::shared_ptr< TFile >	dataFile_
std::string	dataFileName_
std::string	DataHistName_
boost::shared_ptr< TFile >	generatedFile_
std::string	generatedFileName_
std::string	GenHistName_
boost::shared_ptr< TH1 >	MC_distr_
double	Weight3D_ [50][50][50]
std::string	weightFileName_
boost::shared_ptr< TH1 >	weights_

Detailed Description

Definition at line 29 of file Lumi3DReWeighting.h.

Constructor & Destructor Documentation

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_ = boost::shared_ptr<TFile>( new TFile(generatedFileName_.c_str()) ); //MC distribution
    dataFile_      = boost::shared_ptr<TFile>( new TFile(dataFileName_.c_str()) );      //Data distribution

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

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


    MC_distr_ = boost::shared_ptr<TH1>(  (static_cast<TH1*>(generatedFile_->Get( GenHistName_.c_str() )->Clone() )) );
    Data_distr_ = boost::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::Lumi3DReWeighting	(	const std::vector< float > &	MC_distr,
		const std::vector< float > &	Lumi_distr,
		std::string	WeightOutputFile = ""
	)

Definition at line 71 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_ = boost::shared_ptr<TH1> ( new TH1F("MC_distr","MC dist",NMCBins,0., float(NMCBins)) );

  Int_t NDBins = Lumi_distr.size();

  Data_distr_ = boost::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());
  }

}

edm::Lumi3DReWeighting::Lumi3DReWeighting ( )

inline

Definition at line 40 of file Lumi3DReWeighting.h.

{ } ;

Member Function Documentation

double Lumi3DReWeighting::weight3D

(

const edm::EventBase &

e

)

Definition at line 124 of file Lumi3DReWeighting.cc.

References rpcdqm::BX, edm::EventBase::getByLabel(), min(), edm::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];
 
}

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

Definition at line 110 of file Lumi3DReWeighting.cc.

References min(), edm::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];

}

void Lumi3DReWeighting::weight3D_init

(

float

Scale

)

Definition at line 195 of file Lumi3DReWeighting.cc.

References newFWLiteAna::base, gather_cfg::cout, Data_distr_, Exception, create_public_lumi_plots::exp, factorial(), i, j, relval_steps::k, MC_distr_, timingPdfMaker::mean, min(), EdgesToViz::outfile, Weight3D_, weightFileName_, and x().

                                                         { 

  // 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;


}

void Lumi3DReWeighting::weight3D_init

(

std::string

WeightFileName

)

Definition at line 379 of file Lumi3DReWeighting.cc.

References gather_cfg::cout, Exception, i, EdgesToViz::infile, j, relval_steps::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;


}

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

Definition at line 409 of file Lumi3DReWeighting.cc.

References gather_cfg::cout, Exception, i, j, relval_steps::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;


}

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

Definition at line 165 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_ = boost::shared_ptr<TFile>( new TFile(generatedFileName_.c_str()) ); //MC distribution
  dataFile_      = boost::shared_ptr<TFile>( new TFile(dataFileName_.c_str()) );      //Data distribution
  
  boost::shared_ptr<TH1> Data_temp = boost::shared_ptr<TH1>(  (static_cast<TH1*>(dataFile_->Get( DataHistName_.c_str() )->Clone() )) );
  
  boost::shared_ptr<TH1> MC_temp = boost::shared_ptr<TH1>(  (static_cast<TH1*>(generatedFile_->Get( GenHistName_.c_str() )->Clone() )) );
  
  
  MC_distr_ = boost::shared_ptr<TH1>(  (static_cast<TH1*>(generatedFile_->Get( GenHistName_.c_str() )->Clone() )) );
  Data_distr_ = boost::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

boost::shared_ptr<TH1> edm::Lumi3DReWeighting::Data_distr_

protected

Definition at line 69 of file Lumi3DReWeighting.h.

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

boost::shared_ptr<TFile> edm::Lumi3DReWeighting::dataFile_

protected

Definition at line 63 of file Lumi3DReWeighting.h.

Referenced by Lumi3DReWeighting(), and weight3D_set().

std::string edm::Lumi3DReWeighting::dataFileName_

protected

Definition at line 58 of file Lumi3DReWeighting.h.

Referenced by Lumi3DReWeighting(), and weight3D_set().

std::string edm::Lumi3DReWeighting::DataHistName_

protected

Definition at line 60 of file Lumi3DReWeighting.h.

Referenced by Lumi3DReWeighting(), and weight3D_set().

boost::shared_ptr<TFile> edm::Lumi3DReWeighting::generatedFile_

protected

Definition at line 62 of file Lumi3DReWeighting.h.

Referenced by Lumi3DReWeighting(), and weight3D_set().

std::string edm::Lumi3DReWeighting::generatedFileName_

protected

Definition at line 57 of file Lumi3DReWeighting.h.

Referenced by Lumi3DReWeighting(), and weight3D_set().

std::string edm::Lumi3DReWeighting::GenHistName_

protected

Definition at line 59 of file Lumi3DReWeighting.h.

Referenced by Lumi3DReWeighting(), and weight3D_set().

boost::shared_ptr<TH1> edm::Lumi3DReWeighting::MC_distr_

protected

Definition at line 68 of file Lumi3DReWeighting.h.

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

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

protected

Definition at line 72 of file Lumi3DReWeighting.h.

Referenced by weight3D(), and weight3D_init().

std::string edm::Lumi3DReWeighting::weightFileName_

protected

Definition at line 61 of file Lumi3DReWeighting.h.

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

boost::shared_ptr<TH1> edm::Lumi3DReWeighting::weights_

protected

Definition at line 64 of file Lumi3DReWeighting.h.