Breit-Wigner convoluted with a gaussian

From: Manuel Diaz Gomez <Manuel.Maria.Diaz.Gomez_at_cern.ch>
Date: Fri, 29 Apr 2005 15:24:48 +0200

Hi All,

I modified the langaus.C example (under $ROOTSYS/tutorials) so that it does a breit-wigner convoluted with a gaussian (see code bellow). The problem is that I happen to be missing a normalisation factor somewhere, but I just can't see it. Perhaps you can spot it right away.

The following macro defines a gaussian convoluted breit-wigner function, fills a histogram out of it with 1000 entries, and the does the fit. I would expect the par[2] variable to be ~1000 but...rather I get ~3K. I would appretiate any help!

/*--------------------------------------------------------------------*/

Double_t breitgausfun(Double_t *x, Double_t *par)

/*--------------------------------------------------------------------*/

{

//Fit parameters:
//par[0]=Width (scale) Breit-Wigner
//par[1]=Most Probable (MP, location) Breit mean
//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 twoPi = 6.2831853071795;//2Pi

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

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

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

      step = (xupp-xlow) / np;

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

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

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

}

/*--------------------------------------------*/

Double_t Genbreitgaus()

/*--------------------------------------------*/

{
  TH1F *brgauss = new TH1F("breitg","", 131, 0, 130);   TF1 *f = new TF1("f",breitgausfun, 0, 130 ,4);

  Double_t par[4];

  par[0] = 2.495;
  par[1] = 80.0;
  par[2] = 1000.0;
  par[3] = 10.0;
  
  f->SetParameters(par);
  f->FixParameter(0, 2.495);

  f->FixParameter(1, 80.0);
  f->FixParameter(3, 10.0);      

  brgauss->FillRandom("f", 1000);
  brgauss->Fit(f, "RBO");  

}

//////////////////////////////////////////////////////////////////////

Cheers
Manuel Received on Fri Apr 29 2005 - 15:24:54 MEST