Page MenuHomeHEPForge
Files F9501216

No OneTemporary
Actions

View Options

diff --git a/doc/release.notes b/doc/release.notes
index b7cc81d..d9f5cc7 100644
--- a/doc/release.notes
+++ b/doc/release.notes
@@ -1,94 +1,103 @@
///////////////////////////////////////////////////////////////
/// ///
/// This is the History file for the Laura++ package. ///
/// ///
///////////////////////////////////////////////////////////////
+27th January 2014 Daniel Craik
+
+* Added spline interpolation to DP backgrounds
+ - Added Lau2DSplineDPPdf which uses a spline to model a normalised PDF across a DP
+ - Added pABC, Lau2DAbsDPPdf, for Lau2DHistDPPdf and Lau2DSplineDPPdf and moved common
+ code in Lau2DHistDPPdf and Lau2DSplineDPPdf into ABC Lau2DAbsHistDPPdf
+ - LauBkgndModel, modified to use Lau2DAbsDPPdf instead of Lau2DHistDPPdf
+ - setBkgndSpline method added to LauBkgndDPModel to allow use of splines
+
16th January 2014 Thomas Latham
* Add LauRealImagCPCoeffSet, which provides an (x,y), (xbar,ybar) way of
parametrising the complex coefficients.
* Try to improve timing in the *CoeffSet classes by making the complex coeffs
into member variables.
20th December 2013 Daniel Craik
* Added Lau2DCubicSpline which provides cubic spline interpolation of a histogram
- Added Lau2DSplineDP which uses a spline to model variation across a DP (eg efficiency)
- Added pABC, Lau2DAbsDP, for Lau2DHistDP and Lau2DSplineDP and moved common code
in Lau2DHistDP and Lau2DSplineDP into ABC Lau2DAbsHistDP
- LauEffModel, LauDPDepSumPdf and LauDPDepMapPdf modified to use Lau2DAbsDP instead of
Lau2DHistDP
- setEffSpline method added to LauEffModel and useSpline flag added to constructor for
LauDPDepSumPdf to allow use of splines
5th December 2013 Thomas Latham
* Fix small bug in GenFitKpipi example where background asymmetry parameter had
its limits the wrong way around
4th December 2013 Daniel Craik
* Updated 2D chi-squared code to use adaptive binning.
3rd December 2013 Thomas Latham
* Generalise the Makefile in the examples directory
- results in minor changes to the names of 3 of the binaries
29th November 2013 Thomas Latham
* Fixed bug in ResultsExtractor where the output file was not written
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Laura++ v1r1p1
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
22nd November 2013 Thomas Latham
* Fixed bug in LauCPFitModel where values of q = -1 extra PDFs
were used regardless of the event charge.
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Laura++ v1r1
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
20th November 2013 Mark Whitehead
* Changed convention for the barrier factors, swapping from p* to p.
This seems to give more physically reasonable results.
The files updated are:
- src/LauGounarisSakuraiRes.cc
- src/LauRelBreitWignerRes.cc
18th October 2013 Thomas Latham
* Fix dependency problem in Makefile
8th October 2013 Thomas Latham
* Some fixes to yield implementation
* Minor bug fix in DP background histogram class
7th October 2013 Mark Whitehead
* Update to accept the yields and yield asymmetries as LauParameters.
All examples have been updated to reflect this change.
This updated the following files:
- inc/LauCPFitModel.hh
- inc/LauSimpleFitModel.hh
- inc/LauAbsFitModel.hh
- src/LauCPFitModel.cc
- src/LauSimpleFitModel.cc
* Addition of the following particles to src/LauResonanceMaker.cc
Ds*+-, Ds0*(2317)+-, Ds2*(2573)+-, Ds1*(2700)+- and Bs*0
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Laura++ v1r0
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
13th September 2013 Thomas Latham
* Initial import of the package into HEPforge
diff --git a/inc/Lau2DAbsDPPdf.hh b/inc/Lau2DAbsDPPdf.hh
new file mode 100644
index 0000000..0dba64f
--- /dev/null
+++ b/inc/Lau2DAbsDPPdf.hh
@@ -0,0 +1,67 @@
+
+// Copyright University of Warwick 2004 - 2013.
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Authors:
+// Thomas Latham
+// John Back
+// Paul Harrison
+
+/*! \file Lau2DAbsDPPdf.hh
+ \brief File containing declaration of Lau2DAbsDPPdf class.
+*/
+
+/*! \class Lau2DAbsDPPdf
+ \brief Pure abstract base class for defining a variation across a 2D DP.
+
+ Pure abstract base class for defining a normalised variation across a 2D DP.
+*/
+
+#ifndef LAU_2DABS_DP_PDF
+#define LAU_2DABS_DP_PDF
+
+#include "Rtypes.h"
+
+class Lau2DAbsDPPdf {
+
+ public:
+ //! Constructor
+ Lau2DAbsDPPdf() {}
+
+ //! Destructor
+ virtual ~Lau2DAbsDPPdf() {}
+
+ //! Retrieve maximum height
+ virtual Double_t getMaxHeight() const=0;
+
+ //! Retrieve PDF normalisation
+ virtual Double_t getHistNorm() const=0;
+
+ //! Perform the interpolation (unnormalised)
+ /*!
+ \param [in] x the x-axis value
+ \param [in] y the y-aixs value
+ \return the unnormalised PDF value
+ */
+ virtual Double_t interpolateXY(Double_t x, Double_t y) const=0;
+
+ //! Perform the interpolation and divide by the normalisation
+ /*!
+ \param [in] x the x-axis abscissa value
+ \param [in] y the y-axis abscissa value
+ \return the normalised PDF value
+ */
+ virtual Double_t interpolateXYNorm(Double_t x, Double_t y) const=0;
+
+ private:
+ //! Copy constructor - not implemented
+ Lau2DAbsDPPdf( const Lau2DAbsDPPdf& rhs );
+
+ //! Copy assignment operator - not implemented
+ Lau2DAbsDPPdf& operator=(const Lau2DAbsDPPdf& rhs);
+
+ ClassDef(Lau2DAbsDPPdf,0) // Abstract base class for 2D DP variation
+};
+
+#endif
diff --git a/inc/Lau2DAbsHistDPPdf.hh b/inc/Lau2DAbsHistDPPdf.hh
new file mode 100644
index 0000000..d647567
--- /dev/null
+++ b/inc/Lau2DAbsHistDPPdf.hh
@@ -0,0 +1,149 @@
+
+// Copyright University of Warwick 2004 - 2013.
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Authors:
+// Thomas Latham
+// John Back
+// Paul Harrison
+
+/*! \file Lau2DAbsHistDPPdf.hh
+ \brief File containing declaration of Lau2DAbsHistDPPdf class.
+*/
+
+/*! \class Lau2DAbsHistDPPdf
+ \brief Abstract base class for defining a variation across a 2D DP based on a histogram.
+
+ Abstract base class for defining a normalised variation across a 2D DP based on a histogram.
+ The returned values are normalised to the total area.
+ The histogram can be defined in the conventional DP (m13Sq vs m23Sq) or in the square DP and
+ the upper half can be used to describe symmetric DPs.
+*/
+
+#ifndef LAU_2DABSHIST_DP_PDF
+#define LAU_2DABSHIST_DP_PDF
+
+#include "Lau2DAbsDPPdf.hh"
+
+class TH2;
+class LauKinematics;
+class LauVetoes;
+
+class Lau2DAbsHistDPPdf : public Lau2DAbsDPPdf {
+
+ public:
+ //! Constructor
+ /*!
+ \param [in] kinematics the current DP kinematics
+ \param [in] vetoes the vetoes within the DP
+ \param [in] useUpperHalfOnly boolean flag to specify that the supplied histogram contains only the upper half of a symmetric DP
+ \param [in] squareDP boolean flag to specify whether the supplied histogram is in square DP coordinates
+ */
+ Lau2DAbsHistDPPdf(LauKinematics* kinematics, const LauVetoes* vetoes,
+ Bool_t useUpperHalfOnly = kFALSE, Bool_t squareDP = kFALSE);
+
+ //! Destructor
+ virtual ~Lau2DAbsHistDPPdf();
+
+ //! Retrieve maximum height of the PDF
+ /*!
+ \return the maximum height
+ */
+ Double_t getMaxHeight() const {return maxHeight_;}
+
+ //! Retrieve PDF normalisation
+ /*!
+ \return the normalisation factor
+ */
+ virtual Double_t getHistNorm() const=0;
+
+ //! Perform the interpolation (unnormalised)
+ /*!
+ \param [in] x the x-axis value
+ \param [in] y the y-aixs value
+ \return the unnormalised PDF value
+ */
+ virtual Double_t interpolateXY(Double_t x, Double_t y) const=0;
+
+ //! Perform the interpolation and divide by the normalisation
+ /*!
+ \param [in] x the x-axis abscissa value
+ \param [in] y the y-axis abscissa value
+ \return the normalised PDF value
+ */
+ virtual Double_t interpolateXYNorm(Double_t x, Double_t y) const=0;
+
+ protected:
+ //! Get the kinematics object
+ /*!
+ \return the kinematics
+ */
+ const LauKinematics * getKinematics() const {return kinematics_;}
+
+ //! Get the vetoes object
+ /*!
+ \return the vetoes
+ */
+ const LauVetoes * getVetoes() const {return vetoes_;}
+
+ //! Calculate maximum height
+ /*!
+ \param [in,out] hist the histogram
+ */
+ void calcMaxHeight(TH2* hist);
+
+ //! Fluctuate the histogram bin contents in accordance with their errors
+ /*!
+ \param [in,out] hist the histogram
+ */
+ void doBinFluctuation(TH2* hist);
+
+ //! Check whether the given co-oridinates are within the kinematic boundary
+ /*!
+ \param [in] x the x co-ordinate
+ \param [in] y the y co-ordinate
+ \return true if the co-oridinates are within the kinematic boundary, otherwise false
+ */
+ Bool_t withinDPBoundaries(Double_t x, Double_t y) const;
+
+ //! Update the current co-ordinates in the kinematic space
+ /*!
+ \param [in] x the x co-ordinate
+ \param [in] y the y co-ordinate
+ */
+ void updateKinematics(Double_t x, Double_t y) const;
+
+ //! If only using the upper half of the (symmetric) DP then transform into the correct half
+ /*!
+ \param [in,out] x the x co-ordinate
+ \param [in,out] y the y co-ordinate
+ */
+ void getUpperHalf(Double_t& x, Double_t& y) const;
+
+
+ private:
+ //! Copy constructor - not implemented
+ Lau2DAbsHistDPPdf( const Lau2DAbsHistDPPdf& rhs );
+
+ //! Copy assignment operator - not implemented
+ Lau2DAbsHistDPPdf& operator=(const Lau2DAbsHistDPPdf& rhs);
+
+ //! DP kinematics
+ LauKinematics* kinematics_;
+
+ //! Vetos within DP
+ const LauVetoes* vetoes_;
+
+ //! The maximum height of 2D histogram
+ Double_t maxHeight_;
+
+ //! Boolean for using the upper half of DP
+ Bool_t upperHalf_;
+ //! Boolean for using square DP variables
+ Bool_t squareDP_;
+
+ ClassDef(Lau2DAbsHistDPPdf,0) // Abstract base class for 2D DP variations based on a histogram
+};
+
+#endif
diff --git a/inc/Lau2DCubicSpline.hh b/inc/Lau2DCubicSpline.hh
index 87f6be7..268ee6a 100644
--- a/inc/Lau2DCubicSpline.hh
+++ b/inc/Lau2DCubicSpline.hh
@@ -1,201 +1,204 @@
// Copyright University of Warwick 2006 - 2013.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
// Authors:
// Thomas Latham
// John Back
// Paul Harrison
/*! \file Lau2DCubicSpline.hh
\brief File containing declaration of Lau2DCubicSpline class.
*/
/*! \class Lau2DCubicSpline
\brief Class for defining a 2D cubic spline based on an input histogram
Class for defining a 2D cubic spline based on RooBinned2DBicubicBase by
Manuel Tobias Schiller <manuel.schiller@nikhef.nl> (2012-08-29).
Smoothly interpolate between bin midpoints of a 2D histogram with
equal-sized bins
The interpolation function coincides with the saved histogram values in bin
midpoints; a cubic interpolation polynomial is used. Bins must all have the
same size. The interpolation polynomial in a "cell" between four bin
midpoints is:
\f[ p(x, y) = \sum_{i=0}^3 \sum_{j=0}^3 a_{ij} x^i y^j \f]
The coefficients \f$a_{ij}\f$ are determined by the requirement that
function value, first derivatives and the mixed second derivative must agree
with that of the binned function at the four bin midpoints at the edges of
the cell. Symmetric finite differences are used to approximate these
derivatives.
For each cell, the coefficients are determined at constuction time. The
object also keeps a cache of 2D integrals over complete cells, such that 2D
integrations can be done analytically in a reasonably short amount of time.
*/
#ifndef LAU_2DCUBICSPLINE
#define LAU_2DCUBICSPLINE
class TH2;
class Lau2DCubicSpline
{
private:
//! Length of coefficient record in array
enum { NCoeff = 16, CoeffRecLen = 17 };
public:
//! Constructor from histogram
/*!
/param h the histogram
*/
Lau2DCubicSpline(const TH2& h);
//! Destructor
virtual ~Lau2DCubicSpline();
//! Evaluate the function at given point
/*!
\param [in] x the x co-ordinate
\param [in] y the y co-ordinate
*/
virtual Double_t evaluate(Double_t x, Double_t y) const;
//! Evaluate analytical integral in x, y, or x and y
/*!
\param [in] x1 the lower x limit
\param [in] x2 the upper x limit
\param [in] y1 the lower y limit
\param [in] y2 the upper y limit
*/
virtual Double_t analyticalIntegral(Double_t x1, Double_t x2, Double_t y1, Double_t y2) const;
+
+ //! Evaluate analytical integral across the whole function
+ virtual Double_t analyticalIntegral() const;
private:
//! Copy constructor - not implemented
Lau2DCubicSpline( const Lau2DCubicSpline& rhs );
//! Copy assignment operator - not implemented
Lau2DCubicSpline& operator=(const Lau2DCubicSpline& rhs);
//! Number of bins in x axis
Int_t nBinsX;
//! Number of bins in y axis
Int_t nBinsY;
//! Bin size in x
Double_t binSizeX;
//! Bin size in y
Double_t binSizeY;
//! Minimum x value
Double_t xmin;
//! Maximum x value
Double_t xmax;
//! Minimum y value
Double_t ymin;
//! Maximum y value
Double_t ymax;
//! Coefficients of interpolation polynomials
std::vector<Double_t> coeffs;
//! Get contents of a given bin from a histogram
/*!
\param [in] h the histogram
\param [in] xbin the x bin index
\param [in] ybin the y bin index
\return the bin contents
*/
Double_t histcont(const TH2& h, Int_t xbin, Int_t ybin) const;
//! Get d/dx finite difference in a given bin from a histogram
/*!
\param [in] h the histogram
\param [in] xbin the x bin index
\param [in] ybin the y bin index
\return the d/dx finite difference
*/
Double_t dhistdx(const TH2& h, Int_t xbin, Int_t ybin) const;
//! Get d/dy finite difference in a given bin from a histogram
/*!
\param [in] h the histogram
\param [in] xbin the x bin index
\param [in] ybin the y bin index
\return the d/dy finite difference
*/
Double_t dhistdy(const TH2& h, Int_t xbin, Int_t ybin) const;
//! Get d^2/dydx finite difference in a given bin from a histogram
/*!
\param [in] h the histogram
\param [in] xbin the x bin index
\param [in] ybin the y bin index
\return the d^2/dydx finite difference
*/
Double_t d2histdxdy(const TH2& h, Int_t xbin, Int_t ybin) const;
//! Const access to coefficients
/*!
\param [in] binx the x bin index
\param [in] biny the y bin index
\param [in] theCoeff the coefficient index
\return the coefficient
*/
inline const Double_t& coeff(Int_t binx, Int_t biny, Int_t theCoeff) const
{ return coeffs[theCoeff + CoeffRecLen * (binx + nBinsX * biny)]; }
//! Access to coefficients
/*!
\param [in] binx the x bin index
\param [in] biny the y bin index
\param [in] theCoeff the coefficient index
\return the coefficient
*/
inline Double_t& coeff(Int_t binx, Int_t biny, Int_t theCoeff)
{ return coeffs[theCoeff + CoeffRecLen * (binx + nBinsX * biny)]; }
//! Evaluate integral over x at a given y from (x1, y) to (x2, y)
/*!
\param [in] x1 lower x limit
\param [in] x2 upper x limit
\param [in] y y value
\return integral over x
*/
Double_t evalX(Double_t x1, Double_t x2, Double_t y) const;
//! Evaluate integral over y at a given x from (x, y1) to (x, y2)
/*!
\param [in] x x value
\param [in] y1 lower y limit
\param [in] y2 upper y limit
\return integral over y
*/
Double_t evalY(Double_t x, Double_t y1, Double_t y2) const;
//! Evaluate integral over x and y from (x1, y1) to (x2, y2)
/*!
\param [in] x1 lower x limit
\param [in] x2 upper x limit
\param [in] y1 lower y limit
\param [in] y2 upper y limit
\return integral over x and y
*/
Double_t evalXY(Double_t x1, Double_t x2, Double_t y1, Double_t y2) const;
- ClassDef(Lau2DCubicSpline, 0);
+ ClassDef(Lau2DCubicSpline, 0); // Class for defining a 2D cubic spline based on an input histogram
};
#endif
diff --git a/inc/Lau2DHistDPPdf.hh b/inc/Lau2DHistDPPdf.hh
index 72bfd12..4d20b7e 100644
--- a/inc/Lau2DHistDPPdf.hh
+++ b/inc/Lau2DHistDPPdf.hh
@@ -1,149 +1,139 @@
// Copyright University of Warwick 2004 - 2013.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
// Authors:
// Thomas Latham
// John Back
// Paul Harrison
/*! \file Lau2DHistDPPdf.hh
\brief File containing declaration of Lau2DHistDPPdf class.
*/
/*! \class Lau2DHistDPPdf
\brief Class for defining a 2D DP histogram PDF.
Class for defining a 2D DP histogram PDF.
Employs linear interpolation to get the histogram value based on how far away a point in (x,y)
is to nearby bin centres. The returned values are normalised to the total area.
The histogram can be defined in the conventional DP (m13Sq vs m23Sq) or in the square DP.
*/
#ifndef LAU_2DHIST_DP_PDF
#define LAU_2DHIST_DP_PDF
-#include "Rtypes.h"
+#include "Lau2DAbsHistDPPdf.hh"
class TH2;
class LauKinematics;
class LauVetoes;
-class Lau2DHistDPPdf {
+class Lau2DHistDPPdf : public Lau2DAbsHistDPPdf {
public:
//! Constructor
/*!
\param [in] hist the 2D DP histogram
\param [in] kinematics the current DP kinematics
\param [in] vetoes the vetoes within the DP
\param [in] useInterpolation boolean flag to determine whether linear interpolation between bins should be used or simply the raw bin values
\param [in] fluctuateBins boolean flag to determine whether the bin contents should be fluctuated in accordance with their errors (useful for systematic error evaluation)
\param [in] useUpperHalfOnly boolean flag to specify that the supplied histogram contains only the upper half of a symmetric DP
\param [in] squareDP boolean flag to specify whether the supplied histogram is in square DP coordinates
*/
Lau2DHistDPPdf(const TH2* hist, LauKinematics* kinematics, const LauVetoes* vetoes,
Bool_t useInterpolation = kTRUE, Bool_t fluctuateBins = kFALSE,
Bool_t useUpperHalfOnly = kFALSE, Bool_t squareDP = kFALSE);
//! Destructor
virtual ~Lau2DHistDPPdf();
- //! Retrieve maximum height
- Double_t getMaxHeight() const {return maxHeight_;}
-
//! Retrieve PDF normalisation
+ /*!
+ \return the normalisation factor
+ */
Double_t getHistNorm() const {return norm_;}
//! Perform the interpolation (unnormalised)
/*!
\param [in] x the x-axis value
\param [in] y the y-aixs value
\return the unnormalised PDF value
*/
Double_t interpolateXY(Double_t x, Double_t y) const;
//! Perform the interpolation and divide by the normalisation
/*!
\param [in] x the x-axis abscissa value
\param [in] y the y-axis abscissa value
\return the normalised PDF value
*/
Double_t interpolateXYNorm(Double_t x, Double_t y) const;
protected:
//! Get the bin content from the histogram
/*!
\param [in] xBinNo the x-axis bin number
\param [in] yBinNo the y-axis bin number
\return the bin conent
*/
Double_t getBinHistValue(Int_t xBinNo, Int_t yBinNo) const;
- //! Calculate maximum height
- void calcMaxHeight();
-
//! Calculate the PDF normalisation
void calcHistNorm();
-
+
//! Check the normalisation calculation
void checkNormalisation();
- //! Fluctuate the histogram bin contents in accordance with their errors
- void doBinFluctuation();
-
private:
+ //! Copy constructor - not implemented
+ Lau2DHistDPPdf( const Lau2DHistDPPdf& rhs );
+
+ //! Copy assignment operator - not implemented
+ Lau2DHistDPPdf& operator=(const Lau2DHistDPPdf& rhs);
+
//! The underlying histogram
TH2* hist_;
- //! DP kinematics
- LauKinematics* kinematics_;
-
- //! Vetos within DP
- const LauVetoes* vetoes_;
-
//! The histogram x-axis minimum
Double_t minX_;
//! The histogram x-axis maximum
Double_t maxX_;
//! The histogram y-axis minimum
Double_t minY_;
//! The histogram y-axis maximum
Double_t maxY_;
//! The histogram x-axis range
Double_t rangeX_;
//! The histogram y-axis range
Double_t rangeY_;
//! The histogram x-axis bin width
Double_t binXWidth_;
//! The histogram y-axis bin width
Double_t binYWidth_;
//! The histogram x-axis inverse bin width
Double_t invBinXWidth_;
//! The histogram y-axis inverse bin width
Double_t invBinYWidth_;
//! The maximum height of 2D histogram
Double_t maxHeight_;
//! The number of bins on the x-axis of the histogram
Int_t nBinsX_;
//! The number of bins on the y-axis of the histogram
Int_t nBinsY_;
//! The histogram normalisation
Double_t norm_;
//! Control boolean for using the linear interpolation
Bool_t useInterpolation_;
- //! Boolean for using the upper half of DP
- Bool_t upperHalf_;
- //! Boolean for using square DP variables
- Bool_t squareDP_;
ClassDef(Lau2DHistDPPdf,0) // 2D Histogram utility class for DP analyses
};
#endif
diff --git a/inc/Lau2DSplineDPPdf.hh b/inc/Lau2DSplineDPPdf.hh
new file mode 100644
index 0000000..1595154
--- /dev/null
+++ b/inc/Lau2DSplineDPPdf.hh
@@ -0,0 +1,91 @@
+
+// Copyright University of Warwick 2004 - 2013.
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Authors:
+// Thomas Latham
+// John Back
+// Paul Harrison
+
+/*! \file Lau2DSplineDPPdf.hh
+ \brief File containing declaration of Lau2DSplineDPPdf class.
+*/
+
+/*! \class Lau2DSplineDPPdf
+ \brief Class for defining a 2D DP spline PDF.
+
+ Class for defining a 2D DP spline PDF.
+ Employs a 2D cubic spline to get the histogram value based on an input histogram.
+ The returned values are normalised to the total area.
+ The histogram can be defined in the conventional DP (m13Sq vs m23Sq) or in the square DP.
+*/
+
+#ifndef LAU_2DSPLINE_DP_PDF
+#define LAU_2DSPLINE_DP_PDF
+
+#include "Lau2DAbsHistDPPdf.hh"
+
+class TH2;
+class Lau2DCubicSpline;
+class LauKinematics;
+class LauVetoes;
+
+class Lau2DSplineDPPdf: public Lau2DAbsHistDPPdf {
+
+ public:
+ //! Constructor
+ /*!
+ \param [in] hist the 2D DP histogram
+ \param [in] kinematics the current DP kinematics
+ \param [in] vetoes the vetoes within the DP
+ \param [in] fluctuateBins boolean flag to determine whether the bin contents should be fluctuated in accordance with their errors (useful for systematic error evaluation)
+ \param [in] useUpperHalfOnly boolean flag to specify that the supplied histogram contains only the upper half of a symmetric DP
+ \param [in] squareDP boolean flag to specify whether the supplied histogram is in square DP coordinates
+ */
+ Lau2DSplineDPPdf(const TH2* hist, LauKinematics* kinematics, const LauVetoes* vetoes,
+ Bool_t fluctuateBins = kFALSE, Bool_t useUpperHalfOnly = kFALSE, Bool_t squareDP = kFALSE);
+
+ //! Destructor
+ virtual ~Lau2DSplineDPPdf();
+
+ //! Retrieve PDF normalisation
+ /*!
+ \return the normalisation factor
+ */
+ Double_t getHistNorm() const {return norm_;}
+
+ //! Perform the interpolation (unnormalised)
+ /*!
+ \param [in] x the x-axis value
+ \param [in] y the y-aixs value
+ \return the unnormalised PDF value
+ */
+ Double_t interpolateXY(Double_t x, Double_t y) const;
+
+ //! Perform the interpolation and divide by the normalisation
+ /*!
+ \param [in] x the x-axis abscissa value
+ \param [in] y the y-axis abscissa value
+ \return the normalised PDF value
+ */
+ Double_t interpolateXYNorm(Double_t x, Double_t y) const;
+
+ protected:
+ //! Calculate the PDF normalisation
+ void calcHistNorm();
+
+ private:
+ //! The maximum height of 2D histogram
+ Double_t maxHeight_;
+
+ //! The histogram normalisation
+ Double_t norm_;
+
+ //! A 2D cubic spline generated from the histogram
+ Lau2DCubicSpline* spline_;
+
+ ClassDef(Lau2DSplineDPPdf,0) // 2D Spline utility class for DP analyses
+};
+
+#endif
diff --git a/inc/LauBkgndDPModel.hh b/inc/LauBkgndDPModel.hh
index adde5a9..cd947b7 100644
--- a/inc/LauBkgndDPModel.hh
+++ b/inc/LauBkgndDPModel.hh
@@ -1,137 +1,147 @@
// Copyright University of Warwick 2004 - 2013.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
// Authors:
// Thomas Latham
// John Back
// Paul Harrison
/*! \file LauBkgndDPModel.hh
\brief File containing declaration of LauBkgndDPModel class.
*/
/*! \class LauBkgndDPModel
\brief Class for defining a histogram-based background Dalitz plot model
Class for defining a histogram-based background Dalitz plot model
*/
#ifndef LAU_BKGND_DP_MODEL
#define LAU_BKGND_DP_MODEL
#include <vector>
#include "LauAbsBkgndDPModel.hh"
class TH2;
-class Lau2DHistDPPdf;
+class Lau2DAbsDPPdf;
class LauDaughters;
class LauFitDataTree;
class LauVetoes;
class LauBkgndDPModel : public LauAbsBkgndDPModel {
public:
//! Constructor
/*!
\param [in] daughters the daughters in the decay
\param [in] vetoes the vetoes in the Datliz plot
*/
LauBkgndDPModel(LauDaughters* daughters, LauVetoes* vetoes);
//! Destructor
virtual ~LauBkgndDPModel();
//! Set background histogram
/*!
\param [in] histo the 2D histogram describing the DP distribution
\param [in] useInterpolation boolean flag to determine whether linear interpolation between bins should be used or simply the raw bin values
\param [in] fluctuateBins boolean flag to determine whether the bin contents should be fluctuated in accordance with their errors
\param [in] useUpperHalfOnly boolean flag to determine whether, in the case of a symmetric DP, the histogram supplied only includes the upper half of the DP
\param [in] squareDP boolean flag to determine whether the supplied histogram is given in square DP coordinates
*/
void setBkgndHisto(const TH2* histo, Bool_t useInterpolation,
Bool_t fluctuateBins, Bool_t useUpperHalfOnly,
Bool_t squareDP = kFALSE);
+ //! Set the background histogram and generate a spline
+ /*!
+ \param [in] histo the 2D histogram describing the DP distribution
+ \param [in] fluctuateBins boolean flag to determine whether the bin contents should be fluctuated in accordance with their errors
+ \param [in] useUpperHalfOnly boolean flag to determine whether, in the case of a symmetric DP, the histogram supplied only includes the upper half of the DP
+ \param [in] squareDP boolean flag to determine whether the supplied histogram is given in square DP coordinates
+ */
+ void setBkgndSpline(const TH2* histo, Bool_t fluctuateBins,
+ Bool_t useUpperHalfOnly, Bool_t squareDP);
+
//! Initialise the model
virtual void initialise();
//! Generate a toy MC event from the model
/*!
\return success/failure flag
*/
virtual Bool_t generate();
//! Get unnormalised likelihood for a given event
/*!
\param [in] iEvt the event number
\return the unnormalised likelihood value
*/
virtual Double_t getUnNormValue(UInt_t iEvt);
//! Get PDF normalisation constant
/*!
\return the PDF normalisation constant
*/
virtual Double_t getPdfNorm() const {return pdfNorm_;}
//! Get likelihood for a given event
/*!
\param [in] iEvt the event number
\return the likelihood value
*/
virtual Double_t getLikelihood(UInt_t iEvt);
//! Cache the input data and (if appropriate) the per-event likelihood values
/*!
\param [in] fitDataTree the input data
*/
virtual void fillDataTree(const LauFitDataTree& fitDataTree);
protected:
//! Calulate histogram value at a given point
/*!
\param [in] xVal the x-value
\param [in] yVal the y-value
\return the histogram value
*/
Double_t calcHistValue(Double_t xVal, Double_t yVal);
//! Set data event number
/*!
\param [in] iEvt the event number
*/
virtual void setDataEventNo(UInt_t iEvt);
private:
//! Flags whether the DP is symmetrical or not
Bool_t symmetricalDP_;
//! Flags whether or not to work in square DP coordinates
Bool_t squareDP_;
//! PDF of Dalitz plot background, from a 2D histogram
- Lau2DHistDPPdf* bgHistDPPdf_;
+ Lau2DAbsDPPdf* bgHistDPPdf_;
//! Cached histogram values for each event
std::vector<Double_t> bgData_;
//! Histogram value for the current event
Double_t curEvtHistVal_;
//! Maximum height of PDF
Double_t maxPdfHeight_;
//! Normalisation of PDF
Double_t pdfNorm_;
//! Boolean to indicate if the warning that there is no histogram has already been issued
Bool_t doneGenWarning_;
ClassDef(LauBkgndDPModel,0) // DP background model
};
#endif
diff --git a/inc/Laura++_LinkDef.h b/inc/Laura++_LinkDef.h
index fc23f72..fae7ccc 100644
--- a/inc/Laura++_LinkDef.h
+++ b/inc/Laura++_LinkDef.h
@@ -1,93 +1,96 @@
#ifdef __CINT__
#pragma link off all globals;
#pragma link off all classes;
#pragma link off all functions;
#pragma link C++ class Lau1DHistPdf;
#pragma link C++ class Lau2DAbsDP;
+#pragma link C++ class Lau2DAbsDPPdf;
#pragma link C++ class Lau2DAbsHistDP;
+#pragma link C++ class Lau2DAbsHistDPPdf;
#pragma link C++ class Lau2DCubicSpline;
#pragma link C++ class Lau2DHistDP;
#pragma link C++ class Lau2DHistDPPdf;
#pragma link C++ class Lau2DHistPdf;
#pragma link C++ class Lau2DSplineDP;
+#pragma link C++ class Lau2DSplineDPPdf;
#pragma link C++ class LauAbsBkgndDPModel;
#pragma link C++ class LauAbsCoeffSet;
#pragma link C++ class LauAbsDPDynamics;
#pragma link C++ class LauAbsFitModel;
#pragma link C++ class LauAbsPdf;
#pragma link C++ class LauAbsResonance;
#pragma link C++ class LauArgusPdf;
#pragma link C++ class LauAsymmCalc;
#pragma link C++ class LauBelleCPCoeffSet;
#pragma link C++ class LauBelleNR;
#pragma link C++ class LauBelleSymNR;
#pragma link C++ class LauBifurcatedGaussPdf;
#pragma link C++ class LauBkgndDPModel;
#pragma link C++ class LauBreitWignerRes;
#pragma link C++ class LauCacheData;
#pragma link C++ class LauCartesianCPCoeffSet;
#pragma link C++ class LauChebychevPdf;
#pragma link C++ class LauCleoCPCoeffSet;
#pragma link C++ class LauComplex;
#pragma link C++ class LauCPFitModel;
#pragma link C++ class LauCruijffPdf;
#pragma link C++ class LauCrystalBallPdf;
#pragma link C++ class LauDabbaRes;
#pragma link C++ class LauDatabasePDG;
#pragma link C++ class LauDaughters;
#pragma link C++ class LauDPDepBifurGaussPdf;
#pragma link C++ class LauDPDepCruijffPdf;
#pragma link C++ class LauDPDepGaussPdf;
#pragma link C++ class LauDPDepMapPdf;
#pragma link C++ class LauDPDepSumPdf;
#pragma link C++ class LauEffModel;
#pragma link C++ class LauEmbeddedData;
#pragma link C++ class LauExponentialPdf;
#pragma link C++ class LauFitDataTree;
#pragma link C++ class LauFitNtuple;
#pragma link C++ class LauFlatteRes;
#pragma link C++ class LauGaussPdf;
#pragma link C++ class LauGenNtuple;
#pragma link C++ class LauGounarisSakuraiRes;
#pragma link C++ class LauIntegrals;
#pragma link C++ class LauIsobarDynamics;
#pragma link C++ class LauKappaRes;
#pragma link C++ class LauKinematics;
#pragma link C++ class LauKMatrixProdPole;
#pragma link C++ class LauKMatrixProdSVP;
#pragma link C++ class LauKMatrixPropagator;
#pragma link C++ class LauKMatrixPropFactory;
#pragma link C++ class LauLASSBWRes;
#pragma link C++ class LauLASSNRRes;
#pragma link C++ class LauLASSRes;
#pragma link C++ class LauLinearPdf;
#pragma link C++ class LauMagPhaseCoeffSet;
#pragma link C++ class LauMagPhaseCPCoeffSet;
#pragma link C++ class LauNovosibirskPdf;
#pragma link C++ class LauNRAmplitude;
#pragma link C++ class LauParameter;
#pragma link C++ class LauParametricStepFuncPdf;
#pragma link C++ class LauParamFixed;
#pragma link C++ class LauParticlePDG;
#pragma link C++ class LauPrint;
#pragma link C++ class LauRealImagCoeffSet;
#pragma link C++ class LauRealImagCPCoeffSet;
#pragma link C++ class LauRelBreitWignerRes;
#pragma link C++ class LauResonanceInfo;
#pragma link C++ class LauResonanceMaker;
#pragma link C++ class LauScfMap;
#pragma link C++ class LauSigmaRes;
#pragma link C++ class LauSigmoidPdf;
#pragma link C++ class LauSimpleFitModel;
#pragma link C++ class LauSPlot;
#pragma link C++ class LauString;
#pragma link C++ class LauSumPdf;
#pragma link C++ class LauTextFileParser;
#pragma link C++ class LauVetoes;
#pragma link C++ namespace LauConstants;
#pragma link C++ namespace LauFitter;
#pragma link C++ namespace LauRandom;
#endif
diff --git a/src/Lau2DAbsHistDPPdf.cc b/src/Lau2DAbsHistDPPdf.cc
new file mode 100644
index 0000000..ac1f4e0
--- /dev/null
+++ b/src/Lau2DAbsHistDPPdf.cc
@@ -0,0 +1,103 @@
+
+// Copyright University of Warwick 2004 - 2013.
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Authors:
+// Thomas Latham
+// John Back
+// Paul Harrison
+
+/*! \file Lau2DAbsHistDPPdf.cc
+ \brief File containing implementation of Lau2DAbsHistDPPdf class.
+*/
+
+#include <iostream>
+
+#include "TAxis.h"
+#include "TH2.h"
+#include "TRandom.h"
+
+#include "Lau2DAbsHistDPPdf.hh"
+#include "LauKinematics.hh"
+#include "LauRandom.hh"
+#include "LauVetoes.hh"
+
+ClassImp(Lau2DAbsHistDPPdf)
+
+
+Lau2DAbsHistDPPdf::Lau2DAbsHistDPPdf(LauKinematics* kinematics, const LauVetoes* vetoes,
+ Bool_t useUpperHalfOnly, Bool_t squareDP) :
+ kinematics_(kinematics),
+ vetoes_(vetoes),
+ upperHalf_(useUpperHalfOnly),
+ squareDP_(squareDP)
+{
+}
+
+Lau2DAbsHistDPPdf::~Lau2DAbsHistDPPdf()
+{
+}
+
+void Lau2DAbsHistDPPdf::calcMaxHeight(TH2* hist)
+{
+ // Get the maximum height of the 2D histogram
+ maxHeight_ = 1.0;
+ if ( hist ) {
+ Int_t maxBin = hist->GetMaximumBin();
+ maxHeight_ = hist->GetBinContent(maxBin);
+ }
+
+ std::cout << "INFO in Lau2DAbsHistDPPdf::calcMaxHeight : Max height = " << maxHeight_ << std::endl;
+}
+
+void Lau2DAbsHistDPPdf::doBinFluctuation(TH2* hist)
+{
+ if ( !hist ) {
+ return;
+ }
+
+ TRandom* random = LauRandom::zeroSeedRandom();
+
+ Int_t nBinsX = static_cast<Int_t>(hist->GetNbinsX());
+ Int_t nBinsY = static_cast<Int_t>(hist->GetNbinsY());
+
+ for (Int_t i(0); i<nBinsX; i++) {
+ for (Int_t j(0); j<nBinsY; j++) {
+ Double_t currentContent = hist->GetBinContent(i+1,j+1);
+ Double_t currentError = hist->GetBinError(i+1,j+1);
+ Double_t newContent = random->Gaus(currentContent,currentError);
+ if (newContent<0.0) {
+ hist->SetBinContent(i+1,j+1,0.0);
+ } else {
+ hist->SetBinContent(i+1,j+1,newContent);
+ }
+ }
+ }
+}
+
+Bool_t Lau2DAbsHistDPPdf::withinDPBoundaries(Double_t x, Double_t y) const
+{
+ return squareDP_ ? kinematics_->withinSqDPLimits(x,y) : kinematics_->withinDPLimits(x,y);
+}
+
+void Lau2DAbsHistDPPdf::getUpperHalf(Double_t& x, Double_t& y) const
+{
+ if ( upperHalf_ == kTRUE ) {
+ if ( squareDP_ == kFALSE && x > y ) {
+ Double_t temp = y;
+ y = x;
+ x = temp;
+ } else if ( squareDP_ == kTRUE && y > 0.5 ) {
+ y = 1.0 - y;
+ }
+ }
+}
+
+void Lau2DAbsHistDPPdf::updateKinematics(Double_t x, Double_t y) const {
+ if (squareDP_ == kTRUE) {
+ kinematics_->updateSqDPKinematics(x,y);
+ } else {
+ kinematics_->updateKinematics(x,y);
+ }
+}
diff --git a/src/Lau2DCubicSpline.cc b/src/Lau2DCubicSpline.cc
index 6456740..1efb11e 100644
--- a/src/Lau2DCubicSpline.cc
+++ b/src/Lau2DCubicSpline.cc
@@ -1,349 +1,356 @@
// Copyright University of Warwick 2004 - 2013.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
// Authors:
// Thomas Latham
// John Back
// Paul Harrison
/*! \file Lau2DCubicSpline.cc
\brief File containing implementation of Lau2DCubicSpline class.
Class for defining a 2D cubic spline based on RooBinned2DBicubicBase by
Manuel Tobias Schiller <manuel.schiller@nikhef.nl> (2012-08-29).
*/
#include <cmath>
#include <cstdlib>
#include <iostream>
#include <TH2.h>
#include <TSystem.h>
#include "Lau2DCubicSpline.hh"
Lau2DCubicSpline::~Lau2DCubicSpline()
{ }
inline Double_t Lau2DCubicSpline::histcont(
const TH2& h, Int_t xbin, Int_t ybin) const
{
//reflect until we're in range
while(xbin < 0 || xbin >= nBinsX - 1) {
if(xbin < 0) xbin = -xbin-1;
if(xbin >= nBinsX -1) xbin = 2*(nBinsX-1) - xbin - 1;
}
while(ybin < 0 || ybin >= nBinsY - 1) {
if(ybin < 0) ybin = -ybin-1;
if(ybin >= nBinsY -1) ybin = 2*(nBinsY-1) - ybin - 1;
}
return h.GetBinContent(1 + xbin, 1 + ybin);
}
inline Double_t Lau2DCubicSpline::dhistdx(
const TH2& h, Int_t xbin, Int_t ybin) const
{
return 0.5 * (histcont(h, xbin + 1, ybin) -
histcont(h, xbin - 1, ybin));
}
inline Double_t Lau2DCubicSpline::dhistdy(
const TH2& h, Int_t xbin, Int_t ybin) const
{
return 0.5 * (histcont(h, xbin, ybin + 1) -
histcont(h, xbin, ybin - 1));
}
inline Double_t Lau2DCubicSpline::d2histdxdy(
const TH2& h, Int_t xbin, Int_t ybin) const
{
return 0.5 * (histcont(h, xbin - 1, ybin - 1) -
histcont(h, xbin + 1, ybin - 1) +
histcont(h, xbin + 1, ybin + 1) -
histcont(h, xbin - 1, ybin + 1));
}
Lau2DCubicSpline::Lau2DCubicSpline(const TH2& h) :
nBinsX(1 + h.GetNbinsX()), nBinsY(1 + h.GetNbinsY()),
binSizeX(h.GetXaxis()->GetBinWidth(1)),
binSizeY(h.GetYaxis()->GetBinWidth(1)),
xmin(h.GetXaxis()->GetBinCenter(1) - binSizeX),
xmax(h.GetXaxis()->GetBinCenter(nBinsX - 1) + binSizeX),
ymin(h.GetYaxis()->GetBinCenter(1) - binSizeY),
ymax(h.GetYaxis()->GetBinCenter(nBinsY - 1) + binSizeY),
coeffs(CoeffRecLen * nBinsX * nBinsY)
{
TAxis *xaxis = h.GetXaxis(), *yaxis = h.GetYaxis();
// verify that all bins have same size
for (Int_t i = 1; i < nBinsX; ++i) {
- if (std::abs(xaxis->GetBinWidth(i) / binSizeX - 1.) > 1e-9)
+ if (std::abs(xaxis->GetBinWidth(i) / binSizeX - 1.) > 1e-9) {
std::cerr << "ERROR in Lau2DCubicSpline constructor : the histogram has variable bin sizes." << std::endl;
gSystem->Exit(EXIT_FAILURE);
+ }
}
for (Int_t i = 1; i < nBinsY; ++i) {
- if (std::abs(yaxis->GetBinWidth(i) / binSizeY - 1.) > 1e-9)
+ if (std::abs(yaxis->GetBinWidth(i) / binSizeY - 1.) > 1e-9) {
std::cerr << "ERROR in Lau2DCubicSpline constructor : the histogram has variable bin sizes." << std::endl;
gSystem->Exit(EXIT_FAILURE);
+ }
}
// ok, go through histogram to precalculate the interpolation coefficients
// in rectangles between bin centres
//
// for that purpose, we map each of those rectangles to the unit square
for (Int_t j = -1; j < nBinsY - 1; ++j) {
for (Int_t i = -1; i < nBinsX - 1; ++i) {
const Double_t rhs[NCoeff] = {
// function values in bin centres
histcont(h, i, j),
histcont(h, i + 1, j),
histcont(h, i, j + 1),
histcont(h, i + 1, j + 1),
// df/dx in bin centres (finite difference approximation)
dhistdx(h, i, j),
dhistdx(h, i + 1, j),
dhistdx(h, i, j + 1),
dhistdx(h, i + 1, j + 1),
// df/dy in bin centres (finite difference approximation)
dhistdy(h, i, j),
dhistdy(h, i + 1, j),
dhistdy(h, i, j + 1),
dhistdy(h, i + 1, j + 1),
// d^2f/dxdy in bin centres (finite difference approximation)
d2histdxdy(h, i, j),
d2histdxdy(h, i + 1, j),
d2histdxdy(h, i, j + 1),
d2histdxdy(h, i + 1, j + 1)
};
// work out solution - strange array placement is due to the fact
// that terms with x/y to high powers can be small, so they should
// be added up first during evaluation to avoid cancellation
// issues; at the same time you want to access them in order to
// not confuse the CPU cache, so they're stored back to front
//
// a_00 ... a_30
coeff(1 + i, 1 + j, 15) = rhs[0];
coeff(1 + i, 1 + j, 14) = rhs[4];
coeff(1 + i, 1 + j, 13) =
3. * (-rhs[0] + rhs[1]) - 2. * rhs[4] - rhs[5];
coeff(1 + i, 1 + j, 12) =
2. * (rhs[0] - rhs[1]) + rhs[4] + rhs[5];
// a_31 ... a_31
coeff(1 + i, 1 + j, 11) = rhs[8];
coeff(1 + i, 1 + j, 10) = rhs[12];
coeff(1 + i, 1 + j, 9) =
3. * (-rhs[8] + rhs[9]) - 2. * rhs[12] - rhs[13];
coeff(1 + i, 1 + j, 8) =
2. * (rhs[8] - rhs[9]) + rhs[12] + rhs[13];
// a_02 ... a_32
coeff(1 + i, 1 + j, 7) =
3. * (-rhs[0] + rhs[2]) - 2. * rhs[8] - rhs[10];
coeff(1 + i, 1 + j, 6) =
3. * (-rhs[4] + rhs[6]) - 2. * rhs[12] - rhs[14];
coeff(1 + i, 1 + j, 5) =
9. * (rhs[0] - rhs[1] - rhs[2] + rhs[3]) +
6. * (rhs[4] - rhs[6] + rhs[8] - rhs[9]) + 4. * rhs[12] +
3. * (rhs[5] - rhs[7] + rhs[10] - rhs[11]) +
2. * (rhs[13] + rhs[14]) + rhs[15];
coeff(1 + i, 1 + j, 4) =
6. * (-rhs[0] + rhs[1] + rhs[2] - rhs[3]) +
4. * (-rhs[8] + rhs[9]) +
3. * (-rhs[4] - rhs[5] + rhs[6] + rhs[7]) +
2. * (-rhs[10] + rhs[11] - rhs[12] - rhs[13]) -
rhs[14] - rhs[15];
// a_03 ... a_33
coeff(1 + i, 1 + j, 3) =
2. * (rhs[0] - rhs[2]) + rhs[8] + rhs[10];
coeff(1 + i, 1 + j, 2) =
2. * (rhs[4] - rhs[6]) + rhs[12] + rhs[14];
coeff(1 + i, 1 + j, 1) =
6. * (-rhs[0] + rhs[1] + rhs[2] - rhs[3]) +
4. * (-rhs[4] + rhs[6]) +
3. * (-rhs[8] + rhs[9] - rhs[10] + rhs[11]) +
2. * (- rhs[5] + rhs[7] - rhs[12] - rhs[14]) -
rhs[13] - rhs[15];
coeff(1 + i, 1 + j, 0) =
4. * (rhs[0] - rhs[1] - rhs[2] + rhs[3]) +
2. * (rhs[4] + rhs[5] - rhs[6] - rhs[7] +
rhs[8] - rhs[9] + rhs[10] - rhs[11]) +
rhs[12] + rhs[13] + rhs[14] + rhs[15];
// coeff(1 + i, 1 + j, 17) contains integral of function over the
// square in "unit square coordinates", i.e. neglecting bin widths
// this is done to help speed up calculations of 2D integrals
Double_t sum = 0.;
for (Int_t k = 0; k < NCoeff; ++k)
sum += coeff(1 + i, 1 + j, k) /
Double_t((4 - (k % 4)) * (4 - (k / 4)));
coeff(1 + i, 1 + j, NCoeff) = sum;
}
}
}
Double_t Lau2DCubicSpline::evaluate(Double_t x, Double_t y) const
{
// protect against NaN and out of range
if (x <= xmin || x >= xmax || y <= ymin || y >= ymax || x != x || y != y)
return 0.;
// find the bin in question
const Int_t binx = Int_t(Double_t(nBinsX) * (x - xmin) / (xmax - xmin));
const Int_t biny = Int_t(Double_t(nBinsY) * (y - ymin) / (ymax - ymin));
// get low edge of bin
const Double_t xlo = Double_t(nBinsX - binx) / Double_t(nBinsX) * xmin +
Double_t(binx) / Double_t(nBinsX) * xmax;
const Double_t ylo = Double_t(nBinsY - biny) / Double_t(nBinsY) * ymin +
Double_t(biny) / Double_t(nBinsY) * ymax;
// normalise to coordinates in unit sqare
const Double_t hx = (x - xlo) / binSizeX;
const Double_t hy = (y - ylo) / binSizeY;
// monomials
const Double_t hxton[4] = { hx * hx * hx, hx * hx, hx, 1. };
const Double_t hyton[4] = { hy * hy * hy, hy * hy, hy, 1. };
// sum up
Double_t retVal = 0.;
for (Int_t k = 0; k < NCoeff; ++k)
retVal += coeff(binx, biny, k) * hxton[k % 4] * hyton[k / 4];
return retVal;
}
+Double_t Lau2DCubicSpline::analyticalIntegral() const
+{
+ return evalXY(xmin,xmax,ymin,ymax);
+}
+
Double_t Lau2DCubicSpline::analyticalIntegral(Double_t x1, Double_t x2, Double_t y1, Double_t y2) const
{
if(y1==y2) return evalX(x1, x2, y1);
if(x1==x2) return evalY(x1, y1, y2);
return evalXY(x1, x2, y1, y2);
}
Double_t Lau2DCubicSpline::evalX(Double_t x1, Double_t x2, Double_t y) const
{
// protect against NaN
if (x1 != x1 || x2 != x2 || y != y) return 0.;
// find the bin in question
const Int_t biny = Int_t(Double_t(nBinsY) * (y - ymin) / (ymax - ymin));
// get low edge of bin
const Double_t ylo = Double_t(nBinsY - biny) / Double_t(nBinsY) * ymin +
Double_t(biny) / Double_t(nBinsY) * ymax;
// normalise to coordinates in unit sqare
const Double_t hy = (y - ylo) / binSizeY;
// monomials
const Double_t hyton[4] = { hy * hy * hy, hy * hy, hy, 1. };
// integral
Double_t sum = 0.;
for (Int_t binx = 0; binx < nBinsX; ++binx) {
// get low/high edge of bin
const Double_t xhi = Double_t(nBinsX - binx - 1) / Double_t(nBinsX) * xmin +
Double_t(binx + 1) / Double_t(nBinsX) * xmax;
if (xhi < x1) continue;
const Double_t xlo = Double_t(nBinsX - binx) / Double_t(nBinsX) * xmin +
Double_t(binx) / Double_t(nBinsX) * xmax;
if (xlo > x2) break;
// work out integration range
const Double_t a = ((xlo > x1) ? 0. : (x1 - xlo)) / binSizeX;
const Double_t b = ((xhi < x2) ? binSizeX : (x2 - xlo)) / binSizeX;
// integrated monomials
const Double_t hxton[4] = { 0.25 * (b * b * b * b - a * a * a * a),
(b * b * b - a * a * a) / 3., 0.5 * (b * b - a * a), b - a };
Double_t lsum = 0.;
for (Int_t k = 0; k < NCoeff; ++k)
lsum += coeff(binx, biny, k) * hxton[k % 4] * hyton[k / 4];
sum += lsum;
}
// move from unit square coordinates to user coordinates
return sum * binSizeX;
}
Double_t Lau2DCubicSpline::evalY(Double_t x, Double_t y1, Double_t y2) const
{
// protect against NaN
if (x != x || y1 != y1 || y2 != y2) return 0.;
// find the bin in question
const Int_t binx = Int_t(Double_t(nBinsX) * (x - xmin) / (xmax - xmin));
// get low edge of bin
const Double_t xlo = Double_t(nBinsX - binx) / Double_t(nBinsX) * xmin +
Double_t(binx) / Double_t(nBinsX) * xmax;
// normalise to coordinates in unit sqare
const Double_t hx = (x - xlo) / binSizeX;
// monomials
const Double_t hxton[4] = { hx * hx * hx, hx * hx, hx, 1. };
// integral
Double_t sum = 0.;
for (Int_t biny = 0; biny < nBinsY; ++biny) {
// get low/high edge of bin
const Double_t yhi = Double_t(nBinsY - biny - 1) / Double_t(nBinsY) * ymin +
Double_t(biny + 1) / Double_t(nBinsY) * ymax;
if (yhi < y1) continue;
const Double_t ylo = Double_t(nBinsY - biny) / Double_t(nBinsY) * ymin +
Double_t(biny) / Double_t(nBinsY) * ymax;
if (ylo > y2) break;
// work out integration range
const Double_t a = ((ylo > y1) ? 0. : (y1 - ylo)) / binSizeY;
const Double_t b = ((yhi < y2) ? binSizeY : (y2 - ylo)) / binSizeY;
// integrated monomials
const Double_t hyton[4] = { 0.25 * (b * b * b * b - a * a * a * a),
(b * b * b - a * a * a) / 3., 0.5 * (b * b - a * a), b - a };
Double_t lsum = 0.;
for (Int_t k = 0; k < NCoeff; ++k)
lsum += coeff(binx, biny, k) * hxton[k % 4] * hyton[k / 4];
sum += lsum;
}
// move from unit square coordinates to user coordinates
return sum * binSizeY;
}
Double_t Lau2DCubicSpline::evalXY(
Double_t x1, Double_t x2, Double_t y1, Double_t y2) const
{
// protect against NaN
if (x1 != x1 || x2 != x2 || y1 != y1 || y2 != y2) return 0.;
// integral
Double_t sum = 0.;
for (Int_t biny = 0; biny < nBinsY; ++biny) {
// get low/high edge of bin
const Double_t yhi = Double_t(nBinsY - biny - 1) / Double_t(nBinsY) * ymin +
Double_t(biny + 1) / Double_t(nBinsY) * ymax;
if (yhi < y1) continue;
const Double_t ylo = Double_t(nBinsY - biny) / Double_t(nBinsY) * ymin +
Double_t(biny) / Double_t(nBinsY) * ymax;
if (ylo > y2) break;
// work out integration range
const Double_t ay = ((ylo > y1) ? 0. : (y1 - ylo)) / binSizeY;
const Double_t by = ((yhi < y2) ? binSizeY : (y2 - ylo)) / binSizeY;
const Bool_t yFullyContained = std::abs(by - ay - 1.0) < 1e-15;
// integrated monomials
const Double_t hyton[4] = {
0.25 * (by * by * by * by - ay * ay * ay * ay),
(by * by * by - ay * ay * ay) / 3., 0.5 * (by * by - ay * ay),
by - ay };
for (Int_t binx = 0; binx < nBinsX; ++binx) {
// get low/high edge of bin
const Double_t xhi = Double_t(nBinsX - binx - 1) / Double_t(nBinsX) * xmin +
Double_t(binx + 1) / Double_t(nBinsX) * xmax;
if (xhi < x1) continue;
const Double_t xlo = Double_t(nBinsX - binx) / Double_t(nBinsX) * xmin +
Double_t(binx) / Double_t(nBinsX) * xmax;
if (xlo > x2) break;
// work out integration range
const Double_t ax = ((xlo > x1) ? 0. : (x1 - xlo)) / binSizeX;
const Double_t bx = ((xhi < x2) ? binSizeX : (x2 - xlo)) / binSizeX;
const Bool_t xFullyContained = std::abs(bx - ax - 1.0) < 1e-15;
if (xFullyContained && yFullyContained) {
// for fully contained bins, we have cached the integral
sum += coeff(binx, biny, NCoeff);
continue;
}
// integrated monomials
const Double_t hxton[4] = {
0.25 * (bx * bx * bx * bx - ax * ax * ax * ax),
(bx * bx * bx - ax * ax * ax) / 3., 0.5 * (bx * bx - ax * ax),
bx - ax };
// integrate over bin
Double_t lsum = 0.;
for (Int_t k = 0; k < NCoeff; ++k)
lsum += coeff(binx, biny, k) * hxton[k % 4] * hyton[k / 4];
sum += lsum;
}
}
// move from unit square coordinates to user coordinates
return sum * binSizeX * binSizeY;
}
diff --git a/src/Lau2DHistDPPdf.cc b/src/Lau2DHistDPPdf.cc
index a7ddeb3..b546c4f 100644
--- a/src/Lau2DHistDPPdf.cc
+++ b/src/Lau2DHistDPPdf.cc
@@ -1,435 +1,357 @@
// Copyright University of Warwick 2004 - 2013.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
// Authors:
// Thomas Latham
// John Back
// Paul Harrison
/*! \file Lau2DHistDPPdf.cc
\brief File containing implementation of Lau2DHistDPPdf class.
*/
#include <iostream>
#include "TAxis.h"
#include "TH2.h"
#include "TRandom.h"
#include "Lau2DHistDPPdf.hh"
#include "LauKinematics.hh"
#include "LauRandom.hh"
#include "LauVetoes.hh"
ClassImp(Lau2DHistDPPdf)
Lau2DHistDPPdf::Lau2DHistDPPdf(const TH2* hist, LauKinematics* kinematics, const LauVetoes* vetoes,
Bool_t useInterpolation, Bool_t fluctuateBins, Bool_t useUpperHalfOnly, Bool_t squareDP) :
+ Lau2DAbsHistDPPdf(kinematics,vetoes,useUpperHalfOnly,squareDP),
hist_(hist ? dynamic_cast<TH2*>(hist->Clone()) : 0),
- kinematics_(kinematics),
- vetoes_(vetoes),
minX_(0.0), maxX_(0.0), minY_(0.0), maxY_(0.0), rangeX_(0.0), rangeY_(0.0),
binXWidth_(0.0), binYWidth_(0.0),
invBinXWidth_(0.0), invBinYWidth_(0.0),
- maxHeight_(0.0),
nBinsX_(0), nBinsY_(0),
norm_(0.0),
- useInterpolation_(useInterpolation),
- upperHalf_(useUpperHalfOnly),
- squareDP_(squareDP)
+ useInterpolation_(useInterpolation)
{
// For square Dalitz plots, the co-ordinates must be m', theta'
// The input data to the fit is still in m13^2, m23^2 co-ords, and a
// transformation is applied to go from one co-ordinate system to the
// other.
// Save various attributes of the histogram
// (axis ranges, number of bins, areas)
if ( hist_ ) {
TAxis* xAxis = hist_->GetXaxis();
minX_ = xAxis->GetXmin();
maxX_ = xAxis->GetXmax();
TAxis* yAxis = hist_->GetYaxis();
minY_ = yAxis->GetXmin();
maxY_ = yAxis->GetXmax();
nBinsX_ = hist_->GetNbinsX();
nBinsY_ = hist_->GetNbinsY();
} else {
- minX_ = kinematics_->getm13SqMin();
- maxX_ = kinematics_->getm13SqMax();
+ minX_ = getKinematics()->getm13SqMin();
+ maxX_ = getKinematics()->getm13SqMax();
- minY_ = kinematics_->getm23SqMin();
- maxY_ = kinematics_->getm23SqMax();
+ minY_ = getKinematics()->getm23SqMin();
+ maxY_ = getKinematics()->getm23SqMax();
nBinsX_ = 1;
nBinsY_ = 1;
}
rangeX_ = maxX_ - minX_;
rangeY_ = maxY_ - minY_;
if (nBinsX_ > 0) { binXWidth_ = TMath::Abs(rangeX_)/(nBinsX_*1.0); }
if (nBinsY_ > 0) { binYWidth_ = TMath::Abs(rangeY_)/(nBinsY_*1.0); }
if (binXWidth_ > 1e-10) { invBinXWidth_ = 1.0/binXWidth_; }
if (binYWidth_ > 1e-10) { invBinYWidth_ = 1.0/binYWidth_; }
// If the bins are to be fluctuated then do so now before
// calculating anything that depends on the bin content.
if (fluctuateBins) {
- this->doBinFluctuation();
+ this->doBinFluctuation(hist_);
}
// Calculate the PDF normalisation.
this->calcHistNorm();
// And check it is OK.
this->checkNormalisation();
// Also obtain the maximum height
- this->calcMaxHeight();
+ this->calcMaxHeight(hist_);
}
Lau2DHistDPPdf::~Lau2DHistDPPdf()
{
// Destructor
delete hist_; hist_ = 0;
}
-void Lau2DHistDPPdf::calcMaxHeight()
-{
- // Get the maximum height of the 2D histogram
- maxHeight_ = 1.0;
- if ( hist_ ) {
- Int_t maxBin = hist_->GetMaximumBin();
- maxHeight_ = hist_->GetBinContent(maxBin);
- }
-
- std::cout << "INFO in Lau2DHistDPPdf::calcMaxHeight : Max height = " << maxHeight_ << std::endl;
-}
-
void Lau2DHistDPPdf::calcHistNorm()
{
// Calculate the histogram normalisation. We must integrate
// over the allowed kinematic DP region.
// Loop over the total x and y range to get the total area
// under x and y. Just sum the contributions up using 1e-3 increments
// of the range in x and y. Multiply the end result by dx and dy.
const Double_t axisMin = TMath::Min( minX_, minY_ );
const Double_t axisMax = TMath::Max( maxX_, maxY_ );
const Double_t axisRange = axisMax - axisMin;
const Double_t dx(1e-3 * axisRange);
const Double_t dy(dx);
Double_t area(0.0);
Double_t x(axisMin + dx/2.0);
while (x < axisMax) {
Double_t y(axisMin + dy/2.0);
while (y < axisMax) {
area += this->interpolateXY(x,y);
y += dy;
} // y while loop
x += dx;
} // x while loop
norm_ = area*dx*dy;
std::cout << "INFO in Lau2DHistDPPdf::calcHistNorm : Norm = " << norm_ << ", bX*bY = " << binXWidth_ << "*" << binYWidth_ << " = " << binXWidth_*binYWidth_ << std::endl;
}
Double_t Lau2DHistDPPdf::getBinHistValue(Int_t xBinNo, Int_t yBinNo) const
{
if (xBinNo < 0) {
xBinNo = 0;
} else if (xBinNo >= nBinsX_) {
return 0.0;
}
if (yBinNo < 0) {
yBinNo = 0;
} else if (yBinNo >= nBinsY_) {
return 0.0;
}
if ( hist_ == 0 ) {
return 1.0;
}
Double_t value = hist_->GetBinContent(xBinNo+1, yBinNo+1);
return value;
}
Double_t Lau2DHistDPPdf::interpolateXYNorm(Double_t x, Double_t y) const
{
// Get the normalised interpolated value.
Double_t value = this->interpolateXY(x,y);
return value/norm_;
}
Double_t Lau2DHistDPPdf::interpolateXY(Double_t x, Double_t y) const
{
// This function returns the interpolated value of the histogram function
// for the given values of x and y by finding the adjacent bins and extrapolating
// using weights based on the inverse distance of the point from the adajcent
// bin centres.
// Here, x = m13^2, y = m23^2, or m', theta' for square DP co-ordinates
// If we're only using one half then flip co-ordinates
// appropriately for conventional or square DP
- if ( upperHalf_ == kTRUE ) {
- if ( squareDP_ == kFALSE && x > y ) {
- Double_t temp = y;
- y = x;
- x = temp;
- } else if ( squareDP_ == kTRUE && y > 0.5 ) {
- y = 1.0 - y;
- }
- }
+ getUpperHalf(x,y);
// First ask whether the point is inside the kinematic region.
- Bool_t withinDP(kFALSE);
- if (squareDP_ == kTRUE) {
- withinDP = kinematics_->withinSqDPLimits(x,y);
- } else {
- withinDP = kinematics_->withinDPLimits(x,y);
- }
- if (withinDP == kFALSE) {return 0.0;}
+ if (withinDPBoundaries(x,y) == kFALSE) {return 0.0;}
// Update the kinematics to the position of interest.
- if (squareDP_ == kTRUE) {
- kinematics_->updateSqDPKinematics(x,y);
- } else {
- kinematics_->updateKinematics(x,y);
- }
+ updateKinematics(x,y);
// Check that we're not inside a veto
Bool_t vetoOK(kTRUE);
- if (vetoes_) {
- vetoOK = vetoes_->passVeto(kinematics_);
+ if (getVetoes()) {
+ vetoOK = getVetoes()->passVeto(getKinematics());
}
if (vetoOK == kFALSE) {return 0.0;}
// Find the 2D histogram bin for x and y
Int_t i = TMath::FloorNint((x - minX_)*invBinXWidth_);
Int_t j = TMath::FloorNint((y - minY_)*invBinYWidth_);
if (i >= nBinsX_) {i = nBinsX_ - 1;}
if (j >= nBinsY_) {j = nBinsY_ - 1;}
// Ask whether we want to do the interpolation, since this method is
// not reliable for low statistics histograms.
if (useInterpolation_ == kFALSE) {
return this->getBinHistValue(i,j);
}
// Find the bin centres (actual co-ordinate positions, not histogram indices)
Double_t cbinx = (i+0.5)*binXWidth_ + minX_;
Double_t cbiny = (j+0.5)*binYWidth_ + minY_;
// If bin centres are outside kinematic region, do not extrapolate
- if (squareDP_ == kTRUE) {
- withinDP = kinematics_->withinSqDPLimits(cbinx, cbiny);
- } else {
- withinDP = kinematics_->withinDPLimits(cbinx, cbiny);
- }
- if (withinDP == kFALSE) {
+ if (withinDPBoundaries(cbinx, cbiny) == kFALSE) {
return this->getBinHistValue(i,j);
}
// Find the adjacent bins
Double_t deltax = x - cbinx;
Double_t deltay = y - cbiny;
Int_t i_adj(0), j_adj(0);
if (deltax > 0.0) {
i_adj = i + 1;
} else {
i_adj = i - 1;
}
if (deltay > 0.0) {
j_adj = j + 1;
} else {
j_adj = j - 1;
}
Bool_t isXBoundary(kFALSE), isYBoundary(kFALSE);
Double_t value(0.0);
if (i_adj >= nBinsX_ || i_adj < 0) {isYBoundary = kTRUE;}
if (j_adj >= nBinsY_ || j_adj < 0) {isXBoundary = kTRUE;}
// In the corners, do no interpolation. Use entry in bin (i,j)
if (isXBoundary == kTRUE && isYBoundary == kTRUE) {
value = this->getBinHistValue(i,j);
} else if (isXBoundary == kTRUE && isYBoundary == kFALSE) {
// Find the adjacent x bin centre
Double_t cbinx_adj = Double_t(i_adj+0.5)*binXWidth_ + minX_;
- if (squareDP_ == kTRUE) {
- withinDP = kinematics_->withinSqDPLimits(cbinx_adj, y);
- } else {
- withinDP = kinematics_->withinDPLimits(cbinx_adj, y);
- }
-
- if (withinDP == kFALSE) {
+ if (withinDPBoundaries(cbinx_adj, y) == kFALSE) {
// The adjacent bin is outside the DP range. Don't extrapolate.
value = this->getBinHistValue(i,j);
} else {
Double_t dx0 = TMath::Abs(x - cbinx);
Double_t dx1 = TMath::Abs(cbinx_adj - x);
Double_t inter_denom = dx0 + dx1;
Double_t value1 = this->getBinHistValue(i,j);
Double_t value2 = this->getBinHistValue(i_adj,j);
value = (value1*dx1 + value2*dx0)/inter_denom;
}
} else if (isYBoundary == kTRUE && isXBoundary == kFALSE) {
// Find the adjacent y bin centre
Double_t cbiny_adj = Double_t(j_adj+0.5)*binYWidth_ + minY_;
- if (squareDP_ == kTRUE) {
- withinDP = kinematics_->withinSqDPLimits(x, cbiny_adj);
- } else {
- withinDP = kinematics_->withinDPLimits(x, cbiny_adj);
- }
-
- if (withinDP == kFALSE) {
+ if (withinDPBoundaries(x, cbiny_adj) == kFALSE) {
// The adjacent bin is outside the DP range. Don't extrapolate.
value = this->getBinHistValue(i,j);
} else {
Double_t dy0 = TMath::Abs(y - cbiny);
Double_t dy1 = TMath::Abs(cbiny_adj - y);
Double_t inter_denom = dy0 + dy1;
Double_t value1 = this->getBinHistValue(i,j);
Double_t value2 = this->getBinHistValue(i,j_adj);
value = (value1*dy1 + value2*dy0)/inter_denom;
}
} else {
// 2D linear interpolation using inverse distance as weights.
// Find the adjacent x and y bin centres
Double_t cbinx_adj = Double_t(i_adj+0.5)*binXWidth_ + minX_;
Double_t cbiny_adj = Double_t(j_adj+0.5)*binYWidth_ + minY_;
- if (squareDP_ == kTRUE) {
- withinDP = kinematics_->withinSqDPLimits(cbinx_adj, cbiny_adj);
- } else {
- withinDP = kinematics_->withinDPLimits(cbinx_adj, cbiny_adj);
- }
-
- if (withinDP == kFALSE) {
+ if (withinDPBoundaries(cbinx_adj, cbiny_adj) == kFALSE) {
// The adjacent bin is outside the DP range. Don't extrapolate.
value = this->getBinHistValue(i,j);
} else {
Double_t dx0 = TMath::Abs(x - cbinx);
Double_t dx1 = TMath::Abs(cbinx_adj - x);
Double_t dy0 = TMath::Abs(y - cbiny);
Double_t dy1 = TMath::Abs(cbiny_adj - y);
Double_t inter_denom = (dx0 + dx1)*(dy0 + dy1);
Double_t value1 = this->getBinHistValue(i,j);
Double_t value2 = this->getBinHistValue(i_adj,j);
Double_t value3 = this->getBinHistValue(i,j_adj);
Double_t value4 = this->getBinHistValue(i_adj,j_adj);
value = value1*dx1*dy1 + value2*dx0*dy1 + value3*dx1*dy0 + value4*dx0*dy0;
value /= inter_denom;
}
}
return value;
}
void Lau2DHistDPPdf::checkNormalisation()
{
// Loop over the total x and y range to get the total area
// under x and y. Just sum the contributions up using 1e-3 increments
// of the range in x and y. Multiply the end result by dx and dy.
const Double_t axisMin = TMath::Min( minX_, minY_ );
const Double_t axisMax = TMath::Max( maxX_, maxY_ );
const Double_t axisRange = axisMax - axisMin;
const Double_t dx(1e-3 * axisRange);
const Double_t dy(dx);
Double_t area(0.0);
Double_t areaNoNorm(0.0);
// Preserve the value of a variable we change temporarily
Bool_t interpolate = useInterpolation_;
Double_t x(axisMin + dx/2.0);
while (x < axisMax) {
Double_t y(axisMin + dy/2.0);
while (y < axisMax) {
area += this->interpolateXYNorm(x,y);
useInterpolation_ = kFALSE;
areaNoNorm += this->interpolateXY(x,y);
useInterpolation_ = kTRUE;
y += dy;
} // y while loop
x += dx;
} // x while loop
Double_t norm = area*dx*dy;
useInterpolation_ = interpolate; //Restore old value
std::cout << "INFO in Lau2DHistDPPdf::checkNormalisation : Area = " << area << ", dx = " << dx << ", dy = " << dy << ", dx*dy = " << dx*dy << std::endl;
std::cout << " : Area with no norm = " << areaNoNorm << "*dx*dy = " << areaNoNorm*dx*dy << std::endl;
std::cout << " : The total area of the normalised histogram PDF is " << norm << std::endl;
}
-void Lau2DHistDPPdf::doBinFluctuation()
-{
- if ( !hist_ ) {
- return;
- }
-
- TRandom* random = LauRandom::zeroSeedRandom();
- for (Int_t i(0); i<nBinsX_; i++) {
- for (Int_t j(0); j<nBinsY_; j++) {
- Double_t currentContent = hist_->GetBinContent(i+1,j+1);
- Double_t currentError = hist_->GetBinError(i+1,j+1);
- Double_t newContent = random->Gaus(currentContent,currentError);
- if (newContent<0.0) {
- hist_->SetBinContent(i+1,j+1,0.0);
- } else {
- hist_->SetBinContent(i+1,j+1,newContent);
- }
- }
- }
-}
-
diff --git a/src/Lau2DSplineDPPdf.cc b/src/Lau2DSplineDPPdf.cc
new file mode 100644
index 0000000..524db7c
--- /dev/null
+++ b/src/Lau2DSplineDPPdf.cc
@@ -0,0 +1,101 @@
+
+// Copyright University of Warwick 2004 - 2013.
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Authors:
+// Thomas Latham
+// John Back
+// Paul Harrison
+
+/*! \file Lau2DSplineDPPdf.cc
+ \brief File containing implementation of Lau2DSplineDPPdf class.
+*/
+
+#include <iostream>
+
+#include "TAxis.h"
+#include "TH2.h"
+#include "TRandom.h"
+#include "TSystem.h"
+
+#include "Lau2DSplineDPPdf.hh"
+#include "Lau2DCubicSpline.hh"
+#include "LauDaughters.hh"
+#include "LauKinematics.hh"
+#include "LauRandom.hh"
+
+ClassImp(Lau2DSplineDPPdf)
+
+
+Lau2DSplineDPPdf::Lau2DSplineDPPdf(const TH2* hist, LauKinematics* kinematics, const LauVetoes* vetoes,
+ Bool_t fluctuateBins, Bool_t useUpperHalfOnly, Bool_t squareDP) :
+ Lau2DAbsHistDPPdf(kinematics,vetoes,useUpperHalfOnly,squareDP),
+ spline_(0)
+{
+ //We may need to modify the histogram so clone it
+ TH2* tempHist(hist ? dynamic_cast<TH2*>(hist->Clone()) : 0);
+
+ if ( ! tempHist ) {
+ std::cerr << "ERROR in Lau2DSplineDPPdf constructor : the histogram pointer is null." << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+
+ if (fluctuateBins) {
+ this->doBinFluctuation(tempHist);
+ }
+
+ spline_ = new Lau2DCubicSpline(*tempHist);
+
+ // Calculate the PDF normalisation.
+ //this->calcHistNorm();
+
+ // And check it is OK.
+ //this->checkNormalisation();
+
+ // Also obtain the maximum height
+ this->calcMaxHeight(tempHist);
+
+ delete tempHist;
+}
+
+Lau2DSplineDPPdf::~Lau2DSplineDPPdf()
+{
+ delete spline_;
+ spline_ = 0;
+}
+
+Double_t Lau2DSplineDPPdf::interpolateXYNorm(Double_t x, Double_t y) const
+{
+ // Get the normalised interpolated value.
+ Double_t value = this->interpolateXY(x,y);
+ return value/norm_;
+
+}
+
+Double_t Lau2DSplineDPPdf::interpolateXY(Double_t x, Double_t y) const
+{
+ // This function returns the interpolated value of the histogram function
+ // for the given values of x and y by finding the adjacent bins and extrapolating
+ // using weights based on the inverse distance of the point from the adajcent
+ // bin centres.
+ // Here, x = m13^2, y = m23^2, or m', theta' for square DP co-ordinates
+
+ // If we're only using one half then flip co-ordinates
+ // appropriately for conventional or square DP
+ getUpperHalf(x,y);
+
+ // First ask whether the point is inside the kinematic region.
+ if (withinDPBoundaries(x,y) == kFALSE) {
+ std::cerr << "WARNING in Lau2DSplineDPPdf::interpolateXY : Given position is outside the DP boundary, returning 0.0." << std::endl;
+ return 0.0;
+ }
+
+ return spline_->evaluate(x,y);
+
+}
+
+void Lau2DSplineDPPdf::calcHistNorm()
+{
+ norm_ = spline_->analyticalIntegral();
+}
diff --git a/src/LauBkgndDPModel.cc b/src/LauBkgndDPModel.cc
index b537dbb..633eeef 100644
--- a/src/LauBkgndDPModel.cc
+++ b/src/LauBkgndDPModel.cc
@@ -1,255 +1,272 @@
// Copyright University of Warwick 2004 - 2013.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
// Authors:
// Thomas Latham
// John Back
// Paul Harrison
/*! \file LauBkgndDPModel.cc
\brief File containing implementation of LauBkgndDPModel class.
*/
#include <cstdlib>
#include <iostream>
using std::cout;
using std::cerr;
using std::endl;
#include "TRandom.h"
#include "TSystem.h"
#include "Lau2DHistDPPdf.hh"
+#include "Lau2DSplineDPPdf.hh"
#include "LauBkgndDPModel.hh"
#include "LauDaughters.hh"
#include "LauFitDataTree.hh"
#include "LauKinematics.hh"
#include "LauRandom.hh"
#include "LauVetoes.hh"
ClassImp(LauBkgndDPModel)
LauBkgndDPModel::LauBkgndDPModel(LauDaughters* daughters, LauVetoes* vetoes) :
LauAbsBkgndDPModel(daughters, vetoes),
symmetricalDP_(kFALSE),
squareDP_(kFALSE),
bgHistDPPdf_(0),
curEvtHistVal_(0.0),
maxPdfHeight_(1.0),
pdfNorm_(1.0),
doneGenWarning_(kFALSE)
{
if (daughters != 0) {
symmetricalDP_ = daughters->gotSymmetricalDP();
}
}
LauBkgndDPModel::~LauBkgndDPModel()
{
if (bgHistDPPdf_ != 0) {
delete bgHistDPPdf_; bgHistDPPdf_ = 0;
}
}
void LauBkgndDPModel::setBkgndHisto(const TH2* histo, Bool_t useInterpolation, Bool_t fluctuateBins, Bool_t useUpperHalfOnly, Bool_t squareDP)
{
Bool_t upperHalf = kFALSE;
if (symmetricalDP_ == kTRUE && useUpperHalfOnly == kTRUE) {upperHalf = kTRUE;}
cout<<"Bg histogram has upperHalf = "<<static_cast<Int_t>(upperHalf)<<endl;
squareDP_ = squareDP;
LauKinematics* kinematics = this->getKinematics();
const LauVetoes* vetoes = this->getVetoes();
bgHistDPPdf_ = new Lau2DHistDPPdf(histo, kinematics, vetoes, useInterpolation, fluctuateBins, upperHalf, squareDP);
maxPdfHeight_ = bgHistDPPdf_->getMaxHeight();
pdfNorm_ = bgHistDPPdf_->getHistNorm();
}
+void LauBkgndDPModel::setBkgndSpline(const TH2* histo, Bool_t fluctuateBins, Bool_t useUpperHalfOnly, Bool_t squareDP)
+{
+ Bool_t upperHalf = kFALSE;
+ if (symmetricalDP_ == kTRUE && useUpperHalfOnly == kTRUE) {upperHalf = kTRUE;}
+ cout<<"Bg histogram has upperHalf = "<<static_cast<Int_t>(upperHalf)<<endl;
+
+ squareDP_ = squareDP;
+
+ LauKinematics* kinematics = this->getKinematics();
+ const LauVetoes* vetoes = this->getVetoes();
+ bgHistDPPdf_ = new Lau2DSplineDPPdf(histo, kinematics, vetoes, fluctuateBins, upperHalf, squareDP);
+
+ maxPdfHeight_ = bgHistDPPdf_->getMaxHeight();
+ pdfNorm_ = bgHistDPPdf_->getHistNorm();
+}
+
Double_t LauBkgndDPModel::calcHistValue(Double_t xVal, Double_t yVal)
{
// Get the likelihood value of the background in the Dalitz plot.
// Check that we have a valid histogram PDF
if (bgHistDPPdf_ == 0) {
cerr << "WARNING in LauBkgndDPModel::calcHistValue : We don't have a histogram so assuming the likelihood is flat in the Dalitz plot." << endl;
this->setBkgndHisto( 0, kFALSE, kFALSE, kFALSE, kFALSE );
}
// Find out the un-normalised PDF value
Double_t value = bgHistDPPdf_->interpolateXY(xVal, yVal);
LauKinematics* kinematics = this->getKinematics();
// For square DP co-ordinates, need to divide by Jacobian
if (squareDP_ == kTRUE) {
// Make sure that square DP kinematics are correct, then get Jacobian
kinematics->updateSqDPKinematics(xVal, yVal);
Double_t jacobian = kinematics->calcSqDPJacobian();
value /= jacobian;
}
return value;
}
void LauBkgndDPModel::initialise()
{
}
Bool_t LauBkgndDPModel::generate()
{
// Routine to generate the background, using data provided by an
// already defined histogram.
LauKinematics* kinematics = this->getKinematics();
Bool_t gotBG(kFALSE);
while (gotBG == kFALSE) {
if (squareDP_ == kTRUE) {
// Generate a point in m', theta' space. By construction, this point
// is already within the DP region.
Double_t mPrime(0.0), thetaPrime(0.0);
kinematics->genFlatSqDP(mPrime, thetaPrime);
// If we're in a symmetrical DP then we should only generate events in one half
if ( symmetricalDP_ && thetaPrime > 0.5 ) {
thetaPrime = 1.0 - thetaPrime;
}
// Calculate histogram height for DP point and
// compare with the maximum height
if ( bgHistDPPdf_ != 0 ) {
Double_t bgContDP = bgHistDPPdf_->interpolateXY(mPrime, thetaPrime)/maxPdfHeight_;
if (LauRandom::randomFun()->Rndm() < bgContDP) {
kinematics->updateSqDPKinematics(mPrime, thetaPrime);
gotBG = kTRUE;
}
} else {
if ( !doneGenWarning_ ) {
cerr << "WARNING in LauBkgndDPModel::generate : We don't have a histogram so generating flat in the square DP, which won't be flat in the conventional DP!" << endl;
cerr << "WARNING in LauBkgndDPModel::generate : This should never happen!! What have you done?!" << endl;
doneGenWarning_ = kTRUE;
}
kinematics->updateSqDPKinematics(mPrime, thetaPrime);
gotBG = kTRUE;
}
} else {
// Generate a point in the Dalitz plot (phase-space).
Double_t m13Sq(0.0), m23Sq(0.0);
kinematics->genFlatPhaseSpace(m13Sq, m23Sq);
// If we're in a symmetrical DP then we should only generate events in one half
if ( symmetricalDP_ && m13Sq > m23Sq ) {
Double_t tmpSq = m13Sq;
m13Sq = m23Sq;
m23Sq = tmpSq;
}
// Calculate histogram height for DP point and
// compare with the maximum height
if ( bgHistDPPdf_ != 0 ) {
Double_t bgContDP = bgHistDPPdf_->interpolateXY(m13Sq, m23Sq)/maxPdfHeight_;
if (LauRandom::randomFun()->Rndm() < bgContDP) {
kinematics->updateKinematics(m13Sq, m23Sq);
gotBG = kTRUE;
}
} else {
if ( !doneGenWarning_ ) {
cerr << "WARNING in LauBkgndDPModel::generate : We don't have a histogram so generating flat in the DP." << endl;
doneGenWarning_ = kTRUE;
}
kinematics->updateKinematics(m13Sq, m23Sq);
gotBG = kTRUE;
}
}
}
// Implement veto
Bool_t vetoOK(kTRUE);
const LauVetoes* vetoes = this->getVetoes();
if (vetoes) {
vetoOK = vetoes->passVeto(kinematics);
}
// Call this function recusively until we pass the veto.
if (vetoOK == kFALSE) {
this->generate();
}
return kTRUE;
}
void LauBkgndDPModel::fillDataTree(const LauFitDataTree& inputFitTree)
{
// In LauFitDataTree, the first two variables should always be
// m13^2 and m23^2. Other variables follow thus: charge.
Int_t nBranches = inputFitTree.nBranches();
if (nBranches < 2) {
cout<<"Error in LauBkgndDPModel::fillDataTree. Expecting at least 2 variables "
<<"in input data tree, but there are "<<nBranches<<"! Make sure you have "
<<"the right number of variables in your input data file!"<<endl;
return;
}
Double_t m13Sq(0.0), m23Sq(0.0), mPrime(0.0), thetaPrime(0.0);
UInt_t nEvents = inputFitTree.nEvents();
LauKinematics* kinematics = this->getKinematics();
// clear and resize the data vector
bgData_.clear();
bgData_.resize(nEvents);
for (UInt_t iEvt = 0; iEvt < nEvents; iEvt++) {
const LauFitData& dataValues = inputFitTree.getData(iEvt);
LauFitData::const_iterator iter = dataValues.find("m13Sq");
m13Sq = iter->second;
iter = dataValues.find("m23Sq");
m23Sq = iter->second;
// Update the kinematics. This will also update m' and theta' if squareDP = true
kinematics->updateKinematics(m13Sq, m23Sq);
if (squareDP_ == kTRUE) {
mPrime = kinematics->getmPrime();
thetaPrime = kinematics->getThetaPrime();
curEvtHistVal_ = this->calcHistValue(mPrime, thetaPrime);
} else {
curEvtHistVal_ = this->calcHistValue(m13Sq, m23Sq);
}
bgData_[iEvt] = curEvtHistVal_;
}
}
Double_t LauBkgndDPModel::getUnNormValue(UInt_t iEvt)
{
// Retrieve the likelihood for the given event
this->setDataEventNo(iEvt);
return curEvtHistVal_;
}
Double_t LauBkgndDPModel::getLikelihood(UInt_t iEvt)
{
// Retrieve the likelihood for the given event
this->setDataEventNo(iEvt);
Double_t llhd = curEvtHistVal_ / this->getPdfNorm();
return llhd;
}
void LauBkgndDPModel::setDataEventNo(UInt_t iEvt)
{
// Retrieve the data for event iEvt
if (bgData_.size() > iEvt) {
curEvtHistVal_ = bgData_[iEvt];
} else {
cerr<<"ERROR in LauBkgndDPModel::setDataEventNo : Event index too large: "<<iEvt<<" >= "<<bgData_.size()<<"."<<endl;
}
}

File Metadata

Mime Type
text/x-diff
Expires
Sun, Feb 23, 2:04 PM (2 h, 23 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4486450
Default Alt Text
(73 KB)

Event Timeline