chebyshev_fit.C

Go to the documentation of this file.

// Chebyshev polynomial functions
Double_t T0(Double_t x){return 1;}
Double_t T1(Double_t x){return x;}
Double_t T2(Double_t x){return 2*pow(x,2)-1;}
Double_t T3(Double_t x){return 4*pow(x,3)-3*x;}
Double_t T4(Double_t x){return 8*pow(x,4)-8*pow(x,2)+1;}
Double_t T5(Double_t x){return 16*pow(x,5)-20*pow(x,3)+5*x;}
Double_t T6(Double_t x){return 32*pow(x,6)-48*pow(x,4)+18*pow(x,2)-1;}
Double_t T7(Double_t x){return 64*pow(x,7)-112*pow(x,5)+56*pow(x,3)-7*x;}
Double_t T8(Double_t x){return 128*pow(x,8)-256*pow(x,6)+160*pow(x,4)-32*pow(x,2)+1;}
Double_t T9(Double_t x){return 256*pow(x,9)-576*pow(x,7)+432*pow(x,5)-120*pow(x,3)+9*x;}


//---------------------
// chebyshev_fit
//---------------------
void chebyshev_fit(void)
{
   gROOT->Reset();
   TColor::CreateColorWheel();
   
   ChebyshevTestFit();
}

//---------------------
// ChebyshevTestFit
//---------------------
void ChebyshevTestFit()
{
   TCanvas *c1 = new TCanvas("c1");
   c1->SetGrid();
   c1->SetTicks();
   
   TH1D *h = new TH1D("h","A test histogram", 1000, -1.0, 1.0);
   h->FillRandom("gaus", 1000000);
   
   TF1 *f = chebyshev_FindBestFunction(h);

   h->Draw();
   f->SetLineColor(kRed);
   f->Draw("same");
   
   UInt_t order = f->GetNpar() - 1;
   cout<<"Chebyshev polynomial fit to data of order: "<<order<<endl;
}

//---------------------
// chebyshev_Fit
//---------------------
TF1* chebyshev_Fit(TH1D *h, UInt_t order=9)
{
   double xmin = h->GetXaxis()->GetXmin();
   double xmax = h->GetXaxis()->GetXmax();
   TF1 *f = chebyshev_MakeTF1(xmin, xmax, order);

   h->Fit(f, "0Q");
   return f;
}

//---------------------
// chebyshev_FindBestFunction
//---------------------
TF1* chebyshev_FindBestFunction(TH1D *h, UInt_t max_order=9)
{
   // Verify histogram limits
   if(h->GetXaxis()->GetXmin()!=-1.0 || h->GetXaxis()->GetXmax()!=+1.0){
      cerr<<"-- ERROR! The limits of the histo passed to chebyshev_FindBestFunction"<<endl;
      cerr<<"-- ERROR! MUST be -1.0 to +1.0 ("<<h->GetXaxis()->GetXmin()<<" and "<<h->GetXaxis()->GetXmax()<<" passed)"<<endl;
   }

   // Loop over all orders keeping track of the one with the best chisq
   UInt_t best_order = 0;
   Double_t best_chisq_per_dof = 1.0E6;
   for(UInt_t order=0; order<=max_order; order++){
      TF1 *f = chebyshev_MakeTF1(order);
      h->Fit(f, "0Q");
      Double_t chisq_per_dof = f->GetChisquare()/(double)f->GetNDF();
      cout<<"order "<<order<<": chisq="<<f->GetChisquare()<<"  NDF="<<f->GetNDF()<<"  chisq/NDF="<<chisq_per_dof<<endl;
      if(chisq_per_dof<best_chisq_per_dof){
         best_chisq_per_dof = chisq_per_dof;
         best_order = order;
      }
   }
   
   // Regenerate and fit with best function
   TF1 *best_func = chebyshev_MakeTF1(best_order);
   h->Fit(best_func, "0Q");
   
   return best_func;
}

//---------------------
// chebyshev_MakeTF1
//---------------------
TF1* chebyshev_MakeTF1(Double_t xmin, Double_t xmax, UInt_t order)
{
   // ROOT gives a "Too many operators !" error when we try and define
   // past order 9. Truncate it with a warning if too high of an order
   // is specified.
   if(order>9){
      cerr<<"--WARNING! The maximum order that can be used is 9"<<endl;
      cerr<<"--WARNING! order="<<order<<" will be reduced to 9"<<endl;
      order = 9; 
   }

   // Build function string of given order
   stringstream func_str;
   for(UInt_t i=0; i<=order; i++){
      if(i!=0)func_str<<"+";
      func_str<<"["<<i+2<<"]*T"<<i<<"((x-[0])/[1])";
   }

   // Create Function
   stringstream funcname;
   funcname<<"polT"<<order;
   TF1 *myfunc = new TF1(funcname.str().c_str(), func_str.str().c_str(), xmin, xmax);

   // Parmeters used to convert x to -1 to 1 coordinates
   double xmid  = (xmax+xmin)/2.0;
   double xnorm = (xmax-xmin)/2.0;

   // Set initial parameter values and parameter names
   myfunc->SetParName(0, "xmid");
   myfunc->SetParName(1, "xnorm");
   myfunc->FixParameter(0, xmid);
   myfunc->FixParameter(1, xnorm);
   for(UInt_t i=0; i<=order; i++){
      stringstream fname;
      fname<<"T"<<i;
      myfunc->SetParName(i+2, fname.str().c_str());
      myfunc->SetParameter(i+2, 1.0);
   }

   return myfunc;
}