Use log(N) search

kmohrman · Apr 6, 2022 · 267e2db · 267e2db
1 parent 34c9ca5
commit 267e2db
Show file tree

Hide file tree

Showing 2 changed files with 34 additions and 49 deletions.
diff --git a/src/RooParametricHist.cxx b/src/RooParametricHist.cxx
@@ -210,15 +210,16 @@ double RooParametricHist::evaluateMorphFunction(int j) const
 }
 double RooParametricHist::evaluatePartial() const
 {
-  int bin_i;
-  if (x < bins[0]) 	   	return 0;    // should set to 0 instead?
-  else if (x >= bins[N_bins])   return 0;
-
-  else {
-    for(bin_i=0; bin_i<N_bins; bin_i++) {   // faster way to loop through ?
-      if (x>=bins[bin_i] && x < bins[bin_i+1] ) break;
-    }
+  auto it = std::upper_bound(std::begin(bins), std::end(bins), x);
+  if ( it == std::begin(bins) ) {
+    // underflow
+    return 0;
+  }
+  else if ( it == std::end(bins) ) {
+    // overflow
+    return 0;
   }
+  size_t bin_i = std::distance(std::begin(bins), it) - 1;
   RooAbsReal *retVar = (RooAbsReal*)pars.at(bin_i);
 
   double ret = retVar->getVal();

diff --git a/src/RooParametricHist2D.cxx b/src/RooParametricHist2D.cxx
@@ -201,50 +201,34 @@ Double_t RooParametricHist2D::analyticalIntegral(Int_t code, const char* rangeNa
 
 Double_t RooParametricHist2D::evaluate() const 
 { 
+  auto itx = std::upper_bound(std::begin(bins_x), std::end(bins_x), x);
+  if ( itx == std::begin(bins_x) ) {
+    // underflow
+    return 0;
+  }
+  else if ( itx == std::end(bins_x) ) {
+    // overflow
+    return 0;
+  }
+  auto ity = std::upper_bound(std::begin(bins_y), std::end(bins_y), y);
+  if ( ity == std::begin(bins_y) ) {
+    // underflow
+    return 0;
+  }
+  else if ( ity == std::end(bins_y) ) {
+    // overflow
+    return 0;
+  }
+  size_t bin_ix = std::distance(std::begin(bins_x), itx) - 1;
+  size_t bin_iy = std::distance(std::begin(bins_y), ity) - 1;
 
-    int bin_ix;
-    int bin_iy;
-
-    // std::cout << "Evaluating..." << std::endl;
-    // std::cout << "bins_x size = " << bins_x.size() << std::endl;
-    // std::cout << "bins_x type = " << typeid(bins_x).name() << std::endl;
-    // for (auto i = bins_x.begin(); i != bins_x.end(); ++i){
-    //     std::cout << *i << ' ';
-    // }
-    // std::cout << std::endl;
-
-    // std::cout << "bins_y size = " << bins_y.size() << std::endl;
-    // std::cout << "bins_y type = " << typeid(bins_y).name() << std::endl;
-    // for (auto i = bins_y.begin(); i != bins_y.end(); ++i){
-    //     std::cout << *i << ' ';
-    // }
-    // std::cout << std::endl;
-
-    if (x < bins_x[0]) {
-        return 0; // should set to 0 instead?
-    } else if (y < bins_y[0]) {
-        return 0; 
-    } else if (x >= bins_x[N_bins_x]) {
-        return 0;
-    } else if (y >= bins_y[N_bins_y]) {
-        return 0;
-    } else {
-        for(bin_ix=0; bin_ix<N_bins_x; bin_ix++) {   // faster way to loop through ?
-            if (x>=bins_x[bin_ix] && x < bins_x[bin_ix+1] ) break;
-        }
-        for(bin_iy=0; bin_iy<N_bins_y; bin_iy++) {   // faster way to loop through ?
-            if (y>=bins_y[bin_iy] && y < bins_y[bin_iy+1] ) break;
-        }
-    }
-
-    int globalbin = N_bins_x * bin_iy + bin_ix;
+  int globalbin = N_bins_x * bin_iy + bin_ix;
 
-    RooAbsReal *retVar = (RooAbsReal*)pars.at(globalbin);
+  RooAbsReal *retVar = (RooAbsReal*)pars.at(globalbin);
 
-    double ret = retVar->getVal() / (widths_x[bin_ix]*widths_y[bin_iy]);
+  double ret = retVar->getVal() / (widths_x[bin_ix]*widths_y[bin_iy]);
 
-    cval=ret;
-    // std::cout << "Evaluated" << std::endl;
-    return ret; 
+  cval=ret;
+  return ret;
 }