norma2.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 norma2() {

  Int_t nummips = 6;

   TF1* fitFunc = NULL;
   TF1 *fitsnr = NULL;
   TF1* fitsnr2 = NULL;
   Double_t sigRange = 2.0;

   TCanvas *c1 = new TCanvas("c1","c1",10, 100, 900,700);
   TCanvas *c2 = new TCanvas("c2","c2",1000, 100, 900,700);

   TRandom2 r(0);
   gRandom->SetSeed();

   Double_t l = 30.0;
   Double_t aa = 0.4;
   Double_t min = 0.0;
   Double_t max = 40.0;
   Int_t bin = 100;
   Double_t mpv = 1.0;
   Double_t mpv_a = 1.0;
   Double_t a = 1.0;
   Double_t width_a = 0.05;
   width_a = a * width_a;
   Double_t width = 0.05;

   Double_t fr[2] = { 0.0, 0.0 };
   fr[0] = min;
   fr[1] = max;
   Double_t sv[4] = { 0.0, 0.0, 0.0, 0.0 };
   sv[0] = width;
   sv[1] = mpv;
   sv[2] = 1000;
   sv[3] = width;
   Double_t pllo[4] = { 0, 0, 0, 0 };
   Double_t plhi[4] = { 0, 0, 0, 0 };
   pllo[0] = 0.0; pllo[1] = 0.0; pllo[2] = 1.0; pllo[3] = 0.0;
   plhi[0] = 1.0; plhi[1] = 1000000.0; plhi[2] = 1000000.0; plhi[3] = 1.0;
   Double_t fp[4] = { 0.0, 0.0, 0.0, 0.0 };
   Double_t fpe[4] = { 0.0, 0.0, 0.0, 0.0 };
   Double_t chisqr = 0.0;
   Int_t ndf = 0;
   

   TH1F *h = new TH1F( "h", "h", bin, min, max );
   TH1F *h2 = new TH1F( "h2", "h2", bin, min, max );

   fitFunc = new TF1("fitFunc", "landau", mpv - sigRange*width, mpv + sigRange*width);

   for (Int_t i = 0; i < 1000; i++) { 
     Double_t ltemp1 = 0.0;
     Double_t ltemp2 = 0.0;
     for (Int_t jj = 0; jj < nummips; jj++)  {

     Double_t l = 0;
     Double_t l2 = 0;
     l = gRandom -> Landau( mpv, width );
     Double_t xx = r.Uniform( l );
     Double_t yy = r.Uniform( l );

     Double_t tt1 = 2.0 * aa * aa * l * l;
     Double_t tt2 = ( xx - ( 0.5 * l )  ) * ( xx - ( 0.5 * l ) );
     Double_t tt3 = ( yy - ( 0.5 * l )  ) * ( yy - ( 0.5 * l ) );
     Double_t cc = exp( -1.0 * ( ( tt2 / tt1 ) + ( tt3 / tt1 ) ) ); // / mm;

        ltemp1 = ltemp1 + l;
        l2 = l * cc;
        ltemp2 = ltemp2 + l2;

//     l2 = gRandom -> Landau( mpv_a, width_a );
//     if ( l > 5 ) l = 0;
//     if ( ( l2 > 10 ) ) l2 = 0;
}
     h -> Fill( ltemp1 );
     h2 -> Fill( ltemp2 );
   
   }

   Double_t fp2[4] = { 0.0, 0.0, 0.0, 0.0 };
   Double_t fpe2[4] = { 0.0, 0.0, 0.0, 0.0 };
   Double_t chisqr2 = 0.0;
   Int_t ndf2 = 0;
  
   fitsnr2 = langaufit(h,fr,sv,pllo,plhi,fp2,fpe2,&chisqr2,&ndf2);

   Double_t SNRPeak2 = 0;
   Double_t SNRFWHM2 = 0;
   langaupro(fp2,SNRPeak2,SNRFWHM2);

//   h -> Fit( fitFunc );

/*   mpv_fit = fitFunc->GetParameter(1);
   sigma_fit = fitFunc->GetParameter(2);
   Double_t chi2_fit = 0;
   chi2_fit = fitFunc->GetChisquare() / fitFunc->GetNDF();  
 
   cout << " mpv_fit === " << mpv_fit << endl;
   cout << " sigma_fit === " << sigma_fit << endl;
   cout << " fitFunc->GetChisquare() === " << fitFunc->GetChisquare() << endl;
   cout << " fitFunc->GetNDF() === " << fitFunc->GetNDF() << endl;
   cout << " chi2 === " << chi2_fit << endl;
*/
   fitsnr = langaufit(h2,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);
   std::cout << " chisqr === " << chisqr << std::endl;
   std::cout << " ndf === " << ndf << std::endl;

   Double_t SNRPeak = 0;
   Double_t SNRFWHM = 0;
   langaupro(fp,SNRPeak,SNRFWHM);
   std::cout << " SNRPeak / SNRPeak2   " << SNRPeak << " / " << SNRPeak2 << " = " << SNRPeak / SNRPeak2 << std::endl;

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

   c1->cd();
   h -> Draw();
   fitsnr2 -> SetLineColor(kGreen);
   fitsnr2 -> Draw("same");
   c2->cd();
   gStyle->SetOptStat(1111);
   gStyle->SetOptFit(111);
   h2 -> Draw();
   fitsnr -> SetLineColor(kBlue);
   fitsnr -> Draw("same");

}

/*

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");
}
*/