SiStripLorentzAngle/src/SymmetryFit.cc

0001 #include "CalibTracker/SiStripLorentzAngle/interface/SymmetryFit.h"
0002 #include <cmath>
0003 #include <cassert>
0004 #include <memory>
0005
0006 TH1* SymmetryFit::symmetryChi2(std::string basename,
0007                                const std::vector<TH1*>& candidates,
0008                                const std::pair<unsigned, unsigned> range) {
0009   TH1* fake = (TH1*)(candidates[0]->Clone(basename.c_str()));
0010   fake->Reset();
0011   SymmetryFit combined(fake, range);
0012   delete fake;
0013
0014   for (auto const* candidate : candidates) {
0015     SymmetryFit sf(candidate, range);
0016     combined += sf;
0017     delete sf.chi2_;
0018   }
0019
0020   int status = combined.fit();
0021   if (status) {
0022     delete combined.chi2_;
0023     combined.chi2_ = nullptr;
0024   }
0025   return combined.chi2_;
0026 }
0027
0028 TH1* SymmetryFit::symmetryChi2(const TH1* candidate, const std::pair<unsigned, unsigned> range) {
0029   SymmetryFit sf(candidate, range);
0030   int status = sf.fit();
0031   if (status) {
0032     delete sf.chi2_;
0033     sf.chi2_ = nullptr;
0034   }
0035   return sf.chi2_;
0036 }
0037
0038 SymmetryFit::SymmetryFit(const TH1* h, const std::pair<unsigned, unsigned> r)
0039     : symm_candidate_(h),
0040       minDF_(r.second - r.first),
0041       range_(r),
0042       minmaxUsable_(findUsableMinMax()),
0043       ndf_(minmaxUsable_.first < minmaxUsable_.second ? minmaxUsable_.second - minmaxUsable_.first : 0),
0044       chi2_(nullptr) {
0045   makeChi2Histogram();
0046   fillchi2();
0047 }
0048
0049 void SymmetryFit::makeChi2Histogram() {
0050   std::string XXname = name(symm_candidate_->GetName());
0051   unsigned Nbins = 2 * (range_.second - range_.first) + 1;
0052   double binwidth = symm_candidate_->GetBinWidth(1);
0053   double low = symm_candidate_->GetBinCenter(range_.first) - 3 * binwidth / 4;
0054   double up = symm_candidate_->GetBinCenter(range_.second - 1) + 3 * binwidth / 4;
0055   chi2_ = new TH1F(XXname.c_str(), "", Nbins, low, up);
0056 }
0057
0058 std::pair<unsigned, unsigned> SymmetryFit::findUsableMinMax() const {
0059   unsigned NEAR(0), FAR(0);
0060   const std::vector<std::pair<unsigned, unsigned> > cont = continuousRanges();
0061   for (unsigned L = 0; L < cont.size(); L++) {
0062     if (cont[L].first > range_.first)
0063       continue;
0064     for (unsigned R = L; R < cont.size(); R++) {
0065       if (cont[R].second < range_.second)
0066         continue;
0067
0068       const unsigned far = std::min(range_.first - cont[L].first, cont[R].second - range_.second);
0069       const unsigned near = std::max(range_.second < cont[L].second ? 0 : range_.second - cont[L].second,
0070                                      cont[R].first < range_.first ? 0 : cont[R].first - range_.first);
0071
0072       if ((far > near) && (far - near) > (FAR - NEAR)) {
0073         FAR = far;
0074         NEAR = near;
0075       }
0076     }
0077   }
0078   return std::make_pair(NEAR, FAR);
0079 }
0080
0081 std::vector<std::pair<unsigned, unsigned> > SymmetryFit::continuousRanges() const {
0082   std::vector<std::pair<unsigned, unsigned> > ranges;
0083   const unsigned Nbins = symm_candidate_->GetNbinsX();
0084   unsigned i = 0;
0085   while (++i <= Nbins) {
0086     if (symm_candidate_->GetBinError(i)) {
0087       std::pair<unsigned, unsigned> range(i, i + 1);
0088       while (++i <= Nbins && symm_candidate_->GetBinError(i))
0089         range.second++;
0090       ranges.push_back(range);
0091     }
0092   }
0093   return ranges;
0094 }
0095
0096 void SymmetryFit::fillchi2() {
0097   if (ndf_ < minDF_)
0098     return;
0099
0100   for (int i = 1; i <= chi2_->GetNbinsX(); ++i) {
0101     const unsigned L(range_.first - 1 + (i - 1) / 2), R(range_.first + i / 2);
0102     chi2_->SetBinContent(i, chi2(std::make_pair(L, R)));
0103   }
0104 }
0105
0106 float SymmetryFit::chi2(std::pair<unsigned, unsigned> point) {
0107   point.first -= minmaxUsable_.first;
0108   point.second += minmaxUsable_.first;
0109   float XX = 0;
0110   unsigned i = ndf_;
0111   while (i-- > 0) {
0112     XX += chi2_element(point);
0113     point.first--;
0114     point.second++;
0115   }
0116   return XX;
0117 }
0118
0119 float SymmetryFit::chi2_element(std::pair<unsigned, unsigned> range) {
0120   float y0(symm_candidate_->GetBinContent(range.first)), y1(symm_candidate_->GetBinContent(range.second)),
0121       e0(symm_candidate_->GetBinError(range.first)), e1(symm_candidate_->GetBinError(range.second));
0122   assert(e0 && e1);
0123
0124   return pow(y0 - y1, 2) / (e0 * e0 + e1 * e1);
0125 }
0126
0127 int SymmetryFit::fit() {
0128   std::vector<double> p = pol2_from_pol3(chi2_);
0129   if (p.empty() || p[0] < chi2_->GetBinCenter(1) || p[0] > chi2_->GetBinCenter(chi2_->GetNbinsX()))
0130     return 7;
0131
0132   std::unique_ptr<TF1> f(fitfunction());
0133   f->SetParameter(0, p[0]);
0134   f->SetParLimits(0, p[0], p[0]);
0135   f->SetParameter(1, p[1]);
0136   f->SetParLimits(1, p[1], p[1]);
0137   f->SetParameter(2, p[2]);
0138   f->SetParLimits(2, p[2], p[2]);
0139   f->SetParameter(3, ndf_);
0140   f->SetParLimits(3, ndf_, ndf_);  //Fixed
0141   chi2_->Fit(f.get(), "WQ");
0142   return 0;
0143 }
0144
0145 TF1* SymmetryFit::fitfunction() {
0146   TF1* f = new TF1("SymmetryFit", "((x-[0])/[1])**2+[2]+0*[3]");
0147   f->SetParName(0, "middle");
0148   f->SetParName(1, "uncertainty");
0149   f->SetParName(2, "chi2");
0150   f->SetParName(3, "NDF");
0151   return f;
0152 }
0153
0154 std::vector<double> SymmetryFit::pol2_from_pol2(TH1* hist) {
0155   std::vector<double> v;
0156
0157   //Need our own copy for thread safety
0158   TF1 func("mypol2", "pol2");
0159   int status = hist->Fit(&func, "WQ");
0160   if (!status) {
0161     std::vector<double> p;
0162     p.push_back(func.GetParameter(0));
0163     p.push_back(func.GetParameter(1));
0164     p.push_back(func.GetParameter(2));
0165     if (p[2] > 0) {
0166       v.push_back(-0.5 * p[1] / p[2]);
0167       v.push_back(1. / sqrt(p[2]));
0168       v.push_back(p[0] - 0.25 * p[1] * p[1] / p[2]);
0169     }
0170   }
0171   return v;
0172 }
0173
0174 std::vector<double> SymmetryFit::pol2_from_pol3(TH1* hist) {
0175   std::vector<double> v;
0176
0177   auto func = std::make_unique<TF1>("mypol3", "pol3");
0178   int status = hist->Fit(func.get(), "WQ");
0179   if (!status) {
0180     std::vector<double> p;
0181     p.push_back(func->GetParameter(0));
0182     p.push_back(func->GetParameter(1));
0183     p.push_back(func->GetParameter(2));
0184     p.push_back(func->GetParameter(3));
0185     double radical = p[2] * p[2] - 3 * p[1] * p[3];
0186     if (radical > 0) {
0187       double x0 = (-p[2] + sqrt(radical)) / (3 * p[3]);
0188       v.push_back(x0);
0189       v.push_back(pow(radical, -0.25));
0190       v.push_back(p[0] + p[1] * x0 + p[2] * x0 * x0 + p[3] * x0 * x0 * x0);
0191     }
0192   }
0193   return v;
0194 }