MuonResiduals1DOFFitter.cc

Go to the documentation of this file.

#include "Alignment/MuonAlignmentAlgorithms/interface/MuonResiduals1DOFFitter.h"
#include "TMath.h"

static TMinuit *MuonResiduals1DOFFitter_TMinuit;
static double MuonResiduals1DOFFitter_sum_of_weights;
static double MuonResiduals1DOFFitter_number_of_hits;
static bool MuonResiduals1DOFFitter_weightAlignment;

void MuonResiduals1DOFFitter::inform(TMinuit *tMinuit) {
  MuonResiduals1DOFFitter_TMinuit = tMinuit;
}

void MuonResiduals1DOFFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag) {
  MuonResidualsFitterFitInfo *fitinfo = (MuonResidualsFitterFitInfo*)(MuonResiduals1DOFFitter_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)[MuonResiduals1DOFFitter::kResid];
    const double redchi2 = (*resiter)[MuonResiduals1DOFFitter::kRedChi2];

    const double residpeak = par[MuonResiduals1DOFFitter::kAlign];
    const double residsigma = par[MuonResiduals1DOFFitter::kSigma];
    const double residgamma = par[MuonResiduals1DOFFitter::kGamma];

    double weight = (1./redchi2) * MuonResiduals1DOFFitter_number_of_hits / MuonResiduals1DOFFitter_sum_of_weights;
    if (!MuonResiduals1DOFFitter_weightAlignment) weight = 1.;

    if (!MuonResiduals1DOFFitter_weightAlignment  ||  TMath::Prob(redchi2*8, 8) < 0.99) {  // no spikes allowed
      if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian) {
    fval += -weight * MuonResidualsFitter_logPureGaussian(residual, residpeak, residsigma);
      }
      else if (fitter->residualsModel() == MuonResidualsFitter::kPowerLawTails) {
    fval += -weight * MuonResidualsFitter_logPowerLawTails(residual, residpeak, residsigma, residgamma);
      }
      else if (fitter->residualsModel() == MuonResidualsFitter::kROOTVoigt) {
    fval += -weight * MuonResidualsFitter_logROOTVoigt(residual, residpeak, residsigma, residgamma);
      }
      else if (fitter->residualsModel() == MuonResidualsFitter::kGaussPowerTails) {
        fval += -weight * MuonResidualsFitter_logGaussPowerTails(residual, residpeak, residsigma);
      }
      else { assert(false); }

    }
  }
}

double MuonResiduals1DOFFitter::sumofweights() {
  MuonResiduals1DOFFitter_sum_of_weights = 0.;
  MuonResiduals1DOFFitter_number_of_hits = 0.;
  MuonResiduals1DOFFitter_weightAlignment = m_weightAlignment;
  for (std::vector<double*>::const_iterator resiter = residuals_begin();  resiter != residuals_end();  ++resiter) {
    if (m_weightAlignment) {
       double redchi2 = (*resiter)[MuonResiduals1DOFFitter::kRedChi2];
       if (TMath::Prob(redchi2*8, 8) < 0.99) {  // no spikes allowed
      MuonResiduals1DOFFitter_sum_of_weights += 1./redchi2;
      MuonResiduals1DOFFitter_number_of_hits += 1.;
       }
    }
    else {
      MuonResiduals1DOFFitter_sum_of_weights += 1.;
      MuonResiduals1DOFFitter_number_of_hits += 1.;
    }
  }
  return MuonResiduals1DOFFitter_sum_of_weights;
}

bool MuonResiduals1DOFFitter::fit(Alignable *ali) {
  initialize_table();  // if not already initialized
  sumofweights();

  double resid_sum = 0.;
  double resid_sum2 = 0.;
  double resid_N = 0.;
  int N = 0;

  for (std::vector<double*>::const_iterator resiter = residuals_begin();  resiter != residuals_end();  ++resiter) {
    const double residual = (*resiter)[MuonResiduals1DOFFitter::kResid];
    const double redchi2 = (*resiter)[MuonResiduals1DOFFitter::kRedChi2];
    double weight = 1./redchi2;
    if (!m_weightAlignment) weight = 1.;

    if (!m_weightAlignment  ||  TMath::Prob(redchi2*8, 8) < 0.99) {  // no spikes allowed
      if (fabs(residual) < 10.) {   // 10 cm
    resid_sum += weight * residual;
    resid_sum2 += weight * residual * residual;
    resid_N += weight;
    N++;
      }
    }
  }

  double resid_mean = resid_sum/resid_N;
  double resid_stdev = sqrt(resid_sum2/resid_N - pow(resid_sum/resid_N, 2));

  resid_sum = 0.;
  resid_sum2 = 0.;
  resid_N = 0.;

  for (std::vector<double*>::const_iterator resiter = residuals_begin();  resiter != residuals_end();  ++resiter) {
    const double residual = (*resiter)[MuonResiduals1DOFFitter::kResid];
    const double redchi2 = (*resiter)[MuonResiduals1DOFFitter::kRedChi2];
    double weight = 1./redchi2;
    if (!m_weightAlignment) weight = 1.;

    if (!m_weightAlignment  ||  TMath::Prob(redchi2*8, 8) < 0.99) {  // no spikes allowed
      if (fabs(residual - resid_mean) < 2.5*resid_stdev) {
    resid_sum += weight * residual;
    resid_sum2 += weight * residual * residual;
    resid_N += weight;
      }
    }
  }

  resid_mean = resid_sum/resid_N;
  resid_stdev = sqrt(resid_sum2/resid_N - pow(resid_sum/resid_N, 2));

  std::vector<int> num;
  std::vector<std::string> name;
  std::vector<double> start;
  std::vector<double> step;
  std::vector<double> low;
  std::vector<double> high;

  if (fixed(kAlign)) {
  num.push_back(kAlign);          name.push_back(std::string("Align"));          start.push_back(0.);              step.push_back(0.01*resid_stdev);    low.push_back(0.);   high.push_back(0.);
  } else {
  num.push_back(kAlign);          name.push_back(std::string("Align"));          start.push_back(resid_mean);      step.push_back(0.01*resid_stdev);    low.push_back(0.);   high.push_back(0.);
  }
  num.push_back(kSigma);          name.push_back(std::string("Sigma"));          start.push_back(resid_stdev);     step.push_back(0.01*resid_stdev);    low.push_back(0.);   high.push_back(0.);
  if (residualsModel() != kPureGaussian && residualsModel() != kGaussPowerTails) {
  num.push_back(kGamma);          name.push_back(std::string("Gamma"));          start.push_back(0.1*resid_stdev); step.push_back(0.01*resid_stdev);    low.push_back(0.);   high.push_back(0.);
  }

  return dofit(&MuonResiduals1DOFFitter_FCN, num, name, start, step, low, high);
}

