MuonResiduals6DOFFitter.cc File Reference

#include "Alignment/MuonAlignmentAlgorithms/interface/MuonResiduals6DOFFitter.h"
#include "TH2F.h"
#include "TMath.h"
#include "TTree.h"
#include "TFile.h"

Go to the source code of this file.

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

Function Documentation

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

Definition at line 88 of file MuonResiduals6DOFFitter.cc.

References assert(), MuonResidualsFitterFitInfo::fitter(), MuonResidualsFitter::k0010, MuonResidualsFitter::k1010, MuonResidualsFitter::k1100, MuonResidualsFitter::k1110, MuonResidualsFitter::k1111, MuonResiduals6DOFFitter::kAlignPhiX, MuonResiduals6DOFFitter::kAlignPhiY, MuonResiduals6DOFFitter::kAlignPhiZ, MuonResiduals6DOFFitter::kAlignX, MuonResiduals6DOFFitter::kAlignY, MuonResiduals6DOFFitter::kAlignZ, MuonResiduals6DOFFitter::kAlphaX, MuonResiduals6DOFFitter::kAlphaY, MuonResiduals6DOFFitter::kAngleX, MuonResiduals6DOFFitter::kAngleY, MuonResidualsFitter::kGaussPowerTails, MuonResiduals6DOFFitter::kPositionX, MuonResiduals6DOFFitter::kPositionY, MuonResidualsFitter::kPowerLawTails, MuonResidualsFitter::kPureGaussian, MuonResidualsFitter::kPureGaussian2D, MuonResiduals6DOFFitter::kRedChi2, MuonResiduals6DOFFitter::kResidX, MuonResiduals6DOFFitter::kResidXGamma, MuonResiduals6DOFFitter::kResidXSigma, MuonResiduals6DOFFitter::kResidY, MuonResiduals6DOFFitter::kResidYGamma, MuonResiduals6DOFFitter::kResidYSigma, MuonResiduals6DOFFitter::kResSlopeX, MuonResiduals6DOFFitter::kResSlopeXGamma, MuonResiduals6DOFFitter::kResSlopeXSigma, MuonResiduals6DOFFitter::kResSlopeY, MuonResiduals6DOFFitter::kResSlopeYGamma, MuonResiduals6DOFFitter::kResSlopeYSigma, MuonResidualsFitter::kROOTVoigt, MuonResidualsFitter_logGaussPowerTails(), MuonResidualsFitter_logPowerLawTails(), MuonResidualsFitter_logPureGaussian(), MuonResidualsFitter_logPureGaussian2D(), MuonResidualsFitter_logROOTVoigt(), MuonResidualsFitter::residuals_begin(), MuonResidualsFitter::residuals_end(), MuonResidualsFitter::residualsModel(), MuonResidualsFitter::useRes(), and puppiForMET_cff::weight.

