[Rivet-svn] r4124 - in trunk: . include/Rivet/Math

[blackhole at projects.hepforge.org]

Author: dgrell
Date: Tue Feb  5 10:18:28 2013
New Revision: 4124

Log:
Added Breit-Wigner binning helper

Modified:
   trunk/ChangeLog
   trunk/include/Rivet/Math/MathUtils.hh

Modified: trunk/ChangeLog
==============================================================================
--- trunk/ChangeLog	Mon Feb  4 10:13:28 2013	(r4123)
+++ trunk/ChangeLog	Tue Feb  5 10:18:28 2013	(r4124)
@@ -1,3 +1,8 @@
+2013-02-05  David Grellscheid  <David.Grellscheid at durham.ac.uk>
+
+	* include/Rivet/Math/MathUtils.hh: added BWspace bin edge method
+	to give equal-area Breit-Wigner bins
+
 2013-02-01  Andy Buckley  <andy.buckley at cern.ch>
 
 	* Adding an element to the PhiMapping enum and a new mapAngle(angle, mapping) function.

Modified: trunk/include/Rivet/Math/MathUtils.hh
==============================================================================
--- trunk/include/Rivet/Math/MathUtils.hh	Mon Feb  4 10:13:28 2013	(r4123)
+++ trunk/include/Rivet/Math/MathUtils.hh	Tue Feb  5 10:18:28 2013	(r4124)
@@ -248,6 +248,44 @@
     return rtn;
   }
 
+  namespace BWHelpers {
+    /// @brief CDF for the Breit-Wigner distribution
+    inline double CDF(double x, double mu, double gamma) {
+      // normalize to (0;1) distribution
+      const double xn = (x - mu)/gamma;
+      return std::atan(xn)/M_PI + 0.5;
+    }
+
+    /// @brief inverse CDF for the Breit-Wigner distribution
+    inline double antiCDF(double p, double mu, double gamma) {
+      const double xn = std::tan(M_PI*(p-0.5));
+      return gamma*xn + mu;
+    }
+  }
+
+  /// @brief Make a list of @a nbins + 1 values spaced for equal area
+  /// Breit-Wigner binning between @a start and @a end inclusive. @a
+  /// mu and @a gamma are the Breit-Wigner parameters.
+  ///
+  /// NB. The arg ordering and the meaning of the nbins variable is "histogram-like",
+  /// as opposed to the Numpy/Matlab version, and the start and end arguments are expressed
+  /// in "normal" space.
+  inline vector<double> BWspace(size_t nbins, double start, double end, 
+				double mu, double gamma) {
+    assert(end >= start);
+    assert(nbins > 0);
+    const double pmin = BWHelpers::CDF(start,mu,gamma);
+    const double pmax = BWHelpers::CDF(end,  mu,gamma);
+    const vector<double> edges = linspace(nbins, pmin, pmax);
+    assert(edges.size() == nbins+1);
+    vector<double> rtn;
+    foreach (double edge, edges) {
+      rtn.push_back(BWHelpers::antiCDF(edge,mu,gamma));
+    }
+    assert(rtn.size() == nbins+1);
+    return rtn;
+  }
+
 
   /// @brief Return the bin index of the given value, @a val, given a vector of bin edges
   ///