NuHistSmoother Class Reference

#include <NuHistSmoother.h>

List of all members.

Public Member Functions

 NuHistSmoother (const TH1D *hist)
virtual ~NuHistSmoother ()
void SmoothHist ()
void SmoothHist (const TH1D *hist)
void SetHist (const TH1D *hist)
void SetNumIterations (int numit)
void SetNumTwicing (int ntwc)
void SetInterpPower (int ipow)
void SetNdfPerPower (float ndf)
void SetSigmaConstraint (float sig=6)
void SetParameters (int numit=3, int ntwc=2, int ipow=2, float ndf=1, float sig=-1)
double Eval (double x)
void MakeSpline ()
TSpline5 * GetSpline ()
TH1D * GetHist ()

Private Member Functions

void Cleanup ()
void SmoothHistOnce (TH1D *hist)
void SmoothFit (TH1D *hist, int nb, bool twc=false)
double Solve (vector< vector< double > > A, double *b, int dim)

Private Attributes

int fNumIterations
int fNumTwicing
int fInterpPower
float fNdfPerPower
float fSigmaConstraint
TH1D * fOrigHist
TH1D * fSmoothHist
TSpline5 * fSpline

Detailed Description

Definition at line 12 of file NuHistSmoother.h.


Constructor & Destructor Documentation

NuHistSmoother::NuHistSmoother ( const TH1D *  hist  ) 

Definition at line 8 of file NuHistSmoother.cxx.

00008                                               {
00009 
00010   fNumIterations = 3;
00011   fNumTwicing = 2;
00012   fInterpPower = 2;
00013   fNdfPerPower = 1;
00014   fSigmaConstraint = -1;
00015 
00016   fOrigHist = 0;
00017   fSmoothHist = 0;
00018   fSpline = 0;
00019   
00020   this->SetHist(hist);
00021 
00022 }

NuHistSmoother::~NuHistSmoother (  )  [virtual]

Definition at line 24 of file NuHistSmoother.cxx.

References Cleanup().

00024                                {
00025   Cleanup();
00026 }


Member Function Documentation

void NuHistSmoother::Cleanup (  )  [private]

Definition at line 52 of file NuHistSmoother.cxx.

References fOrigHist, fSmoothHist, and fSpline.

Referenced by SetHist(), and ~NuHistSmoother().

00052                             {
00053 
00054   if(fOrigHist){ delete fOrigHist; fOrigHist = 0; }
00055   if(fSmoothHist){ delete fSmoothHist; fSmoothHist = 0; }
00056   if(fSpline){ delete fSpline; fSpline = 0; }
00057 
00058 }

double NuHistSmoother::Eval ( double  x  ) 

Definition at line 64 of file NuHistSmoother.cxx.

References fSpline, and MakeSpline().

00064                                    {
00065 
00066   if(!fSpline) this->MakeSpline();
00067 
00068   return fSpline->Eval(x);
00069 
00070 }

TH1D* NuHistSmoother::GetHist (  )  [inline]

Definition at line 37 of file NuHistSmoother.h.

References fSmoothHist, and SmoothHist().

00037 { if(!fSmoothHist) SmoothHist(); return fSmoothHist; }

TSpline5* NuHistSmoother::GetSpline (  )  [inline]

Definition at line 36 of file NuHistSmoother.h.

References fSpline, and MakeSpline().

00036 { if(!fSpline) MakeSpline(); return fSpline; }

void NuHistSmoother::MakeSpline (  ) 

Definition at line 76 of file NuHistSmoother.cxx.

References fSmoothHist, fSpline, and SmoothHist().

Referenced by Eval(), GetSpline(), and SmoothHist().

