TrackerMaterialAnalysis/plugins/XHistogram.cc

0001 #include <algorithm>
0002 #include <cmath>
0003 #include <vector>
0004
0005 #include "XHistogram.h"
0006
0007 std::vector<XHistogram::position> XHistogram::splitSegment(Range rangeX, Range rangeY) const {
0008   double deltaX = rangeX.second - rangeX.first;
0009   double deltaY = rangeY.second - rangeY.first;
0010   double length = hypot(deltaX, deltaY);
0011   double stepX = (m_xRange.second - m_xRange.first) / m_xBins;
0012   double stepY = (m_yRange.second - m_yRange.first) / m_yBins;
0013
0014   int min_i, max_i, min_j, max_j;
0015   if (rangeX.first < rangeX.second) {
0016     min_i = (int)ceil(rangeX.first / stepX);        // included
0017     max_i = (int)floor(rangeX.second / stepX) + 1;  // excluded
0018   } else {
0019     min_i = (int)ceil(rangeX.second / stepX);
0020     max_i = (int)floor(rangeX.first / stepX) + 1;
0021   }
0022   if (rangeY.first < rangeY.second) {
0023     min_j = (int)ceil(rangeY.first / stepY);
0024     max_j = (int)floor(rangeY.second / stepY) + 1;
0025   } else {
0026     min_j = (int)ceil(rangeY.second / stepY);
0027     max_j = (int)floor(rangeY.first / stepY) + 1;
0028   }
0029
0030   int steps = max_i - min_i + max_j - min_j + 2;
0031   std::vector<position> v;
0032   v.clear();
0033   v.reserve(steps);
0034
0035   v.push_back(position(0., rangeX.first, rangeY.first));
0036   double x, y, f;
0037   for (int i = min_i; i < max_i; ++i) {
0038     x = i * stepX;
0039     y = rangeY.first + (x - rangeX.first) * deltaY / deltaX;
0040     f = std::fabs((x - rangeX.first) / deltaX);
0041     v.push_back(position(f, x, y));
0042   }
0043   for (int i = min_j; i < max_j; ++i) {
0044     y = i * stepY;
0045     x = rangeX.first + (y - rangeY.first) * deltaX / deltaY;
0046     f = std::fabs((y - rangeY.first) / deltaY);
0047     v.push_back(position(f, x, y));
0048   }
0049   v.push_back(position(1., rangeX.second, rangeY.second));
0050
0051   // sort by distance from the start of the segment
0052   std::sort(v.begin(), v.end());
0053
0054   // filter away the fragments shorter than m_minDl, and save the center of each fragment along with its fractionary length
0055   std::vector<position> result;
0056   result.push_back(v.front());
0057   for (int i = 1, s = v.size(); i < s; ++i) {
0058     double mx = (v[i].x + v[i - 1].x) / 2.;
0059     double my = (v[i].y + v[i - 1].y) / 2.;
0060     double df = (v[i].f - v[i - 1].f);
0061     if (df * length < m_minDl)
0062       continue;
0063     result.push_back(position(df, mx, my));
0064   }
0065
0066   return result;
0067 }
0068
0069 /// fill one point
0070 void XHistogram::fill(double x, double y, const std::vector<double>& weight, double norm) {
0071   check_weight(weight);
0072
0073   for (size_t h = 0; h < m_size; ++h)
0074     m_histograms[h]->Fill(x, y, weight[h]);
0075   m_normalization->Fill(x, y, norm);
0076 }
0077
0078 /// fill one point and set its color
0079 void XHistogram::fill(double x, double y, const std::vector<double>& weight, double norm, unsigned int colour) {
0080   check_weight(weight);
0081
0082   for (size_t h = 0; h < m_size; ++h)
0083     m_histograms[h]->Fill(x, y, weight[h]);
0084   m_normalization->Fill(x, y, norm);
0085   m_colormap->SetBinContent(m_colormap->FindBin(x, y), (float)colour);
0086 }
0087
0088 /// fill one segment, normalizing each bin's weight to the fraction of the segment it contains
0089 void XHistogram::fill(const Range& x, const Range& y, const std::vector<double>& weight, double norm) {
0090   check_weight(weight);
0091
0092   std::vector<position> v = splitSegment(x, y);
0093   for (size_t i = 0, s = v.size(); i < s; ++i) {
0094     for (size_t h = 0; h < m_size; ++h)
0095       m_histograms[h]->Fill(v[i].x, v[i].y, v[i].f * weight[h]);
0096     m_normalization->Fill(v[i].x, v[i].y, v[i].f * norm);
0097   }
0098 }
0099
0100 /// fill one segment and set its color, normalizing each bin's weight to the fraction of the segment it contains
0101 void XHistogram::fill(
0102     const Range& x, const Range& y, const std::vector<double>& weight, double norm, unsigned int colour) {
0103   check_weight(weight);
0104
0105   std::vector<position> v = splitSegment(x, y);
0106   for (size_t i = 0, s = v.size(); i < s; ++i) {
0107     for (size_t h = 0; h < m_size; ++h)
0108       m_histograms[h]->Fill(v[i].x, v[i].y, v[i].f * weight[h]);
0109     m_normalization->Fill(v[i].x, v[i].y, v[i].f * norm);
0110     m_colormap->SetBinContent(m_colormap->FindBin(v[i].x, v[i].y), (float)colour);
0111   }
0112 }
0113
0114 /// normalize the histograms
0115 void XHistogram::normalize(void) {
0116   for (int i = 0; i < m_normalization->GetSize(); ++i) {
0117     if ((*m_normalization)[i] > 0.) {
0118       for (size_t h = 0; h < m_size; ++h)
0119         (*m_histograms[h])[i] /= (*m_normalization)[i];
0120       (*m_normalization)[i] = 1.;
0121     }
0122   }
0123 }