#include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsFitter.h"
#include <fstream>
#include "TMath.h"

Functions
double	MuonResidualsFitter_compute_log_convolution (double toversigma, double gammaoversigma, double max, double step, double power)
Double_t	MuonResidualsFitter_GaussPowerTails_TF1 (Double_t xvec, Double_t par)
double	MuonResidualsFitter_integrate_pureGaussian (double low, double high, double center, double sigma)
double	MuonResidualsFitter_logGaussPowerTails (double residual, double center, double sigma)
double	MuonResidualsFitter_logPowerLawTails (double residual, double center, double sigma, double gamma)
double	MuonResidualsFitter_logPureGaussian (double residual, double center, double sigma)
double	MuonResidualsFitter_logROOTVoigt (double residual, double center, double sigma, double gamma)
Double_t	MuonResidualsFitter_powerLawTails_TF1 (Double_t xvec, Double_t par)
Double_t	MuonResidualsFitter_pureGaussian_TF1 (Double_t xvec, Double_t par)
Double_t	MuonResidualsFitter_ROOTVoigt_TF1 (Double_t xvec, Double_t par)
Variables
const double	MuonResidualsFitter_gsbinsize = 0.01
double	MuonResidualsFitter_lookup_table [MuonResidualsFitter_numgsbins][MuonResidualsFitter_numtsbins]
const int	MuonResidualsFitter_numgsbins = 500
const int	MuonResidualsFitter_numtsbins = 500
bool	MuonResidualsFitter_table_initialized = false
static TMinuit *	MuonResidualsFitter_TMinuit
const double	MuonResidualsFitter_tsbinsize = 0.1

Function Documentation

double MuonResidualsFitter_compute_log_convolution	(	double	toversigma,
		double	gammaoversigma,
		double	max,
		double	step,
		double	power
	)

Definition at line 30 of file MuonResidualsFitter.cc.

References funct::exp(), funct::log(), M_PI, max(), funct::pow(), mathSSE::sqrt(), and launcher::step.

Referenced by MuonResidualsFitter::initialize_table().

                                                                                                                                    {
  if (gammaoversigma == 0.) return (-toversigma*toversigma/2.) - log(sqrt(2*M_PI));

  double sum = 0.;
  double uplus = 0.;
  double integrandplus_last = 0.;
  double integrandminus_last = 0.;
  for (double inc = 0.;  uplus < max;  inc += step) {
    double uplus_last = uplus;
    uplus = pow(inc, power);

    double integrandplus = exp(-pow(uplus - toversigma, 2)/2.) / (uplus*uplus/gammaoversigma + gammaoversigma) / M_PI / sqrt(2.*M_PI);
    double integrandminus = exp(-pow(-uplus - toversigma, 2)/2.) / (uplus*uplus/gammaoversigma + gammaoversigma) / M_PI / sqrt(2.*M_PI);

    sum += integrandplus * (uplus - uplus_last);
    sum += 0.5 * fabs(integrandplus_last - integrandplus) * (uplus - uplus_last);

    sum += integrandminus * (uplus - uplus_last);
    sum += 0.5 * fabs(integrandminus_last - integrandminus) * (uplus - uplus_last);

    integrandplus_last = integrandplus;
    integrandminus_last = integrandminus;
  }
  return log(sum);
}

Double_t MuonResidualsFitter_GaussPowerTails_TF1	(	Double_t *	xvec,
		Double_t *	par
	)

Definition at line 112 of file MuonResidualsFitter.cc.

References funct::exp(), and MuonResidualsFitter_logGaussPowerTails().

Referenced by MuonResiduals6DOFrphiFitter::plot(), MuonResiduals6DOFFitter::plot(), MuonResidualsBfieldAngleFitter::plot(), MuonResiduals5DOFFitter::plot(), MuonResidualsPositionFitter::plot(), MuonResiduals1DOFFitter::plot(), and MuonResidualsAngleFitter::plot().

                                                                                {
  return par[0] * exp(MuonResidualsFitter_logGaussPowerTails(xvec[0], par[1], par[2]));
}

