#include <MuonAnalysis/MomentumScaleCalibration/plugins/ErrorsAnalyzer.cc>

Inheritance diagram for ErrorsAnalyzer:

Public Member Functions
	ErrorsAnalyzer (const edm::ParameterSet &)

	~ErrorsAnalyzer ()

Public Member Functions inherited from edm::EDAnalyzer
	EDAnalyzer ()

std::string	workerType () const

virtual	~EDAnalyzer ()

Private Member Functions
virtual void	analyze (const edm::Event &, const edm::EventSetup &)

void	drawHistograms (const TProfile histo, const TProfile histoPlusErr, const TProfile *histoMinusErr, const TString &type)

virtual void	endJob ()

void	fillHistograms ()

void	fillValueError ()

Private Attributes
bool	debug_

std::vector< int >	errorFactors_

std::vector< double >	errors_

int	etaBins_

double	etaMax_

double	etaMin_

uint32_t	maxEvents_

TString	outputFileName_

std::vector< double >	parameters_

int	ptBins_

double	ptMax_

double	ptMin_

int	resolFitType_

TProfile *	sigmaMassVsEta_

TProfile *	sigmaMassVsEtaMinusErr_

TProfile *	sigmaMassVsEtaPlusErr_

TProfile *	sigmaMassVsPt_

TProfile *	sigmaMassVsPtMinusErr_

TProfile *	sigmaMassVsPtPlusErr_

TProfile *	sigmaPtVsEta_

TProfile *	sigmaPtVsEtaMinusErr_

TProfile *	sigmaPtVsEtaPlusErr_

TProfile *	sigmaPtVsPt_

TProfile *	sigmaPtVsPtMinusErr_

TProfile *	sigmaPtVsPtPlusErr_

TString	treeFileName_

std::vector< double >	valueMinusError_

std::vector< double >	valuePlusError_

Additional Inherited Members
Public Types inherited from edm::EDAnalyzer
typedef EDAnalyzer	ModuleType

typedef WorkerT< EDAnalyzer >	WorkerType

Static Public Member Functions inherited from edm::EDAnalyzer
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &)

Protected Member Functions inherited from edm::EDAnalyzer
CurrentProcessingContext const *	currentContext () const

Detailed Description

Description: <one line="" class="" summary>="">

Implementation: <Notes on="" implementation>="">

Definition at line 49 of file ErrorsAnalyzer.h.

Constructor & Destructor Documentation

ErrorsAnalyzer::ErrorsAnalyzer ( const edm::ParameterSet & iConfig )

explicit

Definition at line 6 of file ErrorsAnalyzer.cc.

