Page MenuHomeHEPForge

No OneTemporary

diff --git a/doc/release.notes b/doc/release.notes
index d89e171..4a3747c 100644
--- a/doc/release.notes
+++ b/doc/release.notes
@@ -1,202 +1,206 @@
///////////////////////////////////////////////////////////////
/// ///
/// This is the History file for the Laura++ package. ///
/// ///
///////////////////////////////////////////////////////////////
+28th March 2014 Daniel Craik
+
+* Added support for asymmetric errors to Lau2DHistDP, Lau2DSplineDP and LauEffModel.
+
27th March 2014 Daniel Craik
* Changed histogram classes to use seeded random number generator for fluctuation and raising or lowering of bins and updated doxygen.
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Laura++ v2r0
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
8th March 2014 Thomas Latham
* Some additional functionality for the CoeffSet classes:
- allow the parameter values to be set (optionally setting the initial and generated values as well)
- allow the parameters to be set to float or to be fixed in the fit
These are needed when cloning but wanting some of the parameters to have different values and/or floating behaviour from the cloned set.
* Improve the printout of the setting of the coefficient values in the fit models and the creation of resonances
* Add LauFlatNR for the unform NR model - ends its rather strange special treatment
* Fix bug setting resAmpInt to 0 for LauPolNR
* Many other improvements to the info/warning/error printouts
* Modify GenFitBelleCPKpipi example to demonstrate cloning in action
* Add -Werror to compiler flags (treats warnings as errors)
5th March 2014 Thomas Latham
* Some improvements to LauPolNR to speed up amplitude calculation
2nd March 2014 Thomas Latham
* A number of updates to the CoeffSet classes:
- allow specification of the basename just after construction (before being given to the fit model)
- allow configuration of the parameter fit ranges (through static methods of base class)
- more adaptable cloning of the parameters (e.g. can only clone phase but allow magnitude to float independently)
- allow CP-violating parameters to be second-stage in Belle and CLEO formats
* Some improvements to the Doxygen and runtime information printouts
20th February 2014 Louis Henry
* Add LauPolNR - class for modelling the nonresonant contribution based on BaBar 3K model (arXiv:1201.5897)
6th February 2014 Thomas Latham
* Correct helicity convention information in Doxygen for LauKinematics
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Laura++ v1r2
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
5th February 2014 Thomas Latham
* Add rule to the Makefile that creates a rootmap file for the library
4th February 2014 Daniel Craik
* Fixed bug in Lau2DSplineDPPdf - normalisation was not being calculated
* Added out-of-range warning in LauBkgndDPModel and supressed excessive warnings
3rd February 2014 Mark Whitehead
(in branch for release in v2r0)
* Added a new class to allow parameters to be a function of other parameters
- inc/LauFormulaPar.hh
- src/LauFormulaPar.cc
28th January 2014 Daniel Craik
* Improved out-of-range efficiency warnings in LauEffModel and supressed excessive errors
* Modified LauIsobarDynamics to allow LASS parameters to be configured for LauLASSBWRes and
LauLASSNRRes
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
- LauBkgndDPModel, modified to use Lau2DAbsDPPdf instead of Lau2DHistDPPdf
- setBkgndSpline method added to LauBkgndDPModel to allow use of splines
22nd January 2014 Thomas Latham
* Improve some error checks and corresponding warning messages in
LauCPFitModel::setSignalDPParameters
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
18th December 2013 Mark Whitehead
(in branch for release in v2r0)
* Added functionality to include Gaussian constraints on floated
parameters in the fit.
The files updated are:
- inc/LauAbsFitModel.hh
- inc/LauParameter.hh
- src/LauAbsFitModel.cc
- src/LauParameter.cc
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
3rd December 2013 Thomas Latham
(in branch for release in v2r0)
* Have the master save an ntuple with all fitter parameters and the full correlation matrix information.
29th November 2013 Thomas Latham
* Fixed bug in ResultsExtractor where the output file was not written
29th November 2013 Thomas Latham
(in branch for release in v2r0)
* Allow the slave ntuples to store the partial covariance matrices in the simultaneous fitting
26th November 2013 Thomas Latham
(in branch for release in v2r0)
* Added first version of the simultaneous fitting framework
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
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/Lau2DAbsHistDP.hh b/inc/Lau2DAbsHistDP.hh
index 76b4953..6962a32 100644
--- a/inc/Lau2DAbsHistDP.hh
+++ b/inc/Lau2DAbsHistDP.hh
@@ -1,117 +1,125 @@
// 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 Lau2DAbsHistDP.hh
\brief File containing declaration of Lau2DAbsHistDP class.
*/
/*! \class Lau2DAbsHistDP
\brief Abstract base class for defining a variation across a 2D DP based on a histogram.
Abstract base class for defining an unnormalised variation across a 2D DP based on a histogram.
Contains helper methods to vary bin contents in the input histogram within uncertainties and
scale the input histogram to match a given average value.
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
#define LAU_2DABSHIST_DP
#include "Lau2DAbsDP.hh"
class TH2;
class LauDaughters;
class LauKinematics;
class Lau2DAbsHistDP : public Lau2DAbsDP {
public:
//! Constructor
/*!
\param [in] daughters the daughter particles
\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
*/
Lau2DAbsHistDP(const LauDaughters* daughters, Bool_t useUpperHalfOnly = kFALSE, Bool_t squareDP = kFALSE);
//! Copy constructor
//Lau2DAbsHistDP( const Lau2DAbsHistDP& rhs );
//! Destructor
virtual ~Lau2DAbsHistDP();
//! Perform the interpolation
/*!
\param [in] x the x-axis value
\param [in] y the y-axis value
\return the interpolated histogram value
*/
virtual Double_t interpolateXY(Double_t x, Double_t y) const=0;
protected:
//! Fluctuate the contents of each histogram bin independently, in accordance with their errors
/*!
\param [in,out] hist the histogram
*/
void doBinFluctuation(TH2* hist);
+ //! Fluctuate the contents of each histogram bin independently, in accordance with their errors
+ /*!
+ \param [in,out] hist the histogram
+ \param [in] errorHi the histogram containing the upper uncertainties on the bins
+ \param [in] errorLo the histogram containing the lower uncertainties on the bins
+ */
+ void doBinFluctuation(TH2* hist, TH2* errorHi, TH2* errorLo);
+
//! Rescale the histogram bin contents based on the desired average efficiency and its uncertainty
/*!
The new average is sampled from a Gaussian distribution G(x;avEff,avEffError).
The histogram is then rescaled according to newAvg/oldAvg.
\param [in,out] hist the histogram
\param [in] avEff the desired average efficiency
\param [in] avEffError the error on that efficiency
*/
void raiseOrLowerBins(TH2* hist, Double_t avEff, Double_t avEffError);
//! Compute the average bin content for bins within the kinematic boundary
/*!
This method just uses the raw bin contents with no interpolation
\param [in,out] hist the histogram
\return the average value over the DP
*/
Double_t computeAverageContents(TH2 const * const hist) const;
//! Check whether the given co-ordinates are within the kinematic boundary
/*!
\param [in] x the x co-ordinate
\param [in] y the y co-ordinate
\return true if the co-ordinates are within the kinematic boundary, otherwise false
*/
Bool_t withinDPBoundaries(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
Lau2DAbsHistDP( const Lau2DAbsHistDP& rhs );
//! Copy assignment operator - not implemented
Lau2DAbsHistDP& operator=(const Lau2DAbsHistDP& rhs);
//! Kinematics used to check events are within DP boundary
const LauKinematics* kinematics_;
//! Boolean for using the upper half of DP
Bool_t upperHalf_;
//! Boolean for using square DP variables
Bool_t squareDP_;
ClassDef(Lau2DAbsHistDP,0) // Abstract base class for 2D DP variations based on a histogram
};
#endif
diff --git a/inc/Lau2DHistDP.hh b/inc/Lau2DHistDP.hh
index 8bb82fc..66931ac 100644
--- a/inc/Lau2DHistDP.hh
+++ b/inc/Lau2DHistDP.hh
@@ -1,113 +1,136 @@
// 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 Lau2DHistDP.hh
\brief File containing declaration of Lau2DHistDP class.
*/
/*! \class Lau2DHistDP
\brief Class for defining a 2D DP histogram.
Class for defining a 2D DP histogram.
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 not normalised to the total histogram area
(useful for efficiency histograms for example).
The histogram can be defined in the conventional DP (m13Sq vs m23Sq) or in the square DP.
*/
#ifndef LAU_2DHIST_DP
#define LAU_2DHIST_DP
#include "Lau2DAbsHistDP.hh"
class TH2;
class LauDaughters;
class LauKinematics;
class Lau2DHistDP : public Lau2DAbsHistDP {
public:
//! Constructor
/*!
\param [in] hist the 2D DP histogram
\param [in] daughters the daughter particles
\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).
The seed for the random number generator used to fluctuate the bins should first be set using LauRandom::setSeed.
\param [in] avEff the desired average efficiency - see Lau2DHistDP::raiseOrLowerBins.
The seed for the random number generator used to raise or lower the bins should first be set using LauRandom::setSeed.
\param [in] avEffError the error on that efficiency - see Lau2DHistDP::raiseOrLowerBins
\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
*/
Lau2DHistDP(const TH2* hist, const LauDaughters* daughters,
Bool_t useInterpolation = kTRUE, Bool_t fluctuateBins = kFALSE,
Double_t avEff = -1.0, Double_t avEffError = -1.0, Bool_t useUpperHalfOnly = kFALSE, Bool_t squareDP = kFALSE);
+ //! Constructor
+ /*!
+ \param [in] hist the 2D DP histogram
+ \param [in] errorHi the 2D DP histogram containing the upper uncertainty
+ \param [in] errorLo the 2D DP histogram containing the lower uncertainty
+ \param [in] daughters the daughter particles
+ \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).
+ The seed for the random number generator used to fluctuate the bins should first be set using LauRandom::setSeed.
+ \param [in] avEff the desired average efficiency - see Lau2DHistDP::raiseOrLowerBins.
+ The seed for the random number generator used to raise or lower the bins should first be set using LauRandom::setSeed.
+ \param [in] avEffError the error on that efficiency - see Lau2DHistDP::raiseOrLowerBins
+ \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
+ */
+ Lau2DHistDP(const TH2* hist, const TH2* errorHi, const TH2* errorLo, const LauDaughters* daughters,
+ Bool_t useInterpolation = kTRUE, Bool_t fluctuateBins = kFALSE,
+ Double_t avEff = -1.0, Double_t avEffError = -1.0, Bool_t useUpperHalfOnly = kFALSE, Bool_t squareDP = kFALSE);
+
//! Destructor
virtual ~Lau2DHistDP();
//! Perform the interpolation
/*!
\param [in] x the x-axis value
\param [in] y the y-aixs value
\return the interpolated histogram value
*/
Double_t interpolateXY(Double_t x, Double_t y) const;
protected:
//! Get the raw 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;
private:
//! Copy constructor - not implemented
Lau2DHistDP( const Lau2DHistDP& rhs );
//! Copy assignment operator - not implemented
Lau2DHistDP& operator=(const Lau2DHistDP& rhs);
//! The underlying histogram
TH2* hist_;
+ //! The histogram containing the upper errors
+ TH2* errorHi_;
+ //! The histogram containing the lower errors
+ TH2* errorLo_;
//! 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 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_;
//! Control boolean for using the linear interpolation
Bool_t useInterpolation_;
ClassDef(Lau2DHistDP,0) // 2D Histogram utility class for DP analyses
};
#endif
diff --git a/inc/Lau2DSplineDP.hh b/inc/Lau2DSplineDP.hh
index 2e52e68..132e527 100644
--- a/inc/Lau2DSplineDP.hh
+++ b/inc/Lau2DSplineDP.hh
@@ -1,78 +1,96 @@
// 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 Lau2DSplineDP.hh
\brief File containing declaration of Lau2DSplineDP class.
*/
/*! \class Lau2DSplineDP
\brief Class for defining variations across a 2D DP using a spline.
Class for defining variations across a 2D DP using a spline.
Employs a 2D cubic spline to get the histogram value based on an input histogram.
The returned values are not normalised to the total histogram area
(useful for efficiency histograms for example).
The histogram can be defined in the conventional DP (m13Sq vs m23Sq) or in the square DP.
*/
#ifndef LAU_2DSPLINE_DP
#define LAU_2DSPLINE_DP
#include "Lau2DAbsHistDP.hh"
class TH2;
class Lau2DCubicSpline;
class LauDaughters;
class LauKinematics;
class Lau2DSplineDP : public Lau2DAbsHistDP {
public:
//! Constructor
/*!
\param [in] hist the 2D DP histogram
\param [in] daughters the daughter particles
\param [in] fluctuateBins boolean flag to determine whether the bin contents should be fluctuated in accordance with their errors (useful for systematic error evaluation).
The seed for the random number generator used to fluctuate the bins should first be set using LauRandom::setSeed.
\param [in] avEff the desired average efficiency - see Lau2DSplineDP::raiseOrLowerBins.
The seed for the random number generator used to raise or lower the bins should first be set using LauRandom::setSeed.
\param [in] avEffError the error on that efficiency - see Lau2DSplineDP::raiseOrLowerBins
\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
*/
Lau2DSplineDP(const TH2* hist, const LauDaughters* daughters,
Bool_t fluctuateBins = kFALSE, Double_t avEff = -1.0,
Double_t avEffError = -1.0, Bool_t useUpperHalfOnly = kFALSE, Bool_t squareDP = kFALSE);
+ //! Constructor
+ /*!
+ \param [in] hist the 2D DP histogram
+ \param [in] errorHi the 2D DP histogram containing the upper uncertainty
+ \param [in] errorLo the 2D DP histogram containing the lower uncertainty
+ \param [in] daughters the daughter particles
+ \param [in] fluctuateBins boolean flag to determine whether the bin contents should be fluctuated in accordance with their errors (useful for systematic error evaluation).
+ The seed for the random number generator used to fluctuate the bins should first be set using LauRandom::setSeed.
+ \param [in] avEff the desired average efficiency - see Lau2DSplineDP::raiseOrLowerBins.
+ The seed for the random number generator used to raise or lower the bins should first be set using LauRandom::setSeed.
+ \param [in] avEffError the error on that efficiency - see Lau2DSplineDP::raiseOrLowerBins
+ \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
+ */
+ Lau2DSplineDP(const TH2* hist, const TH2* errorHi, const TH2* errorLo, const LauDaughters* daughters,
+ Bool_t fluctuateBins = kFALSE, Double_t avEff = -1.0,
+ Double_t avEffError = -1.0, Bool_t useUpperHalfOnly = kFALSE, Bool_t squareDP = kFALSE);
+
//! Destructor
virtual ~Lau2DSplineDP();
//! Perform the interpolation
/*!
\param [in] x the x-axis value
\param [in] y the y-aixs value
\return the interpolated histogram value
*/
Double_t interpolateXY(Double_t x, Double_t y) const;
private:
//! Copy constructor - not implemented
Lau2DSplineDP( const Lau2DSplineDP& rhs );
//! Copy assignment operator - not implemented
Lau2DSplineDP& operator=(const Lau2DSplineDP& rhs);
//! A 2D cubic spline generated from the histogram
Lau2DCubicSpline* spline_;
ClassDef(Lau2DSplineDP,0) // 2D Spline utility class for DP analyses
};
#endif
diff --git a/inc/LauEffModel.hh b/inc/LauEffModel.hh
index 8ef388d..a33c09f 100644
--- a/inc/LauEffModel.hh
+++ b/inc/LauEffModel.hh
@@ -1,141 +1,182 @@
// 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 LauEffModel.hh
\brief File containing declaration of LauEffModel class.
*/
/*! \class LauEffModel
\brief Class that implements the efficiency description across the signal Dalitz plot.
Class that defines the efficiency model variation across the signal Dalitz plot.
The method is based in a predetermined two-dimensional histogram to characterize the phase space acceptance.
The efficiency variation is defined in terms of x = m_13^2, y = m_23^2 for the Dalitz plot (default)
or x = m', y = theta' for the square Dalitz plot
*/
#ifndef LAUEFFMODEL
#define LAUEFFMODEL
#include "Rtypes.h"
class TH2;
class LauDaughters;
class LauKinematics;
class LauVetoes;
class Lau2DAbsDP;
class LauEffModel {
public:
//! Constructor
/*!
\param [in] daughters the daughters particles of the Dalitz plot model
\param [in] vetoes the object describing the vetoes applied in the phase space
*/
LauEffModel(const LauDaughters* daughters, const LauVetoes* vetoes);
//! Destructor
virtual ~LauEffModel();
//! Set the efficiency variation across the phase space using a predetermined 2D histogram.
/*!
The efficiency is defined in terms of x = m_13^2, y = m_23^2 or x = m', y = theta' for the square Dalitz plot
\param [in] effHisto the 2-dimensional histogram that describes the efficiency variation
\param [in] useInterpolation boolean flag decision to switch on/off 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.
The seed for the random number generator used to fluctuate the bins should first be set using LauRandom::setSeed.
\param [in] avEff the desired average efficiency - see Lau2DHistDP::raiseOrLowerBins, values less than zero switch off this behaviour.
The seed for the random number generator used to raise or lower the bins should first be set using LauRandom::setSeed.
\param [in] absError the error on that efficiency - see Lau2DHistDP::raiseOrLowerBins, values less than zero switch off this behaviour
\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 setEffHisto(const TH2* effHisto,
Bool_t useInterpolation = kTRUE, Bool_t fluctuateBins = kFALSE,
Double_t avEff = -1.0, Double_t absError = -1.0,
Bool_t useUpperHalfOnly = kFALSE, Bool_t squareDP = kFALSE);
+ //! Set the efficiency variation across the phase space using a predetermined 2D histogram.
+ /*!
+ The efficiency is defined in terms of x = m_13^2, y = m_23^2 or x = m', y = theta' for the square Dalitz plot
+
+ \param [in] effHisto the 2-dimensional histogram that describes the efficiency variation
+ \param [in] errorHi the 2-dimensional histogram that describes the upper uncertainty on the efficiency variation
+ \param [in] errorLo the 2-dimensional histogram that describes the lower uncertainty on the efficiency variation
+ \param [in] useInterpolation boolean flag decision to switch on/off 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.
+ The seed for the random number generator used to fluctuate the bins should first be set using LauRandom::setSeed.
+ \param [in] avEff the desired average efficiency - see Lau2DHistDP::raiseOrLowerBins, values less than zero switch off this behaviour.
+ The seed for the random number generator used to raise or lower the bins should first be set using LauRandom::setSeed.
+ \param [in] absError the error on that efficiency - see Lau2DHistDP::raiseOrLowerBins, values less than zero switch off this behaviour
+ \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 setEffHisto(const TH2* effHisto, const TH2* errorHi, const TH2* errorLo,
+ Bool_t useInterpolation = kTRUE, Bool_t fluctuateBins = kFALSE,
+ Double_t avEff = -1.0, Double_t absError = -1.0,
+ Bool_t useUpperHalfOnly = kFALSE, Bool_t squareDP = kFALSE);
+
//! Set the efficiency variation across the phase space using a spline based on a predetermined 2D histogram.
/*!
The efficiency is defined in terms of x = m_13^2, y = m_23^2 or x = m', y = theta' for the square Dalitz plot
\param [in] effHisto the 2-dimensional histogram that describes the efficiency variation
\param [in] fluctuateBins boolean flag to determine whether the bin contents should be fluctuated in accordance with their errors.
The seed for the random number generator used to fluctuate the bins should first be set using LauRandom::setSeed.
\param [in] avEff the desired average efficiency - see Lau2DHistDP::raiseOrLowerBins, values less than zero switch off this behaviour.
The seed for the random number generator used to raise or lower the bins should first be set using LauRandom::setSeed.
\param [in] absError the error on that efficiency - see Lau2DHistDP::raiseOrLowerBins, values less than zero switch off this behaviour
\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 setEffSpline(const TH2* effHisto,
Bool_t fluctuateBins = kFALSE,
Double_t avEff = -1.0, Double_t absError = -1.0,
Bool_t useUpperHalfOnly = kFALSE, Bool_t squareDP = kFALSE);
+ //! Set the efficiency variation across the phase space using a spline based on a predetermined 2D histogram.
+ /*!
+ The efficiency is defined in terms of x = m_13^2, y = m_23^2 or x = m', y = theta' for the square Dalitz plot
+
+ \param [in] effHisto the 2-dimensional histogram that describes the efficiency variation
+ \param [in] errorHi the 2-dimensional histogram that describes the upper uncertainty on the efficiency variation
+ \param [in] errorLo the 2-dimensional histogram that describes the lower uncertainty on the efficiency variation
+ \param [in] fluctuateBins boolean flag to determine whether the bin contents should be fluctuated in accordance with their errors.
+ The seed for the random number generator used to fluctuate the bins should first be set using LauRandom::setSeed.
+ \param [in] avEff the desired average efficiency - see Lau2DHistDP::raiseOrLowerBins, values less than zero switch off this behaviour.
+ The seed for the random number generator used to raise or lower the bins should first be set using LauRandom::setSeed.
+ \param [in] absError the error on that efficiency - see Lau2DHistDP::raiseOrLowerBins, values less than zero switch off this behaviour
+ \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 setEffSpline(const TH2* effHisto, const TH2* errorHi, const TH2* errorLo,
+ Bool_t fluctuateBins = kFALSE,
+ Double_t avEff = -1.0, Double_t absError = -1.0,
+ Bool_t useUpperHalfOnly = kFALSE, Bool_t squareDP = kFALSE);
+
//! Determine the efficiency for a given point in the Dalitz plot.
/*!
The method uses the 2D histogram set by the setEffHisto() function and the vetoes information.
\param [in] kinematics the object that defines the DP position
\return the efficiency value at the given point in the DP
*/
Double_t calcEfficiency( const LauKinematics* kinematics ) const;
//! Determine whether the given DP position is outside the vetoes
/*!
\param [in] kinematics the object that defines the DP position
\return kTRUE if the DP position is outside all veto regions, kFALSE otherwise
*/
Bool_t passVeto( const LauKinematics* kinematics ) const;
//! Determine whether the efficiency histogram has had its bins fluctuated within their errors
Bool_t fluctuateEffHisto() const {return fluctuateEffHisto_;}
private:
//! Copy constructor - not implemented
LauEffModel( const LauEffModel& rhs );
//! Get the efficiency from a two-dimensional histogram by interpolating in x and y
/*!
\param [in] xVal the value to be interpolated in x
\param [in] yVal the value to be interpolated in y
*/
Double_t getEffHistValue(Double_t xVal, Double_t yVal) const;
//! The daughters object
const LauDaughters* daughters_;
//! The vetoes object
const LauVetoes* vetoes_;
//! The efficiency histogram object
Lau2DAbsDP* effHisto_;
//! Use of the square Dalitz plot
Bool_t squareDP_;
//! Fluctuate histogram within the error
Bool_t fluctuateEffHisto_;
//! Flag to track whether a warning has been issued for bin values less than zero
mutable Bool_t lowBinWarningIssued_;
//! Flag to track whether a warning has been issued for bin values greater than one
mutable Bool_t highBinWarningIssued_;
ClassDef(LauEffModel, 0) // Implement the signal efficiency across the DP
};
#endif
diff --git a/src/Lau2DAbsHistDP.cc b/src/Lau2DAbsHistDP.cc
index 93226a8..8d1254d 100644
--- a/src/Lau2DAbsHistDP.cc
+++ b/src/Lau2DAbsHistDP.cc
@@ -1,123 +1,150 @@
// 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 Lau2DAbsHistDP.cc
\brief File containing implementation of Lau2DAbsHistDP class.
*/
#include <iostream>
#include "TAxis.h"
#include "TH2.h"
#include "TRandom.h"
#include "TSystem.h"
#include "Lau2DAbsHistDP.hh"
#include "LauDaughters.hh"
#include "LauKinematics.hh"
#include "LauRandom.hh"
ClassImp(Lau2DAbsHistDP)
Lau2DAbsHistDP::Lau2DAbsHistDP(const LauDaughters* daughters, Bool_t useUpperHalfOnly, Bool_t squareDP) :
kinematics_( (daughters!=0) ? daughters->getKinematics() : 0 ),
upperHalf_(useUpperHalfOnly),
squareDP_(squareDP)
{
}
Lau2DAbsHistDP::~Lau2DAbsHistDP()
{
}
void Lau2DAbsHistDP::doBinFluctuation(TH2* hist)
{
TRandom* random = LauRandom::randomFun();
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);
}
}
}
}
+void Lau2DAbsHistDP::doBinFluctuation(TH2* hist, TH2* errorHi, TH2* errorLo)
+{
+ TRandom* random = LauRandom::randomFun();
+
+ 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 currentErrorHi = errorHi->GetBinContent(i+1,j+1);
+ Double_t currentErrorLo = errorLo->GetBinContent(i+1,j+1);
+ Double_t shift = random->Gaus(0.,1.);
+ Double_t newContent(currentContent);
+
+ if(shift>0) newContent += shift*currentErrorHi;
+ else newContent += shift*currentErrorLo;
+
+ if (newContent<0.0) {
+ hist->SetBinContent(i+1,j+1,0.0);
+ } else {
+ hist->SetBinContent(i+1,j+1,newContent);
+ }
+ }
+ }
+}
+
void Lau2DAbsHistDP::raiseOrLowerBins(TH2* hist, Double_t avEff, Double_t avEffError)
{
TRandom* random = LauRandom::randomFun();
Double_t curAvg = this->computeAverageContents(hist);
Double_t newAvg = random->Gaus(avEff,avEffError);
hist->Scale( newAvg / curAvg );
}
Double_t Lau2DAbsHistDP::computeAverageContents(TH2 const * const hist) const
{
Double_t totalContent(0.0);
Int_t binsWithinDPBoundary(0);
Int_t nBinsX = static_cast<Int_t>(hist->GetNbinsX());
Int_t nBinsY = static_cast<Int_t>(hist->GetNbinsY());
// Loop through the bins and include any that have their centre or any
// of the four corners within the kinematic boundary
for ( Int_t i(0); i<nBinsX; ++i ) {
Double_t binXCentre = hist->GetXaxis()->GetBinCenter(i+1);
Double_t binXLowerEdge = hist->GetXaxis()->GetBinLowEdge(i+1);
Double_t binXUpperEdge = hist->GetXaxis()->GetBinUpEdge(i+1);
for ( Int_t j(0); j<nBinsY; ++j ) {
Double_t binYCentre = hist->GetYaxis()->GetBinCenter(i+1);
Double_t binYLowerEdge = hist->GetYaxis()->GetBinLowEdge(i+1);
Double_t binYUpperEdge = hist->GetYaxis()->GetBinUpEdge(i+1);
if ( withinDPBoundaries( binXCentre, binYCentre ) ||
withinDPBoundaries( binXLowerEdge, binYLowerEdge ) ||
withinDPBoundaries( binXUpperEdge, binYUpperEdge ) ||
withinDPBoundaries( binXLowerEdge, binYUpperEdge ) ||
withinDPBoundaries( binXUpperEdge, binYLowerEdge ) ) {
totalContent += hist->GetBinContent(i+1, j+1);
++binsWithinDPBoundary;
}
}
}
return totalContent/binsWithinDPBoundary;
}
Bool_t Lau2DAbsHistDP::withinDPBoundaries(Double_t x, Double_t y) const
{
return squareDP_ ? kinematics_->withinSqDPLimits(x,y) : kinematics_->withinDPLimits(x,y);
}
void Lau2DAbsHistDP::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;
}
}
}
diff --git a/src/Lau2DHistDP.cc b/src/Lau2DHistDP.cc
index d414aa9..3755221 100644
--- a/src/Lau2DHistDP.cc
+++ b/src/Lau2DHistDP.cc
@@ -1,242 +1,346 @@
// 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 Lau2DHistDP.cc
\brief File containing implementation of Lau2DHistDP class.
*/
#include <iostream>
#include "TAxis.h"
#include "TH2.h"
#include "TRandom.h"
#include "TSystem.h"
#include "Lau2DHistDP.hh"
#include "LauDaughters.hh"
#include "LauKinematics.hh"
#include "LauRandom.hh"
ClassImp(Lau2DHistDP)
Lau2DHistDP::Lau2DHistDP(const TH2* hist, const LauDaughters* daughters,
Bool_t useInterpolation, Bool_t fluctuateBins,
Double_t avEff, Double_t avEffError, Bool_t useUpperHalfOnly, Bool_t squareDP) :
Lau2DAbsHistDP(daughters,useUpperHalfOnly,squareDP),
hist_(hist ? dynamic_cast<TH2*>(hist->Clone()) : 0),
+ errorHi_(0), errorLo_(0),
minX_(0.0), maxX_(0.0),
minY_(0.0), maxY_(0.0),
rangeX_(0.0), rangeY_(0.0),
binXWidth_(0.0), binYWidth_(0.0),
nBinsX_(0), nBinsY_(0),
useInterpolation_(useInterpolation)
{
if ( ! hist_ ) {
std::cerr << "ERROR in Lau2DHistDP constructor : the histogram pointer is null." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
// Save various attributes of the histogram
// (axis ranges, number of bins, areas)
TAxis* xAxis = hist_->GetXaxis();
minX_ = static_cast<Double_t>(xAxis->GetXmin());
maxX_ = static_cast<Double_t>(xAxis->GetXmax());
rangeX_ = maxX_ - minX_;
TAxis* yAxis = hist_->GetYaxis();
minY_ = static_cast<Double_t>(yAxis->GetXmin());
maxY_ = static_cast<Double_t>(yAxis->GetXmax());
rangeY_ = maxY_ - minY_;
nBinsX_ = static_cast<Int_t>(hist_->GetNbinsX());
nBinsY_ = static_cast<Int_t>(hist_->GetNbinsY());
binXWidth_ = static_cast<Double_t>(TMath::Abs(rangeX_)/(nBinsX_*1.0));
binYWidth_ = static_cast<Double_t>(TMath::Abs(rangeY_)/(nBinsY_*1.0));
if (fluctuateBins) {
this->doBinFluctuation(hist_);
}
if (avEff > 0.0 && avEffError > 0.0) {
this->raiseOrLowerBins(hist_,avEff,avEffError);
}
}
+Lau2DHistDP::Lau2DHistDP(const TH2* hist, const TH2* errorHi, const TH2* errorLo, const LauDaughters* daughters,
+ Bool_t useInterpolation, Bool_t fluctuateBins,
+ Double_t avEff, Double_t avEffError, Bool_t useUpperHalfOnly, Bool_t squareDP) :
+ Lau2DAbsHistDP(daughters,useUpperHalfOnly,squareDP),
+ hist_(hist ? dynamic_cast<TH2*>(hist->Clone()) : 0),
+ errorHi_(errorHi ? dynamic_cast<TH2*>(errorHi->Clone()) : 0),
+ errorLo_(errorLo ? dynamic_cast<TH2*>(errorLo->Clone()) : 0),
+ minX_(0.0), maxX_(0.0),
+ minY_(0.0), maxY_(0.0),
+ rangeX_(0.0), rangeY_(0.0),
+ binXWidth_(0.0), binYWidth_(0.0),
+ nBinsX_(0), nBinsY_(0),
+ useInterpolation_(useInterpolation)
+{
+ if ( ! hist_ ) {
+ std::cerr << "ERROR in Lau2DHistDP constructor : the histogram pointer is null." << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+ if ( ! errorHi_ ) {
+ std::cerr << "ERROR in Lau2DHistDP constructor : the upper error histogram pointer is null." << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+ if ( ! errorLo_ ) {
+ std::cerr << "ERROR in Lau2DHistDP constructor : the lower error histogram pointer is null." << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+
+ // Save various attributes of the histogram
+ // (axis ranges, number of bins, areas)
+ TAxis* xAxis = hist_->GetXaxis();
+ minX_ = static_cast<Double_t>(xAxis->GetXmin());
+ maxX_ = static_cast<Double_t>(xAxis->GetXmax());
+ rangeX_ = maxX_ - minX_;
+
+ TAxis* yAxis = hist_->GetYaxis();
+ minY_ = static_cast<Double_t>(yAxis->GetXmin());
+ maxY_ = static_cast<Double_t>(yAxis->GetXmax());
+ rangeY_ = maxY_ - minY_;
+
+ nBinsX_ = static_cast<Int_t>(hist_->GetNbinsX());
+ nBinsY_ = static_cast<Int_t>(hist_->GetNbinsY());
+
+ binXWidth_ = static_cast<Double_t>(TMath::Abs(rangeX_)/(nBinsX_*1.0));
+ binYWidth_ = static_cast<Double_t>(TMath::Abs(rangeY_)/(nBinsY_*1.0));
+
+ if(static_cast<Int_t>(errorLo_->GetNbinsX()) != nBinsX_ ||
+ static_cast<Int_t>(errorLo_->GetNbinsY()) != nBinsY_) {
+ std::cerr << "ERROR in Lau2DHistDP constructor : the lower error histogram has a different number of bins to the main histogram." << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+
+ if(static_cast<Int_t>(errorHi_->GetNbinsX()) != nBinsX_ ||
+ static_cast<Int_t>(errorHi_->GetNbinsY()) != nBinsY_) {
+ std::cerr << "ERROR in Lau2DHistDP constructor : the upper error histogram has a different number of bins to the main histogram." << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+
+ xAxis = errorLo_->GetXaxis();
+ yAxis = errorLo_->GetYaxis();
+
+ if(static_cast<Double_t>(xAxis->GetXmin()) != minX_ ||
+ static_cast<Double_t>(xAxis->GetXmax()) != maxX_) {
+ std::cerr << "ERROR in Lau2DHistDP constructor : the lower error histogram has a different x range to the main histogram." << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+
+ if(static_cast<Double_t>(yAxis->GetXmin()) != minY_ ||
+ static_cast<Double_t>(yAxis->GetXmax()) != maxY_) {
+ std::cerr << "ERROR in Lau2DHistDP constructor : the lower error histogram has a different y range to the main histogram." << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+
+ xAxis = errorHi_->GetXaxis();
+ yAxis = errorHi_->GetYaxis();
+
+ if(static_cast<Double_t>(xAxis->GetXmin()) != minX_ ||
+ static_cast<Double_t>(xAxis->GetXmax()) != maxX_) {
+ std::cerr << "ERROR in Lau2DHistDP constructor : the upper error histogram has a different x range to the main histogram." << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+
+ if(static_cast<Double_t>(yAxis->GetXmin()) != minY_ ||
+ static_cast<Double_t>(yAxis->GetXmax()) != maxY_) {
+ std::cerr << "ERROR in Lau2DHistDP constructor : the upper error histogram has a different y range to the main histogram." << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+
+ if (fluctuateBins) {
+ this->doBinFluctuation(hist_,errorHi_,errorLo_);
+ }
+ if (avEff > 0.0 && avEffError > 0.0) {
+ this->raiseOrLowerBins(hist_,avEff,avEffError);
+ }
+}
+
Lau2DHistDP::~Lau2DHistDP()
{
delete hist_;
hist_ = 0;
+ if(errorHi_) {
+ delete errorHi_;
+ errorHi_=0;
+ }
+ if(errorLo_) {
+ delete errorLo_;
+ errorLo_=0;
+ }
}
Double_t Lau2DHistDP::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;
}
Double_t value = hist_->GetBinContent(xBinNo+1, yBinNo+1);
return value;
}
Double_t Lau2DHistDP::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 Lau2DHistDP::interpolateXY : Given position is outside the DP boundary, returning 0.0." << std::endl;
return 0.0;
}
// Find the 2D histogram bin for x and y
Int_t i = Int_t((x - minX_)*nBinsX_/rangeX_);
Int_t j = Int_t((y - minY_)*nBinsY_/rangeY_);
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 getBinHistValue(i,j);}
// Find the bin centres (actual co-ordinate positions, not histogram indices)
Double_t cbinx = Double_t(i+0.5)*rangeX_/nBinsX_ + minX_;
Double_t cbiny = Double_t(j+0.5)*rangeY_/nBinsY_ + minY_;
// If bin centres are outside kinematic region, do not extrapolate
if (withinDPBoundaries(cbinx,cbiny) == kFALSE) {return 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 = 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)*rangeX_/nBinsX_ + minX_;
if (withinDPBoundaries(cbinx_adj, y) == kFALSE) {
// The adjacent bin is outside the DP range. Don't extrapolate.
value = 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 = getBinHistValue(i,j);
Double_t value2 = 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)*rangeY_/nBinsY_ + minY_;
if (withinDPBoundaries(x, cbiny_adj) == kFALSE) {
// The adjacent bin is outside the DP range. Don't extrapolate.
value = 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 = getBinHistValue(i,j);
Double_t value2 = 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)*rangeX_/nBinsX_ + minX_;
Double_t cbiny_adj = Double_t(j_adj+0.5)*rangeY_/nBinsY_ + minY_;
if (withinDPBoundaries(cbinx_adj, cbiny_adj) == kFALSE) {
// The adjacent bin is outside the DP range. Don't extrapolate.
value = 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 = getBinHistValue(i,j);
Double_t value2 = getBinHistValue(i_adj,j);
Double_t value3 = getBinHistValue(i,j_adj);
Double_t value4 = getBinHistValue(i_adj,j_adj);
value = value1*dx1*dy1 + value2*dx0*dy1 + value3*dx1*dy0 + value4*dx0*dy0;
value /= inter_denom;
}
}
return value;
}
diff --git a/src/Lau2DSplineDP.cc b/src/Lau2DSplineDP.cc
index ba0f302..2250b9a 100644
--- a/src/Lau2DSplineDP.cc
+++ b/src/Lau2DSplineDP.cc
@@ -1,83 +1,175 @@
// 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 Lau2DSplineDP.cc
\brief File containing implementation of Lau2DSplineDP class.
*/
#include <iostream>
#include "TAxis.h"
#include "TH2.h"
#include "TRandom.h"
#include "TSystem.h"
#include "Lau2DSplineDP.hh"
#include "Lau2DCubicSpline.hh"
#include "LauDaughters.hh"
#include "LauKinematics.hh"
#include "LauRandom.hh"
ClassImp(Lau2DSplineDP)
Lau2DSplineDP::Lau2DSplineDP(const TH2* hist, const LauDaughters* daughters,
Bool_t fluctuateBins, Double_t avEff, Double_t avEffError,
Bool_t useUpperHalfOnly, Bool_t squareDP) :
Lau2DAbsHistDP(daughters,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 Lau2DSplineDP constructor : the histogram pointer is null." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
if (fluctuateBins) {
this->doBinFluctuation(tempHist);
}
if (avEff > 0.0 && avEffError > 0.0) {
this->raiseOrLowerBins(tempHist,avEff,avEffError);
}
spline_ = new Lau2DCubicSpline(*tempHist);
delete tempHist;
}
+Lau2DSplineDP::Lau2DSplineDP(const TH2* hist, const TH2* errorHi, const TH2* errorLo, const LauDaughters* daughters,
+ Bool_t fluctuateBins, Double_t avEff, Double_t avEffError,
+ Bool_t useUpperHalfOnly, Bool_t squareDP) :
+ Lau2DAbsHistDP(daughters,useUpperHalfOnly,squareDP),
+ spline_(0)
+{
+ //We may need to modify the histogram so clone it
+ TH2* tempHist(hist ? dynamic_cast<TH2*>(hist->Clone()) : 0);
+ TH2* tempErrorHi(errorHi ? dynamic_cast<TH2*>(errorHi->Clone()) : 0);
+ TH2* tempErrorLo(errorLo ? dynamic_cast<TH2*>(errorLo->Clone()) : 0);
+
+ if ( ! tempHist ) {
+ std::cerr << "ERROR in Lau2DSplineDP constructor : the histogram pointer is null." << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+ if ( ! tempErrorHi ) {
+ std::cerr << "ERROR in Lau2DHistDP constructor : the upper error histogram pointer is null." << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+ if ( ! tempErrorLo ) {
+ std::cerr << "ERROR in Lau2DHistDP constructor : the lower error histogram pointer is null." << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+
+ TAxis* xAxis = tempHist->GetXaxis();
+ Double_t minX = static_cast<Double_t>(xAxis->GetXmin());
+ Double_t maxX = static_cast<Double_t>(xAxis->GetXmax());
+
+ TAxis* yAxis = tempHist->GetYaxis();
+ Double_t minY = static_cast<Double_t>(yAxis->GetXmin());
+ Double_t maxY = static_cast<Double_t>(yAxis->GetXmax());
+
+ Int_t nBinsX = static_cast<Int_t>(tempHist->GetNbinsX());
+ Int_t nBinsY = static_cast<Int_t>(tempHist->GetNbinsY());
+
+ if(static_cast<Int_t>(tempErrorLo->GetNbinsX()) != nBinsX ||
+ static_cast<Int_t>(tempErrorLo->GetNbinsY()) != nBinsY) {
+ std::cerr << "ERROR in Lau2DHistDP constructor : the lower error histogram has a different number of bins to the main histogram." << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+
+ if(static_cast<Int_t>(tempErrorHi->GetNbinsX()) != nBinsX ||
+ static_cast<Int_t>(tempErrorHi->GetNbinsY()) != nBinsY) {
+ std::cerr << "ERROR in Lau2DHistDP constructor : the upper error histogram has a different number of bins to the main histogram." << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+
+ xAxis = tempErrorLo->GetXaxis();
+ yAxis = tempErrorLo->GetYaxis();
+
+ if(static_cast<Double_t>(xAxis->GetXmin()) != minX ||
+ static_cast<Double_t>(xAxis->GetXmax()) != maxX) {
+ std::cerr << "ERROR in Lau2DHistDP constructor : the lower error histogram has a different x range to the main histogram." << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+
+ if(static_cast<Double_t>(yAxis->GetXmin()) != minY ||
+ static_cast<Double_t>(yAxis->GetXmax()) != maxY) {
+ std::cerr << "ERROR in Lau2DHistDP constructor : the lower error histogram has a different y range to the main histogram." << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+
+ xAxis = tempErrorHi->GetXaxis();
+ yAxis = tempErrorHi->GetYaxis();
+
+ if(static_cast<Double_t>(xAxis->GetXmin()) != minX ||
+ static_cast<Double_t>(xAxis->GetXmax()) != maxX) {
+ std::cerr << "ERROR in Lau2DHistDP constructor : the upper error histogram has a different x range to the main histogram." << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+
+ if(static_cast<Double_t>(yAxis->GetXmin()) != minY ||
+ static_cast<Double_t>(yAxis->GetXmax()) != maxY) {
+ std::cerr << "ERROR in Lau2DHistDP constructor : the upper error histogram has a different y range to the main histogram." << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+
+
+ if (fluctuateBins) {
+ this->doBinFluctuation(tempHist,tempErrorHi,tempErrorLo);
+ }
+ if (avEff > 0.0 && avEffError > 0.0) {
+ this->raiseOrLowerBins(tempHist,avEff,avEffError);
+ }
+
+ spline_ = new Lau2DCubicSpline(*tempHist);
+
+ delete tempHist;
+ delete tempErrorHi;
+ delete tempErrorLo;
+}
+
Lau2DSplineDP::~Lau2DSplineDP()
{
delete spline_;
spline_ = 0;
}
Double_t Lau2DSplineDP::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 Lau2DSplineDP::interpolateXY : Given position is outside the DP boundary, returning 0.0." << std::endl;
return 0.0;
}
return spline_->evaluate(x,y);
}
diff --git a/src/LauEffModel.cc b/src/LauEffModel.cc
index 464d101..c938ab1 100644
--- a/src/LauEffModel.cc
+++ b/src/LauEffModel.cc
@@ -1,191 +1,246 @@
// 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 LauEffModel.cc
\brief File containing implementation of LauEffModel class.
*/
#include <cstdlib>
#include <iostream>
using std::cout;
using std::cerr;
using std::endl;
#include "TSystem.h"
#include "Lau2DHistDP.hh"
#include "Lau2DSplineDP.hh"
#include "LauDaughters.hh"
#include "LauEffModel.hh"
#include "LauKinematics.hh"
#include "LauVetoes.hh"
ClassImp(LauEffModel)
LauEffModel::LauEffModel(const LauDaughters* daughters, const LauVetoes* vetoes) :
daughters_( daughters ),
vetoes_( vetoes ),
effHisto_( 0 ),
squareDP_( kFALSE ),
fluctuateEffHisto_( kFALSE ),
lowBinWarningIssued_( kFALSE ),
highBinWarningIssued_( kFALSE )
{
if ( daughters_ == 0 ) {
cerr << "ERROR in LauEffModel Constructor : invalid pointer to daughters object supplied." << endl;
gSystem->Exit(EXIT_FAILURE);
}
}
/*LauEffModel::LauEffModel( const LauEffModel& rhs ) :
daughters_( rhs.daughters_ ),
vetoes_( rhs.vetoes_ ),
effHisto_( rhs.effHisto_ ? new Lau2DHistDP( *rhs.effHisto_ ) : 0 ),
squareDP_( rhs.squareDP_ ),
fluctuateEffHisto_( rhs.fluctuateEffHisto_ )
{
}*/
LauEffModel::~LauEffModel()
{
if (effHisto_ != 0) {
delete effHisto_; effHisto_ = 0;
}
}
void LauEffModel::setEffHisto(const TH2* effHisto, Bool_t useInterpolation,
Bool_t fluctuateBins, Double_t avEff, Double_t absError,
Bool_t useUpperHalfOnly, Bool_t squareDP)
{
// Set the efficiency across the Dalitz plot using a predetermined 2D histogram
// with x = m_13^2, y = m_23^2.
Bool_t upperHalf( kFALSE );
if ( daughters_->gotSymmetricalDP() && useUpperHalfOnly == kTRUE) {upperHalf = kTRUE;}
cout<<"Efficiency histogram has upperHalf = "<<static_cast<Int_t>(upperHalf)<<endl;
squareDP_ = squareDP;
if(effHisto_) {
delete effHisto_;
effHisto_=0;
}
// Copy the histogram.
effHisto_ = new Lau2DHistDP(effHisto, daughters_,
useInterpolation, fluctuateBins,
avEff, absError, upperHalf, squareDP);
fluctuateEffHisto_ = fluctuateBins;
if (avEff > 0.0 && absError > 0.0) {
fluctuateEffHisto_ = kTRUE;
}
}
+void LauEffModel::setEffHisto(const TH2* effHisto, const TH2* errorHi, const TH2* errorLo, Bool_t useInterpolation,
+ Bool_t fluctuateBins, Double_t avEff, Double_t absError,
+ Bool_t useUpperHalfOnly, Bool_t squareDP)
+{
+ // Set the efficiency across the Dalitz plot using a predetermined 2D histogram
+ // with x = m_13^2, y = m_23^2.
+ Bool_t upperHalf( kFALSE );
+ if ( daughters_->gotSymmetricalDP() && useUpperHalfOnly == kTRUE) {upperHalf = kTRUE;}
+ cout<<"Efficiency histogram has upperHalf = "<<static_cast<Int_t>(upperHalf)<<endl;
+
+ squareDP_ = squareDP;
+
+ if(effHisto_) {
+ delete effHisto_;
+ effHisto_=0;
+ }
+
+ // Copy the histogram.
+ effHisto_ = new Lau2DHistDP(effHisto, errorHi, errorLo, daughters_,
+ useInterpolation, fluctuateBins,
+ avEff, absError, upperHalf, squareDP);
+ fluctuateEffHisto_ = fluctuateBins;
+
+ if (avEff > 0.0 && absError > 0.0) {
+ fluctuateEffHisto_ = kTRUE;
+ }
+}
+
void LauEffModel::setEffSpline(const TH2* effHisto,
Bool_t fluctuateBins, Double_t avEff, Double_t absError,
Bool_t useUpperHalfOnly, Bool_t squareDP)
{
// Set the efficiency across the Dalitz plot using a predetermined 2D histogram
// with x = m_13^2, y = m_23^2.
Bool_t upperHalf( kFALSE );
if ( daughters_->gotSymmetricalDP() && useUpperHalfOnly == kTRUE) {upperHalf = kTRUE;}
cout<<"Efficiency histogram has upperHalf = "<<static_cast<Int_t>(upperHalf)<<endl;
squareDP_ = squareDP;
if(effHisto_) {
delete effHisto_;
effHisto_=0;
}
// Copy the histogram.
effHisto_ = new Lau2DSplineDP(effHisto, daughters_,
fluctuateBins, avEff, absError, upperHalf, squareDP);
fluctuateEffHisto_ = fluctuateBins;
if (avEff > 0.0 && absError > 0.0) {
fluctuateEffHisto_ = kTRUE;
}
}
+void LauEffModel::setEffSpline(const TH2* effHisto, const TH2* errorHi, const TH2* errorLo,
+ Bool_t fluctuateBins, Double_t avEff, Double_t absError,
+ Bool_t useUpperHalfOnly, Bool_t squareDP)
+{
+ // Set the efficiency across the Dalitz plot using a predetermined 2D histogram
+ // with x = m_13^2, y = m_23^2.
+ Bool_t upperHalf( kFALSE );
+ if ( daughters_->gotSymmetricalDP() && useUpperHalfOnly == kTRUE) {upperHalf = kTRUE;}
+ cout<<"Efficiency histogram has upperHalf = "<<static_cast<Int_t>(upperHalf)<<endl;
+
+ squareDP_ = squareDP;
+
+ if(effHisto_) {
+ delete effHisto_;
+ effHisto_=0;
+ }
+
+ // Copy the histogram.
+ effHisto_ = new Lau2DSplineDP(effHisto, errorHi, errorLo, daughters_,
+ fluctuateBins, avEff, absError, upperHalf, squareDP);
+ fluctuateEffHisto_ = fluctuateBins;
+
+ if (avEff > 0.0 && absError > 0.0) {
+ fluctuateEffHisto_ = kTRUE;
+ }
+}
+
Double_t LauEffModel::getEffHistValue(Double_t xVal, Double_t yVal) const
{
// Get the efficiency from the 2D histo.
Double_t eff(1.0);
if (effHisto_ != 0) {
eff = effHisto_->interpolateXY(xVal, yVal);
}
return eff;
}
Double_t LauEffModel::calcEfficiency( const LauKinematics* kinematics ) const
{
// Routine to calculate the efficiency for the given event/point in
// the Dalitz plot. This routine uses the 2D histogram set by the
// setEffHisto() funciton.
Double_t eff(1.0);
// Check for vetoes
Bool_t vetoOK(kTRUE);
if ( vetoes_ != 0 ) {
vetoOK = vetoes_->passVeto(kinematics);
}
if (vetoOK == kFALSE) {
// We failed the veto.
eff = 0.0;
} else {
// We are using a 2D histogram representation of the efficiency across the Dalitz plot.
// First find out which bin we are in given the x and y Dalitz plot mass values
// Here, x = m13^2, y = m23^2 or if using square DP x = m', y = theta'.
if (squareDP_ == kTRUE) {
eff = this->getEffHistValue(kinematics->getmPrime(), kinematics->getThetaPrime());
} else {
eff = this->getEffHistValue(kinematics->getm13Sq(), kinematics->getm23Sq());
}
}
// Check that the efficiency is in the allowed range (0-1)
// If we're using a spline then out-of-range efficiencies can be caused by adjacent bins that all contain a value of either zero or one.
// The spline requires the efficiency, its first derivatives and the mixed second derivative to be continuous and to match the input histogram
// at the bin centres. Derivatives are calculated using a finite difference approximation taking the difference between the neighbouring bins.
// If two bins are zero but the third is not then the second bin will have a positive first derivative causing the spline to dip below zero
// between the two zero bins to remain smooth. The analogous case with adjacent maximised bins will cause peaks above one. Such dips are
// unavoidable but are correctly removed here.
if ( eff < 0.0 ) {
if(!lowBinWarningIssued_) {
std::cerr << "WARNING in LauEffModel::calcEfficiency : Efficiency " << eff << " is less than 0 - setting to 0. You may want to check your histogram!" << std::endl
<< " : If you are using a spline then this could be caused by adjacent empty bins. Further warnings will be suppressed." << std::endl;
lowBinWarningIssued_=kTRUE;
}
eff = 0.0;
} else if ( eff > 1.0 ) {
if(!highBinWarningIssued_) {
std::cerr << "WARNING in LauEffModel::calcEfficiency : Efficiency " << eff << " is greater than 1 - setting to 1. You may want to check your histogram!" << std::endl
<< " : If you are using a spline then this could be caused by adjacent full bins. Further warnings will be suppressed." << std::endl;
highBinWarningIssued_=kTRUE;
}
eff = 1.0;
}
return eff;
}
Bool_t LauEffModel::passVeto( const LauKinematics* kinematics ) const
{
Bool_t pass = kTRUE;
if ( vetoes_ != 0 ) {
pass = vetoes_->passVeto( kinematics );
}
return pass;
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 8:34 PM (1 d, 3 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3806063
Default Alt Text
(62 KB)

Event Timeline