double MuonResidualsFitter_integrate_pureGaussian	(	double	low,
		double	high,
		double	center,
		double	sigma
	)

Definition at line 116 of file MuonResidualsFitter.cc.

References funct::exp(), and mathSSE::sqrt().

                                                                                                        {
  return (erf((high + center) / sqrt(2.) / sigma) - erf((low + center) / sqrt(2.) / sigma)) * exp(0.5/sigma/sigma) / 2.;
}

double MuonResidualsFitter_logGaussPowerTails	(	double	residual,
		double	center,
		double	sigma
	)

Definition at line 101 of file MuonResidualsFitter.cc.

References a, funct::exp(), funct::log(), m, M_PI, n, funct::pow(), asciidump::s, mathSSE::sqrt(), and x.

Referenced by MuonResiduals1DOFFitter_FCN(), MuonResiduals5DOFFitter_FCN(), MuonResiduals6DOFFitter_FCN(), MuonResiduals6DOFrphiFitter_FCN(), MuonResidualsAngleFitter_FCN(), MuonResidualsBfieldAngleFitter_FCN(), MuonResidualsFitter_GaussPowerTails_TF1(), and MuonResidualsPositionFitter_FCN().

                                                                                            {
  double x = residual-center;
  double s = fabs(sigma);
  double m = 2*s;
  double a = pow(m,4)*exp(-2);
  double n = sqrt(2*M_PI)*s*erf(sqrt(2))+2*a*pow(m,-3)/3;

  if (fabs(x)<m) return -x*x/(2*s*s) - log(n);
  else return log(a) -4*log(fabs(x)) - log(n);
}

double MuonResidualsFitter_logPowerLawTails	(	double	residual,
		double	center,
		double	sigma,
		double	gamma
	)

Definition at line 57 of file MuonResidualsFitter.cc.

References funct::log(), M_PI, MuonResidualsFitter_gsbinsize, MuonResidualsFitter_lookup_table, MuonResidualsFitter_numgsbins, MuonResidualsFitter_numtsbins, and MuonResidualsFitter_tsbinsize.

Referenced by MuonResiduals1DOFFitter_FCN(), MuonResiduals5DOFFitter_FCN(), MuonResiduals6DOFFitter_FCN(), MuonResiduals6DOFrphiFitter_FCN(), MuonResidualsAngleFitter_FCN(), MuonResidualsBfieldAngleFitter_FCN(), MuonResidualsFitter_powerLawTails_TF1(), and MuonResidualsPositionFitter_FCN().

                                                                                                        {
  sigma = fabs(sigma);
  double gammaoversigma = fabs(gamma / sigma);
  double toversigma = fabs((residual - center) / sigma);

  int gsbin0 = int(floor(gammaoversigma / MuonResidualsFitter_gsbinsize));
  int gsbin1 = int(ceil(gammaoversigma / MuonResidualsFitter_gsbinsize));
  int tsbin0 = int(floor(toversigma / MuonResidualsFitter_tsbinsize));
  int tsbin1 = int(ceil(toversigma / MuonResidualsFitter_tsbinsize));

  bool gsisbad = (gsbin0 >= MuonResidualsFitter_numgsbins  ||  gsbin1 >= MuonResidualsFitter_numgsbins);
  bool tsisbad = (tsbin0 >= MuonResidualsFitter_numtsbins  ||  tsbin1 >= MuonResidualsFitter_numtsbins);

  if (gsisbad  ||  tsisbad) {
    return log(gammaoversigma/M_PI) - log(toversigma*toversigma + gammaoversigma*gammaoversigma) - log(sigma);
  }
  else {
    double val00 = MuonResidualsFitter_lookup_table[gsbin0][tsbin0];
    double val01 = MuonResidualsFitter_lookup_table[gsbin0][tsbin1];
    double val10 = MuonResidualsFitter_lookup_table[gsbin1][tsbin0];
    double val11 = MuonResidualsFitter_lookup_table[gsbin1][tsbin1];

    double val0 = val00 + ((toversigma / MuonResidualsFitter_tsbinsize) - tsbin0) * (val01 - val00);
    double val1 = val10 + ((toversigma / MuonResidualsFitter_tsbinsize) - tsbin0) * (val11 - val10);
    
    double val = val0 + ((gammaoversigma / MuonResidualsFitter_gsbinsize) - gsbin0) * (val1 - val0);
    
    return val - log(sigma);
  }
}

