RecoParticleFlow/src/TH2Analyzer.cc

0001 #include "Validation/RecoParticleFlow/interface/TH2Analyzer.h"
0002
0003 #include "TF1.h"
0004 #include "TH1D.h"
0005 #include "TH2D.h"
0006 #include "TPad.h"
0007
0008 #include <algorithm>
0009 #include <iostream>
0010 #include <sstream>
0011 #include <string>
0012
0013 using namespace std;
0014
0015 // remove?
0016 void TH2Analyzer::Eval(int rebinFactor) {
0017   // std::cout << "Eval!" << std::endl;
0018
0019   Reset();
0020
0021   const string bname = hist2D_->GetName();
0022   const string rebinName = bname + "_rebin";
0023   rebinnedHist2D_ = (TH2D *)hist2D_->Clone(rebinName.c_str());
0024   rebinnedHist2D_->RebinX(rebinFactor);
0025
0026   const string averageName = bname + "_average";
0027   average_ = new TH1D(averageName.c_str(),
0028                       "arithmetic average",
0029                       rebinnedHist2D_->GetNbinsX(),
0030                       rebinnedHist2D_->GetXaxis()->GetXmin(),
0031                       rebinnedHist2D_->GetXaxis()->GetXmax());
0032
0033   const string rmsName = bname + "_RMS";
0034   RMS_ = new TH1D(rmsName.c_str(),
0035                   "RMS",
0036                   rebinnedHist2D_->GetNbinsX(),
0037                   rebinnedHist2D_->GetXaxis()->GetXmin(),
0038                   rebinnedHist2D_->GetXaxis()->GetXmax());
0039
0040   const string sigmaGaussName = bname + "_sigmaGauss";
0041   sigmaGauss_ = new TH1D(sigmaGaussName.c_str(),
0042                          "sigmaGauss",
0043                          rebinnedHist2D_->GetNbinsX(),
0044                          rebinnedHist2D_->GetXaxis()->GetXmin(),
0045                          rebinnedHist2D_->GetXaxis()->GetXmax());
0046
0047   const string meanXName = bname + "_meanX";
0048   meanXslice_ = new TH1D(meanXName.c_str(),
0049                          "meanX",
0050                          rebinnedHist2D_->GetNbinsX(),
0051                          rebinnedHist2D_->GetXaxis()->GetXmin(),
0052                          rebinnedHist2D_->GetXaxis()->GetXmax());
0053
0054   ProcessSlices(rebinnedHist2D_);
0055 }
0056
0057 void TH2Analyzer::Reset() {
0058   if (rebinnedHist2D_)
0059     delete rebinnedHist2D_;
0060   if (average_)
0061     delete average_;
0062   if (RMS_)
0063     delete RMS_;
0064   if (sigmaGauss_)
0065     delete sigmaGauss_;
0066   if (meanXslice_)
0067     delete meanXslice_;
0068
0069   // for( unsigned i=0; i<parameters_.size(); ++i) {
0070   //  delete parameters_[i];
0071   //}
0072
0073   // parameters_.clear();
0074 }
0075
0076 void TH2Analyzer::Eval(const int rebinFactor, const int binxmin, const int binxmax, const bool cst_binning) {
0077   Reset();
0078   const string bname = hist2D_->GetName();
0079   const string rebinName = bname + "_rebin";
0080
0081   if (binxmax > hist2D_->GetNbinsX()) {
0082     std::cout << "Error: TH2Analyzer.cc: binxmax>hist2D_->GetNbinsX()" << std::endl;
0083     return;
0084   }
0085
0086   if (cst_binning) {
0087     // std::cout << "hist2D_->GetXaxis()->GetBinLowEdge(" << binxmin << ") = "
0088     // << hist2D_->GetXaxis()->GetBinLowEdge(binxmin) << std::endl; std::cout <<
0089     // "hist2D_->GetXaxis()->GetBinUpEdge(" << binxmax << ") = " <<
0090     // hist2D_->GetXaxis()->GetBinUpEdge(binxmax) << std::endl; std::cout <<
0091     // "hist2D_->GetNbinsY() = " << hist2D_->GetNbinsY() << std::endl; std::cout
0092     // << "hist2D_->GetYaxis()->GetXmin() = " << hist2D_->GetYaxis()->GetXmin()
0093     // << std::endl; std::cout << "hist2D_->GetYaxis()->GetXmax() = " <<
0094     // hist2D_->GetYaxis()->GetXmax() << std::endl;
0095     rebinnedHist2D_ = new TH2D(rebinName.c_str(),
0096                                "rebinned histo",
0097                                binxmax - binxmin + 1,
0098                                hist2D_->GetXaxis()->GetBinLowEdge(binxmin),
0099                                hist2D_->GetXaxis()->GetBinUpEdge(binxmax),
0100                                hist2D_->GetNbinsY(),
0101                                hist2D_->GetYaxis()->GetXmin(),
0102                                hist2D_->GetYaxis()->GetXmax());
0103     for (int binyc = 1; binyc < hist2D_->GetNbinsY() + 1; ++binyc) {
0104       for (int binxc = binxmin; binxc < binxmax + 1; ++binxc) {
0105         // std::cout << "hist2D_->GetBinContent(" << binxc << "," << binyc << ")
0106         // = " << hist2D_->GetBinContent(binxc,binyc) << std::endl; std::cout <<
0107         // "hist2D_->GetBinError(" << binxc << "," << binyc << ") = " <<
0108         // hist2D_->GetBinError(binxc,binyc) << std::endl; std::cout <<
0109         // "binxc-binxmin+1 = " << binxc-binxmin+1 << std::endl;
0110         rebinnedHist2D_->SetBinContent(binxc - binxmin + 1, binyc, hist2D_->GetBinContent(binxc, binyc));
0111         rebinnedHist2D_->SetBinError(binxc - binxmin + 1, binyc, hist2D_->GetBinError(binxc, binyc));
0112       }
0113     }
0114     rebinnedHist2D_->RebinX(rebinFactor);
0115
0116     const string averageName = bname + "_average";
0117     average_ = new TH1D(averageName.c_str(),
0118                         "arithmetic average",
0119                         rebinnedHist2D_->GetNbinsX(),
0120                         rebinnedHist2D_->GetXaxis()->GetXmin(),
0121                         rebinnedHist2D_->GetXaxis()->GetXmax());
0122
0123     const string rmsName = bname + "_RMS";
0124     RMS_ = new TH1D(rmsName.c_str(),
0125                     "RMS",
0126                     rebinnedHist2D_->GetNbinsX(),
0127                     rebinnedHist2D_->GetXaxis()->GetXmin(),
0128                     rebinnedHist2D_->GetXaxis()->GetXmax());
0129
0130     const string sigmaGaussName = bname + "_sigmaGauss";
0131     sigmaGauss_ = new TH1D(sigmaGaussName.c_str(),
0132                            "sigmaGauss",
0133                            rebinnedHist2D_->GetNbinsX(),
0134                            rebinnedHist2D_->GetXaxis()->GetXmin(),
0135                            rebinnedHist2D_->GetXaxis()->GetXmax());
0136
0137     const string meanXName = bname + "_meanX";
0138     meanXslice_ = new TH1D(meanXName.c_str(),
0139                            "meanX",
0140                            rebinnedHist2D_->GetNbinsX(),
0141                            rebinnedHist2D_->GetXaxis()->GetXmin(),
0142                            rebinnedHist2D_->GetXaxis()->GetXmax());
0143   } else {
0144     // binning is not constant, but made to obtain almost the same number of
0145     // events in each bin
0146
0147     // std::cout << "binxmax-binxmin+1 = " << binxmax-binxmin+1 << std::endl;
0148     // std::cout << "rebinFactor = " << rebinFactor << std::endl;
0149     // std::cout << "(binxmax-binxmin+1)/rebinFactor = " <<
0150     // (binxmax-binxmin+1)/rebinFactor << std::endl; std::cout <<
0151     // "((binxmax-binxmin+1)%rebinFactor) = " <<
0152     // ((binxmax-binxmin+1)%rebinFactor) << std::endl; std::cout <<
0153     // "abs((binxmax-binxmin+1)/rebinFactor) = " <<
0154     // std::abs((binxmax-binxmin+1)/rebinFactor) << std::endl;
0155
0156     unsigned int nbin = 0;
0157     if (((binxmax - binxmin + 1) % rebinFactor) != 0.0) {
0158       nbin = std::abs((binxmax - binxmin + 1) / rebinFactor) + 1;
0159     } else
0160       nbin = (binxmax - binxmin + 1) / rebinFactor;
0161
0162     double *xlow = new double[nbin + 1];
0163     int *binlow = new int[nbin + 1];
0164
0165     TH1D *h0_slice1 = hist2D_->ProjectionY("h0_slice1", binxmin, binxmax, "");
0166     const unsigned int totalNumberOfEvents = static_cast<unsigned int>(h0_slice1->GetEntries());
0167     // std::cout << "totalNumberOfEvents = " << totalNumberOfEvents <<
0168     // std::endl;
0169     delete h0_slice1;
0170
0171     unsigned int neventsc = 0;
0172     unsigned int binXmaxc = binxmin + 1;
0173     xlow[0] = hist2D_->GetXaxis()->GetBinLowEdge(binxmin);
0174     binlow[0] = binxmin;
0175     for (unsigned int binc = 1; binc < nbin; ++binc) {
0176       while (neventsc < binc * totalNumberOfEvents / nbin) {
0177         TH1D *h0_slice1c = hist2D_->ProjectionY("h0_slice1", binxmin, binXmaxc, "");
0178         neventsc = static_cast<unsigned int>(h0_slice1c->GetEntries());
0179         //        //std::cout << "FL : neventsc = " << neventsc << std::endl;
0180         //        //std::cout << "FL : binXmaxc = " << binXmaxc << std::endl;
0181         ++binXmaxc;
0182         delete h0_slice1c;
0183       }
0184       // std::cout << "binXmaxc-1 = " << binXmaxc-1 << std::endl;
0185       binlow[binc] = binXmaxc - 1;
0186       xlow[binc] = hist2D_->GetXaxis()->GetBinLowEdge(binXmaxc - 1);
0187     }
0188     xlow[nbin] = hist2D_->GetXaxis()->GetBinUpEdge(binxmax);
0189     binlow[nbin] = binxmax;
0190
0191     // for (unsigned int binc=0;binc<nbin+1;++binc)
0192     //{
0193     //  std::cout << "xlow[" << binc << "] = " << xlow[binc] << std::endl;
0194     //}
0195
0196     rebinnedHist2D_ = new TH2D(rebinName.c_str(),
0197                                "rebinned histo",
0198                                nbin,
0199                                xlow,
0200                                hist2D_->GetNbinsY(),
0201                                hist2D_->GetYaxis()->GetXmin(),
0202                                hist2D_->GetYaxis()->GetXmax());
0203     for (int binyc = 1; binyc < hist2D_->GetNbinsY() + 1; ++binyc) {
0204       for (unsigned int binxc = 1; binxc < nbin + 1; ++binxc) {
0205         double sum = 0.0;
0206         for (int binxh2c = binlow[binxc - 1]; binxh2c < binlow[binxc]; ++binxh2c) {
0207           sum += hist2D_->GetBinContent(binxh2c, binyc);
0208         }
0209         rebinnedHist2D_->SetBinContent(binxc, binyc, sum);
0210       }
0211     }
0212
0213     const string averageName = bname + "_average";
0214     average_ = new TH1D(averageName.c_str(), "arithmetic average", nbin, xlow);
0215
0216     const string rmsName = bname + "_RMS";
0217     RMS_ = new TH1D(rmsName.c_str(), "RMS", nbin, xlow);
0218
0219     const string sigmaGaussName = bname + "_sigmaGauss";
0220     sigmaGauss_ = new TH1D(sigmaGaussName.c_str(), "sigmaGauss", nbin, xlow);
0221
0222     const string meanXName = bname + "_meanX";
0223     meanXslice_ = new TH1D(meanXName.c_str(), "meanX", nbin, xlow);
0224     delete[] xlow;
0225     delete[] binlow;
0226   }
0227   ProcessSlices(rebinnedHist2D_);
0228 }
0229
0230 void TH2Analyzer::ProcessSlices(const TH2D *histo) {
0231   // std::cout << "ProcessSlices!" << std::endl;
0232
0233   TH1::AddDirectory(false);
0234
0235   for (int i = 1; i <= histo->GetNbinsX(); ++i) {
0236     TH1D *proj = histo->ProjectionY("toto", i, i);
0237     const double mean = proj->GetMean();
0238     const double rms = proj->GetRMS();
0239     // std::cout << "mean = " << mean << std::endl;
0240     // std::cout << "rms = " << rms << std::endl;
0241     average_->SetBinContent(i, mean);
0242     average_->SetBinError(i, proj->GetMeanError());
0243     RMS_->SetBinContent(i, rms);
0244     RMS_->SetBinError(i, proj->GetRMSError());
0245
0246     const double error = (histo->GetXaxis()->GetBinUpEdge(i) - histo->GetXaxis()->GetBinLowEdge(i)) / 2.0;
0247     // std::cout << "error = " << error << std::endl;
0248     meanXslice_->SetBinContent(i, histo->GetXaxis()->GetBinLowEdge(i) + error);
0249     meanXslice_->SetBinError(i, error);
0250     // std::cout << "histo->GetXaxis()->GetBinLowEdge(" << i << ") = "
0251     //          << histo->GetXaxis()->GetBinLowEdge(i) << std::endl;
0252     // std::cout << "meanXslice_->GetBinError(" << i << ") = "
0253     //        << meanXslice_->GetBinError(i) << std::endl;
0254     ProcessSlice(i, proj);
0255     delete proj;
0256   }
0257
0258   TH1::AddDirectory(true);
0259 }
0260
0261 void TH2Analyzer::ProcessSlice(const int i, TH1D *proj) const {
0262   // const double mean = proj->GetMean();
0263   const double rms = proj->GetRMS();
0264   // std::cout << "FL: mean = " << mean << std::endl;
0265   // std::cout << "FL: rms = " << rms << std::endl;
0266
0267   if (rms != 0.0) {
0268     const double fitmin = proj->GetMean() - proj->GetRMS();
0269     const double fitmax = proj->GetMean() + proj->GetRMS();
0270
0271     // std::cout << "i = " << i << std::endl;
0272     // std::cout << "fitmin = " << fitmin << std::endl;
0273     // std::cout << "fitmax = " << fitmax << std::endl;
0274
0275     // proj->Draw();
0276     TF1 *f1 = new TF1("f1", "gaus", fitmin, fitmax);
0277     f1->SetParameters(proj->GetRMS(), proj->GetMean(), proj->GetBinContent(proj->GetMaximumBin()));
0278     proj->Fit(f1, "R", "", proj->GetXaxis()->GetXmin(), proj->GetXaxis()->GetXmax());
0279
0280     // std::ostringstream oss;
0281     // oss << i;
0282     // const std::string plotfitname="Plots/fitbin_"+oss.str()+".eps";
0283     // gPad->SaveAs( plotfitname.c_str() );
0284     // std::cout << "param1 = " << f1->GetParameter(1) << std::endl;
0285     // std::cout << "param2 = " << f1->GetParameter(2) << std::endl;
0286     // std::cout << "paramError2 = " << f1->GetParError(2) << std::endl;
0287
0288     sigmaGauss_->SetBinContent(i, f1->GetParameter(2));
0289     sigmaGauss_->SetBinError(i, f1->GetParError(2));
0290     delete f1;
0291   } else {
0292     sigmaGauss_->SetBinContent(i, 0.0);
0293     sigmaGauss_->SetBinError(i, 0.0);
0294   }
0295 }