References gather_cfg::cout, errorFactors_, errors_, etaBins_, etaMax_, etaMin_, cmsRelvalreport::exit, fillValueError(), edm::ParameterSet::getParameter(), parameters_, ptBins_, ptMax_, ptMin_, sigmaMassVsEta_, sigmaMassVsEtaMinusErr_, sigmaMassVsEtaPlusErr_, sigmaMassVsPt_, sigmaMassVsPtMinusErr_, sigmaMassVsPtPlusErr_, sigmaPtVsEta_, sigmaPtVsEtaMinusErr_, sigmaPtVsEtaPlusErr_, sigmaPtVsPt_, sigmaPtVsPtMinusErr_, and sigmaPtVsPtPlusErr_.

                                                              :
   treeFileName_( iConfig.getParameter<std::string>("InputFileName") ),
   resolFitType_( iConfig.getParameter<int>("ResolFitType") ),
   maxEvents_( iConfig.getParameter<int>("MaxEvents") ),
   outputFileName_( iConfig.getParameter<std::string>("OutputFileName") ),
   ptBins_( iConfig.getParameter<int>("PtBins") ),
   ptMin_( iConfig.getParameter<double>("PtMin") ),
   ptMax_( iConfig.getParameter<double>("PtMax") ),
   etaBins_( iConfig.getParameter<int>("EtaBins") ),
   etaMin_( iConfig.getParameter<double>("EtaMin") ),
   etaMax_( iConfig.getParameter<double>("EtaMax") ),
   debug_( iConfig.getParameter<bool>("Debug") )
 {
   parameters_ = iConfig.getParameter<std::vector<double> >("Parameters");
   errors_ = iConfig.getParameter<std::vector<double> >("Errors");
   errorFactors_ = iConfig.getParameter<std::vector<int> >("ErrorFactors");
 
   if( (parameters_.size() != errors_.size()) || (parameters_.size() != errorFactors_.size()) ) {
     std::cout << "Error: parameters and errors have different number of values" << std::endl;
     exit(1);
   }
 
   fillValueError();
 
   sigmaPtVsPt_ = new TProfile("sigmaPtVsPtProfile", "sigmaPtVsPt", ptBins_, ptMin_, ptMax_);
   sigmaPtVsPtPlusErr_ = new TProfile("sigmaPtVsPtPlusErrProfile", "sigmaPtVsPtPlusErr", ptBins_, ptMin_, ptMax_);
   sigmaPtVsPtMinusErr_ = new TProfile("sigmaPtVsPtMinusErrProfile", "sigmaPtVsPtMinusErr", ptBins_, ptMin_, ptMax_);
 
   sigmaPtVsEta_ = new TProfile("sigmaPtVsEtaProfile", "sigmaPtVsEta", etaBins_, etaMin_, etaMax_);
   sigmaPtVsEtaPlusErr_ = new TProfile("sigmaPtVsEtaPlusErrProfile", "sigmaPtVsEtaPlusErr", etaBins_, etaMin_, etaMax_);
   sigmaPtVsEtaMinusErr_ = new TProfile("sigmaPtVsEtaMinusErrProfile", "sigmaPtVsEtaMinusErr", etaBins_, etaMin_, etaMax_);
 
   sigmaMassVsPt_ = new TProfile("sigmaMassVsPtProfile", "sigmaMassVsPt", ptBins_, ptMin_, ptMax_);
   sigmaMassVsPtPlusErr_ = new TProfile("sigmaMassVsPtPlusErrProfile", "sigmaMassVsPtPlusErr", ptBins_, ptMin_, ptMax_);
   sigmaMassVsPtMinusErr_ = new TProfile("sigmaMassVsPtMinusErrProfile", "sigmaMassVsPtMinusErr", ptBins_, ptMin_, ptMax_);
 
   sigmaMassVsEta_ = new TProfile("sigmaMassVsEtaProfile", "sigmaMassVsEta", etaBins_, etaMin_, etaMax_);
   sigmaMassVsEtaPlusErr_ = new TProfile("sigmaMassVsEtaPlusErrProfile", "sigmaMassVsEtaPlusErr", etaBins_, etaMin_, etaMax_);
   sigmaMassVsEtaMinusErr_ = new TProfile("sigmaMassVsEtaMinusErrProfile", "sigmaMassVsEtaMinusErr", etaBins_, etaMin_, etaMax_);
 }

ErrorsAnalyzer::~ErrorsAnalyzer ( )

Definition at line 62 of file ErrorsAnalyzer.cc.

References drawHistograms(), fillHistograms(), download_sqlite_cfg::outputFile, outputFileName_, sigmaMassVsEta_, sigmaMassVsEtaMinusErr_, sigmaMassVsEtaPlusErr_, sigmaMassVsPt_, sigmaMassVsPtMinusErr_, sigmaMassVsPtPlusErr_, sigmaPtVsEta_, sigmaPtVsEtaMinusErr_, sigmaPtVsEtaPlusErr_, sigmaPtVsPt_, sigmaPtVsPtMinusErr_, and sigmaPtVsPtPlusErr_.

 {
   gROOT->SetStyle("Plain");
 
   fillHistograms();
 
   TFile * outputFile = new TFile(outputFileName_, "RECREATE");
 
   drawHistograms(sigmaPtVsEta_, sigmaPtVsEtaPlusErr_, sigmaPtVsEtaMinusErr_, "sigmaPtVsEta");
   drawHistograms(sigmaPtVsPt_, sigmaPtVsPtPlusErr_, sigmaPtVsPtMinusErr_, "sigmaPtVsPt");
 
   drawHistograms(sigmaMassVsEta_, sigmaMassVsEtaPlusErr_, sigmaMassVsEtaMinusErr_, "sigmaMassVsEta");
   drawHistograms(sigmaMassVsPt_, sigmaMassVsPtPlusErr_, sigmaMassVsPtMinusErr_, "sigmaMassVsPt");
 
   outputFile->Write();
   outputFile->Close();
 }

