[ROOT] Breit-Wigner convoluted with a gaussian

Discussion:

Manuel Diaz Gomez

2005-04-29 13:24:48 UTC

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

Rene Brun

2005-04-30 05:04:32 UTC

Permalink

Hi Manuel,

Change the convolution range from .4 sigma to 4 (as in langaus.C)
ie:

Double_t sc = 4; // convolution extends to +-sc Gaussian sigmas

Rene Brun

On Fri,

Post by Manuel Diaz Gomez
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)
/*--------------------------------------------------------------------*/
{
//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

Jan Friedrich

2005-05-02 07:45:45 UTC

Permalink

Hi Manuel,

if you need that specific function (and don't want to do numerical
integration on your own), TMath::Voigt() may be an option for you.

Greetings Jan

Post by Rene Brun
Hi Manuel,
Change the convolution range from .4 sigma to 4 (as in langaus.C)
Double_t sc = 4; // convolution extends to +-sc Gaussian sigmas
Rene Brun
On Fri,

Manuel Diaz Gomez

2005-05-02 09:28:08 UTC

Permalink

Thanks Rene,

Changing the convolution range solved the problem.

Cheers
Manuel

-----Original Message-----
From: Rene Brun [mailto:***@pcroot.cern.ch]
Sent: Sat 4/30/2005 7:04 AM
To: Manuel Diaz Gomez
Cc: ***@pcroot.cern.ch
Subject: Re: [ROOT] Breit-Wigner convoluted with a gaussian
Hi Manuel,

Change the convolution range from .4 sigma to 4 (as in langaus.C)
ie:

Double_t sc = 4; // convolution extends to +-sc Gaussian sigmas

Rene Brun

On Fri,