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 <boost/algorithm/string/erase.hpp>
0006 #include <TF1.h>
0007
0008 void LA_Filler_Fitter::fit_width_profile(Book& book) {
0009 for (Book::iterator it = book.begin(".*" + method(WIDTH)); it != book.end(); ++it) {
0010 it->second->SetTitle("Mean Cluster Width;tan#theta_{t}");
0011 TH1* const p = it->second;
0012 if (p->GetEntries() < 400) {
0013 delete p;
0014 book[it->first] = nullptr;
0015 continue;
0016 }
0017 p->SetTitle(";tan#theta_{t};");
0018 const float min = p->GetMinimum();
0019 const float max = p->GetMaximum();
0020 float xofmin = p->GetBinCenter(p->GetMinimumBin());
0021 if (xofmin > 0.0 || xofmin < -0.15)
0022 xofmin = -0.05;
0023 const float xofmax = p->GetBinCenter(p->GetMaximumBin());
0024
0025 TF1* const fit = new TF1("LA_profile_fit", "[2]*(TMath::Abs(x-[0]))+[1]", -1, 1);
0026 fit->SetParLimits(0, -0.15, 0.01);
0027 fit->SetParLimits(1, 0.6 * min, 1.25 * min);
0028 fit->SetParLimits(2, 0.1, 10);
0029 fit->SetParameters(xofmin, min, (max - min) / fabs(xofmax - xofmin));
0030
0031 int badfit = p->Fit(fit, "IEQ", "", -.5, .3);
0032 if (badfit)
0033 badfit = p->Fit(fit, "IEQ", "", -.46, .26);
0034 if (badfit) {
0035 book.erase(it);
0036 }
0037 }
0038 }
0039
0040 void LA_Filler_Fitter::make_and_fit_symmchi2(Book& book) {
0041 for (Book::iterator it = book.begin(".*_all"); it != book.end(); ++it) {
0042 const std::string base = boost::erase_all_copy(it->first, "_all");
0043
0044 std::vector<Book::iterator> rebin_hists;
0045 const Book::iterator& all = it;
0046 rebin_hists.push_back(all);
0047 Book::iterator w1 = book.find(base + "_w1");
0048 rebin_hists.push_back(w1);
0049 Book::iterator var_w2 = book.find(base + method(AVGV2, false));
0050 rebin_hists.push_back(var_w2);
0051 Book::iterator var_w3 = book.find(base + method(AVGV3, false));
0052 rebin_hists.push_back(var_w3);
0053
0054 const unsigned rebin = std::max(var_w2 == book.end() ? 0 : find_rebin(var_w2->second),
0055 var_w3 == book.end() ? 0 : find_rebin(var_w3->second));
0056 for (const auto& it : rebin_hists)
0057 if (it != book.end())
0058 it->second->Rebin(rebin > 1 ? rebin < 7 ? rebin : 6 : 1);
0059
0060 TH1* const prob_w1 =
0061 w1 == book.end() ? nullptr : subset_probability(base + method(PROB1, false), w1->second, all->second);
0062 TH1* const rmsv_w2 =
0063 var_w2 == book.end() ? nullptr : rms_profile(base + method(RMSV2, false), (TProfile* const)var_w2->second);
0064 TH1* const rmsv_w3 =
0065 var_w3 == book.end() ? nullptr : rms_profile(base + method(RMSV3, false), (TProfile* const)var_w3->second);
0066
0067 std::vector<TH1*> fit_hists;
0068 if (prob_w1) {
0069 book.book(base + method(PROB1, false), prob_w1);
0070 fit_hists.push_back(prob_w1);
0071 prob_w1->SetTitle("Width==1 Probability;tan#theta_{t}-(dx/dz)_{reco}");
0072 }
0073 if (var_w2 != book.end()) {
0074 book.book(base + method(RMSV2, false), rmsv_w2);
0075 fit_hists.push_back(var_w2->second);
0076 var_w2->second->SetTitle("Width==2 Mean Variance;tan#theta_{t}-(dx/dz)_{reco}");
0077 fit_hists.push_back(rmsv_w2);
0078 rmsv_w2->SetTitle("Width==2 RMS Variance;tan#theta_{t}-(dx/dz)_{reco}");
0079 }
0080 if (var_w3 != book.end()) {
0081 book.book(base + method(RMSV3, false), rmsv_w3);
0082 fit_hists.push_back(var_w3->second);
0083 var_w3->second->SetTitle("Width==3 Mean Variance;tan#theta_{t}-(dx/dz)_{reco}");
0084 fit_hists.push_back(rmsv_w3);
0085 rmsv_w3->SetTitle("Width==3 RMS Variance;tan#theta_{t}-(dx/dz)_{reco}");
0086 }
0087
0088 if (fit_hists.empty())
0089 continue;
0090 const unsigned bins = fit_hists[0]->GetNbinsX();
0091 const unsigned guess = fit_hists[0]->FindBin(0);
0092 const std::pair<unsigned, unsigned> range(guess - bins / 30, guess + bins / 30);
0093
0094 for (auto const& hist : fit_hists) {
0095 TH1* const chi2 = SymmetryFit::symmetryChi2(hist, range);
0096 if (chi2) {
0097 book.book(chi2->GetName(), chi2);
0098 chi2->SetTitle("Symmetry #chi^{2};tan#theta_{t}-(dx/dz)_{reco}");
0099 }
0100 }
0101 }
0102 }
0103
0104 unsigned LA_Filler_Fitter::find_rebin(const TH1* const hist) {
0105 const double mean = hist->GetMean();
0106 const double rms = hist->GetRMS();
0107 const int begin = std::min(1, hist->GetXaxis()->FindFixBin(mean - rms));
0108 const int end = std::max(hist->GetNbinsX(), hist->GetXaxis()->FindFixBin(mean + rms)) + 1;
0109 unsigned current_hole(0), max_hole(0);
0110 for (int i = begin; i < end; i++) {
0111 if (!hist->GetBinError(i))
0112 current_hole++;
0113 else if (current_hole) {
0114 max_hole = std::max(current_hole, max_hole);
0115 current_hole = 0;
0116 }
0117 }
0118 return max_hole + 1;
0119 }
0120
0121 TH1* LA_Filler_Fitter::rms_profile(const std::string name, const TProfile* const prof) {
0122 const int bins = prof->GetNbinsX();
0123 TH1* const rms =
0124 new TH1F(name.c_str(), "", bins, prof->GetBinLowEdge(1), prof->GetBinLowEdge(bins) + prof->GetBinWidth(bins));
0125 for (int i = 1; i <= bins; i++) {
0126 const double Me = prof->GetBinError(i);
0127 const double neff = prof->GetBinEntries(
0128 i);
0129 rms->SetBinContent(i, Me * sqrt(neff));
0130 rms->SetBinError(i, Me / sqrt(2));
0131 }
0132 return rms;
0133 }
0134
0135 TH1* LA_Filler_Fitter::subset_probability(const std::string name, const TH1* const subset, const TH1* const total) {
0136 const int bins = subset->GetNbinsX();
0137 TH1* const prob = new TH1F(
0138 name.c_str(), "", bins, subset->GetBinLowEdge(1), subset->GetBinLowEdge(bins) + subset->GetBinWidth(bins));
0139 for (int i = 1; i <= bins; i++) {
0140 const double s = subset->GetBinContent(i);
0141 const double T = total->GetBinContent(i);
0142 const double B = T - s;
0143
0144 const double p = T ? s / T : 0;
0145 const double perr = T ? ((s && B) ? sqrt(s * s * B + B * B * s) / (T * T) : 1 / T) : 0;
0146
0147 prob->SetBinContent(i, p);
0148 prob->SetBinError(i, perr);
0149 }
0150 return prob;
0151 }
