test/Macros/Smoother.C

0001 #include <iostream>
0002 #include "TH1F.h"
0003 #include "TString.h"
0004
0005 /**
0006  * This class receives a TH1 histogram and applies a smoothing.
0007  * The smoothing can be done in two ways:
0008  * - with the applySmooth method:
0009  * -- for each bin take the two bins on the left and the right
0010  * -- compute the average of the sidepoints and compare with the middle
0011  * -- if the difference is more than N% (configurable) replace the point with the average
0012  * - with the applySmoothSingle method:
0013  * -- the same but using only the closest bin on left and right.
0014  */
0015
0016 class Smoother {
0017  public:
0018   Smoother( const int N = 0.5 ) : ratio_(N) {}
0019
0020   TH1F * smooth( TH1F * histo, const int iterations, const bool * single, const double & newHistoXmin = 71.1876, const double & newHistoXmax = 111.1876, const int newHistoBins = 1001 );
0021
0022  protected:
0023   TH1F * applySmooth( TH1F * histo, const double & newHistoXmin, const double & newHistoXmax, const int newHistoBins );
0024   TH1F * applySmoothSingle( TH1F * histo, const double & newHistoXmin, const double & newHistoXmax, const int newHistoBins );
0025   int ratio_;
0026 };
0027
0028 TH1F * Smoother::applySmooth( TH1F * histo, const double & newHistoXmin, const double & newHistoXmax, const int newHistoBins )
0029 {
0030   TH1F * smoothedHisto = (TH1F*)histo->Clone();
0031   smoothedHisto->Reset();
0032
0033   histo->SetAxisRange(newHistoXmin, newHistoXmax);
0034   int xBinMin = histo->FindBin(newHistoXmin);
0035   int xBinMax = histo->FindBin(newHistoXmax);
0036   std::cout << "xBinMin = " << xBinMin << std::endl;
0037   std::cout << "xBinMax = " << xBinMax << std::endl;
0038
0039   double xMin = histo->GetXaxis()->GetXmin();
0040   double xMax = histo->GetXaxis()->GetXmax();
0041   int nBins = histo->GetNbinsX();
0042   std::cout << "xMin = " << xMin << std::endl;
0043   std::cout << "xMax = " << xMax << std::endl;
0044   std::cout << "nBins = " << nBins << std::endl;
0045
0046   // Start from the third bin up to the second to last bin
0047   int smoothedBin = 2;
0048   smoothedHisto->SetBinContent(0, histo->GetBinContent(0));
0049   smoothedHisto->SetBinContent(1, histo->GetBinContent(1));
0050   for( int i=2+xBinMin; i<xBinMax-2; ++i, ++smoothedBin ) {
0051     double min2 = histo->GetBinContent(i-2);
0052     double min1 = histo->GetBinContent(i-1);
0053     double med = histo->GetBinContent(i);
0054     double plus1 = histo->GetBinContent(i+1);
0055     double plus2 = histo->GetBinContent(i+2);
0056     // If the slope is the same before and after (med is not a min/max)
0057     if( (min1 - min2)*(plus2 - plus1) > 0 ) {
0058       // compare the med with the mean of the four points
0059       double newMed = ((min1+min2+plus1+plus2)/4);
0060       if( fabs(med/newMed - 1) > ratio_ ) {
0061     std::cout << "Replacing value for bin " << i << " from " << med << " to " << newMed << std::endl;
0062     // smoothedHisto->SetBinContent(smoothedBin, newMed);
0063     smoothedHisto->SetBinContent(i, newMed);
0064     continue;
0065       }
0066     }
0067     // smoothedHisto->SetBinContent(smoothedBin, med);
0068     smoothedHisto->SetBinContent(i, med);
0069     std::cout << "bin["<<i<<"] = " << histo->GetBinContent(i) << std::endl;
0070   }
0071   return smoothedHisto;
0072 }
0073
0074 TH1F * Smoother::applySmoothSingle( TH1F * histo, const double & newHistoXmin, const double & newHistoXmax, const int newHistoBins )
0075 {
0076   TH1F * smoothedHisto = (TH1F*)histo->Clone();
0077   smoothedHisto->Reset();
0078
0079   histo->SetAxisRange(newHistoXmin, newHistoXmax);
0080   int xBinMin = histo->FindBin(newHistoXmin);
0081   int xBinMax = histo->FindBin(newHistoXmax);
0082   std::cout << "xBinMin = " << xBinMin << std::endl;
0083   std::cout << "xBinMax = " << xBinMax << std::endl;
0084
0085   double xMin = histo->GetXaxis()->GetXmin();
0086   double xMax = histo->GetXaxis()->GetXmax();
0087   int nBins = histo->GetNbinsX();
0088   std::cout << "xMin = " << xMin << std::endl;
0089   std::cout << "xMax = " << xMax << std::endl;
0090   std::cout << "nBins = " << nBins << std::endl;
0091
0092   // Start from the third bin up to the second to last bin
0093   int smoothedBin = 2;
0094   smoothedHisto->SetBinContent(0, histo->GetBinContent(0));
0095   smoothedHisto->SetBinContent(1, histo->GetBinContent(1));
0096   for( int i=2+xBinMin; i<xBinMax-2; ++i, ++smoothedBin ) {
0097     double min1 = histo->GetBinContent(i-1);
0098     double med = histo->GetBinContent(i);
0099     double plus1 = histo->GetBinContent(i+1);
0100     // compare the med with the mean of the four points
0101     double newMed = ((min1+plus1)/2);
0102     if( fabs(med/newMed - 1) > ratio_ ) {
0103       std::cout << "Replacing value for bin " << i << " from " << med << " to " << newMed << std::endl;
0104       smoothedHisto->SetBinContent(i, newMed);
0105       continue;
0106     }
0107     smoothedHisto->SetBinContent(i, med);
0108     std::cout << "bin["<<i<<"] = " << histo->GetBinContent(i) << std::endl;
0109   }
0110   return smoothedHisto;
0111 }
0112
0113
0114 TH1F * Smoother::smooth( TH1F * histo, const int iterations, const bool * single, const double & newHistoXmin, const double & newHistoXmax, const int newHistoBins )
0115 {
0116   histo->Draw();
0117   TH1F * smoothedHisto = histo;
0118   for( int i=0; i<iterations; ++i ) {
0119     if( single[i] ) {
0120       smoothedHisto = applySmoothSingle(smoothedHisto, newHistoXmin, newHistoXmax, newHistoBins);
0121     }
0122     else {
0123       smoothedHisto = applySmooth(smoothedHisto, newHistoXmin, newHistoXmax, newHistoBins);
0124     }
0125 //     smoothedHisto->Draw("same");
0126 //     smoothedHisto->SetLineColor(2+i);
0127   }
0128   smoothedHisto->Draw("same");
0129   smoothedHisto->SetLineColor(kRed);
0130   return smoothedHisto;
0131 }