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