MuonResidualsFitter.h File Reference

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"
#include "Alignment/CommonAlignment/interface/Alignable.h"
#include "TMinuit.h"
#include "TH1F.h"
#include "TProfile.h"
#include "TF1.h"
#include "TMath.h"
#include "TRandom3.h"
#include "TMatrixDSym.h"
#include <cstdio>
#include <iostream>
#include <string>
#include <sstream>
#include <map>

Go to the source code of this file.

Classes
struct	MuonResidualsFitter::MuonAlignmentTreeRow
class	MuonResidualsFitter
class	MuonResidualsFitterFitInfo

Functions
void	MuonResidualsAngleFitter_FCN (int &npar, double *gin, double &fval, double *par, int iflag)
double	MuonResidualsFitter_compute_log_convolution (double toversigma, double gammaoversigma, double max=1000., double step=0.001, double power=4.)
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_logPureGaussian2D (double x, double y, double x0, double y0, double sx, double sy, double r)
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)
void	MuonResidualsPositionFitter_FCN (int &npar, double *gin, double &fval, double *par, int iflag)

Function Documentation

void MuonResidualsAngleFitter_FCN	(	int &	npar,
		double *	gin,
		double &	fval,
		double *	par,
		int	iflag
	)

Definition at line 9 of file MuonResidualsAngleFitter.cc.

References MuonResidualsFitterFitInfo::fitter(), MuonResidualsAngleFitter::kAngle, MuonResidualsAngleFitter::kGamma, MuonResidualsFitter::kGaussPowerTails, MuonResidualsFitter::kPowerLawTails, MuonResidualsFitter::kPureGaussian, MuonResidualsAngleFitter::kResidual, MuonResidualsFitter::kROOTVoigt, MuonResidualsAngleFitter::kSigma, MuonResidualsAngleFitter::kXAngle, MuonResidualsAngleFitter::kXControl, MuonResidualsAngleFitter::kYAngle, MuonResidualsAngleFitter::kYControl, MuonResidualsAngleFitter_TMinuit, MuonResidualsFitter_logGaussPowerTails(), MuonResidualsFitter_logPowerLawTails(), MuonResidualsFitter_logPureGaussian(), MuonResidualsFitter_logROOTVoigt(), MuonResidualsFitter::residuals_begin(), MuonResidualsFitter::residuals_end(), MuonResidualsFitter::residualsModel(), and base_cff::xangle.

Referenced by MuonResidualsAngleFitter::fit().

                                                                                                {
  MuonResidualsFitterFitInfo *fitinfo = (MuonResidualsFitterFitInfo*)(MuonResidualsAngleFitter_TMinuit->GetObjectFit());
  MuonResidualsFitter *fitter = fitinfo->fitter();

  fval = 0.;
  for (std::vector<double*>::const_iterator resiter = fitter->residuals_begin();  resiter != fitter->residuals_end();  ++resiter) {
    const double residual = (*resiter)[MuonResidualsAngleFitter::kResidual];
    const double xangle = (*resiter)[MuonResidualsAngleFitter::kXAngle];
    const double yangle = (*resiter)[MuonResidualsAngleFitter::kYAngle];
    
    double center = 0.;
    center += par[MuonResidualsAngleFitter::kAngle];
    center += par[MuonResidualsAngleFitter::kXControl] * xangle;
    center += par[MuonResidualsAngleFitter::kYControl] * yangle;
    
    if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian) {
      fval += -MuonResidualsFitter_logPureGaussian(residual, center, par[MuonResidualsAngleFitter::kSigma]);
    }
    else if (fitter->residualsModel() == MuonResidualsFitter::kPowerLawTails) {
      fval += -MuonResidualsFitter_logPowerLawTails(residual, center, par[MuonResidualsAngleFitter::kSigma], par[MuonResidualsAngleFitter::kGamma]);
    }
    else if (fitter->residualsModel() == MuonResidualsFitter::kROOTVoigt) {
      fval += -MuonResidualsFitter_logROOTVoigt(residual, center, par[MuonResidualsAngleFitter::kSigma], par[MuonResidualsAngleFitter::kGamma]);
    }
    else if (fitter->residualsModel() == MuonResidualsFitter::kGaussPowerTails) {
      fval += -MuonResidualsFitter_logGaussPowerTails(residual, center, par[MuonResidualsAngleFitter::kSigma]);
    }
    else { assert(false); }
  }
}

