Invariant-Mass Binning

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

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

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}
+#   )

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));
+}
+
+}

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