diff --git a/include/InvariantMassBinning.h b/include/InvariantMassBinning.h
new file mode 100644
index 00000000..3e63dd29
--- /dev/null
+++ b/include/InvariantMassBinning.h
@@ -0,0 +1,98 @@
+/*  YAP - Yet another PWA toolkit
+    Copyright 2015, Technische Universitaet Muenchen,
+    Authors: Daniel Greenwald, Johannes Rauch
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/// \file
+
+
+#ifndef  yap_InvariantMassBinning_h
+#define  yap_InvariantMassBinning_h
+
+#include "fwd/InvariantMassBinning.h"
+
+#include "fwd/CachedValue.h"
+#include "fwd/DataPoint.h"
+#include "fwd/Model.h"
+#include "fwd/StatusManager.h"
+
+#include "StaticDataAccessor.h"
+
+#include <memory>
+#include <vector>
+
+namespace yap {
+
+/// \brief Partitions the invariant-mass range in \f$n\f$ bins, whose (constant) lower
+/// edges are stored in `BinLowEdges_`.
+/// \details `BinLowEdges_` is a \f$(n+1)\f$-element vector whose entries partition the
+/// axis in an closed-open fashion as follows:
+/// \f[
+/// (-\infty, m_0) \cup [m_0,m_1) \cup [m_1,m_2) \dots [m_{n-1}, m_n) \cup [m_n, +\infty).
+/// \f]
+/// The `0`-th value of `BinLowEdges_` corresponds to \f$m_0\f$; the `n`-th to \f$m_n\f$.
+/// InvariantMassBinning has a function to calculate which bin the invariant mass of a
+/// ParticleCombination lies in:
+///  * In case of underflow (i.e. \f$m < m_0\f$), the bin will be set to `-1`.
+///  * In case of overflow (i.e. \f$m > m_n\f$), the bin will be set to `n`.
+///
+/// \author Paolo Di Giglio.
+class InvariantMassBinning : public StaticDataAccessor
+{
+public:
+    /// Constructor.
+    /// \param m         The owning Model.
+    /// \param low_edges Low edges of the bins; the last element is the upper edge of the last bin.
+    explicit InvariantMassBinning(Model& m, const std::vector<double>& low_edges);
+
+    /// \brief Calculate which bin the invatiant mass of the
+    /// ParticleCombination's belong to.
+    /// \details The class description explains which bin value is used
+    /// in case of overflow or underflow.
+    /// \param d  The DataPoint to calculate into.
+    /// \param sm The StatusManager to update.
+    virtual void calculate(DataPoint& d, StatusManager& sm) const override;
+
+    /// Access the partition.
+    const std::vector<double>& lowEdges() const
+    { return BinLowEdges_; }
+
+    /// Return the bin number.
+    double bin(const DataPoint& d, const std::shared_ptr<ParticleCombination>& pc) const;
+
+private:
+
+    /// Lower edges of the bins.
+    const std::vector<double> BinLowEdges_;
+
+    /// The bin the masses of the ParticleCombination's belong to.
+    std::shared_ptr<RealCachedValue> Bin_;
+};
+
+/// Check if n corresponds to an underflow value.
+/// \param n The index to check agains the bin underflow value.
+inline const bool is_underflow(int n)
+{ return n < 0; }
+
+/// Check if n corresponds to an overflow value.
+/// \param n The index to check agains the bin overflow value.
+/// \param B The partition the bin overflow value has to be evaluated from.
+inline const bool is_overflow(int n, const InvariantMassBinning& B)
+{ return n >= static_cast<int>(B.lowEdges().size()) - 1; }
+
+}
+
+#endif
diff --git a/include/fwd/InvariantMassBinning.h b/include/fwd/InvariantMassBinning.h
new file mode 100644
index 00000000..e299f4ba
--- /dev/null
+++ b/include/fwd/InvariantMassBinning.h
@@ -0,0 +1,32 @@
+/*  YAP - Yet another PWA toolkit
+    Copyright 2015, Technische Universitaet Muenchen,
+    Authors: Daniel Greenwald, Johannes Rauch
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/// \file
+/// \author Paolo Di Giglio.
+/// Contains forward declarations only.
+
+#ifndef  yap_InvariantMassBinningFwd_h
+#define  yap_InvariantMassBinningFwd_h
+
+namespace yap {
+
+class InvariantMassBinning;
+
+}
+
+#endif
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 40258c21..b6b2b610 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,52 +1,53 @@
 set(INCLUDE_DIRECTORIES
-	${YAP_SOURCE_DIR}/include)
+    ${YAP_SOURCE_DIR}/include)
 include_directories(${INCLUDE_DIRECTORIES})
 
 set(YAP_SOURCES
-	BlattWeisskopf.cxx
-	BreitWigner.cxx
-	CachedValue.cxx
-	ClebschGordan.cxx
-	DataAccessor.cxx
-	DataPartition.cxx
-	DataPoint.cxx
-	DataSet.cxx
-	DecayChannel.cxx
-	DecayingParticle.cxx
-	DecayTree.cxx
-	DecayTreeVectorIntegral.cxx
-  Filters.cxx
-	FinalStateParticle.cxx
-	Flatte.cxx
-	FreeAmplitude.cxx
-	FourMomenta.cxx
-	HelicityAngles.cxx
-	HelicityFormalism.cxx
-	ImportanceSampler.cxx
-	Integrator.cxx
-	MassRange.cxx
-	MassShape.cxx
-	MassShapeWithNominalMass.cxx
-	MeasuredBreakupMomenta.cxx
-	Model.cxx
-	ModelIntegral.cxx
-	Parameter.cxx
-	ParticleCombination.cxx
-	ParticleCombinationCache.cxx
-	Particle.cxx
-	ParticleFactory.cxx
-	PDL.cxx
-	PoleMass.cxx
-	QuantumNumbers.cxx
-	RecalculableDataAccessor.cxx
-	RelativisticBreitWigner.cxx
-	Resonance.cxx
-	SpinAmplitude.cxx
-  SpinAmplitudeCache.cxx
-  StaticDataAccessor.cxx
-	StatusManager.cxx
-	WignerD.cxx
-	ZemachFormalism.cxx
+    BlattWeisskopf.cxx
+    BreitWigner.cxx
+    CachedValue.cxx
+    ClebschGordan.cxx
+    DataAccessor.cxx
+    DataPartition.cxx
+    DataPoint.cxx
+    DataSet.cxx
+    DecayChannel.cxx
+    DecayTree.cxx
+    DecayTreeVectorIntegral.cxx
+    DecayingParticle.cxx
+    Filters.cxx
+    FinalStateParticle.cxx
+    Flatte.cxx
+    FourMomenta.cxx
+    FreeAmplitude.cxx
+    HelicityAngles.cxx
+    HelicityFormalism.cxx
+    ImportanceSampler.cxx
+    Integrator.cxx
+    InvariantMassBinning.cxx
+    MassRange.cxx
+    MassShape.cxx
+    MassShapeWithNominalMass.cxx
+    MeasuredBreakupMomenta.cxx
+    Model.cxx
+    ModelIntegral.cxx
+    PDL.cxx
+    Parameter.cxx
+    Particle.cxx
+    ParticleCombination.cxx
+    ParticleCombinationCache.cxx
+    ParticleFactory.cxx
+    PoleMass.cxx
+    QuantumNumbers.cxx
+    RecalculableDataAccessor.cxx
+    RelativisticBreitWigner.cxx
+    Resonance.cxx
+    SpinAmplitude.cxx
+    SpinAmplitudeCache.cxx
+    StaticDataAccessor.cxx
+    StatusManager.cxx
+    WignerD.cxx
+    ZemachFormalism.cxx
 )
 
 add_library(YAP SHARED ${YAP_SOURCES})
@@ -55,21 +56,21 @@ add_library(YAP SHARED ${YAP_SOURCES})
 # cmake -DLIBRARY_OUTPUT_DIRECTORY:PATH=<lib_path> <path_to_CMakeLists.tex>
 # otherwise, default LD_LIBRARY_PATH
 if(NOT DEFINED LIBRARY_OUTPUT_DIRECTORY)
-	set(LIBRARY_OUTPUT_DIRECTORY ${YAP_SOURCE_DIR}/lib/${CMAKE_BUILD_TYPE})
+    set(LIBRARY_OUTPUT_DIRECTORY ${YAP_SOURCE_DIR}/lib/${CMAKE_BUILD_TYPE})
 endif()
 
 if(NOT DEFINED INCLUDE_OUTPUT_DIRECTORY)
-	set(INCLUDE_OUTPUT_DIRECTORY ${YAP_SOURCE_DIR}/include/YAP)
+    set(INCLUDE_OUTPUT_DIRECTORY ${YAP_SOURCE_DIR}/include/YAP)
 endif()
 
 install(TARGETS YAP LIBRARY DESTINATION ${LIBRARY_OUTPUT_DIRECTORY})
 
 # Matches all the headers in ${YAPDID}/include and its subdirs
 file(GLOB_RECURSE
-	 INSTALL_INCLUDES ${YAP_SOURCE_DIR}/include/*.h)
+     INSTALL_INCLUDES ${YAP_SOURCE_DIR}/include/*.h)
 
 #message(STATUS "${INSTALL_INCLUDES}")
 
 #install(FILE       ${INSTALL_INCLUDES}
-#	    DESTINATION ${INCLUDE_OUTPUT_DIRECTORY}
-#	)
+#       DESTINATION ${INCLUDE_OUTPUT_DIRECTORY}
+#   )
diff --git a/src/InvariantMassBinning.cxx b/src/InvariantMassBinning.cxx
new file mode 100644
index 00000000..89e93c6e
--- /dev/null
+++ b/src/InvariantMassBinning.cxx
@@ -0,0 +1,66 @@
+#include "InvariantMassBinning.h"
+
+#include "CachedValue.h"
+#include "CalculationStatus.h"
+#include "DataPoint.h"
+#include "Exceptions.h"
+#include "FourMomenta.h"
+#include "Model.h"
+#include "ParticleCombination.h"
+#include "StatusManager.h"
+
+#include <algorithm>
+#include <functional>
+#include <vector>
+
+namespace yap {
+
+//-------------------------
+InvariantMassBinning::InvariantMassBinning(Model& m, const std::vector<double>& bins) :
+    StaticDataAccessor(m, equal_by_orderless_content),
+    BinLowEdges_(bins),
+    Bin_(RealCachedValue::create(*this))
+{
+    if (BinLowEdges_.size() < 2)
+        throw exceptions::Exception("At least two bin edges must be specified",
+                                    "InvariantMassBinning::InvariantMassBinning");
+
+    // use std::less_equal so that if two elements are the same, std::is_sorted
+    // will return false
+    if (!std::is_sorted(BinLowEdges_.cbegin(), BinLowEdges_.cend(), std::less_equal<double>()))
+        throw exceptions::Exception("Elements of the vector are not monotonically increasing",
+                                    "InvariantMassBinning::InvariantMassBinning");
+
+    registerWithModel();
+}
+
+//-------------------------
+void InvariantMassBinning::calculate(DataPoint& d, StatusManager& sm) const
+{
+    // set all bins to uncalculated
+    sm.set(*this, CalculationStatus::uncalculated);
+
+    // loop over particle combinations -> indices
+    for (auto& pc_i : symmetrizationIndices()) {
+        // check if calculation is necessary
+        if (sm.status(*Bin_, pc_i.second ) == CalculationStatus::uncalculated) {
+            auto invariant_mass = model()->fourMomenta()->m(d, pc_i.first);
+
+            // get the first lower bound that is greater than the invariant mass
+            auto bin = std::upper_bound(BinLowEdges_.cbegin(), BinLowEdges_.cend(), invariant_mass);
+
+            // set the value into the DataPoint
+            // NOTE: std::distance() - 1 accounts for the fact that the actual bin value
+            // is the one preceding the one that is found above
+            Bin_->setValue(std::distance(BinLowEdges_.cbegin(), bin) - 1, d, pc_i.second, sm);
+        }
+    }
+}
+
+//-------------------------
+double InvariantMassBinning::bin(const DataPoint& d, const std::shared_ptr<ParticleCombination>& pc) const
+{
+    return Bin_->value(d, symmetrizationIndex(pc));
+}
+
+}
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 9b608244..c8d74103 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -25,6 +25,7 @@ set(YAP_TEST_SOURCES
   test_HelicityAngles.cxx
   test_HelicityAngles_boostRotate.cxx
   test_integration.cxx
+  test_InvariantMassBinning.cxx
   test_Matrix.cxx
   test_swapDalitzAxes.cxx
   test_swapFinalStates.cxx
diff --git a/test/test_InvariantMassBinning.cxx b/test/test_InvariantMassBinning.cxx
new file mode 100644
index 00000000..f7cbc880
--- /dev/null
+++ b/test/test_InvariantMassBinning.cxx
@@ -0,0 +1,115 @@
+/// \file test_InvariantMassBinning.cxx
+/// \brief Test the binning of the mass axis.
+/// \author Paolo Di Giglio
+
+#include <catch.hpp>
+#include <catch_capprox.hpp>
+
+#include <logging.h>
+
+#include <DataSet.h>
+#include <DataPoint.h>
+#include <InvariantMassBinning.h>
+#include <FourMomenta.h>
+#include <MassRange.h>
+#include <Model.h>
+#include <ParticleCombination.h>
+
+#include "helperFunctions.h"
+#include <iostream>
+#include <memory>
+#include <vector>
+
+const std::vector<double> monospaced_partition(const yap::MassRange& range, std::size_t bins)
+{
+    const double bin_width = (range[1] - range[0]) / bins;
+    // Reserve the space for the partition
+    std::vector<double> p;
+    p.reserve(bins + 1);
+
+    for (std::size_t i = 0; i < bins + 1; ++i)
+        p.push_back(i * bin_width);
+
+    return p;
+}
+
+inline const std::vector<yap::MassRange> model_mass_range(std::shared_ptr<yap::Model> model)
+{
+    return yap::mass_range(
+               model->massAxes(),
+               model->initialStateParticles().begin()->first,
+               model->finalStateParticles());
+}
+
+inline const bool mass_belongs_to_bin(double mass, int bin, const std::vector<double>& bin_low_edges)
+{
+    if (bin == static_cast<int>(bin_low_edges.size()) - 1)
+        return bin_low_edges[bin] <= mass;
+    if (bin == -1)
+        return mass < bin_low_edges[0];
+
+    return bin_low_edges[bin] <= mass && mass < bin_low_edges[bin + 1];
+}
+
+class InvariantMassBinningTest : public yap::InvariantMassBinning
+{
+public:
+    InvariantMassBinningTest(yap::Model& m, const std::vector<double>& bins) :
+        yap::InvariantMassBinning(m, bins)
+    {}
+
+    void addPC(std::shared_ptr<yap::ParticleCombination> pc)
+    { yap::InvariantMassBinning::addParticleCombination(pc); }
+};
+
+TEST_CASE( "InvariantMassBinning" )
+{
+    // disable debug logs in test
+    yap::disableLogs(el::Level::Debug);
+    //yap::plainLogs(el::Level::Debug);
+
+    // Create the model
+    auto model = create_model();
+    REQUIRE(model->consistent());
+
+    // Get the allowed mass range for the model mass axes
+    auto m_range = model_mass_range(model);
+//    for (const auto& it : m_range)
+//        std::cerr << it[0] << " " << it[1] << std::endl;
+
+    // Number of bins
+    const unsigned int nBins = 10;
+    // Vector of the low edges of the partition
+    auto bin_low_edges = monospaced_partition(*m_range.begin(), nBins);
+    // Create the binning
+    InvariantMassBinningTest binning(*model, bin_low_edges);
+
+    // Loop over particle combinations -> indices
+    for (const auto& pc_i : model->fourMomenta()->symmetrizationIndices()) {
+        // Add particle combinations to the binning
+        binning.addPC(pc_i.first);
+    }
+
+    // Generate the DataSet
+    const unsigned int nPoints = 1000;
+    auto data_set = generate_data(model, nPoints);
+
+    // Calculate and check the bins
+    for (const auto& pc_i : model->fourMomenta()->symmetrizationIndices()) {
+        for (auto& dp : data_set) {
+            // Calculate the bin
+            binning.calculate(dp, data_set);
+            // Check the bin
+            auto mass = model->fourMomenta()->m(dp, pc_i.first);
+            auto bin  = binning.bin(dp, pc_i.first);
+            auto bin_is_correct = mass_belongs_to_bin(mass, bin, bin_low_edges);
+            if (!bin_is_correct)
+                std::cout << "bin = " << bin << " "
+                          << bin_low_edges[bin] << " <= "
+                          << mass << " < "
+                          << bin_low_edges[bin + 1] << std::endl;
+
+            REQUIRE(bin_is_correct);
+        }
+    }
+}