00076                                {
00077 
00078   if(!fSmoothHist){
00079     this->SmoothHist(); 
00080     return;
00081   }
00082   
00083   if(fSpline){ delete fSpline; fSpline = 0; }
00084 
00085   // Determine the interpolation nodes
00086   int nbins = fSmoothHist->GetNbinsX();
00087 
00088   vector<double> hx;
00089   vector<double> hy;
00090   for(int i=1; i<=nbins; ++i){
00091     if(fSmoothHist->GetBinError(i)==0) continue;
00092     hx.push_back(fSmoothHist->GetBinCenter(i));
00093     hy.push_back(fSmoothHist->GetBinContent(i));
00094   }
00095 
00096   int npx = hx.size();
00097    
00098   // Interpolate nodes
00099   fSpline = new TSpline5("",&hx[0],&hy[0],npx,"b1e1b2e2",0,0,0,0);
00100 
00101 }

void NuHistSmoother::SetHist ( const TH1D *  hist  ) 

Definition at line 43 of file NuHistSmoother.cxx.

References Cleanup(), and fOrigHist.

Referenced by SmoothHist().

00043                                             {
00044 
00045   this->Cleanup();
00046   fOrigHist = (TH1D*)hist->Clone("");
00047 
00048 }

void NuHistSmoother::SetInterpPower ( int  ipow  )  [inline]

Definition at line 25 of file NuHistSmoother.h.

References fInterpPower.

00025 { fInterpPower = ipow; }

void NuHistSmoother::SetNdfPerPower ( float  ndf  )  [inline]

Definition at line 26 of file NuHistSmoother.h.

References fNdfPerPower.

00026 { fNdfPerPower = ndf; }

void NuHistSmoother::SetNumIterations ( int  numit  )  [inline]

Definition at line 23 of file NuHistSmoother.h.

References fNumIterations.

00023 { fNumIterations = numit; }

void NuHistSmoother::SetNumTwicing ( int  ntwc  )  [inline]

Definition at line 24 of file NuHistSmoother.h.

References fNumTwicing.

00024 { fNumTwicing = ntwc; }

void NuHistSmoother::SetParameters ( int  numit = 3,
int  ntwc = 2,
int  ipow = 2,
float  ndf = 1,
float  sig = -1 
)

Definition at line 30 of file NuHistSmoother.cxx.

References fInterpPower, fNdfPerPower, fNumIterations, fNumTwicing, and fSigmaConstraint.

00031                                                         {
00032 
00033   fNumIterations = numit;
00034   fNumTwicing = ntwc;
00035   fInterpPower = ipow;
00036   fNdfPerPower = ndf;
00037   fSigmaConstraint = sig;
00038 
00039 }

void NuHistSmoother::SetSigmaConstraint ( float  sig = 6  )  [inline]

Definition at line 27 of file NuHistSmoother.h.

References fSigmaConstraint.

00027 { fSigmaConstraint = sig; }

void NuHistSmoother::SmoothFit ( TH1D *  hist,
int  nb,
bool  twc = false 
) [private]

Definition at line 203 of file NuHistSmoother.cxx.

References MinosMaterial::A(), fInterpPower, fOrigHist, fSigmaConstraint, init(), and Solve().

Referenced by SmoothHistOnce().

