HDQMUtil Namespace Reference

Functions
double	Gauss (double *x, double *par)
double	langaufun (double *x, double *par)
int32_t	langaupro (double *params, double &maxx, double &FWHM)

Detailed Description

@ class HDQMfitUtilities @ fit Landau distributions to historic monitoring elements @ fits from Susy's analysis (DQM/SiStripHistoricInfoClient/test/TrendsWithFits)

Function Documentation

double HDQMUtil::Gauss	(	double *	x,
		double *	par
	)

Definition at line 159 of file HDQMfitUtilities.cc.

Referenced by HDQMfitUtilities::doGaussFit().

                                       {
    // The noise function: a gaussian

    double arg = 0;
    if (par[2])
      arg = (x[0] - par[1]) / par[2];

    double noise = par[0] * TMath::Exp(-0.5 * arg * arg);
    return noise;
  }

double HDQMUtil::langaufun	(	double *	x,
		double *	par
	)

Definition at line 5 of file HDQMfitUtilities.cc.

References mps_fire::i, np, SimDataFormats::CaloAnalysis::sc, and geometryCSVtoXML::xx.

Referenced by HDQMfitUtilities::doLanGaussFit(), and langaupro().

                                           {
    // 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 invsq2pi = 0.3989422804014;  // (2 pi)^(-1/2)
    double mpshift = -0.22278298;       // Landau maximum location

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

    // Variables
    double xx;
    double mpc;
    double fland;
    double sum = 0.0;
    double xlow, xupp;
    double step;
    double 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;

    // Landau Distribution Production
    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]);
  }

int32_t HDQMUtil::langaupro	(	double *	params,
		double &	maxx,
		double &	FWHM
	)

Definition at line 58 of file HDQMfitUtilities.cc.

References Abs(), MillePedeFileConverter_cfg::e, mps_fire::i, checklumidiff::l, langaufun(), and AlCaHLTBitMon_ParallelJobs::p.

Referenced by MuonTestSummary::doEnergyTests(), and HDQMfitUtilities::doLanGaussFit().

                                                                {
    edm::LogInfo("fitUtility") << "inside langaupro " << std::endl;
    // Seaches for the location (x value) at the maximum of the
    // Landau and its full width at half-maximum.
    //
    // The search is probably not very efficient, but it's a first try.

    double p, x, fy, fxr, fxl;
    double step;
    double l, lold, dl;
    int32_t i = 0;
    const int32_t MAXCALLS = 10000;
    const double dlStop = 1e-3;  // relative change < .001

    // Search for maximum
    p = params[1] - 0.1 * params[0];
    step = 0.05 * params[0];
    lold = -2.0;
    l = -1.0;

    dl = (l - lold) / lold;  // FIXME catch divide by zero
    while ((TMath::Abs(dl) > dlStop) && (i < MAXCALLS)) {
      i++;

      lold = l;
      x = p + step;
      l = langaufun(&x, params);
      dl = (l - lold) / lold;  // FIXME catch divide by zero

      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;

    dl = (l - lold) / lold;  // FIXME catch divide by zero
    while ((TMath::Abs(dl) > dlStop) && (i < MAXCALLS)) {
      i++;

      lold = l;
      x = p + step;
      l = TMath::Abs(langaufun(&x, params) - fy);
      dl = (l - lold) / lold;  // FIXME catch divide by zero

      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;

    dl = (l - lold) / lold;  // FIXME catch divide by zero
    while ((TMath::Abs(dl) > dlStop) && (i < MAXCALLS)) {
      i++;

      lold = l;
      x = p + step;
      l = TMath::Abs(langaufun(&x, params) - fy);
      dl = (l - lold) / lold;  // FIXME catch divide by zero

      if (l > lold)
        step = -step / 10;

      p += step;
    }

    if (i == MAXCALLS)
      return (-3);

    fxl = x;

    FWHM = fxr - fxl;
    return (0);
  }