double MuonResidualsFitter_compute_log_convolution	(	double	toversigma,
		double	gammaoversigma,
		double	max = 1000.,
		double	step = 0.001,
		double	power = 4.
	)

Definition at line 57 of file MuonResidualsFitter.cc.

References JetChargeProducer_cfi::exp, cmsBatch::log, M_PI, SiStripPI::max, funct::pow(), and mathSSE::sqrt().

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 153 of file MuonResidualsFitter.cc.

References JetChargeProducer_cfi::exp, and MuonResidualsFitter_logGaussPowerTails().

Referenced by MuonResiduals1DOFFitter::plot(), MuonResidualsAngleFitter::plot(), MuonResidualsBfieldAngleFitter::plot(), MuonResidualsPositionFitter::plot(), MuonResiduals5DOFFitter::plot(), MuonResiduals6DOFrphiFitter::plot(), and MuonResiduals6DOFFitter::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 159 of file MuonResidualsFitter.cc.

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

{
  static const double isqr2 = 1./sqrt(2.);
  return (erf((high + center) * isqr2 / sigma) - erf((low + center) * isqr2 / sigma)) * exp(0.5/sigma/sigma) * 0.5;
}

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

Definition at line 140 of file MuonResidualsFitter.cc.

References a, JetChargeProducer_cfi::exp, cmsBatch::log, funct::m, M_PI, gen::n, funct::pow(), alignCSCRings::s, mathSSE::sqrt(), and vertices_cff::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;
  static const double e2 = exp(-2.), sqerf = sqrt(2.*M_PI)*erf(sqrt(2.));
  double a = pow(m,4)*e2;
  double n = sqerf*s + 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 86 of file MuonResidualsFitter.cc.

References createfilelist::int, cmsBatch::log, M_PI, MuonResidualsFitter_gsbinsize, MuonResidualsFitter_lookup_table, MuonResidualsFitter_numgsbins, MuonResidualsFitter_numtsbins, MuonResidualsFitter_tsbinsize, and heppy_batch::val.

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 31 of file MuonResidualsFitter.cc.

References cmsBatch::log, M_PI, and funct::pow().

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

{
  sigma = fabs(sigma);
  static const double cgaus = 0.5 * log( 2.*M_PI );
  return (-pow(residual - center, 2) *0.5 / sigma / sigma) - cgaus - log(sigma);
}

double MuonResidualsFitter_logPureGaussian2D	(	double	x,
		double	y,
		double	x0,
		double	y0,
		double	sx,
		double	sy,
		double	r
	)

Definition at line 46 of file MuonResidualsFitter.cc.

References PVValHelper::dx, PVValHelper::dy, cmsBatch::log, M_PI, funct::pow(), alignCSCRings::r, mathSSE::sqrt(), and fftjetcommon_cfi::sy.

Referenced by MuonResiduals5DOFFitter_FCN(), MuonResiduals6DOFFitter_FCN(), and MuonResiduals6DOFrphiFitter_FCN().

{
  sx = fabs(sx);  sy = fabs(sy);
  static const double c2gaus = log( 2.*M_PI );
  double one_r2 = 1. - r*r;
  double dx = x - x0;
  double dy = y - y0;
  return -0.5/one_r2 * ( pow(dx/sx, 2) + pow(dy/sy, 2) - 2.*r*dx/sx*dy/sy ) - c2gaus - log(sx*sy*sqrt(one_r2));
}

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

Definition at line 128 of file MuonResidualsFitter.cc.

References cmsBatch::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 122 of file MuonResidualsFitter.cc.

References JetChargeProducer_cfi::exp, and MuonResidualsFitter_logPowerLawTails().

Referenced by MuonResiduals1DOFFitter::plot(), MuonResidualsAngleFitter::plot(), MuonResidualsBfieldAngleFitter::plot(), MuonResidualsPositionFitter::plot(), MuonResiduals5DOFFitter::plot(), MuonResiduals6DOFrphiFitter::plot(), and MuonResiduals6DOFFitter::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 39 of file MuonResidualsFitter.cc.

