ROOT logo

From $ROOTSYS/tutorials/fit/langaus.C

//-----------------------------------------------------------------------
//
// Convoluted Landau and Gaussian Fitting Function
//         (using ROOT's Landau and Gauss functions)
//
//  Based on a Fortran code by R.Fruehwirth (fruhwirth@hephy.oeaw.ac.at)
//  Adapted for C++/ROOT by H.Pernegger (Heinz.Pernegger@cern.ch) and
//   Markus Friedl (Markus.Friedl@cern.ch)
//
//  to execute this example, do:
//  root > .x langaus.C
// or
//  root > .x langaus.C++
//
//-----------------------------------------------------------------------

#include "TH1.h"
#include "TF1.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TMath.h"

Double_t langaufun(Double_t *x, Double_t *par) {

   //Fit parameters:
   //par[0]=Width (scale) parameter of Landau density
   //par[1]=Most Probable (MP, location) parameter of Landau density
   //par[2]=Total area (integral -inf to inf, normalization constant)
   //par[3]=Width (sigma) of convoluted Gaussian function
   //
   //In the Landau distribution (represented by the CERNLIB approximation), 
   //the maximum is located at x=-0.22278298 with the location parameter=0.
   //This shift is corrected within this function, so that the actual
   //maximum is identical to the MP parameter.

      // Numeric constants
      Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
      Double_t mpshift  = -0.22278298;       // Landau maximum location

      // Control constants
      Double_t np = 100.0;      // number of convolution steps
      Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas

      // Variables
      Double_t xx;
      Double_t mpc;
      Double_t fland;
      Double_t sum = 0.0;
      Double_t xlow,xupp;
      Double_t step;
      Double_t i;


      // MP shift correction
      mpc = par[1] - mpshift * par[0]; 

      // Range of convolution integral
      xlow = x[0] - sc * par[3];
      xupp = x[0] + sc * par[3];

      step = (xupp-xlow) / np;

      // Convolution integral of Landau and Gaussian by sum
      for(i=1.0; i<=np/2; i++) {
         xx = xlow + (i-.5) * step;
         fland = TMath::Landau(xx,mpc,par[0]) / par[0];
         sum += fland * TMath::Gaus(x[0],xx,par[3]);

         xx = xupp - (i-.5) * step;
         fland = TMath::Landau(xx,mpc,par[0]) / par[0];
         sum += fland * TMath::Gaus(x[0],xx,par[3]);
      }

      return (par[2] * step * sum * invsq2pi / par[3]);
}



TF1 *langaufit(TH1F *his, Double_t *fitrange, Double_t *startvalues, Double_t *parlimitslo, Double_t *parlimitshi, Double_t *fitparams, Double_t *fiterrors, Double_t *ChiSqr, Int_t *NDF)
{
   // Once again, here are the Landau * Gaussian parameters:
   //   par[0]=Width (scale) parameter of Landau density
   //   par[1]=Most Probable (MP, location) parameter of Landau density
   //   par[2]=Total area (integral -inf to inf, normalization constant)
   //   par[3]=Width (sigma) of convoluted Gaussian function
   //
   // Variables for langaufit call:
   //   his             histogram to fit
   //   fitrange[2]     lo and hi boundaries of fit range
   //   startvalues[4]  reasonable start values for the fit
   //   parlimitslo[4]  lower parameter limits
   //   parlimitshi[4]  upper parameter limits
   //   fitparams[4]    returns the final fit parameters
   //   fiterrors[4]    returns the final fit errors
   //   ChiSqr          returns the chi square
   //   NDF             returns ndf

   Int_t i;
   Char_t FunName[100];

   sprintf(FunName,"Fitfcn_%s",his->GetName());

   TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName);
   if (ffitold) delete ffitold;

   TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4);
   ffit->SetParameters(startvalues);
   ffit->SetParNames("Width","MP","Area","GSigma");
   
   for (i=0; i<4; i++) {
      ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
   }

   his->Fit(FunName,"RB0");   // fit within specified range, use ParLimits, do not plot

   ffit->GetParameters(fitparams);    // obtain fit parameters
   for (i=0; i<4; i++) {
      fiterrors[i] = ffit->GetParError(i);     // obtain fit parameter errors
   }
   ChiSqr[0] = ffit->GetChisquare();  // obtain chi^2
   NDF[0] = ffit->GetNDF();           // obtain ndf

   return (ffit);              // return fit function

}