Member Function Documentation

void ErrorsAnalyzer::analyze	(	const edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)

privatevirtual

Implements edm::EDAnalyzer.

Definition at line 174 of file ErrorsAnalyzer.cc.

175 {

176 }

void ErrorsAnalyzer::drawHistograms	(	const TProfile *	histo,
		const TProfile *	histoPlusErr,
		const TProfile *	histoMinusErr,
		const TString &	type
	)

private

Definition at line 80 of file ErrorsAnalyzer.cc.

References svgfig::canvas(), i, sigmaMassVsEta_, sigmaMassVsEtaMinusErr_, sigmaMassVsEtaPlusErr_, sigmaPtVsEta_, sigmaPtVsEtaMinusErr_, sigmaPtVsEtaPlusErr_, indexGen::title, and makeHLTPrescaleTable::values.

Referenced by ~ErrorsAnalyzer().

 {
   TH1D * sigmaPtVsEtaTH1D = new TH1D(type, type, histo->GetNbinsX(),
                      histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
   TH1D * sigmaPtVsEtaPlusErrTH1D = new TH1D(type+"PlusErr", type+"PlusErr", histo->GetNbinsX(),
                         histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
   TH1D * sigmaPtVsEtaMinusErrTH1D = new TH1D(type+"MinusErr", type+"MinusErr", histo->GetNbinsX(),
                          histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
 
   TH1D * sigmaMassVsEtaTH1D = new TH1D(type, type, histo->GetNbinsX(),
                        histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
   TH1D * sigmaMassVsEtaPlusErrTH1D = new TH1D(type+"PlusErr", type+"PlusErr", histo->GetNbinsX(),
                           histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
   TH1D * sigmaMassVsEtaMinusErrTH1D = new TH1D(type+"MinusErr", type+"MinusErr", histo->GetNbinsX(),
                            histo->GetXaxis()->GetXmin(), histo->GetXaxis()->GetXmax());
 
   TCanvas * canvas = new TCanvas("canvas"+type, "canvas"+type, 1000, 800);
   for( int iBin = 1; iBin <= histo->GetNbinsX(); ++iBin ) {
     sigmaPtVsEtaTH1D->SetBinContent(iBin, histo->GetBinContent(iBin));
     sigmaPtVsEtaPlusErrTH1D->SetBinContent(iBin, histoPlusErr->GetBinContent(iBin));
     sigmaPtVsEtaMinusErrTH1D->SetBinContent(iBin, histoMinusErr->GetBinContent(iBin));
   }
   int numBins = sigmaPtVsEtaTH1D->GetNbinsX();
   // Draw TGraph with asymmetric errors
   Double_t * values = sigmaPtVsEtaTH1D->GetArray();
   Double_t * valuesPlus = sigmaPtVsEtaPlusErrTH1D->GetArray();
   Double_t * valuesMinus = sigmaPtVsEtaMinusErrTH1D->GetArray();
   double * posErrors = new double[numBins];
   double * negErrors = new double[numBins];
 
   TGraphAsymmErrors * graphAsymmErrors = new TGraphAsymmErrors(sigmaPtVsEtaTH1D);
   TGraph * graph = new TGraph(sigmaPtVsEtaTH1D);
 
   for( int i=1; i<=numBins; ++i ) {
     // std::cout << "filling " << i << std::endl;
     posErrors[i-1] = valuesPlus[i] - values[i];
     if( valuesMinus[i] < 0 ) negErrors[i-1] = values[i];
     else negErrors[i-1] = values[i] - valuesMinus[i];
 
     graphAsymmErrors->SetPointEYlow(i-1, negErrors[i-1]);
     graphAsymmErrors->SetPointEYhigh(i-1, posErrors[i-1]);
   }
 
   canvas->Draw();
 
   graphAsymmErrors->SetTitle("");
   graphAsymmErrors->SetFillColor(kGray);
   graphAsymmErrors->Draw("A2");
   TString title(type);
   if( type == "Eta" ) title = "#eta";
   graphAsymmErrors->GetXaxis()->SetTitle(title);
   graphAsymmErrors->GetYaxis()->SetTitle("#sigmaPt/Pt");
   graph->Draw();
 
   //  graph->SetLineColor(kGray);
   //  graph->SetMarkerColor(kBlack);
   //  graph->Draw("AP");
 
 //   if( debug_ ) {
 //     sigmaPtVsEtaPlusErrTH1D->SetLineColor(kRed);
 //     sigmaPtVsEtaMinusErrTH1D->SetLineColor(kBlue);
 //   }
 //   else {
 //     sigmaPtVsEtaPlusErrTH1D->SetFillColor(kGray);
 //     sigmaPtVsEtaPlusErrTH1D->SetLineColor(kWhite);
 //     sigmaPtVsEtaMinusErrTH1D->SetFillColor(kWhite);
 //     sigmaPtVsEtaMinusErrTH1D->SetLineColor(kWhite);
 //   }
 //   sigmaPtVsEtaPlusErrTH1D->Draw();
 //   sigmaPtVsEtaTH1D->Draw("SAME");
 //   sigmaPtVsEtaMinusErrTH1D->Draw("SAME");
 
   sigmaPtVsEtaPlusErrTH1D->Write();
   sigmaPtVsEtaTH1D->Write();
   sigmaPtVsEtaMinusErrTH1D->Write();
 
   sigmaPtVsEtaPlusErr_->Write();
   sigmaPtVsEta_->Write();
   sigmaPtVsEtaMinusErr_->Write();
 
   // Mass
   sigmaMassVsEtaPlusErrTH1D->Write();
   sigmaMassVsEtaTH1D->Write();
   sigmaMassVsEtaMinusErrTH1D->Write();
 
   sigmaMassVsEtaPlusErr_->Write();
   sigmaMassVsEta_->Write();
   sigmaMassVsEtaMinusErr_->Write();
 
   canvas->Write();
 }

virtual void ErrorsAnalyzer::endJob ( void )

inlineprivatevirtual

Reimplemented from edm::EDAnalyzer.

Definition at line 60 of file ErrorsAnalyzer.h.

60 {};

void ErrorsAnalyzer::fillHistograms ( void )

private

Definition at line 178 of file ErrorsAnalyzer.cc.

References gather_cfg::cout, MuScleFitUtils::debug, debug_, MuScleFitUtils::debugMassResol_, i, MuScleFitUtils::massResolution(), maxEvents_, parameters_, RootTreeHandler::readTree(), MuScleFitUtils::resfind, resolFitType_, MuScleFitUtils::resolutionFunction, resolutionFunctionService(), resolutionFunctionVecService(), sigmaMassVsEta_, sigmaMassVsEtaMinusErr_, sigmaMassVsEtaPlusErr_, sigmaMassVsPt_, sigmaMassVsPtMinusErr_, sigmaMassVsPtPlusErr_, resolutionFunctionBase< T >::sigmaPt(), sigmaPtVsEta_, sigmaPtVsEtaMinusErr_, sigmaPtVsEtaPlusErr_, sigmaPtVsPt_, sigmaPtVsPtMinusErr_, sigmaPtVsPtPlusErr_, treeFileName_, valueMinusError_, and valuePlusError_.

Referenced by ~ErrorsAnalyzer().

 {
   std::cout << "Reading muon pairs from Root Tree in " << treeFileName_ << std::endl;
   RootTreeHandler rootTreeHandler;
 
   typedef std::vector<std::pair<lorentzVector,lorentzVector> > MuonPairVector;
   MuonPairVector savedPair;
   std::vector<std::pair<int, int> > evtRun;
   rootTreeHandler.readTree(maxEvents_, treeFileName_, &savedPair, 0, &evtRun);
   // rootTreeHandler.readTree(maxEvents, inputRootTreeFileName_, &savedPair, &(MuScleFitUtils::genPair));
 
   resolutionFunctionBase<std::vector<double> > * resolutionFunctionForVec = resolutionFunctionVecService( resolFitType_ );
   MuScleFitUtils::resolutionFunction = resolutionFunctionService( resolFitType_ );
   MuScleFitUtils::debugMassResol_ = false;
   MuScleFitUtils::debug = 0;
   MuScleFitUtils::resfind = std::vector<int>(6, 0);
 
   // Loop on all the pairs
   unsigned int i = 0;
   MuonPairVector::iterator it = savedPair.begin();
   std::cout << "Starting loop on " << savedPair.size() << " muons" << std::endl;
   for( ; it != savedPair.end(); ++it, ++i ) {
     double pt1 = it->first.pt();
     double eta1 = it->first.eta();
     double pt2 = it->second.pt();
     double eta2 = it->second.eta();
 
     if( debug_ ) {
       std::cout << "pt1 = " << pt1 << ", eta1 = " << eta1 << ", pt2 = " << pt2 << ", eta2 = " << eta2 << std::endl;
     }
 
     if( pt1 == 0 && pt2 == 0 && eta1 == 0 && eta2 == 0 ) continue;
 
     double sigmaPt1 = resolutionFunctionForVec->sigmaPt( pt1,eta1,parameters_ );
     double sigmaPt2 = resolutionFunctionForVec->sigmaPt( pt2,eta2,parameters_ );
     double sigmaPtPlusErr1 = resolutionFunctionForVec->sigmaPt( pt1,eta1,valuePlusError_ );
     double sigmaPtPlusErr2 = resolutionFunctionForVec->sigmaPt( pt2,eta2,valuePlusError_ );
     double sigmaPtMinusErr1 = resolutionFunctionForVec->sigmaPt( pt1,eta1,valueMinusError_ );
     double sigmaPtMinusErr2 = resolutionFunctionForVec->sigmaPt( pt2,eta2,valueMinusError_ );
 
     double sigmaMass = MuScleFitUtils::massResolution( it->first, it->second, parameters_ );
     double sigmaMassPlusErr = MuScleFitUtils::massResolution( it->first, it->second, valuePlusError_ );
     double sigmaMassMinusErr = MuScleFitUtils::massResolution( it->first, it->second, valueMinusError_ );
 
     if( debug_ ) {
       std::cout << "sigmaPt1 = " << sigmaPt1 << " + " << sigmaPtPlusErr1 << " - " << sigmaPtMinusErr1 << std::endl;
       std::cout << "sigmaPt2 = " << sigmaPt2 << " + " << sigmaPtPlusErr2 << " - " << sigmaPtMinusErr2 << std::endl;
       std::cout << "sigmaMass = " << sigmaMass << " + " << sigmaMassPlusErr << " - " << sigmaMassMinusErr << std::endl;
     }
 
     // Protections from nans
     if( pt1 != pt1 ) continue;
     if( pt2 != pt2 ) continue;
     if( eta1 != eta1 ) continue;
     if( eta2 != eta2 ) continue;
     if( sigmaPt1 != sigmaPt1 ) continue;
     if( sigmaPt2 != sigmaPt2 ) continue;
     if( sigmaPtPlusErr1 != sigmaPtPlusErr1 ) continue;
     if( sigmaPtPlusErr2 != sigmaPtPlusErr2 ) continue;
     if( sigmaPtMinusErr1 != sigmaPtMinusErr1 ) continue;
     if( sigmaPtMinusErr2 != sigmaPtMinusErr2 ) continue;
     if( sigmaMass != sigmaMass ) continue;
     if( sigmaMassPlusErr != sigmaMassPlusErr ) continue;
     if( sigmaMassMinusErr != sigmaMassMinusErr ) continue;
 
     std::cout << "Muon pair number " << i << std::endl;
 
     // std::cout << "sigmaPt1 = " << sigmaPt1 << ", sigmaPt2 = " << sigmaPt2 << std::endl;
     // std::cout << "sigmaPtPlusErr1 = " << sigmaPtPlusErr1 << ", sigmaPtMinusErr2 = " << sigmaPtMinusErr2 << std::endl;
     // std::cout << "sigmaPtMinusErr1 = " << sigmaPtPlusErr1 << ", sigmaPtMinusErr2 = " << sigmaPtMinusErr2 << std::endl;
 
     sigmaPtVsPt_->Fill(pt1, sigmaPt1);
     sigmaPtVsPt_->Fill(pt2, sigmaPt2);
     sigmaPtVsPtPlusErr_->Fill(pt1, sigmaPtPlusErr1);
     sigmaPtVsPtPlusErr_->Fill(pt2, sigmaPtPlusErr2);
     sigmaPtVsPtMinusErr_->Fill(pt1, sigmaPtMinusErr1);
     sigmaPtVsPtMinusErr_->Fill(pt2, sigmaPtMinusErr2);
 
     sigmaPtVsEta_->Fill(eta1, sigmaPt1);
     sigmaPtVsEta_->Fill(eta2, sigmaPt2);
     sigmaPtVsEtaPlusErr_->Fill(eta1, sigmaPtPlusErr1);
     sigmaPtVsEtaPlusErr_->Fill(eta2, sigmaPtPlusErr2);
     sigmaPtVsEtaMinusErr_->Fill(eta1, sigmaPtMinusErr1);
     sigmaPtVsEtaMinusErr_->Fill(eta2, sigmaPtMinusErr2);
 
 
     sigmaMassVsPt_->Fill(pt1, sigmaMass);
     sigmaMassVsPt_->Fill(pt2, sigmaMass);
     sigmaMassVsPtPlusErr_->Fill(pt1, sigmaMassPlusErr);
     sigmaMassVsPtPlusErr_->Fill(pt2, sigmaMassPlusErr);
     sigmaMassVsPtMinusErr_->Fill(pt1, sigmaMassMinusErr);
     sigmaMassVsPtMinusErr_->Fill(pt2, sigmaMassMinusErr);
 
     sigmaMassVsEta_->Fill(eta1, sigmaMass);
     sigmaMassVsEta_->Fill(eta2, sigmaMass);
     sigmaMassVsEtaPlusErr_->Fill(eta1, sigmaMassPlusErr);
     sigmaMassVsEtaPlusErr_->Fill(eta2, sigmaMassPlusErr);
     sigmaMassVsEtaMinusErr_->Fill(eta1, sigmaMassMinusErr);
     sigmaMassVsEtaMinusErr_->Fill(eta2, sigmaMassMinusErr);
   }
 }

void ErrorsAnalyzer::fillValueError ( )

private

Definition at line 47 of file ErrorsAnalyzer.cc.

References errorFactors_, errors_, i, parameters_, valueMinusError_, and valuePlusError_.

Referenced by ErrorsAnalyzer().

 {
   valuePlusError_.resize(parameters_.size());
   valueMinusError_.resize(parameters_.size());
 
   std::vector<double>::const_iterator parIt = parameters_.begin();
   std::vector<double>::const_iterator errIt = errors_.begin();
   std::vector<int>::const_iterator errFactorIt = errorFactors_.begin();
   int i=0;
   for( ; parIt != parameters_.end(); ++parIt, ++errIt, ++errFactorIt, ++i ) {
     valuePlusError_[i] = *parIt + (*errIt)*(*errFactorIt);
     valueMinusError_[i] = *parIt - (*errIt)*(*errFactorIt);
   }
 }