double MuonResiduals1DOFFitter::plot(std::string name, TFileDirectory *dir, Alignable *ali) {
  sumofweights();

  std::stringstream name_residual, name_residual_raw;
  name_residual << name << "_residual";
  name_residual_raw << name << "_residual_raw";

  double min_residual = -100.;     double max_residual = 100.;
  TH1F *hist_residual = dir->make<TH1F>(name_residual.str().c_str(), "", 100, min_residual, max_residual);
  TH1F *hist_residual_raw = dir->make<TH1F>(name_residual_raw.str().c_str(), "", 100, min_residual, max_residual);

  name_residual << "_fit";
  TF1 *fit_residual = nullptr;
  if (residualsModel() == kPureGaussian) {
    fit_residual = new TF1(name_residual.str().c_str(), MuonResidualsFitter_pureGaussian_TF1, min_residual, max_residual, 3);
    fit_residual->SetParameters(MuonResiduals1DOFFitter_sum_of_weights * (max_residual - min_residual)/100., 10.*value(kAlign), 10.*value(kSigma));
  }
  else if (residualsModel() == kPowerLawTails) {
    fit_residual = new TF1(name_residual.str().c_str(), MuonResidualsFitter_powerLawTails_TF1, min_residual, max_residual, 4);
    fit_residual->SetParameters(MuonResiduals1DOFFitter_sum_of_weights * (max_residual - min_residual)/100., 10.*value(kAlign), 10.*value(kSigma), 10.*value(kGamma));
  }
  else if (residualsModel() == kROOTVoigt) {
    fit_residual = new TF1(name_residual.str().c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_residual, max_residual, 4);
    fit_residual->SetParameters(MuonResiduals1DOFFitter_sum_of_weights * (max_residual - min_residual)/100., 10.*value(kAlign), 10.*value(kSigma), 10.*value(kGamma));
  }
  else if (residualsModel() == kGaussPowerTails) {
    fit_residual = new TF1(name_residual.str().c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_residual, max_residual, 3);
    fit_residual->SetParameters(MuonResiduals1DOFFitter_sum_of_weights * (max_residual - min_residual)/100., 10.*value(kAlign), 10.*value(kSigma));
  }
  else { assert(false); }

  fit_residual->SetLineColor(2);  fit_residual->SetLineWidth(2);
  fit_residual->Write();

  for (std::vector<double*>::const_iterator resiter = residuals_begin();  resiter != residuals_end();  ++resiter) {
    const double resid = (*resiter)[MuonResiduals1DOFFitter::kResid];
    const double redchi2 = (*resiter)[MuonResiduals1DOFFitter::kRedChi2];
    double weight = 1./redchi2;
    if (!m_weightAlignment) weight = 1.;

    if (!m_weightAlignment  ||  TMath::Prob(redchi2*8, 8) < 0.99) {  // no spikes allowed
      hist_residual->Fill(10.*(resid + value(kAlign)), weight);
    }

    hist_residual_raw->Fill(10.*resid);
  }

  double chi2 = 0.;
  double ndof = 0.;
  for (int i = 1;  i <= hist_residual->GetNbinsX();  i++) {
    double xi = hist_residual->GetBinCenter(i);
    double yi = hist_residual->GetBinContent(i);
    double yerri = hist_residual->GetBinError(i);
    double yth = fit_residual->Eval(xi);
    if (yerri > 0.) {
      chi2 += pow((yth - yi)/yerri, 2);
      ndof += 1.;
    }
  }
  ndof -= npar();

  return (ndof > 0. ? chi2 / ndof : -1.);
}