Chi2VertexFitter.cc

Go to the documentation of this file.

/* From SimpleFits Package
 * Designed an written by
 * author: Ian M. Nugent
 * Humboldt Foundations
 */
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "RecoTauTag/ImpactParameter/interface/Chi2VertexFitter.h"
#include "RecoTauTag/ImpactParameter/interface/ChiSquareFunctionUpdator.h"
#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/MnUserParameters.h"
#include "Minuit2/MnPrint.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/MnSimplex.h"
#include "Minuit2/CombinedMinimizer.h"
#include "Minuit2/MnMinimize.h"
#include "Minuit2/MnMinos.h"
#include "Minuit2/MnHesse.h"
#include "Minuit2/MnContours.h"
#include "Minuit2/MnPlot.h"
#include "Minuit2/MinosError.h"
#include "Minuit2/ContoursError.h"
#include <iostream>
 
using namespace tauImpactParameter;
 
bool Chi2VertexFitter::fit() {
  if (isFit_ == true)
    return true;  // do not refit
  if (!isConfigured_)
    return false;  // do not fit if configuration failed
  ChiSquareFunctionUpdator updator(this);
  ROOT::Minuit2::MnUserParameters MnPar;
  for (int i = 0; i < par_.GetNrows(); i++) {
    TString name = freeParName(i);
    // if not limited (vhigh <= vlow)
    MnPar.Add(name.Data(),
              par_(i),
              sqrt(fabs(parcov_(i, i))),
              par_(i) - nsigma_ * sqrt(fabs(parcov_(i, i))),
              par_(i) + nsigma_ * sqrt(fabs(parcov_(i, i))));
  }
 
  unsigned int max = 10;
  int numberofcalls = 200 + par_.GetNrows() * 100 + par_.GetNrows() * par_.GetNrows() * 5;
  double tolerance(0.01);
  double edmMin(0.001 * updator.Up() * tolerance);
 
  ROOT::Minuit2::MnMinimize minimize(updator, MnPar);
  ROOT::Minuit2::FunctionMinimum min = minimize(numberofcalls, tolerance);
  for (unsigned int i = 0; i <= max && min.Edm() > edmMin; i++) {
    if (i == max)
      return false;
    min = minimize(i * numberofcalls, tolerance);
  }
  // give return flag based on status
  if (min.IsAboveMaxEdm()) {
    edm::LogWarning("Chi2VertexFitter::Fit") << "Found Vertex that is above EDM " << std::endl;
    return false;
  }
  if (!min.IsValid()) {
    edm::LogWarning("Chi2VertexFitter::Fit") << "Failed min.IsValid()" << std::endl;
    if (!min.HasValidParameters()) {
      edm::LogWarning("Chi2VertexFitter::Fit") << "Failed min.HasValidParameters()" << std::endl;
    }
    if (!min.HasValidCovariance()) {
      edm::LogWarning("Chi2VertexFitter::Fit") << "Failed min.HasValidCovariance()" << std::endl;
    }
    if (!min.HesseFailed()) {
      edm::LogWarning("Chi2VertexFitter::Fit") << "Failed min.HesseFailed()" << std::endl;
    }
    if (!min.HasReachedCallLimit()) {
      edm::LogWarning("Chi2VertexFitter::Fit") << "Failed min.HasReachedCallLimit()" << std::endl;
    }
    return false;
  }
  chi2_ = min.Fval();
  // Get output parameters
  for (int i = 0; i < par_.GetNrows(); i++) {
    par_(i) = min.UserParameters().Value(i);
  }
  // Get output covariance
  for (int i = 0; i < par_.GetNrows(); i++) {
    for (int j = 0; j < par_.GetNrows(); j++) {
      parcov_(i, j) = min.UserCovariance()(i, j);
    }
  }
 
  isFit_ = true;
  return isFit_;
}