/data/refman/pasoursint/CMSSW_6_1_1/src/DQM/SiStripHistoricInfoClient/interface/fitUtilities.h File Reference

#include <vector>
#include <cstring>
#include <iostream>
#include <sstream>
#include <string>
#include "TH1.h"
#include "TF1.h"
#include "TMath.h"
#include "DQMServices/Core/interface/MonitorElement.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

Go to the source code of this file.

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

---

Function Documentation

double Gauss	(	double *	x,
		double *	par
	)

Definition at line 225 of file fitUtilities.h.

References cmsCodeRulesChecker::arg.

Referenced by fitUtilities::doGaussFit(), and StabilizedGauss::StabilizedGauss().

{
  // 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 langaufun	(	double *	x,
		double *	par
	)

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

Definition at line 66 of file fitUtilities.h.

References i, np, and relval_parameters_module::step.

                                        {
  //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 langaupro	(	double *	params,
		double &	maxx,
		double &	FWHM
	)

Definition at line 120 of file fitUtilities.h.

References alignCSCRings::e, i, prof2calltree::l, langaufun(), AlCaHLTBitMon_ParallelJobs::p, relval_parameters_module::step, and x.

Referenced by MuonTestSummary::doEnergyTests(), HDQMfitUtilities::doLanGaussFit(), and fitUtilities::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);
}