double MuonResidualsFitter_logPureGaussian	(	double	residual,
		double	center,
		double	sigma
	)

Definition at line 20 of file MuonResidualsFitter.cc.

References funct::log(), M_PI, funct::pow(), and mathSSE::sqrt().

Referenced by MuonResiduals1DOFFitter_FCN(), MuonResiduals5DOFFitter_FCN(), MuonResiduals6DOFFitter_FCN(), MuonResiduals6DOFrphiFitter_FCN(), MuonResidualsAngleFitter_FCN(), MuonResidualsBfieldAngleFitter_FCN(), MuonResidualsFitter_pureGaussian_TF1(), and MuonResidualsPositionFitter_FCN().

                                                                                         {
  sigma = fabs(sigma);
  return (-pow(residual - center, 2) / 2. / sigma / sigma) + log(1. / sqrt(2.*M_PI) / sigma);
}

double MuonResidualsFitter_logROOTVoigt	(	double	residual,
		double	center,
		double	sigma,
		double	gamma
	)

Definition at line 93 of file MuonResidualsFitter.cc.

References funct::log().

Referenced by MuonResiduals1DOFFitter_FCN(), MuonResiduals5DOFFitter_FCN(), MuonResiduals6DOFFitter_FCN(), MuonResiduals6DOFrphiFitter_FCN(), MuonResidualsAngleFitter_FCN(), MuonResidualsBfieldAngleFitter_FCN(), and MuonResidualsPositionFitter_FCN().

                                                                                                    {
  return log(TMath::Voigt(residual - center, fabs(sigma), fabs(gamma)*2.));
}

Double_t MuonResidualsFitter_powerLawTails_TF1	(	Double_t *	xvec,
		Double_t *	par
	)

Definition at line 89 of file MuonResidualsFitter.cc.

References funct::exp(), and MuonResidualsFitter_logPowerLawTails().

Referenced by MuonResiduals6DOFrphiFitter::plot(), MuonResiduals6DOFFitter::plot(), MuonResidualsBfieldAngleFitter::plot(), MuonResiduals5DOFFitter::plot(), MuonResidualsPositionFitter::plot(), MuonResiduals1DOFFitter::plot(), and MuonResidualsAngleFitter::plot().

                                                                              {
  return par[0] * exp(MuonResidualsFitter_logPowerLawTails(xvec[0], par[1], par[2], par[3]));
}

Double_t MuonResidualsFitter_pureGaussian_TF1	(	Double_t *	xvec,
		Double_t *	par
	)

Definition at line 26 of file MuonResidualsFitter.cc.

References funct::exp(), and MuonResidualsFitter_logPureGaussian().

Referenced by MuonResiduals6DOFrphiFitter::plot(), MuonResiduals6DOFFitter::plot(), MuonResidualsBfieldAngleFitter::plot(), MuonResiduals5DOFFitter::plot(), MuonResidualsPositionFitter::plot(), MuonResiduals1DOFFitter::plot(), and MuonResidualsAngleFitter::plot().

                                                                             {
  return par[0] * exp(MuonResidualsFitter_logPureGaussian(xvec[0], par[1], par[2]));
}

Double_t MuonResidualsFitter_ROOTVoigt_TF1	(	Double_t *	xvec,
		Double_t *	par
	)

Definition at line 97 of file MuonResidualsFitter.cc.

Referenced by MuonResiduals6DOFrphiFitter::plot(), MuonResiduals6DOFFitter::plot(), MuonResidualsBfieldAngleFitter::plot(), MuonResiduals5DOFFitter::plot(), MuonResidualsPositionFitter::plot(), MuonResiduals1DOFFitter::plot(), and MuonResidualsAngleFitter::plot().

                                                                          {
  return par[0] * TMath::Voigt(xvec[0] - par[1], fabs(par[2]), fabs(par[3])*2.);
}