Utilities/src/SimpleHistogramGenerator.cc

0001 #include "FastSimulation/Utilities/interface/SimpleHistogramGenerator.h"
0002 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0003
0004 #include <cmath>
0005 #include "TH1.h"
0006 // #include <iostream>
0007
0008 SimpleHistogramGenerator::SimpleHistogramGenerator(TH1* histo)
0009     :  //myHisto(histo),
0010       //theXaxis(histo->GetXaxis()),
0011       nBins(histo->GetXaxis()->GetNbins()),
0012       xMin(histo->GetXaxis()->GetXmin()),
0013       xMax(histo->GetXaxis()->GetXmax()),
0014       binWidth((xMax - xMin) / (float)nBins) {
0015   integral.reserve(nBins + 2);
0016   integral.push_back(0.);
0017   for (int i = 1; i <= nBins; ++i)
0018     integral.push_back(integral[i - 1] + histo->GetBinContent(i));
0019   integral.push_back(integral[nBins]);
0020   nEntries = integral[nBins + 1];
0021   for (int i = 1; i <= nBins; ++i)
0022     integral[i] /= nEntries;
0023 }
0024
0025 double SimpleHistogramGenerator::generate(RandomEngineAndDistribution const* random) const {
0026   // return a random number distributed according the histogram bin contents.
0027   // NB Only valid for 1-d histograms, with fixed bin width.
0028
0029   double r1 = random->flatShoot();
0030   int ibin = binarySearch(nBins, integral, r1);
0031   double x = xMin + (double)(ibin)*binWidth;
0032   if (r1 > integral[ibin])
0033     x += binWidth * (r1 - integral[ibin]) / (integral[ibin + 1] - integral[ibin]);
0034   return x;
0035 }
0036
0037 int SimpleHistogramGenerator::binarySearch(const int& n, const std::vector<float>& array, const double& value) const {
0038   // Binary search in an array of n values to locate value.
0039   //
0040   // Array is supposed  to be sorted prior to this call.
0041   // If match is found, function returns position of element.
0042   // If no match found, function gives nearest element smaller than value.
0043
0044   int nabove, nbelow, middle;
0045   nabove = n + 1;
0046   nbelow = 0;
0047   while (nabove - nbelow > 1) {
0048     middle = (nabove + nbelow) / 2;
0049     if (value == array[middle - 1])
0050       return middle - 1;
0051     if (value < array[middle - 1])
0052       nabove = middle;
0053     else
0054       nbelow = middle;
0055   }
0056   return nbelow - 1;
0057 }