Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:14:57

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);  //Should be prof->GetBinEffectiveEntries(i);, not availible this version ROOT.  This is only ok for unweighted fills
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 }