00203                                                           {
00204 
00205   vector<int> hidx;
00206   vector<double> hx;
00207   vector<double> hy;
00208   vector<double> herr;
00209 
00210   // Store non-zero entries
00211   for(int i=1; i<=hist->GetNbinsX(); ++i){
00212     if(hist->GetBinError(i)==0) continue;
00213     hidx.push_back(i);
00214     hx.push_back(hist->GetBinCenter(i));
00215     hy.push_back(hist->GetBinContent(i));
00216     herr.push_back(pow(hist->GetBinError(i),2));
00217   }
00218 
00219   // Number of non-zero entries
00220   int nbins = hx.size();
00221 
00222   // Make sure we have enough neighbors
00223   if(nbins<2*nb+1) nb = (nbins-1)/2;
00224 
00225   // Declare matrix A
00226   vector< vector<double> > A(fInterpPower+1, vector<double>(fInterpPower+1,0));
00227 
00228   // Declare arrays
00229   double *x = new double[2*fInterpPower+1];
00230   double *yx = new double[fInterpPower+1];
00231 
00232   // Loop over non-zero bins
00233   for(int i=0; i<nbins; ++i){
00234 
00235     // Reset arrays
00236     for(int k=0;k<2*fInterpPower+1;++k) x[k]=0;
00237     for(int k=0;k<fInterpPower+1;++k)  yx[k]=0;
00238 
00239     // Define first bin
00240     int init = i - nb;
00241     if(init<0) init = 0;
00242     if(init+2*nb>=nbins) init = nbins - 2*nb -1;
00243 
00244     // Fill x and yx arrays
00245     for(int j=init; j<=init+2*nb; ++j){
00246       for(int k=0;k<2*fInterpPower+1;++k) x[k] += pow(hx[j]-hx[i],k) / herr[j];
00247       for(int k=0;k<fInterpPower+1;++k)  yx[k] += hy[j] * pow(hx[j]-hx[i],k) / herr[j];
00248     }
00249 
00250     // Define matrix A
00251     for(int k=0;k<fInterpPower+1;++k){
00252       for(int l=0;l<fInterpPower+1;++l){
00253         A[k][l] = x[k+l];
00254       }
00255     }
00256 
00257     // Solve linear system
00258     double ynew = Solve(A,yx,fInterpPower+1);        
00259 
00260     // Set limit on deviations to N sigma
00261     if(fSigmaConstraint>0){
00262 
00263       double yorig   = 0;
00264       if(!twc) yorig = fOrigHist->GetBinContent(hidx[i]);
00265 
00266       double yerr  = fOrigHist->GetBinError(hidx[i]);
00267 
00268       if(ynew-yorig > fSigmaConstraint * yerr){
00269         ynew = yorig + fSigmaConstraint * yerr;
00270       }
00271       else if(yorig-ynew > fSigmaConstraint * yerr){
00272         ynew = yorig - fSigmaConstraint * yerr;
00273       }
00274 
00275     }
00276     
00277     // Change histogram
00278     hist->SetBinContent(hidx[i], ynew);
00279 
00280   }// End of loop over non-zero bins
00281 
00282   // Delete arrays
00283   delete x;
00284   delete yx;
00285 
00286 }

void NuHistSmoother::SmoothHist ( const TH1D *  hist  ) 

Definition at line 107 of file NuHistSmoother.cxx.

References SetHist(), and SmoothHist().

00107                                                {
00108 
00109   this->SetHist(hist);
00110   this->SmoothHist();  
00111 
00112 }

void NuHistSmoother::SmoothHist (  ) 

Definition at line 118 of file NuHistSmoother.cxx.

References fNumIterations, fOrigHist, fSmoothHist, MakeSpline(), and SmoothHistOnce().

Referenced by GetHist(), MakeSpline(), and SmoothHist().

00118                                {
00119 
00120   assert(fOrigHist);
00121 
00122   if(fSmoothHist){ delete fSmoothHist; fSmoothHist = 0; }
00123 
00124   fSmoothHist = (TH1D*)fOrigHist->Clone("");
00125   
00126   // Smooth N times
00127   for(int i=0; i<fNumIterations; ++i){
00128     SmoothHistOnce(fSmoothHist);
00129   }
00130 
00131   // Generate the spline interpolation
00132   this->MakeSpline();
00133 
00134 }

void NuHistSmoother::SmoothHistOnce ( TH1D *  hist  )  [private]

Definition at line 141 of file NuHistSmoother.cxx.

References fInterpPower, fNdfPerPower, fNumTwicing, and SmoothFit().

Referenced by SmoothHist().

