SymmetryFit.cc

Go to the documentation of this file.

#include "CalibTracker/SiStripLorentzAngle/interface/SymmetryFit.h"
#include <cmath>
#include <cassert>
#include <memory>

TH1* SymmetryFit::symmetryChi2(std::string basename,
                               const std::vector<TH1*>& candidates,
                               const std::pair<unsigned, unsigned> range) {
  TH1* fake = (TH1*)(candidates[0]->Clone(basename.c_str()));
  fake->Reset();
  SymmetryFit combined(fake, range);
  delete fake;

  for (auto const* candidate : candidates) {
    SymmetryFit sf(candidate, range);
    combined += sf;
    delete sf.chi2_;
  }

  int status = combined.fit();
  if (status) {
    delete combined.chi2_;
    combined.chi2_ = nullptr;
  }
  return combined.chi2_;
}

TH1* SymmetryFit::symmetryChi2(const TH1* candidate, const std::pair<unsigned, unsigned> range) {
  SymmetryFit sf(candidate, range);
  int status = sf.fit();
  if (status) {
    delete sf.chi2_;
    sf.chi2_ = nullptr;
  }
  return sf.chi2_;
}

SymmetryFit::SymmetryFit(const TH1* h, const std::pair<unsigned, unsigned> r)
    : symm_candidate_(h),
      minDF_(r.second - r.first),
      range_(r),
      minmaxUsable_(findUsableMinMax()),
      ndf_(minmaxUsable_.first < minmaxUsable_.second ? minmaxUsable_.second - minmaxUsable_.first : 0),
      chi2_(nullptr) {
  makeChi2Histogram();
  fillchi2();
}

void SymmetryFit::makeChi2Histogram() {
  std::string XXname = name(symm_candidate_->GetName());
  unsigned Nbins = 2 * (range_.second - range_.first) + 1;
  double binwidth = symm_candidate_->GetBinWidth(1);
  double low = symm_candidate_->GetBinCenter(range_.first) - 3 * binwidth / 4;
  double up = symm_candidate_->GetBinCenter(range_.second - 1) + 3 * binwidth / 4;
  chi2_ = new TH1F(XXname.c_str(), "", Nbins, low, up);
}

std::pair<unsigned, unsigned> SymmetryFit::findUsableMinMax() const {
  unsigned NEAR(0), FAR(0);
  const std::vector<std::pair<unsigned, unsigned> > cont = continuousRanges();
  for (unsigned L = 0; L < cont.size(); L++) {
    if (cont[L].first > range_.first)
      continue;
    for (unsigned R = L; R < cont.size(); R++) {
      if (cont[R].second < range_.second)
        continue;

      const unsigned far = std::min(range_.first - cont[L].first, cont[R].second - range_.second);
      const unsigned near = std::max(range_.second < cont[L].second ? 0 : range_.second - cont[L].second,
                                     cont[R].first < range_.first ? 0 : cont[R].first - range_.first);

      if ((far > near) && (far - near) > (FAR - NEAR)) {
        FAR = far;
        NEAR = near;
      }
    }
  }
  return std::make_pair(NEAR, FAR);
}

std::vector<std::pair<unsigned, unsigned> > SymmetryFit::continuousRanges() const {
  std::vector<std::pair<unsigned, unsigned> > ranges;
  const unsigned Nbins = symm_candidate_->GetNbinsX();
  unsigned i = 0;
  while (++i <= Nbins) {
    if (symm_candidate_->GetBinError(i)) {
      std::pair<unsigned, unsigned> range(i, i + 1);
      while (++i <= Nbins && symm_candidate_->GetBinError(i))
        range.second++;
      ranges.push_back(range);
    }
  }
  return ranges;
}

void SymmetryFit::fillchi2() {
  if (ndf_ < minDF_)
    return;

  for (int i = 1; i <= chi2_->GetNbinsX(); ++i) {
    const unsigned L(range_.first - 1 + (i - 1) / 2), R(range_.first + i / 2);
    chi2_->SetBinContent(i, chi2(std::make_pair(L, R)));
  }
}

float SymmetryFit::chi2(std::pair<unsigned, unsigned> point) {
  point.first -= minmaxUsable_.first;
  point.second += minmaxUsable_.first;
  float XX = 0;
  unsigned i = ndf_;
  while (i-- > 0) {
    XX += chi2_element(point);
    point.first--;
    point.second++;
  }
  return XX;
}

float SymmetryFit::chi2_element(std::pair<unsigned, unsigned> range) {
  float y0(symm_candidate_->GetBinContent(range.first)), y1(symm_candidate_->GetBinContent(range.second)),
      e0(symm_candidate_->GetBinError(range.first)), e1(symm_candidate_->GetBinError(range.second));
  assert(e0 && e1);

  return pow(y0 - y1, 2) / (e0 * e0 + e1 * e1);
}

int SymmetryFit::fit() {
  std::vector<double> p = pol2_from_pol3(chi2_);
  if (p.empty() || p[0] < chi2_->GetBinCenter(1) || p[0] > chi2_->GetBinCenter(chi2_->GetNbinsX()))
    return 7;

  std::unique_ptr<TF1> f(fitfunction());
  f->SetParameter(0, p[0]);
  f->SetParLimits(0, p[0], p[0]);
  f->SetParameter(1, p[1]);
  f->SetParLimits(1, p[1], p[1]);
  f->SetParameter(2, p[2]);
  f->SetParLimits(2, p[2], p[2]);
  f->SetParameter(3, ndf_);
  f->SetParLimits(3, ndf_, ndf_);  //Fixed
  chi2_->Fit(f.get(), "WQ");
  return 0;
}

TF1* SymmetryFit::fitfunction() {
  TF1* f = new TF1("SymmetryFit", "((x-[0])/[1])**2+[2]+0*[3]");
  f->SetParName(0, "middle");
  f->SetParName(1, "uncertainty");
  f->SetParName(2, "chi2");
  f->SetParName(3, "NDF");
  return f;
}

std::vector<double> SymmetryFit::pol2_from_pol2(TH1* hist) {
  std::vector<double> v;

  //Need our own copy for thread safety
  TF1 func("mypol2", "pol2");
  int status = hist->Fit(&func, "WQ");
  if (!status) {
    std::vector<double> p;
    p.push_back(func.GetParameter(0));
    p.push_back(func.GetParameter(1));
    p.push_back(func.GetParameter(2));
    if (p[2] > 0) {
      v.push_back(-0.5 * p[1] / p[2]);
      v.push_back(1. / sqrt(p[2]));
      v.push_back(p[0] - 0.25 * p[1] * p[1] / p[2]);
    }
  }
  return v;
}

std::vector<double> SymmetryFit::pol2_from_pol3(TH1* hist) {
  std::vector<double> v;

  auto func = std::make_unique<TF1>("mypol3", "pol3");
  int status = hist->Fit(func.get(), "WQ");
  if (!status) {
    std::vector<double> p;
    p.push_back(func->GetParameter(0));
    p.push_back(func->GetParameter(1));
    p.push_back(func->GetParameter(2));
    p.push_back(func->GetParameter(3));
    double radical = p[2] * p[2] - 3 * p[1] * p[3];
    if (radical > 0) {
      double x0 = (-p[2] + sqrt(radical)) / (3 * p[3]);
      v.push_back(x0);
      v.push_back(pow(radical, -0.25));
      v.push_back(p[0] + p[1] * x0 + p[2] * x0 * x0 + p[3] * x0 * x0 * x0);
    }
  }
  return v;
}