From de285759e05828cc14fa943ea3bb000fa44bbb4b Mon Sep 17 00:00:00 2001 From: Colin Bernet Date: Thu, 12 Nov 2009 14:01:17 +0000 Subject: [PATCH] added comments for Florent --- DQMOffline/PFTau/interface/METBenchmark.h | 1 - .../PFTau/interface/MatchMETBenchmark.h | 21 +- DQMOffline/PFTau/src/MatchMETBenchmark.cc | 5 +- .../RecoParticleFlow/interface/TH2Analyzer.h | 14 + .../RecoParticleFlow/src/TH2Analyzer.cc | 287 +++++++++--------- 5 files changed, 181 insertions(+), 147 deletions(-) diff --git a/DQMOffline/PFTau/interface/METBenchmark.h b/DQMOffline/PFTau/interface/METBenchmark.h index 4a99b83d9b028..53cbfec516d48 100644 --- a/DQMOffline/PFTau/interface/METBenchmark.h +++ b/DQMOffline/PFTau/interface/METBenchmark.h @@ -26,7 +26,6 @@ class METBenchmark : public Benchmark { protected: - TH1F* pt_; TH1F* phi_; TH1F* sumEt_; diff --git a/DQMOffline/PFTau/interface/MatchMETBenchmark.h b/DQMOffline/PFTau/interface/MatchMETBenchmark.h index 310cf9498a14b..a31cba52667e4 100644 --- a/DQMOffline/PFTau/interface/MatchMETBenchmark.h +++ b/DQMOffline/PFTau/interface/MatchMETBenchmark.h @@ -8,6 +8,15 @@ #include "DataFormats/METReco/interface/METFwd.h" +// is this include necessary? +// check all includes + +// integrate and check your benchmarks in PFRootEvent (take PFCandidateManager as an example) + +// integrate and check your benchmarks Validation/RecoParticleFlow (take PFCandidateManager as an example) + +// remove the old benchmarks from these 2 packages (python files, C++ code, ...) + #include class MatchMETBenchmark : public Benchmark { @@ -26,15 +35,23 @@ class MatchMETBenchmark : public Benchmark { protected: - - TH2F* delta_et_Over_et_VS_et_; + // next 3: add to MatchCandidateBenchmark? + + // (rec - true) / true = rec/true - 1 TH2F* delta_et_VS_et_; + TH2F* delta_et_Over_et_VS_et_; + TH2F* delta_phi_VS_et_; + TH1F* delta_ex_; + + // True and Rec: remove. remove the following histo? TH2F* RecEt_VS_TrueEt_; TH2F* delta_set_VS_set_; TH2F* delta_set_Over_set_VS_set_; TH2F* delta_ex_VS_set_; + + // remove the following histo? TH2F* RecSet_Over_TrueSet_VS_TrueSet_; }; diff --git a/DQMOffline/PFTau/src/MatchMETBenchmark.cc b/DQMOffline/PFTau/src/MatchMETBenchmark.cc index 3038a37cbb2e3..ded57ebddcd25 100644 --- a/DQMOffline/PFTau/src/MatchMETBenchmark.cc +++ b/DQMOffline/PFTau/src/MatchMETBenchmark.cc @@ -31,7 +31,7 @@ void MatchMETBenchmark::setup() { dphiPS = PhaseSpace( 100, -3.2, 3.2); dptPS = PhaseSpace( 100, -500, 500); setPS = PhaseSpace( 300, 0.0, 3000); - dsetPS = PhaseSpace( 300, 0.-1000, 1000); + dsetPS = PhaseSpace( 200, 0.-1000, 1000); setOvsetPS = PhaseSpace( 500,0., 2.); break; case DQMOFFLINE: @@ -46,6 +46,8 @@ void MatchMETBenchmark::setup() { break; } + // variable bins to be done here, as they will save a lot of memory. + //float ptBins[11] = {0, 1, 2, 5, 10, 20, 50, 100, 200, 400, 1000}; delta_et_Over_et_VS_et_ = book2D("delta_et_Over_et_VS_et_", @@ -98,6 +100,7 @@ void MatchMETBenchmark::fillOne(const reco::MET& cand, if( !isInRange(cand.pt(), cand.eta(), cand.phi() ) ) return; + // in case it happens, print an error message LogWarning if ( matchedCand.pt()>0.001 ) delta_et_Over_et_VS_et_->Fill( matchedCand.pt(), (cand.pt() - matchedCand.pt())/matchedCand.pt() ); delta_et_VS_et_->Fill( matchedCand.pt(), cand.pt() - matchedCand.pt() ); delta_phi_VS_et_->Fill( matchedCand.pt(), cand.phi() - matchedCand.phi() ); diff --git a/Validation/RecoParticleFlow/interface/TH2Analyzer.h b/Validation/RecoParticleFlow/interface/TH2Analyzer.h index 682e3093269d2..ef2f4bc313730 100644 --- a/Validation/RecoParticleFlow/interface/TH2Analyzer.h +++ b/Validation/RecoParticleFlow/interface/TH2Analyzer.h @@ -9,6 +9,14 @@ class TH2; class TH1D; class TH2D; +// EN FAIT NON ? +// check why white window in colin's case +// are you making copies of histograms without changing the name? +// names could be handled in the following way: +// - give a name to each instance of TH2Analyzer in constructor +// - all histograms of the corresponding TH2Analyzer are created with a name (key) +// which starts with name_RMS + class TH2Analyzer : public TObject { public: @@ -48,9 +56,15 @@ class TH2Analyzer : public TObject { TH1D* SigmaGauss() { return sigmaGauss_; } TH1D* MeanX() { return meanXslice_; } + // add an histo for chi2 / ndof + // add a function FitSlice(int i) + // not now: work along Y + private: void ProcessSlices( const TH2D* histo ); + + // no need for const, because i is copied void ProcessSlice(const int i, TH1D* histo ) const; const TH2* hist2D_; diff --git a/Validation/RecoParticleFlow/src/TH2Analyzer.cc b/Validation/RecoParticleFlow/src/TH2Analyzer.cc index 765ee767f0e4b..09cc268da22fa 100644 --- a/Validation/RecoParticleFlow/src/TH2Analyzer.cc +++ b/Validation/RecoParticleFlow/src/TH2Analyzer.cc @@ -12,6 +12,7 @@ using namespace std; +// remove? void TH2Analyzer::Eval( int rebinFactor ) { //std::cout << "Eval!" << std::endl; @@ -71,145 +72,145 @@ void TH2Analyzer::Eval(const int rebinFactor, const int binxmin, const string rebinName = bname + "_rebin"; if (binxmax>hist2D_->GetNbinsX()) - { - std::cout << "Error: TH2Analyzer.cc: binxmax>hist2D_->GetNbinsX()" << std::endl; - return; - } + { + std::cout << "Error: TH2Analyzer.cc: binxmax>hist2D_->GetNbinsX()" << std::endl; + return; + } if (cst_binning) - { - //std::cout << "hist2D_->GetXaxis()->GetBinLowEdge(" << binxmin << ") = " << hist2D_->GetXaxis()->GetBinLowEdge(binxmin) << std::endl; - //std::cout << "hist2D_->GetXaxis()->GetBinUpEdge(" << binxmax << ") = " << hist2D_->GetXaxis()->GetBinUpEdge(binxmax) << std::endl; - //std::cout << "hist2D_->GetNbinsY() = " << hist2D_->GetNbinsY() << std::endl; - //std::cout << "hist2D_->GetYaxis()->GetXmin() = " << hist2D_->GetYaxis()->GetXmin() << std::endl; - //std::cout << "hist2D_->GetYaxis()->GetXmax() = " << hist2D_->GetYaxis()->GetXmax() << std::endl; - rebinnedHist2D_ = new TH2D(rebinName.c_str(),"rebinned histo", - binxmax-binxmin+1, - hist2D_->GetXaxis()->GetBinLowEdge(binxmin), - hist2D_->GetXaxis()->GetBinUpEdge(binxmax), - hist2D_->GetNbinsY(), - hist2D_->GetYaxis()->GetXmin(), - hist2D_->GetYaxis()->GetXmax() ); - for (int binyc=1;binycGetNbinsY()+1;++binyc) { - for (int binxc=binxmin;binxcGetBinContent(" << binxc << "," << binyc << ") = " << hist2D_->GetBinContent(binxc,binyc) << std::endl; - //std::cout << "hist2D_->GetBinError(" << binxc << "," << binyc << ") = " << hist2D_->GetBinError(binxc,binyc) << std::endl; - //std::cout << "binxc-binxmin+1 = " << binxc-binxmin+1 << std::endl; - rebinnedHist2D_->SetBinContent(binxc-binxmin+1,binyc,hist2D_->GetBinContent(binxc,binyc)); - rebinnedHist2D_->SetBinError(binxc-binxmin+1,binyc,hist2D_->GetBinError(binxc,binyc)); - } - } - rebinnedHist2D_->RebinX( rebinFactor ); + //std::cout << "hist2D_->GetXaxis()->GetBinLowEdge(" << binxmin << ") = " << hist2D_->GetXaxis()->GetBinLowEdge(binxmin) << std::endl; + //std::cout << "hist2D_->GetXaxis()->GetBinUpEdge(" << binxmax << ") = " << hist2D_->GetXaxis()->GetBinUpEdge(binxmax) << std::endl; + //std::cout << "hist2D_->GetNbinsY() = " << hist2D_->GetNbinsY() << std::endl; + //std::cout << "hist2D_->GetYaxis()->GetXmin() = " << hist2D_->GetYaxis()->GetXmin() << std::endl; + //std::cout << "hist2D_->GetYaxis()->GetXmax() = " << hist2D_->GetYaxis()->GetXmax() << std::endl; + rebinnedHist2D_ = new TH2D(rebinName.c_str(),"rebinned histo", + binxmax-binxmin+1, + hist2D_->GetXaxis()->GetBinLowEdge(binxmin), + hist2D_->GetXaxis()->GetBinUpEdge(binxmax), + hist2D_->GetNbinsY(), + hist2D_->GetYaxis()->GetXmin(), + hist2D_->GetYaxis()->GetXmax() ); + for (int binyc=1;binycGetNbinsY()+1;++binyc) + { + for (int binxc=binxmin;binxcGetBinContent(" << binxc << "," << binyc << ") = " << hist2D_->GetBinContent(binxc,binyc) << std::endl; + //std::cout << "hist2D_->GetBinError(" << binxc << "," << binyc << ") = " << hist2D_->GetBinError(binxc,binyc) << std::endl; + //std::cout << "binxc-binxmin+1 = " << binxc-binxmin+1 << std::endl; + rebinnedHist2D_->SetBinContent(binxc-binxmin+1,binyc,hist2D_->GetBinContent(binxc,binyc)); + rebinnedHist2D_->SetBinError(binxc-binxmin+1,binyc,hist2D_->GetBinError(binxc,binyc)); + } + } + rebinnedHist2D_->RebinX( rebinFactor ); - const string averageName = bname + "_average"; - average_ = new TH1D( averageName.c_str(),"arithmetic average", - rebinnedHist2D_->GetNbinsX(), - rebinnedHist2D_->GetXaxis()->GetXmin(), - rebinnedHist2D_->GetXaxis()->GetXmax() ); + const string averageName = bname + "_average"; + average_ = new TH1D( averageName.c_str(),"arithmetic average", + rebinnedHist2D_->GetNbinsX(), + rebinnedHist2D_->GetXaxis()->GetXmin(), + rebinnedHist2D_->GetXaxis()->GetXmax() ); - const string rmsName = bname + "_RMS"; - RMS_ = new TH1D( rmsName.c_str(), "RMS", - rebinnedHist2D_->GetNbinsX(), - rebinnedHist2D_->GetXaxis()->GetXmin(), - rebinnedHist2D_->GetXaxis()->GetXmax() ); + const string rmsName = bname + "_RMS"; + RMS_ = new TH1D( rmsName.c_str(), "RMS", + rebinnedHist2D_->GetNbinsX(), + rebinnedHist2D_->GetXaxis()->GetXmin(), + rebinnedHist2D_->GetXaxis()->GetXmax() ); - const string sigmaGaussName = bname + "_sigmaGauss"; - sigmaGauss_ = new TH1D(sigmaGaussName.c_str(), "sigmaGauss", - rebinnedHist2D_->GetNbinsX(), - rebinnedHist2D_->GetXaxis()->GetXmin(), - rebinnedHist2D_->GetXaxis()->GetXmax() ); + const string sigmaGaussName = bname + "_sigmaGauss"; + sigmaGauss_ = new TH1D(sigmaGaussName.c_str(), "sigmaGauss", + rebinnedHist2D_->GetNbinsX(), + rebinnedHist2D_->GetXaxis()->GetXmin(), + rebinnedHist2D_->GetXaxis()->GetXmax() ); - const string meanXName = bname + "_meanX"; - meanXslice_ = new TH1D(meanXName.c_str(), "meanX", - rebinnedHist2D_->GetNbinsX(), - rebinnedHist2D_->GetXaxis()->GetXmin(), - rebinnedHist2D_->GetXaxis()->GetXmax() ); - } + const string meanXName = bname + "_meanX"; + meanXslice_ = new TH1D(meanXName.c_str(), "meanX", + rebinnedHist2D_->GetNbinsX(), + rebinnedHist2D_->GetXaxis()->GetXmin(), + rebinnedHist2D_->GetXaxis()->GetXmax() ); + } else - { - // binning is not constant, but made to obtain almost the same number of events in each bin + { + // binning is not constant, but made to obtain almost the same number of events in each bin - //std::cout << "binxmax-binxmin+1 = " << binxmax-binxmin+1 << std::endl; - //std::cout << "rebinFactor = " << rebinFactor << std::endl; - //std::cout << "(binxmax-binxmin+1)/rebinFactor = " << (binxmax-binxmin+1)/rebinFactor << std::endl; - //std::cout << "((binxmax-binxmin+1)%rebinFactor) = " << ((binxmax-binxmin+1)%rebinFactor) << std::endl; - //std::cout << "abs((binxmax-binxmin+1)/rebinFactor) = " << abs((binxmax-binxmin+1)/rebinFactor) << std::endl; + //std::cout << "binxmax-binxmin+1 = " << binxmax-binxmin+1 << std::endl; + //std::cout << "rebinFactor = " << rebinFactor << std::endl; + //std::cout << "(binxmax-binxmin+1)/rebinFactor = " << (binxmax-binxmin+1)/rebinFactor << std::endl; + //std::cout << "((binxmax-binxmin+1)%rebinFactor) = " << ((binxmax-binxmin+1)%rebinFactor) << std::endl; + //std::cout << "abs((binxmax-binxmin+1)/rebinFactor) = " << abs((binxmax-binxmin+1)/rebinFactor) << std::endl; - unsigned int nbin=0; - if (((binxmax-binxmin+1)%rebinFactor)!=0.0) - { - nbin=abs((binxmax-binxmin+1)/rebinFactor)+1; - } - else nbin=(binxmax-binxmin+1)/rebinFactor; + unsigned int nbin=0; + if (((binxmax-binxmin+1)%rebinFactor)!=0.0) + { + nbin=abs((binxmax-binxmin+1)/rebinFactor)+1; + } + else nbin=(binxmax-binxmin+1)/rebinFactor; - double *xlow = new double[nbin+1]; - int *binlow = new int[nbin+1]; + double *xlow = new double[nbin+1]; + int *binlow = new int[nbin+1]; - TH1D* h0_slice1 = hist2D_->ProjectionY("h0_slice1", binxmin, binxmax, ""); - const unsigned int totalNumberOfEvents=static_cast(h0_slice1->GetEntries()); - //std::cout << "totalNumberOfEvents = " << totalNumberOfEvents << std::endl; - delete h0_slice1; + TH1D* h0_slice1 = hist2D_->ProjectionY("h0_slice1", binxmin, binxmax, ""); + const unsigned int totalNumberOfEvents=static_cast(h0_slice1->GetEntries()); + //std::cout << "totalNumberOfEvents = " << totalNumberOfEvents << std::endl; + delete h0_slice1; - unsigned int neventsc=0; - unsigned int binXmaxc=binxmin+1; - xlow[0]=hist2D_->GetXaxis()->GetBinLowEdge(binxmin); - binlow[0]=binxmin; - for (unsigned int binc=1;bincProjectionY("h0_slice1",binxmin, binXmaxc, ""); - neventsc=static_cast(h0_slice1c->GetEntries()); -// //std::cout << "FL : neventsc = " << neventsc << std::endl; -// //std::cout << "FL : binXmaxc = " << binXmaxc << std::endl; - ++binXmaxc; - delete h0_slice1c; - } - //std::cout << "binXmaxc-1 = " << binXmaxc-1 << std::endl; - binlow[binc]=binXmaxc-1; - xlow[binc]=hist2D_->GetXaxis()->GetBinLowEdge(binXmaxc-1); - } - xlow[nbin]=hist2D_->GetXaxis()->GetBinUpEdge(binxmax); - binlow[nbin]=binxmax; + unsigned int neventsc=0; + unsigned int binXmaxc=binxmin+1; + xlow[0]=hist2D_->GetXaxis()->GetBinLowEdge(binxmin); + binlow[0]=binxmin; + for (unsigned int binc=1;bincProjectionY("h0_slice1",binxmin, binXmaxc, ""); + neventsc=static_cast(h0_slice1c->GetEntries()); + // //std::cout << "FL : neventsc = " << neventsc << std::endl; + // //std::cout << "FL : binXmaxc = " << binXmaxc << std::endl; + ++binXmaxc; + delete h0_slice1c; + } + //std::cout << "binXmaxc-1 = " << binXmaxc-1 << std::endl; + binlow[binc]=binXmaxc-1; + xlow[binc]=hist2D_->GetXaxis()->GetBinLowEdge(binXmaxc-1); + } + xlow[nbin]=hist2D_->GetXaxis()->GetBinUpEdge(binxmax); + binlow[nbin]=binxmax; - //for (unsigned int binc=0;bincGetNbinsY(), - hist2D_->GetYaxis()->GetXmin(), - hist2D_->GetYaxis()->GetXmax() ); - for (int binyc=1;binycGetNbinsY()+1;++binyc) - { - for (unsigned int binxc=1;binxcGetNbinsY(), + hist2D_->GetYaxis()->GetXmin(), + hist2D_->GetYaxis()->GetXmax() ); + for (int binyc=1;binycGetNbinsY()+1;++binyc) { - sum+=hist2D_->GetBinContent(binxh2c,binyc); + for (unsigned int binxc=1;binxcGetBinContent(binxh2c,binyc); + } + rebinnedHist2D_->SetBinContent(binxc,binyc,sum); + } } - rebinnedHist2D_->SetBinContent(binxc,binyc,sum); - } - } - const string averageName = bname + "_average"; - average_ = new TH1D( averageName.c_str(),"arithmetic average", nbin, xlow); + const string averageName = bname + "_average"; + average_ = new TH1D( averageName.c_str(),"arithmetic average", nbin, xlow); - const string rmsName = bname + "_RMS"; - RMS_ = new TH1D( rmsName.c_str(), "RMS", nbin, xlow); + const string rmsName = bname + "_RMS"; + RMS_ = new TH1D( rmsName.c_str(), "RMS", nbin, xlow); - const string sigmaGaussName = bname + "_sigmaGauss"; - sigmaGauss_ = new TH1D(sigmaGaussName.c_str(), "sigmaGauss", nbin, xlow); + const string sigmaGaussName = bname + "_sigmaGauss"; + sigmaGauss_ = new TH1D(sigmaGaussName.c_str(), "sigmaGauss", nbin, xlow); - const string meanXName = bname + "_meanX"; - meanXslice_ = new TH1D(meanXName.c_str(), "meanX", nbin, xlow); - delete [] xlow; - delete [] binlow; - } + const string meanXName = bname + "_meanX"; + meanXslice_ = new TH1D(meanXName.c_str(), "meanX", nbin, xlow); + delete [] xlow; + delete [] binlow; + } ProcessSlices( rebinnedHist2D_ ); } @@ -254,34 +255,34 @@ void TH2Analyzer::ProcessSlice(const int i, TH1D* proj ) const { //std::cout << "FL: rms = " << rms << std::endl; if (rms!=0.0) - { - const double fitmin=proj->GetMean()-proj->GetRMS(); - const double fitmax=proj->GetMean()+proj->GetRMS(); + { + const double fitmin=proj->GetMean()-proj->GetRMS(); + const double fitmax=proj->GetMean()+proj->GetRMS(); - //std::cout << "i = " << i << std::endl; - //std::cout << "fitmin = " << fitmin << std::endl; - //std::cout << "fitmax = " << fitmax << std::endl; + //std::cout << "i = " << i << std::endl; + //std::cout << "fitmin = " << fitmin << std::endl; + //std::cout << "fitmax = " << fitmax << std::endl; - //proj->Draw(); - TF1* f1= new TF1("f1", "gaus", fitmin, fitmax); - f1->SetParameters(proj->GetRMS(),proj->GetMean(),proj->GetBinContent(proj->GetMaximumBin())); - proj->Fit("f1","R", "", proj->GetXaxis()->GetXmin(), proj->GetXaxis()->GetXmax()); + //proj->Draw(); + TF1* f1= new TF1("f1", "gaus", fitmin, fitmax); + f1->SetParameters(proj->GetRMS(),proj->GetMean(),proj->GetBinContent(proj->GetMaximumBin())); + proj->Fit("f1","R", "", proj->GetXaxis()->GetXmin(), proj->GetXaxis()->GetXmax()); - //std::ostringstream oss; - //oss << i; - //const std::string plotfitname="Plots/fitbin_"+oss.str()+".eps"; - //gPad->SaveAs( plotfitname.c_str() ); - //std::cout << "param1 = " << f1->GetParameter(1) << std::endl; - //std::cout << "param2 = " << f1->GetParameter(2) << std::endl; - //std::cout << "paramError2 = " << f1->GetParError(2) << std::endl; + //std::ostringstream oss; + //oss << i; + //const std::string plotfitname="Plots/fitbin_"+oss.str()+".eps"; + //gPad->SaveAs( plotfitname.c_str() ); + //std::cout << "param1 = " << f1->GetParameter(1) << std::endl; + //std::cout << "param2 = " << f1->GetParameter(2) << std::endl; + //std::cout << "paramError2 = " << f1->GetParError(2) << std::endl; - sigmaGauss_->SetBinContent(i, f1->GetParameter(2)); - sigmaGauss_->SetBinError(i, f1->GetParError(2)); - delete f1; - } + sigmaGauss_->SetBinContent(i, f1->GetParameter(2)); + sigmaGauss_->SetBinError(i, f1->GetParError(2)); + delete f1; + } else - { - sigmaGauss_->SetBinContent(i, 0.0); - sigmaGauss_->SetBinError(i, 0.0); - } + { + sigmaGauss_->SetBinContent(i, 0.0); + sigmaGauss_->SetBinError(i, 0.0); + } }