Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-09-09 02:18:24

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: /*case MULTI:*/ {
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     //Need our own copy for thread safety
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     //Need our own copy for thread safety
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 }