File indexing completed on 2024-04-06 11:59:59
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_);
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
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 }