Langauss with a landau(-x)

[NikkieD]

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?

Thanks in advance,
Nikkie

[NikkieD]

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");

}