File indexing completed on 2023-03-17 10:44:41
0001 #include "CalibTracker/SiStripLorentzAngle/interface/LA_Filler_Fitter.h"
0002 #include "FWCore/Utilities/interface/EDMException.h"
0003
0004 #include <cmath>
0005 #include <regex>
0006 #include <boost/algorithm/string/erase.hpp>
0007 #include <TF1.h>
0008 #include <TGraphErrors.h>
0009 #include <TProfile.h>
0010
0011 LA_Filler_Fitter::Result LA_Filler_Fitter::result(Method m, const std::string name, const Book& book) {
0012 Result p;
0013 const std::string base = boost::erase_all_copy(name, method(m));
0014
0015 const TH1* const h = book[name];
0016 const TH1* const reco = book[base + "_reconstruction"];
0017 const TH1* const field = book[base + "_field"];
0018
0019 if (reco)
0020 p.reco = std::make_pair<float, float>(reco->GetMean(), reco->GetMeanError());
0021 if (field)
0022 p.field = field->GetMean();
0023 if (h) {
0024 switch (m) {
0025 case WIDTH: {
0026 const TF1* const f = h->GetFunction("LA_profile_fit");
0027 if (!f)
0028 break;
0029 p.measured = std::make_pair<float, float>(f->GetParameter(0), f->GetParError(0));
0030 p.chi2 = f->GetChisquare();
0031 p.ndof = f->GetNDF();
0032 p.entries = (unsigned)(h->GetEntries());
0033 break;
0034 }
0035 case PROB1:
0036 case AVGV2:
0037 case AVGV3:
0038 case RMSV2:
0039 case RMSV3: {
0040 const TF1* const f = h->GetFunction("SymmetryFit");
0041 if (!f)
0042 break;
0043 p.measured = std::make_pair<float, float>(p.reco.first + f->GetParameter(0), f->GetParameter(1));
0044 p.chi2 = f->GetParameter(2);
0045 p.ndof = (unsigned)(f->GetParameter(3));
0046 p.entries = (m & PROB1) ? (unsigned)book[base + "_w1"]->GetEntries()
0047 : (m & (AVGV2 | RMSV2)) ? (unsigned)book[base + method(AVGV2, false)]->GetEntries()
0048 : (m & (AVGV3 | RMSV3)) ? (unsigned)book[base + method(AVGV3, false)]->GetEntries()
0049 : 0;
0050 break;
0051 }
0052 default:
0053 break;
0054 }
0055 }
0056 return p;
0057 }
0058
0059 std::map<uint32_t, LA_Filler_Fitter::Result> LA_Filler_Fitter::module_results(const Book& book, const Method m) {
0060 std::map<uint32_t, Result> results;
0061 for (Book::const_iterator it = book.begin(".*_module\\d*" + method(m)); it != book.end(); ++it) {
0062 const uint32_t detid =
0063 std::stoul(std::regex_replace(it->first, std::regex(".*_module(\\d*)_.*"), std::string("\\1")));
0064 results[detid] = result(m, it->first, book);
0065 }
0066 return results;
0067 }
0068
0069 std::map<std::string, LA_Filler_Fitter::Result> LA_Filler_Fitter::layer_results(const Book& book, const Method m) {
0070 std::map<std::string, Result> results;
0071 for (Book::const_iterator it = book.begin(".*layer\\d.*" + method(m)); it != book.end(); ++it) {
0072 const std::string name = boost::erase_all_copy(it->first, method(m));
0073 results[name] = result(m, it->first, book);
0074 }
0075 return results;
0076 }
0077
0078 std::map<std::string, std::vector<LA_Filler_Fitter::Result> > LA_Filler_Fitter::ensemble_results(const Book& book,
0079 const Method m) {
0080 std::map<std::string, std::vector<Result> > results;
0081 for (Book::const_iterator it = book.begin(".*_sample.*" + method(m)); it != book.end(); ++it) {
0082 const std::string name = std::regex_replace(it->first, std::regex("sample\\d*_"), "");
0083 results[name].push_back(result(m, it->first, book));
0084 }
0085 return results;
0086 }
0087
0088 void LA_Filler_Fitter::summarize_ensembles(Book& book) const {
0089 typedef std::map<std::string, std::vector<Result> > results_t;
0090 results_t results;
0091 for (int m = FIRST_METHOD; m <= LAST_METHOD; m <<= 1)
0092 if (methods_ & m) {
0093 results_t g = ensemble_results(book, (Method)(m));
0094 results.insert(g.begin(), g.end());
0095 }
0096
0097 for (auto const& group : results) {
0098 const std::string name = group.first;
0099 for (auto const& p : group.second) {
0100 const float pad = (ensembleUp_ - ensembleLow_) / 10;
0101 book.fill(p.reco.first, name + "_ensembleReco", 12 * ensembleBins_, ensembleLow_ - pad, ensembleUp_ + pad);
0102 book.fill(p.measured.first, name + "_measure", 12 * ensembleBins_, ensembleLow_ - pad, ensembleUp_ + pad);
0103 book.fill(p.measured.second, name + "_merr", 500, 0, 0.01);
0104 book.fill((p.measured.first - p.reco.first) / p.measured.second, name + "_pull", 500, -10, 10);
0105 }
0106
0107 TF1 gaus("mygaus", "gaus");
0108 book[name + "_measure"]->Fit(&gaus, "LLQ");
0109 book[name + "_merr"]->Fit(&gaus, "LLQ");
0110 book[name + "_pull"]->Fit(&gaus, "LLQ");
0111 }
0112 }
0113
0114 std::map<std::string, std::vector<LA_Filler_Fitter::EnsembleSummary> > LA_Filler_Fitter::ensemble_summary(
0115 const Book& book) {
0116 std::map<std::string, std::vector<EnsembleSummary> > summary;
0117 for (Book::const_iterator it = book.begin(".*_ensembleReco"); it != book.end(); ++it) {
0118 const std::string base = boost::erase_all_copy(it->first, "_ensembleReco");
0119
0120 const TH1* const reco = it->second;
0121 const TH1* const measure = book[base + "_measure"];
0122 const TH1* const merr = book[base + "_merr"];
0123 const TH1* const pull = book[base + "_pull"];
0124
0125 EnsembleSummary s;
0126 s.samples = reco->GetEntries();
0127 s.truth = reco->GetMean();
0128 s.meanMeasured = std::make_pair<float, float>(measure->GetFunction("gaus")->GetParameter(1),
0129 measure->GetFunction("gaus")->GetParError(1));
0130 s.sigmaMeasured = std::make_pair<float, float>(measure->GetFunction("gaus")->GetParameter(2),
0131 measure->GetFunction("gaus")->GetParError(2));
0132 s.meanUncertainty = std::make_pair<float, float>(merr->GetFunction("gaus")->GetParameter(1),
0133 merr->GetFunction("gaus")->GetParError(1));
0134 s.pull = std::make_pair<float, float>(pull->GetFunction("gaus")->GetParameter(2),
0135 pull->GetFunction("gaus")->GetParError(2));
0136
0137 const std::string name = std::regex_replace(base, std::regex("ensembleBin\\d*_"), "");
0138 summary[name].push_back(s);
0139 }
0140 return summary;
0141 }
0142
0143 std::pair<std::pair<float, float>, std::pair<float, float> > LA_Filler_Fitter::offset_slope(
0144 const std::vector<LA_Filler_Fitter::EnsembleSummary>& ensembles) {
0145 try {
0146 std::vector<float> x, y, xerr, yerr;
0147 for (auto const& ensemble : ensembles) {
0148 x.push_back(ensemble.truth);
0149 xerr.push_back(0);
0150 y.push_back(ensemble.meanMeasured.first);
0151 yerr.push_back(ensemble.meanMeasured.second);
0152 }
0153 TGraphErrors graph(x.size(), &(x[0]), &(y[0]), &(xerr[0]), &(yerr[0]));
0154
0155 TF1 pol1("mypol1", "pol1");
0156 graph.Fit(&pol1);
0157 const TF1* const fit = graph.GetFunction("pol1");
0158
0159 return std::make_pair(std::make_pair(fit->GetParameter(0), fit->GetParError(0)),
0160 std::make_pair(fit->GetParameter(1), fit->GetParError(1)));
0161 } catch (edm::Exception const& e) {
0162 std::cerr << "Fitting Line Failed " << std::endl << e << std::endl;
0163 return std::make_pair(std::make_pair(0, 0), std::make_pair(0, 0));
0164 }
0165 }
0166
0167 float LA_Filler_Fitter::pull(const std::vector<LA_Filler_Fitter::EnsembleSummary>& ensembles) {
0168 float p(0), w(0);
0169 for (auto const& ensemble : ensembles) {
0170 const float unc2 = pow(ensemble.pull.second, 2);
0171 p += ensemble.pull.first / unc2;
0172 w += 1 / unc2;
0173 }
0174 return p / w;
0175 }
0176
0177 std::ostream& operator<<(std::ostream& strm, const LA_Filler_Fitter::Result& r) {
0178 return strm << r.reco.first << "\t" << r.reco.second << "\t" << r.measured.first << "\t" << r.measured.second << "\t"
0179 << r.calMeasured.first << "\t" << r.calMeasured.second << "\t" << r.field << "\t" << r.chi2 << "\t"
0180 << r.ndof << "\t" << r.entries;
0181 }
0182
0183 std::ostream& operator<<(std::ostream& strm, const LA_Filler_Fitter::EnsembleSummary& e) {
0184 return strm << e.truth << "\t" << e.meanMeasured.first << "\t" << e.meanMeasured.second << "\t"
0185 << e.sigmaMeasured.first << "\t" << e.sigmaMeasured.second << "\t" << e.meanUncertainty.first << "\t"
0186 << e.meanUncertainty.second << "\t" << e.pull.first << "\t" << e.pull.second << "\t" << e.samples;
0187 }
