# Langauss with a landau(-x)

Dear All,

I would like to have the cumulative distribution of a convolution of a Landau(-x) with a gaussian, which I want to use later for a fit. Following the langaus root example, I added another loop around the convolution loop, to do the integration to obtain the cumulative distribution.

My first problem is, that whatever I try to make the landau(-x) it crashes. Second of all, even though I tried making the cumulative distribution, it still looks like the landau*gauss.

The code can be found here:

``````double lanGausCummulatif(double *x, double *par) {

// Fit parameters:
// par[0] = Width (scale) parameter of Landau density
// par[1] = Most Probable Value 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 invSqrd2pi = 0.3989422804014;   // (2 pi)^(-1/2)
double mpshift    = -0.22278298;       // Landau maximum location

// Control constants
double nCsteps     = 100.0;      // number of convolution steps
double sigmasGauss =   5.0;      // convolution extends to +-sc Gaussian sigmas

// Variables
double xx;
double mpc;
double fLandau;
double sum = 0.0;
double step;

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

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

step = (xupp-xlow) / nCsteps;

double integral = 0.0;

// Convolution integral of Landau and Gaussian by sum
for ( double iVcth = 0.0; iVcth < x[0]; iVcth++ ) {
double differential = 0.0;
for( double iConv = 1.0; iConv <= nCsteps/2; iConv++ ) {

xx = xlow + (iConv-.5) * step;
fLandau = TMath::Landau(xx, mpc, par[0]) / par[0];
//sum += fLandau * TMath::Gaus(x[0], xx, par[3]);
differential += fLandau * TMath::Gaus(iVcth, xx, par[3]);

xx = xupp - (iConv-.5) * step;
fLandau = TMath::Landau(xx, mpc, par[0]) / par[0];
//sum += fLandau * TMath::Gaus(x[0], xx, par[3]);
differential += fLandau * TMath::Gaus(iVcth, xx, par[3]);
}
integral += (par[2] * step * differential * invSqrd2pi / par[3]);
//std::cout << integral << std::endl;
}

//return - (par[2] * step * sum * invSqrd2pi / par[3]);
return integral;

}

TF1 * lanGaus() {
// Currently using a Landau in Vcth units, but could use fC and then scale the PreampShaper output
TF1* f = new TF1("f", "lanGausCummulatif", -600, 600, 4);
// Initial parameter values to be checked!
// I set Landau "width" to 9 (is this the same as sigma?)
// MVP is 144
// Total area of the convolution is 1 for now, this is a normalization factor
// Sigma of the gaussian is 6 (CBC3 noise is 6 Vcth units??? To be checked!!!)
f->SetNpx(5000);
f->SetParameters(9, 440, 1, 6);
f->SetParNames("Width","MP","Area","GSigma");
TCanvas *canvas = new TCanvas();
f->Draw();

return f;
}``````

Could someone please give me some hints on this issue?

Nikkie

I managed to get the integration part working again. I’m almost there, just the fact that the Landau’s tail is not on the right side.

I thought that I could fix the tail by making the width negative (-9 in my case).
However, after the fit the width is positive (~39) and it looks like the landau tail is still on the wrong side…

``````#include "./../src/Template.cc"

int runNos[32]       = { 64, 63, 62, 61, 60, 59, 56, 55, 54, 53, 52, 51, 74, 49, 48, 47, 45, 43, 42, 73, 317,318,320,321,322,323,324,326,327,328,329,330 };
double thrNos[32]    = { 380,390,400,410,420,430,440,450,460,470,480,490,500,510,520,530,540,550,560,570,370,360,350,340,330,320,310,300,280,260,240,220 };
double thrNosErr[32] = { 1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1   };
//380,390,400,410,420,430,440,450,460,470,480,490,500,510,520,530,540,550,560,570,370,360,350,340,330,320,310,300,280,260,240,220
//64,63,62,61,60,59,56,55,54,53,52,51,74,49,48,47,45,43,42,73,317,318,320,321,322,323,324,326,327,328,329,330

double lanGausCummulatif(double *x, double *par) {

// Fit parameters:
// par[0] = Width (scale) parameter of Landau density
// par[1] = Most Probable Value 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 invSqrd2pi = 0.3989422804014;   // (2 pi)^(-1/2)
double mpshift    = -0.22278298;       // Landau maximum location

// Control constants
double nCsteps     = 100.0;      // number of convolution steps
double sigmasGauss =   5.0;      // convolution extends to +-sc Gaussian sigmas

// Variables
double xx;
double mpc;
double fLandau;
double sum = 0.0;
double step;

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

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

step = (xupp-xlow) / nCsteps;

double integral = 0.0;

// Convolution integral of Landau and Gaussian by sum
for ( double iVcth = 0.0; iVcth < x[0]; iVcth++ ) {
double differential = 0.0;
for( double iConv = 1.0; iConv <= nCsteps/2; iConv++ ) {

xx = (xlow + (iConv-.5) * step);
fLandau = TMath::Landau(xx, mpc, par[0]) / par[0];
//sum += fLandau * TMath::Gaus(x[0], xx, par[3]);
differential += fLandau * TMath::Gaus(iVcth, xx, par[3]);

xx = (xupp - (iConv-.5) * step);
fLandau = TMath::Landau(xx, mpc, par[0]) / par[0];
//sum += fLandau * TMath::Gaus(x[0], xx, par[3]);
differential += fLandau * TMath::Gaus(iVcth, xx, par[3]);
}
integral += (par[2] * step * differential * invSqrd2pi / par[3]);
//std::cout << integral << std::endl;
}

//return - (par[2] * step * sum * invSqrd2pi / par[3]);
return integral;

}

TF1 * lanGaus() {

TF1* f = new TF1("f", "lanGausCummulatif", 200, 600, 4);
// Initial parameter values to be checked!
// I set Landau "width" to 9 (is this the same as sigma?)
// MVP is 144
// Total area of the convolution is 1 for now, this is a normalization factor
// Sigma of the gaussian is 6 (CBC3 noise is 6 Vcth units??? To be checked!!!)
f->SetParameters(-9, 440, 1, 6);
f->SetParNames("Width","MP","Area","GSigma");
//f->FixParameter(1, 440);
//f->SetParLimits(0, -20, 0);
f->SetParLimits(1, 420, 460);
//f->SetParLimits(3, 1, 12);

f->SetNpx(500);
TCanvas *canv = new TCanvas();
f->Draw();

return f;
}

void runner() {

//test(141);

double thisEfficiency[32] = {0};
double thisEfficiencyError[32] = {0};

//  subtactResidualBackground_Cluster (pRunNumber, 0, -1, 1, 2.0);

for (unsigned int iRun = 0; iRun < 32; iRun++) {

thisEfficiency [iRun] = subtactResidualBackground_Cluster ( runNos[iRun], 0, 1, 1, 1, 2.0 );
thisEfficiencyError[iRun] = subtactResidualBackground_ClusterError ( runNos[iRun], 0, 1, 1, 1, 2.0 );

}//end for iRun

TCanvas *c0 = new TCanvas();

TGraphErrors *eff = new TGraphErrors ( 32, thrNos, thisEfficiency, thrNosErr, thisEfficiencyError );
//eff->SetMarkerStyle(21);
eff->SetMarkerColor(2);
eff->GetXaxis()->SetTitle("threshold [Vcth]");
eff->GetYaxis()->SetTitle("efficiency [%]");
eff->Draw("ap");

TF1 *fitFunction = lanGaus();
eff->Fit(fitFunction, "R");
c0->SaveAs("langausFitTest.pdf");

}``````

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.