Referenced by MuonResiduals6DOFFitter::fit().

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

  fval = 0.;
  for (std::vector<double*>::const_iterator resiter = fitter->residuals_begin();  resiter != fitter->residuals_end();  ++resiter)
  {
    const double residX = (*resiter)[MuonResiduals6DOFFitter::kResidX];
    const double residY = (*resiter)[MuonResiduals6DOFFitter::kResidY];
    const double resslopeX = (*resiter)[MuonResiduals6DOFFitter::kResSlopeX];
    const double resslopeY = (*resiter)[MuonResiduals6DOFFitter::kResSlopeY];
    const double positionX = (*resiter)[MuonResiduals6DOFFitter::kPositionX];
    const double positionY = (*resiter)[MuonResiduals6DOFFitter::kPositionY];
    const double angleX = (*resiter)[MuonResiduals6DOFFitter::kAngleX];
    const double angleY = (*resiter)[MuonResiduals6DOFFitter::kAngleY];
    const double redchi2 = (*resiter)[MuonResiduals6DOFFitter::kRedChi2];

    const double alignx = par[MuonResiduals6DOFFitter::kAlignX];
    const double aligny = par[MuonResiduals6DOFFitter::kAlignY];
    const double alignz = par[MuonResiduals6DOFFitter::kAlignZ];
    const double alignphix = par[MuonResiduals6DOFFitter::kAlignPhiX];
    const double alignphiy = par[MuonResiduals6DOFFitter::kAlignPhiY];
    const double alignphiz = par[MuonResiduals6DOFFitter::kAlignPhiZ];
    const double resXsigma = par[MuonResiduals6DOFFitter::kResidXSigma];
    const double resYsigma = par[MuonResiduals6DOFFitter::kResidYSigma];
    const double slopeXsigma = par[MuonResiduals6DOFFitter::kResSlopeXSigma];
    const double slopeYsigma = par[MuonResiduals6DOFFitter::kResSlopeYSigma];
    const double alphax = par[MuonResiduals6DOFFitter::kAlphaX];
    const double alphay = par[MuonResiduals6DOFFitter::kAlphaY];
    const double resXgamma = par[MuonResiduals6DOFFitter::kResidXGamma];
    const double resYgamma = par[MuonResiduals6DOFFitter::kResidYGamma];
    const double slopeXgamma = par[MuonResiduals6DOFFitter::kResSlopeXGamma];
    const double slopeYgamma = par[MuonResiduals6DOFFitter::kResSlopeYGamma];

    double coefX = alphax, coefY = alphay;
    if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian ||
        fitter->residualsModel() == MuonResidualsFitter::kPureGaussian2D) coefX = coefY = 0.;
    double residXpeak = residual_x(alignx, aligny, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY, coefX, resslopeX);
    double residYpeak = residual_y(alignx, aligny, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY, coefY, resslopeY);
    double slopeXpeak = residual_dxdz(alignx, aligny, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY);
    double slopeYpeak = residual_dydz(alignx, aligny, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY);

    double weight = (1./redchi2) * number_of_hits / sum_of_weights;
    if (!weight_alignment) weight = 1.;

    if (!weight_alignment  ||  TMath::Prob(redchi2*12, 12) < 0.99)  // no spikes allowed
    {
      if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian) {
        if (fitter->useRes() == MuonResidualsFitter::k1111) {
          fval += -weight * MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma);
          fval += -weight * MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma);
          fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma);
          fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeY, slopeYpeak, slopeYsigma);
        }
        else if (fitter->useRes() == MuonResidualsFitter::k1110) {
          fval += -weight * MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma);
          fval += -weight * MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma);
          fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma);
        }
        else if (fitter->useRes() == MuonResidualsFitter::k1100) {
          fval += -weight * MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma);
          fval += -weight * MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma);
        }
        else if (fitter->useRes() == MuonResidualsFitter::k1010) {
          fval += -weight * MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma);
          fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma);
        }
        else if (fitter->useRes() == MuonResidualsFitter::k0010) {
          fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma);
        }
        //std::cout<<"FCNx("<<residX<<","<<residXpeak<<","<<resXsigma<<") = "<<MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma)<<std::endl;
        //std::cout<<"FCNy("<<residY<<","<<residYpeak<<","<<resYsigma<<") = "<<MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma)<<std::endl;
        //std::cout<<"FCNsx("<<resslopeX<<","<<slopeXpeak<<","<<slopeXsigma<<") = "<<MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma)<<std::endl;
        //std::cout<<"FCNsy("<<resslopeY<<","<<slopeYpeak<<","<<slopeYsigma<<") = "<<MuonResidualsFitter_logPureGaussian(resslopeY, slopeYpeak, slopeYsigma)<<std::endl;
      }
      else if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian2D) {
        if (fitter->useRes() == MuonResidualsFitter::k1111) {
          fval += -weight * MuonResidualsFitter_logPureGaussian2D(residX, resslopeX, residXpeak, slopeXpeak, resXsigma, slopeXsigma, alphax);
          fval += -weight * MuonResidualsFitter_logPureGaussian2D(residY, resslopeY, residYpeak, slopeYpeak, resYsigma, slopeYsigma, alphay);
        }
        else if (fitter->useRes() == MuonResidualsFitter::k1110) {
          fval += -weight * MuonResidualsFitter_logPureGaussian2D(residX, resslopeX, residXpeak, slopeXpeak, resXsigma, slopeXsigma, alphax);
          fval += -weight * MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma);
        }
        else if (fitter->useRes() == MuonResidualsFitter::k1100) {
          fval += -weight * MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma);
          fval += -weight * MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma);
        }
        else if (fitter->useRes() == MuonResidualsFitter::k1010) {
          fval += -weight * MuonResidualsFitter_logPureGaussian2D(residX, resslopeX, residXpeak, slopeXpeak, resXsigma, slopeXsigma, alphax);
        }
        else if (fitter->useRes() == MuonResidualsFitter::k0010) {
                fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma);
        }
        //std::cout<<"FCNx("<<residX<<","<<resslopeX<<","<<residXpeak<<","<<slopeXpeak<<","<<resXsigma<<","<<slopeXsigma<<","<<alphax<<") = "<<MuonResidualsFitter_logPureGaussian2D(residX, resslopeX, residXpeak, slopeXpeak, resXsigma, slopeXsigma, alphax)<<std::endl;
        //std::cout<<"FCNy("<<residY<<","<<resslopeY<<","<<residYpeak<<","<<slopeYpeak<<","<<resYsigma<<","<<slopeYsigma<<","<<alphay<<") = "<<MuonResidualsFitter_logPureGaussian2D(residY, resslopeY, residYpeak, slopeYpeak, resYsigma, slopeYsigma, alphay)<<std::endl;
      }
      else if (fitter->residualsModel() == MuonResidualsFitter::kPowerLawTails) {
        fval += -weight * MuonResidualsFitter_logPowerLawTails(residX, residXpeak, resXsigma, resXgamma);
        fval += -weight * MuonResidualsFitter_logPowerLawTails(residY, residYpeak, resYsigma, resYgamma);
        fval += -weight * MuonResidualsFitter_logPowerLawTails(resslopeX, slopeXpeak, slopeXsigma, slopeXgamma);
        fval += -weight * MuonResidualsFitter_logPowerLawTails(resslopeY, slopeYpeak, slopeYsigma, slopeYgamma);
      }
      else if (fitter->residualsModel() == MuonResidualsFitter::kROOTVoigt) {
        fval += -weight * MuonResidualsFitter_logROOTVoigt(residX, residXpeak, resXsigma, resXgamma);
        fval += -weight * MuonResidualsFitter_logROOTVoigt(residY, residYpeak, resYsigma, resYgamma);
        fval += -weight * MuonResidualsFitter_logROOTVoigt(resslopeX, slopeXpeak, slopeXsigma, slopeXgamma);
        fval += -weight * MuonResidualsFitter_logROOTVoigt(resslopeY, slopeYpeak, slopeYsigma, slopeYgamma);
      }
      else if (fitter->residualsModel() == MuonResidualsFitter::kGaussPowerTails) {
        fval += -weight * MuonResidualsFitter_logGaussPowerTails(residX, residXpeak, resXsigma);
        fval += -weight * MuonResidualsFitter_logGaussPowerTails(residY, residYpeak, resYsigma);
        fval += -weight * MuonResidualsFitter_logGaussPowerTails(resslopeX, slopeXpeak, slopeXsigma);
        fval += -weight * MuonResidualsFitter_logGaussPowerTails(resslopeY, slopeYpeak, slopeYsigma);
      }
      else { assert(false); }
    }
  }
/*
  const double resXsigma = par[MuonResiduals6DOFFitter::kResidXSigma];
  const double resYsigma = par[MuonResiduals6DOFFitter::kResidYSigma];
  const double slopeXsigma = par[MuonResiduals6DOFFitter::kResSlopeXSigma];
  const double slopeYsigma = par[MuonResiduals6DOFFitter::kResSlopeYSigma];
  const double alphax = par[MuonResiduals6DOFFitter::kAlphaX];
  const double alphay = par[MuonResiduals6DOFFitter::kAlphaY];
  std::cout<<"fval="<<fval<<","<<resXsigma<<","<<slopeXsigma<<","<<alphax<<","<<resYsigma<<","<<slopeYsigma<<","<<alphay<<std::endl;
*/
}