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

From: Robert Feuerbach (feuerbac@ernest.phys.cmu.edu)
Date: Wed Mar 07 2001 - 19:55:34 MET


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