Int_t langaupro(Double_t *params, Double_t &maxx, Double_t &FWHM) {

   // Seaches for the location (x value) at the maximum of the 
   // Landau-Gaussian convolute and its full width at half-maximum.
   //
   // The search is probably not very efficient, but it's a first try.

   Double_t p,x,fy,fxr,fxl;
   Double_t step;
   Double_t l,lold;
   Int_t i = 0;
   Int_t MAXCALLS = 10000;


   // Search for maximum

   p = params[1] - 0.1 * params[0];
   step = 0.05 * params[0];
   lold = -2.0;
   l    = -1.0;


   while ( (l != lold) && (i < MAXCALLS) ) {
      i++;

      lold = l;
      x = p + step;
      l = langaufun(&x,params);
 
      if (l < lold)
         step = -step/10;
 
      p += step;
   }

   if (i == MAXCALLS)
      return (-1);

   maxx = x;

   fy = l/2;


   // Search for right x location of fy

   p = maxx + params[0];
   step = params[0];
   lold = -2.0;
   l    = -1e300;
   i    = 0;


   while ( (l != lold) && (i < MAXCALLS) ) {
      i++;

      lold = l;
      x = p + step;
      l = TMath::Abs(langaufun(&x,params) - fy);
 
      if (l > lold)
         step = -step/10;
 
      p += step;
   }

   if (i == MAXCALLS)
      return (-2);

   fxr = x;


   // Search for left x location of fy

   p = maxx - 0.5 * params[0];
   step = -params[0];
   lold = -2.0;
   l    = -1e300;
   i    = 0;

   while ( (l != lold) && (i < MAXCALLS) ) {
      i++;

      lold = l;
      x = p + step;
      l = TMath::Abs(langaufun(&x,params) - fy);
 
      if (l > lold)
         step = -step/10;
 
      p += step;
   }

   if (i == MAXCALLS)
      return (-3);


   fxl = x;

   FWHM = fxr - fxl;
   return (0);
}

