OfflineValidation/src/PVValidationHelpers.cc

0001 #include "Alignment/OfflineValidation/interface/PVValidationHelpers.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include <cassert>
0004 #include <numeric>
0005 #include <algorithm>
0006 #include <iostream>
0007 #include "TMath.h"
0008 #include "TH1.h"
0009 #include "TF1.h"
0010
0011 //*************************************************************
0012 void PVValHelper::add(std::map<std::string, TH1*>& h, TH1* hist)
0013 //*************************************************************
0014 {
0015   h[hist->GetName()] = hist;
0016   hist->StatOverflows(kTRUE);
0017 }
0018
0019 //*************************************************************
0020 void PVValHelper::fill(std::map<std::string, TH1*>& h, const std::string& s, double x)
0021 //*************************************************************
0022 {
0023   if (h.count(s) == 0) {
0024     edm::LogWarning("PVValidationHelpers") << "Trying to fill non-existing Histogram named " << s << std::endl;
0025     return;
0026   }
0027   h[s]->Fill(x);
0028 }
0029
0030 //*************************************************************
0031 void PVValHelper::fill(std::map<std::string, TH1*>& h, const std::string& s, double x, double y)
0032 //*************************************************************
0033 {
0034   if (h.count(s) == 0) {
0035     edm::LogWarning("PVValidationHelpers") << "Trying to fill non-existing Histogram named " << s << std::endl;
0036     return;
0037   }
0038   h[s]->Fill(x, y);
0039 }
0040
0041 //*************************************************************
0042 void PVValHelper::fillByIndex(std::vector<TH1F*>& h, unsigned int index, double x, std::string tag)
0043 //*************************************************************
0044 {
0045   assert(!h.empty());
0046   if (index < h.size()) {
0047     h[index]->Fill(x);
0048   } else {
0049     edm::LogWarning("PVValidationHelpers") << "Trying to fill non-existing Histogram with index " << index
0050                                            << " for array with size: " << h.size() << " tag: " << tag << std::endl;
0051     return;
0052   }
0053 }
0054
0055 //*************************************************************
0056 void PVValHelper::shrinkHistVectorToFit(std::vector<TH1F*>& h, unsigned int desired_size)
0057 //*************************************************************
0058 {
0059   h.erase(h.begin() + desired_size, h.end());
0060 }
0061
0062 //*************************************************************
0063 PVValHelper::plotLabels PVValHelper::getTypeString(PVValHelper::residualType type)
0064 //*************************************************************
0065 {
0066   PVValHelper::plotLabels returnType;
0067   switch (type) {
0068       // absoulte
0069
0070     case PVValHelper::dxy:
0071       returnType = std::make_tuple("dxy", "d_{xy}", "[#mum]");
0072       break;
0073     case PVValHelper::dx:
0074       returnType = std::make_tuple("dx", "d_{x}", "[#mum]");
0075       break;
0076     case PVValHelper::dy:
0077       returnType = std::make_tuple("dy", "d_{y}", "[#mum]");
0078       break;
0079     case PVValHelper::dz:
0080       returnType = std::make_tuple("dz", "d_{z}", "[#mum]");
0081       break;
0082     case PVValHelper::IP2D:
0083       returnType = std::make_tuple("IP2D", "IP_{2D}", "[#mum]");
0084       break;
0085     case PVValHelper::resz:
0086       returnType = std::make_tuple("resz", "z_{trk}-z_{vtx}", "[#mum]");
0087       break;
0088     case PVValHelper::IP3D:
0089       returnType = std::make_tuple("IP3D", "IP_{3D}", "[#mum]");
0090       break;
0091     case PVValHelper::d3D:
0092       returnType = std::make_tuple("d3D", "d_{3D}", "[#mum]");
0093       break;
0094
0095       // normalized
0096
0097     case PVValHelper::norm_dxy:
0098       returnType = std::make_tuple("norm_dxy", "d_{xy}/#sigma_{d_{xy}}", "");
0099       break;
0100     case PVValHelper::norm_dx:
0101       returnType = std::make_tuple("norm_dx", "d_{x}/#sigma_{d_{x}}", "");
0102       break;
0103     case PVValHelper::norm_dy:
0104       returnType = std::make_tuple("norm_dy", "d_{y}/#sigma_{d_{y}}", "");
0105       break;
0106     case PVValHelper::norm_dz:
0107       returnType = std::make_tuple("norm_dz", "d_{z}/#sigma_{d_{z}}", "");
0108       break;
0109     case PVValHelper::norm_IP2D:
0110       returnType = std::make_tuple("norm_IP2D", "IP_{2D}/#sigma_{IP_{2D}}", "");
0111       break;
0112     case PVValHelper::norm_resz:
0113       returnType = std::make_tuple("norm_resz", "z_{trk}-z_{vtx}/#sigma_{res_{z}}", "");
0114       break;
0115     case PVValHelper::norm_IP3D:
0116       returnType = std::make_tuple("norm_IP3D", "IP_{3D}/#sigma_{IP_{3D}}", "");
0117       break;
0118     case PVValHelper::norm_d3D:
0119       returnType = std::make_tuple("norm_d3D", "d_{3D}/#sigma_{d_{3D}}", "");
0120       break;
0121
0122     default:
0123       edm::LogWarning("PVValidationHelpers") << " getTypeString() unknown residual type: " << type << std::endl;
0124   }
0125
0126   return returnType;
0127 }
0128
0129 //*************************************************************
0130 PVValHelper::plotLabels PVValHelper::getVarString(PVValHelper::plotVariable var)
0131 //*************************************************************
0132 {
0133   PVValHelper::plotLabels returnVar;
0134   switch (var) {
0135     case PVValHelper::phi:
0136       returnVar = std::make_tuple("phi", "#phi", "[rad]");
0137       break;
0138     case PVValHelper::eta:
0139       returnVar = std::make_tuple("eta", "#eta", "");
0140       break;
0141     case PVValHelper::pT:
0142       returnVar = std::make_tuple("pT", "p_{T}", "[GeV]");
0143       break;
0144     case PVValHelper::pTCentral:
0145       returnVar = std::make_tuple("pTCentral", "p_{T} |#eta|<1.", "[GeV]");
0146       break;
0147     case PVValHelper::ladder:
0148       returnVar = std::make_tuple("ladder", "ladder number", "");
0149       break;
0150     case PVValHelper::modZ:
0151       returnVar = std::make_tuple("modZ", "module number", "");
0152       break;
0153     default:
0154       edm::LogWarning("PVValidationHelpers") << " getVarString() unknown plot variable: " << var << std::endl;
0155   }
0156
0157   return returnVar;
0158 }
0159
0160 //*************************************************************
0161 std::vector<float> PVValHelper::generateBins(int n, float start, float range)
0162 //*************************************************************
0163 {
0164   std::vector<float> v(n);
0165   float interval = range / (n - 1);
0166   std::iota(v.begin(), v.end(), 1.);
0167
0168   /*
0169     std::cout<<" interval:"<<interval<<std::endl;
0170     for(float &a : v) { std::cout<< a << " ";  }
0171     std::cout<< "\n";
0172   */
0173
0174   std::for_each(begin(v), end(v), [&](float& a) { a = start + ((a - 1) * interval); });
0175
0176   return v;
0177 }
0178
0179 //*************************************************************
0180 Measurement1D PVValHelper::getMedian(TH1F* histo)
0181 //*************************************************************
0182 {
0183   double median = 0.;
0184   double q = 0.5;  // 0.5 quantile for "median"
0185   // protect against empty histograms
0186   if (histo->Integral() != 0) {
0187     histo->GetQuantiles(1, &median, &q);
0188   }
0189
0190   Measurement1D result(median, median / TMath::Sqrt(histo->GetEntries()));
0191
0192   return result;
0193 }
0194
0195 //*************************************************************
0196 Measurement1D PVValHelper::getMAD(TH1F* histo)
0197 //*************************************************************
0198 {
0199   int nbins = histo->GetNbinsX();
0200   double median = getMedian(histo).value();
0201   double x_lastBin = histo->GetBinLowEdge(nbins + 1);
0202   const char* HistoName = histo->GetName();
0203   std::string Finalname = "resMed_";
0204   Finalname.append(HistoName);
0205   TH1F* newHisto = new TH1F(Finalname.c_str(), Finalname.c_str(), nbins, 0., x_lastBin);
0206   double* residuals = new double[nbins];
0207   const float* weights = histo->GetArray();
0208
0209   for (int j = 0; j < nbins; j++) {
0210     residuals[j] = std::abs(median - histo->GetBinCenter(j + 1));
0211     newHisto->Fill(residuals[j], weights[j]);
0212   }
0213
0214   double theMAD = (PVValHelper::getMedian(newHisto).value()) * 1.4826;
0215
0216   delete[] residuals;
0217   residuals = nullptr;
0218   newHisto->Delete("");
0219
0220   Measurement1D result(theMAD, theMAD / histo->GetEntries());
0221   return result;
0222 }
0223
0224 //*************************************************************
0225 std::pair<Measurement1D, Measurement1D> PVValHelper::fitResiduals(TH1* hist)
0226 //*************************************************************
0227 {
0228   //float fitResult(9999);
0229   if (hist->GetEntries() < 1) {
0230     return std::make_pair(Measurement1D(0., 0.), Measurement1D(0., 0.));
0231   };
0232
0233   float mean = hist->GetMean();
0234   float sigma = hist->GetRMS();
0235
0236   TF1 func("tmp", "gaus", mean - 1.5 * sigma, mean + 1.5 * sigma);
0237   if (0 == hist->Fit(&func, "QNR")) {  // N: do not blow up file by storing fit!
0238     mean = func.GetParameter(1);
0239     sigma = func.GetParameter(2);
0240     // second fit: three sigma of first fit around mean of first fit
0241     func.SetRange(mean - 2 * sigma, mean + 2 * sigma);
0242     // I: integral gives more correct results if binning is too wide
0243     // L: Likelihood can treat empty bins correctly (if hist not weighted...)
0244     if (0 == hist->Fit(&func, "Q0LR")) {
0245       if (hist->GetFunction(func.GetName())) {  // Take care that it is later on drawn:
0246         hist->GetFunction(func.GetName())->ResetBit(TF1::kNotDraw);
0247       }
0248     }
0249   }
0250
0251   float res_mean = func.GetParameter(1);
0252   float res_width = func.GetParameter(2);
0253
0254   float res_mean_err = func.GetParError(1);
0255   float res_width_err = func.GetParError(2);
0256
0257   Measurement1D resultM(res_mean, res_mean_err);
0258   Measurement1D resultW(res_width, res_width_err);
0259
0260   std::pair<Measurement1D, Measurement1D> result;
0261
0262   result = std::make_pair(resultM, resultW);
0263   return result;
0264 }