Muon/test/langauFit.C

0001 //-----------------------------------------------------------------------
0002 //
0003 //  Convoluted Landau and Gaussian Fitting Function
0004 //         (using ROOT's Landau and Gauss functions)
0005 //
0006 //  Based on a Fortran code by R.Fruehwirth (fruhwirth@hephy.oeaw.ac.at)
0007 //  Adapted for C++/ROOT by H.Pernegger (Heinz.Pernegger@cern.ch) and
0008 //   Markus Friedl (Markus.Friedl@cern.ch)
0009 //
0010 //  to execute this example, do:
0011 //  root > .x langaus.C
0012 // or
0013 //  root > .x langaus.C++
0014 //
0015 //-----------------------------------------------------------------------
0016
0017 #include "TH1.h"
0018 #include "TF1.h"
0019 #include "TROOT.h"
0020 #include "TStyle.h"
0021 #include "TMath.h"
0022
0023 Double_t langaufun(Double_t *x, Double_t *par) {
0024   //Fit parameters:
0025   //par[0]=Width (scale) parameter of Landau density
0026   //par[1]=Most Probable (MP, location) parameter of Landau density
0027   //par[2]=Total area (integral -inf to inf, normalization constant)
0028   //par[3]=Width (sigma) of convoluted Gaussian function
0029   //
0030   //In the Landau distribution (represented by the CERNLIB approximation),
0031   //the maximum is located at x=-0.22278298 with the location parameter=0.
0032   //This shift is corrected within this function, so that the actual
0033   //maximum is identical to the MP parameter.
0034
0035   // Numeric constants
0036   Double_t invsq2pi = 0.3989422804014;  // (2 pi)^(-1/2)
0037   Double_t mpshift = -0.22278298;       // Landau maximum location
0038
0039   // Control constants
0040   Double_t np = 100.0;  // number of convolution steps
0041   Double_t sc = 5.0;    // convolution extends to +-sc Gaussian sigmas
0042
0043   // Variables
0044   Double_t xx;
0045   Double_t mpc;
0046   Double_t fland;
0047   Double_t sum = 0.0;
0048   Double_t xlow, xupp;
0049   Double_t step;
0050   Double_t i;
0051
0052   // MP shift correction
0053   mpc = par[1] - mpshift * par[0];
0054
0055   // Range of convolution integral
0056   xlow = x[0] - sc * par[3];
0057   xupp = x[0] + sc * par[3];
0058
0059   step = (xupp - xlow) / np;
0060
0061   // Convolution integral of Landau and Gaussian by sum
0062   for (i = 1.0; i <= np / 2; i++) {
0063     xx = xlow + (i - .5) * step;
0064     fland = TMath::Landau(xx, mpc, par[0]) / par[0];
0065     sum += fland * TMath::Gaus(x[0], xx, par[3]);
0066
0067     xx = xupp - (i - .5) * step;
0068     fland = TMath::Landau(xx, mpc, par[0]) / par[0];
0069     sum += fland * TMath::Gaus(x[0], xx, par[3]);
0070   }
0071
0072   return (par[2] * step * sum * invsq2pi / par[3]);
0073 }
0074
0075 TF1 *langaufit(TH1F *his,
0076                Double_t *fitrange,
0077                Double_t *startvalues,
0078                Double_t *parlimitslo,
0079                Double_t *parlimitshi,
0080                Double_t *fitparams,
0081                Double_t *fiterrors,
0082                Double_t *ChiSqr,
0083                Int_t *NDF) {
0084   // Once again, here are the Landau * Gaussian parameters:
0085   //   par[0]=Width (scale) parameter of Landau density
0086   //   par[1]=Most Probable (MP, location) parameter of Landau density
0087   //   par[2]=Total area (integral -inf to inf, normalization constant)
0088   //   par[3]=Width (sigma) of convoluted Gaussian function
0089   //
0090   // Variables for langaufit call:
0091   //   his             histogram to fit
0092   //   fitrange[2]     lo and hi boundaries of fit range
0093   //   startvalues[4]  reasonable start values for the fit
0094   //   parlimitslo[4]  lower parameter limits
0095   //   parlimitshi[4]  upper parameter limits
0096   //   fitparams[4]    returns the final fit parameters
0097   //   fiterrors[4]    returns the final fit errors
0098   //   ChiSqr          returns the chi square
0099   //   NDF             returns ndf
0100
0101   Int_t i;
0102   Char_t FunName[100];
0103
0104   sprintf(FunName, "Fitfcn_%s", his->GetName());
0105
0106   TF1 *ffitold = (TF1 *)gROOT->GetListOfFunctions()->FindObject(FunName);
0107   if (ffitold)
0108     delete ffitold;
0109
0110   TF1 *ffit = new TF1(FunName, langaufun, fitrange[0], fitrange[1], 4);
0111   ffit->SetParameters(startvalues);
0112   ffit->SetParNames("Width", "MP", "Area", "GSigma");
0113
0114   for (i = 0; i < 4; i++) {
0115     ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
0116   }
0117
0118   try {
0119     his->Fit(FunName, "RB0Q SERIAL");  // fit within specified range, use ParLimits, do not plot
0120   } catch (...) {
0121   }
0122
0123   if (ffit) {
0124     ffit->GetParameters(fitparams);  // obtain fit parameters
0125     for (i = 0; i < 4; i++) {
0126       fiterrors[i] = ffit->GetParError(i);  // obtain fit parameter errors
0127     }
0128     ChiSqr[0] = ffit->GetChisquare();  // obtain chi^2
0129     NDF[0] = ffit->GetNDF();           // obtain ndf
0130   }
0131
0132   return (ffit);  // return fit function
0133 }
0134
0135 Int_t langaupro(Double_t *params, Double_t &maxx, Double_t &FWHM) {
0136   // Seaches for the location (x value) at the maximum of the
0137   // Landau-Gaussian convolute and its full width at half-maximum.
0138   //
0139   // The search is probably not very efficient, but it's a first try.
0140
0141   Double_t p, x, fy, fxr, fxl;
0142   Double_t step;
0143   Double_t l, lold;
0144   Int_t i = 0;
0145   Int_t MAXCALLS = 10000;
0146
0147   // Search for maximum
0148
0149   p = params[1] - 0.1 * params[0];
0150   step = 0.05 * params[0];
0151   lold = -2.0;
0152   l = -1.0;
0153
0154   while ((l != lold) && (i < MAXCALLS)) {
0155     i++;
0156
0157     lold = l;
0158     x = p + step;
0159     l = langaufun(&x, params);
0160
0161     if (l < lold)
0162       step = -step / 10;
0163
0164     p += step;
0165   }
0166
0167   if (i == MAXCALLS)
0168     return (-1);
0169
0170   maxx = x;
0171
0172   fy = l / 2;
0173
0174   // Search for right x location of fy
0175
0176   p = maxx + params[0];
0177   step = params[0];
0178   lold = -2.0;
0179   l = -1e300;
0180   i = 0;
0181
0182   while ((l != lold) && (i < MAXCALLS)) {
0183     i++;
0184
0185     lold = l;
0186     x = p + step;
0187     l = TMath::Abs(langaufun(&x, params) - fy);
0188
0189     if (l > lold)
0190       step = -step / 10;
0191
0192     p += step;
0193   }
0194
0195   if (i == MAXCALLS)
0196     return (-2);
0197
0198   fxr = x;
0199
0200   // Search for left x location of fy
0201
0202   p = maxx - 0.5 * params[0];
0203   step = -params[0];
0204   lold = -2.0;
0205   l = -1e300;
0206   i = 0;
0207
0208   while ((l != lold) && (i < MAXCALLS)) {
0209     i++;
0210
0211     lold = l;
0212     x = p + step;
0213     l = TMath::Abs(langaufun(&x, params) - fy);
0214
0215     if (l > lold)
0216       step = -step / 10;
0217
0218     p += step;
0219   }
0220
0221   if (i == MAXCALLS)
0222     return (-3);
0223
0224   fxl = x;
0225
0226   FWHM = fxr - fxl;
0227   return (0);
0228 }