void langaus() {
   // Fill Histogram
   Int_t data[100] = {0,0,0,0,0,0,2,6,11,18,18,55,90,141,255,323,454,563,681,
                    737,821,796,832,720,637,558,519,460,357,291,279,241,212,
                    153,164,139,106,95,91,76,80,80,59,58,51,30,49,23,35,28,23,
                    22,27,27,24,20,16,17,14,20,12,12,13,10,17,7,6,12,6,12,4,
                    9,9,10,3,4,5,2,4,1,5,5,1,7,1,6,3,3,3,4,5,4,4,2,2,7,2,4};
   TH1F *hSNR = new TH1F("snr","Signal-to-noise",400,0,400);

   for (Int_t i=0; i<100; i++) hSNR->Fill(i,data[i]);

   // Fitting SNR histo
   printf("Fitting...\n");

   // Setting fit range and start values
   Double_t fr[2];
   Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4];
   fr[0]=0.3*hSNR->GetMean();
   fr[1]=3.0*hSNR->GetMean();

   pllo[0]=0.5; pllo[1]=5.0; pllo[2]=1.0; pllo[3]=0.4;
   plhi[0]=5.0; plhi[1]=50.0; plhi[2]=1000000.0; plhi[3]=5.0;
   sv[0]=1.8; sv[1]=20.0; sv[2]=50000.0; sv[3]=3.0;

   Double_t chisqr;
   Int_t    ndf;
   TF1 *fitsnr = langaufit(hSNR,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);
   
   Double_t SNRPeak, SNRFWHM;
   langaupro(fp,SNRPeak,SNRFWHM);

   printf("Fitting done\nPlotting results...\n");

   // Global style settings
   gStyle->SetOptStat(1111);
   gStyle->SetOptFit(111);
   gStyle->SetLabelSize(0.03,"x");
   gStyle->SetLabelSize(0.03,"y");

   hSNR->GetXaxis()->SetRange(0,70);
   hSNR->Draw();
   fitsnr->Draw("lsame");
}

 langaus.C:1
 langaus.C:2
 langaus.C:3
 langaus.C:4
 langaus.C:5
 langaus.C:6
 langaus.C:7
 langaus.C:8
 langaus.C:9
 langaus.C:10
 langaus.C:11
 langaus.C:12
 langaus.C:13
 langaus.C:14
 langaus.C:15
 langaus.C:16
 langaus.C:17
 langaus.C:18
 langaus.C:19
 langaus.C:20
 langaus.C:21
 langaus.C:22
 langaus.C:23
 langaus.C:24
 langaus.C:25
 langaus.C:26
 langaus.C:27
 langaus.C:28
 langaus.C:29
 langaus.C:30
 langaus.C:31
 langaus.C:32
 langaus.C:33
 langaus.C:34
 langaus.C:35
 langaus.C:36
 langaus.C:37
 langaus.C:38
 langaus.C:39
 langaus.C:40
 langaus.C:41
 langaus.C:42
 langaus.C:43
 langaus.C:44
 langaus.C:45
 langaus.C:46
 langaus.C:47
 langaus.C:48
 langaus.C:49
 langaus.C:50
 langaus.C:51
 langaus.C:52
 langaus.C:53
 langaus.C:54
 langaus.C:55
 langaus.C:56
 langaus.C:57
 langaus.C:58
 langaus.C:59
 langaus.C:60
 langaus.C:61
 langaus.C:62
 langaus.C:63
 langaus.C:64
 langaus.C:65
 langaus.C:66
 langaus.C:67
 langaus.C:68
 langaus.C:69
 langaus.C:70
 langaus.C:71
 langaus.C:72
 langaus.C:73
 langaus.C:74
 langaus.C:75
 langaus.C:76
 langaus.C:77
 langaus.C:78
 langaus.C:79
 langaus.C:80
 langaus.C:81
 langaus.C:82
 langaus.C:83
 langaus.C:84
 langaus.C:85
 langaus.C:86
 langaus.C:87
 langaus.C:88
 langaus.C:89
 langaus.C:90
 langaus.C:91
 langaus.C:92
 langaus.C:93
 langaus.C:94
 langaus.C:95
 langaus.C:96
 langaus.C:97
 langaus.C:98
 langaus.C:99
 langaus.C:100
 langaus.C:101
 langaus.C:102
 langaus.C:103
 langaus.C:104
 langaus.C:105
 langaus.C:106
 langaus.C:107
 langaus.C:108
 langaus.C:109
 langaus.C:110
 langaus.C:111
 langaus.C:112
 langaus.C:113
 langaus.C:114
 langaus.C:115
 langaus.C:116
 langaus.C:117
 langaus.C:118
 langaus.C:119
 langaus.C:120
 langaus.C:121
 langaus.C:122
 langaus.C:123
 langaus.C:124
 langaus.C:125
 langaus.C:126
 langaus.C:127
 langaus.C:128
 langaus.C:129
 langaus.C:130
 langaus.C:131
 langaus.C:132
 langaus.C:133
 langaus.C:134
 langaus.C:135
 langaus.C:136
 langaus.C:137
 langaus.C:138
 langaus.C:139
 langaus.C:140
 langaus.C:141
 langaus.C:142
 langaus.C:143
 langaus.C:144
 langaus.C:145
 langaus.C:146
 langaus.C:147
 langaus.C:148
 langaus.C:149
 langaus.C:150
 langaus.C:151
 langaus.C:152
 langaus.C:153
 langaus.C:154
 langaus.C:155
 langaus.C:156
 langaus.C:157
 langaus.C:158
 langaus.C:159
 langaus.C:160
 langaus.C:161
 langaus.C:162
 langaus.C:163
 langaus.C:164
 langaus.C:165
 langaus.C:166
 langaus.C:167
 langaus.C:168
 langaus.C:169
 langaus.C:170
 langaus.C:171
 langaus.C:172
 langaus.C:173
 langaus.C:174
 langaus.C:175
 langaus.C:176
 langaus.C:177
 langaus.C:178
 langaus.C:179
 langaus.C:180
 langaus.C:181
 langaus.C:182
 langaus.C:183
 langaus.C:184
 langaus.C:185
 langaus.C:186
 langaus.C:187
 langaus.C:188
 langaus.C:189
 langaus.C:190
 langaus.C:191
 langaus.C:192
 langaus.C:193
 langaus.C:194
 langaus.C:195
 langaus.C:196
 langaus.C:197
 langaus.C:198
 langaus.C:199
 langaus.C:200
 langaus.C:201
 langaus.C:202
 langaus.C:203
 langaus.C:204
 langaus.C:205
 langaus.C:206
 langaus.C:207
 langaus.C:208
 langaus.C:209
 langaus.C:210
 langaus.C:211
 langaus.C:212
 langaus.C:213
 langaus.C:214
 langaus.C:215
 langaus.C:216
 langaus.C:217
 langaus.C:218
 langaus.C:219
 langaus.C:220
 langaus.C:221
 langaus.C:222
 langaus.C:223
 langaus.C:224
 langaus.C:225
 langaus.C:226
 langaus.C:227
 langaus.C:228
 langaus.C:229
 langaus.C:230
 langaus.C:231
 langaus.C:232
 langaus.C:233
 langaus.C:234
 langaus.C:235
 langaus.C:236
 langaus.C:237
 langaus.C:238
 langaus.C:239
 langaus.C:240
 langaus.C:241
 langaus.C:242
 langaus.C:243
 langaus.C:244
 langaus.C:245
 langaus.C:246
 langaus.C:247
 langaus.C:248
 langaus.C:249
 langaus.C:250
 langaus.C:251
 langaus.C:252
 langaus.C:253
 langaus.C:254
 langaus.C:255
 langaus.C:256
 langaus.C:257
 langaus.C:258
 langaus.C:259
 langaus.C:260
 langaus.C:261
 langaus.C:262
 langaus.C:263
 langaus.C:264
 langaus.C:265
 langaus.C:266
 langaus.C:267
 langaus.C:268
 langaus.C:269
 langaus.C:270
 langaus.C:271
 langaus.C:272
 langaus.C:273
 langaus.C:274