00141                                              {
00142 
00143   TH1D *hBuffer = (TH1D*)hist->Clone("");
00144 
00145   // Count non-empty bins
00146   int nbins = 0;
00147   for(int i=1;i<=hist->GetNbinsX();i++){
00148     if(hist->GetBinError(i)>0) nbins++;
00149   }
00150 
00151   // Calculate ndf per fit
00152   int ndf = ceil(fInterpPower * fNdfPerPower);
00153 
00154   if(ndf<1) ndf = 1;
00155 
00156   // Calculate number of neighbors
00157   int minnb = (fInterpPower + 1 + ndf)/2;
00158 
00159   if(2*minnb+3>nbins) minnb = (nbins - 3)/2;
00160 
00161   // Smooth 3 times with (n,n+1,n) structure
00162   SmoothFit(hBuffer,minnb);
00163   SmoothFit(hBuffer,minnb+1);
00164   SmoothFit(hBuffer,minnb);
00165 
00166   // Improve matching via twicing N times
00167   for(int twc=0;twc<fNumTwicing;++twc){
00168 
00169     // Reset errors
00170     for(int i=1;i<=hist->GetNbinsX();++i){
00171       hBuffer->SetBinError(i,0);
00172     }
00173     
00174     // Compute residuals
00175     TH1D* hResidual = (TH1D*)hist->Clone("");
00176     hResidual->Add(hBuffer,-1);
00177     
00178     // Smooth residuals
00179     SmoothFit(hResidual, minnb, true);
00180     SmoothFit(hResidual, minnb+1, true);
00181     SmoothFit(hResidual, minnb, true);
00182 
00183     // Add residuals back
00184     hBuffer->Add(hResidual);
00185 
00186     delete hResidual;
00187 
00188   }
00189 
00190   // Set actual histogram bin contents
00191   for(int i=1;i<=hist->GetNbinsX();++i){
00192     hist->SetBinContent(i,hBuffer->GetBinContent(i));
00193   }
00194   
00195   delete hBuffer;
00196 
00197 }

double NuHistSmoother::Solve ( vector< vector< double > >  A,
double *  b,
int  dim 
) [private]

Definition at line 293 of file NuHistSmoother.cxx.

Referenced by SmoothFit().

00293                                                                         {
00294 
00295   for(int col=dim-1; col>0; col--){
00296     for(int row=0; row<col; ++row){
00297       if(A[row][col]==0) continue;
00298       b[row] = b[row]*A[col][col]/A[row][col] - b[col];
00299       for(int coli=0; coli<col; ++coli){
00300         A[row][coli] = A[row][coli]*A[col][col]/A[row][col] - A[col][coli];
00301       }
00302     }
00303   }
00304 
00305   if(A[0][0]==0) return 0;
00306   else           return b[0] / A[0][0];
00307 
00308 }


Member Data Documentation

Definition at line 49 of file NuHistSmoother.h.

Referenced by SetInterpPower(), SetParameters(), SmoothFit(), and SmoothHistOnce().

Definition at line 50 of file NuHistSmoother.h.

Referenced by SetNdfPerPower(), SetParameters(), and SmoothHistOnce().

Definition at line 47 of file NuHistSmoother.h.

Referenced by SetNumIterations(), SetParameters(), and SmoothHist().

Definition at line 48 of file NuHistSmoother.h.

Referenced by SetNumTwicing(), SetParameters(), and SmoothHistOnce().

TH1D* NuHistSmoother::fOrigHist [private]

Definition at line 53 of file NuHistSmoother.h.

Referenced by Cleanup(), SetHist(), SmoothFit(), and SmoothHist().

Definition at line 51 of file NuHistSmoother.h.

Referenced by SetParameters(), SetSigmaConstraint(), and SmoothFit().

TH1D* NuHistSmoother::fSmoothHist [private]

Definition at line 54 of file NuHistSmoother.h.

Referenced by Cleanup(), GetHist(), MakeSpline(), and SmoothHist().

TSpline5* NuHistSmoother::fSpline [private]

Definition at line 56 of file NuHistSmoother.h.

Referenced by Cleanup(), Eval(), GetSpline(), and MakeSpline().


The documentation for this class was generated from the following files:

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1