References JetChargeProducer_cfi::exp, and MuonResidualsFitter_logPureGaussian().

Referenced by MuonResiduals1DOFFitter::plot(), MuonResidualsAngleFitter::plot(), MuonResidualsBfieldAngleFitter::plot(), MuonResidualsPositionFitter::plot(), MuonResiduals5DOFFitter::plot(), MuonResiduals6DOFrphiFitter::plot(), and MuonResiduals6DOFFitter::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 134 of file MuonResidualsFitter.cc.

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

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

void MuonResidualsPositionFitter_FCN	(	int &	npar,
		double *	gin,
		double &	fval,
		double *	par,
		int	iflag
	)

Definition at line 9 of file MuonResidualsPositionFitter.cc.

References MuonResidualsFitterFitInfo::fitter(), MuonResidualsPositionFitter::kAngleError, MuonResidualsPositionFitter::kGamma, MuonResidualsFitter::kGaussPowerTails, MuonResidualsPositionFitter::kPhiz, MuonResidualsPositionFitter::kPosition, MuonResidualsFitter::kPowerLawTails, MuonResidualsFitter::kPureGaussian, MuonResidualsPositionFitter::kResidual, MuonResidualsFitter::kROOTVoigt, MuonResidualsPositionFitter::kScattering, MuonResidualsPositionFitter::kSigma, MuonResidualsPositionFitter::kTrackAngle, MuonResidualsPositionFitter::kTrackPosition, MuonResidualsPositionFitter::kZpos, MuonResidualsFitter_logGaussPowerTails(), MuonResidualsFitter_logPowerLawTails(), MuonResidualsFitter_logPureGaussian(), MuonResidualsFitter_logROOTVoigt(), MuonResidualsPositionFitter_TMinuit, MuonResidualsFitter::residuals_begin(), MuonResidualsFitter::residuals_end(), and MuonResidualsFitter::residualsModel().

Referenced by MuonResidualsPositionFitter::fit().

                                                                                                   {
  MuonResidualsFitterFitInfo *fitinfo = (MuonResidualsFitterFitInfo*)(MuonResidualsPositionFitter_TMinuit->GetObjectFit());
  MuonResidualsFitter *fitter = fitinfo->fitter();

  fval = 0.;
  for (std::vector<double*>::const_iterator resiter = fitter->residuals_begin();  resiter != fitter->residuals_end();  ++resiter) {
    const double residual = (*resiter)[MuonResidualsPositionFitter::kResidual];
    const double angleerror = (*resiter)[MuonResidualsPositionFitter::kAngleError];
    const double trackangle = (*resiter)[MuonResidualsPositionFitter::kTrackAngle];
    const double trackposition = (*resiter)[MuonResidualsPositionFitter::kTrackPosition];
      
    double center = 0.;
    center += par[MuonResidualsPositionFitter::kPosition];
    center += par[MuonResidualsPositionFitter::kZpos] * trackangle;
    center += par[MuonResidualsPositionFitter::kPhiz] * trackposition;
    center += par[MuonResidualsPositionFitter::kScattering] * angleerror;

    if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian) {
      fval += -MuonResidualsFitter_logPureGaussian(residual, center, par[MuonResidualsPositionFitter::kSigma]);
    }
    else if (fitter->residualsModel() == MuonResidualsFitter::kPowerLawTails) {
      fval += -MuonResidualsFitter_logPowerLawTails(residual, center, par[MuonResidualsPositionFitter::kSigma], par[MuonResidualsPositionFitter::kGamma]);
    }
    else if (fitter->residualsModel() == MuonResidualsFitter::kROOTVoigt) {
      fval += -MuonResidualsFitter_logROOTVoigt(residual, center, par[MuonResidualsPositionFitter::kSigma], par[MuonResidualsPositionFitter::kGamma]);
    }
    else if (fitter->residualsModel() == MuonResidualsFitter::kGaussPowerTails) {
      fval += -MuonResidualsFitter_logGaussPowerTails(residual, center, par[MuonResidualsPositionFitter::kSigma]);
    }
    else { assert(false); }
  }
}