Re: [ROOT] Fitting a TH1+TF1 to a TH1

From: Rene Brun (Rene.Brun@cern.ch)
Date: Thu Mar 08 2001 - 12:59:52 MET


Hi Robert,
I have added the two following inline functions:

 Double_t  *TFormula::GetParameters() {return fParams;}
 Double_t  *TF1::GetParErrors() {return fParErrors;}

Rene Brun

Robert Feuerbach wrote:
> 
> On Tue, 6 Mar 2001, Stephen Bailey wrote:
> > Hello.
> >
> > I have a histogram which is comprised of signal and background.
> > I have a template histogram for the background and a function
> > shape for the signal.  I would like to fit the template histogram
> > and the signal function to the data histogram to get their
> > relative normalizations.
> >
> > Is there a way to do this easily in ROOT?  Thanks.
> >
> > Stephen
> >
> 
> Hi Stephen,
> 
> I'm actually doing the same thing -- fitting a TH1 to another TH1 + some
> arbitrary function.
> 
> I ended up creating a TF1 that accesses a C-function (fitfunc), and
> fitfunc then looks up the histogram and function, and adds them together.
> Unfortunately, you have to use some sort of global variable pointer to
> point to the appropriate template histogram, since the permitted
> parameters for fitfunc are rather tightly specified by the TF1 contructor
> you need to use:
> 
> TF1 TF1(const char* name, Double_t (*)(Double_t*, Double_t*) fcn, Double_t
> xmin = 0, Double_t xmax = 1, Int_t npar = 0)
> 
> To give you an idea of what I do, an excerpt of my code is below. It isn't
> the prettiest code you'll find, but it does the job.
> 
> Also, I'd like to petition for a
>    TF1::GetParErrors(Double_t *eparams)
> method,  parallel to
>    TF1::GetParameters(Double_t *params)
> 
> Good luck,
> Rob Feuerbach
> 
> //---------------------------------------------------------
> Double_t linear_bck(Double_t *x, Double_t parms[]);
> 
> Double_t mmfitfunc(Double_t *x, Double_t parms[]);
> 
> // create a namespace to hold all the "global variables" you'll need
> // to keep track of which template histogram you are using.
> // I'm fitting to 2(lambda_hist and sigma0_hist) + bckgrnd function
> 
> namespace MMWQTemplate {
>   TDirectory *directory;   // directory template histograms are in
> 
>   TH1 *lambda_hist, *sigma0_hist;  // template histos
>   Double_t (*backgrnd_func)(Double_t *x, Double_t parms[]);
>   Int_t nbck_params;
> 
>   TF1 *fitfunc=NULL;
>   TF1 *bckfunc=NULL;
>   Double_t xmin, xmax;
>   Double_t binsize;
> 
> 
>   void MMFitWQSetup(TDirectory *template_dir) {
>     // Setup some of the "global's" in this namespace
>     directory = template_dir;
> 
>     xmin = 1.05;  xmax = 1.3;  // fit over this range
> 
>     backgrnd_func = linear_bck;
>     nbck_params = 2;
> 
>     // define background function
>     if (!bckfunc) {
>       bckfunc = new TF1("bckfunc",backgrnd_func,xmin,xmax,nbck_params);
>     }
> 
>     // define total fit function
>     if (!fitfunc) {
>       fitfunc = new TF1("mmfitwq",mmfitfunc,xmin,xmax,2+nbck_params);
>       fitfunc->SetParNames("Lambda","Sigma0","Bck Cnts","Bck Slope");
>     }
>     binsize = -1.;
>   }
> 
>   TH1 * GetHypeTemplate(TH1 * hist, int hype) {
>     char temp[30];
>     int Wid, Qid;
> 
>     sscanf(hist->GetName(),"%[^0-9]%d_%d%s",temp,&Wid,&Qid,temp);
>     sprintf(temp,"WQ%d_%d_%dMM",hype,Wid,Qid);
>     return (TH1*)directory->Get(temp);
>   }
> 
>   void MMPrepareFitWQHist(TH1 *hist) {
>     // prepare to fit this histogram, get the template histograms setup,
>     // setup starting values for parameters
>     char temp[30];
>     int Wid,Qid;
> 
>     // get appropriate template histograms : Lambda
>     lambda_hist = GetHypeTemplate(hist,0);
> 
>     // get appropriate template histograms : Sigma0
>     sigma0_hist = GetHypeTemplate(hist,1);
> 
>     // setup initial parameters
>     Double_t parms[10];
>     parms[0] = hist->Integral(hist->FindBin(1.09),hist->FindBin(1.15));
>     parms[1] = hist->Integral(hist->FindBin(1.09),hist->FindBin(1.15));
>     parms[2] = 10.;
>     parms[3] = 0.;
> 
>     fitfunc->SetParameters(parms);
> 
>   }
> 
> }
> 
> Double_t mmfitfunc(Double_t *x, Double_t parms[]) {
>   // parms are set up as:
>   // parms[0] = lambda scale factor
>   // parms[1] = sigma0 scale factor
>   // parms[2...] = background terms
> 
>   Double_t lambda_contrib=0., sigma0_contrib=0., back_contrib=0.;
>   TH1 *lambda_hist = MMWQTemplate::lambda_hist;
>   TH1 *sigma0_hist = MMWQTemplate::sigma0_hist;
> 
>   if (lambda_hist) {
>     Int_t bin = lambda_hist->FindBin(*x);
>     Double_t cont = lambda_hist->GetBinContent(bin);
>     lambda_contrib = parms[0]*cont;
>   }
>   if (sigma0_hist) {
>     Int_t bin = sigma0_hist->FindBin(*x);
>     Double_t cont = sigma0_hist->GetBinContent(bin);
>     sigma0_contrib = parms[1]*cont;
>   }
>   back_contrib = MMWQTemplate::backgrnd_func(x,&(parms[2]));
>   return lambda_contrib+sigma0_contrib+back_contrib;
> }
> 
> Double_t linear_bck(Double_t *x, Double_t parms[]) {
>   // 1rst parameter is always the area of the background (A) -- must be
>   // positive!
>   // 2nd parameter is the slope (b)
>   // f = a+b*x becomes  A/d + b*(x-d/2)
>   Double_t d1 = MMWQTemplate::xmax - MMWQTemplate::xmin;
>   Double_t d2 = MMWQTemplate::xmax + MMWQTemplate::xmin;
>   static Double_t w = MMWQTemplate::binsize;
>   return parms[0]*w/(d1) + parms[1]*w*(*x - d2/2);
> }
> 
> // I can't get ROOT to interactively look inside the namespace, so
> // some wrapper functions were necessary
> 
> void MMFitWQSetup(TDirectory *template_dir) {
>   MMWQTemplate::MMFitWQSetup(template_dir);
> }
> 
> void MMPrepareFitWQHist(TH1 *hist) {
>   MMWQTemplate::MMFitWQHist(hist);
> }
> 
> // Then, to fit, you do:
> 
> void fit_hist(TH1 *hist, TDirectory *template_dir) {
>   // example of how to fit the histogram "hist"
> 
>   // only have to do this Setup once
>   MMFitWQSetup(template_dir);
> 
>   // for each histogram: prepare by selecting the template histos
>   MMPrepareFitWQHist(hist);
> 
>   // do the fit
>   hist->Fit("mmfitwq","RQ");
> 
>   Double_t parms[10];
>   hist->GetFunction("mmfitwq")->GetParameters(parms);
> 
>   // process on the fit results....
> }
> --
> /***************************************************
>  * Robert Feuerbach   feuerbac@ernest.phys.cmu.edu *
>  * CMU Physics Department           (412) 268-2772 *
>  * Pittsburgh, PA 15213        FAX: (412) 681-0648 *
>  ***************************************************/



This archive was generated by hypermail 2b29 : Tue Jan 01 2002 - 17:50:39 MET