Page MenuHomeHEPForge
Files F9501496

No OneTemporary
Actions

View Options

diff --git a/doc/release.notes b/doc/release.notes
index 3edbdb3..d89e171 100644
--- a/doc/release.notes
+++ b/doc/release.notes
@@ -1,198 +1,202 @@
///////////////////////////////////////////////////////////////
/// ///
/// This is the History file for the Laura++ package. ///
/// ///
///////////////////////////////////////////////////////////////
+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/Lau1DHistPdf.hh b/inc/Lau1DHistPdf.hh
index 92571bd..84995b5 100644
--- a/inc/Lau1DHistPdf.hh
+++ b/inc/Lau1DHistPdf.hh
@@ -1,125 +1,126 @@
// 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 Lau1DHistPdf.hh
\brief File containing declaration of Lau1DHistPdf class.
*/
/*! \class Lau1DHistPdf
\brief Class for defining a 1D histogram PDF.
Class for defining a 1D histogram PDF.
Employs linear interpolation to get the PDF value based on how far away a
point is from nearby bin centres. The returned values are normalised to
the total area.
*/
#ifndef LAU_1DHIST_PDF
#define LAU_1DHIST_PDF
#include "LauAbsPdf.hh"
class TH1;
class Lau1DHistPdf : public LauAbsPdf {
public:
//! Constructor
/*!
\param [in] theVarName the name of the abscissa variable
\param [in] hist the 1D histogram from which the PDF should be constructed
\param [in] minAbscissa the minimum value of the abscissa
\param [in] maxAbscissa the maximum value of the abscissa
\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] 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.
*/
Lau1DHistPdf(const TString& theVarName, const TH1* hist, Double_t minAbscissa, Double_t maxAbscissa,
Bool_t useInterpolation = kTRUE, Bool_t fluctuateBins = kFALSE);
//! Destructor
virtual ~Lau1DHistPdf();
//! Copy constructor
Lau1DHistPdf(const Lau1DHistPdf& other);
//! Calculate the likelihood (and intermediate info) for a given abscissa
/*!
\param [in] abscissas the values of the abscissa(s)
*/
virtual void calcLikelihoodInfo(const LauAbscissas& abscissas);
using LauAbsPdf::calcLikelihoodInfo;
//! Check that PDF is positive
/*!
Check that PDF is positive.
Dealt with in getBinHistValue, so nothing to do here.
*/
virtual void checkPositiveness() {};
//! Calculate the normalisation
virtual void calcNorm();
//! Calculate the PDF height
/*!
\param [in] kinematics the current DP kinematics (not strictly required in this case since PDF has no DP dependence)
*/
virtual void calcPDFHeight( const LauKinematics* kinematics );
protected:
//! Fluctuate the histogram bin contents in accorance with their errors
void doBinFluctuation();
//! Check the normalisation calculation
void checkNormalisation();
//! Get the bin content from the histogram
/*!
\param [in] bin the bin number
\return the bin content
*/
Double_t getBinHistValue(Int_t bin) const;
//! Perform the interpolation (unnormalised)
/*!
\param [in] x the abscissa value
\return the unnormalised PDF value
*/
Double_t interpolate(Double_t x) const;
//! Perform the interpolation and divide by the normalisation
/*!
\param [in] x the abscissa value
\return the normalised PDF value
*/
Double_t interpolateNorm(Double_t x) const;
private:
//! The underlying histogram
TH1* hist_;
//! Control boolean for using the linear interpolation
Bool_t useInterpolation_;
//! Control boolean for performing the fluctuation of the histogram bin contents
Bool_t fluctuateBins_;
//! The number of bins in the histogram
Int_t nBins_;
//! The histogram axis minimum
Double_t axisMin_;
//! The histogram axis maximum
Double_t axisMax_;
//! The histogram axis range
Double_t axisRange_;
ClassDef(Lau1DHistPdf,0) // 1D histogram pdf class
};
#endif
diff --git a/inc/Lau2DHistDP.hh b/inc/Lau2DHistDP.hh
index 64134cd..8bb82fc 100644
--- a/inc/Lau2DHistDP.hh
+++ b/inc/Lau2DHistDP.hh
@@ -1,111 +1,113 @@
// 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)
- \param [in] avEff the desired average efficiency - see Lau2DHistDP::raiseOrLowerBins
+ \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);
//! 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 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/Lau2DHistDPPdf.hh b/inc/Lau2DHistDPPdf.hh
index 4d20b7e..adf1109 100644
--- a/inc/Lau2DHistDPPdf.hh
+++ b/inc/Lau2DHistDPPdf.hh
@@ -1,139 +1,140 @@
// 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 "Lau2DAbsHistDPPdf.hh"
class TH2;
class LauKinematics;
class LauVetoes;
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] 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] 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 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 the PDF normalisation
void calcHistNorm();
//! Check the normalisation calculation
void checkNormalisation();
private:
//! Copy constructor - not implemented
Lau2DHistDPPdf( const Lau2DHistDPPdf& rhs );
//! Copy assignment operator - not implemented
Lau2DHistDPPdf& operator=(const Lau2DHistDPPdf& rhs);
//! The underlying histogram
TH2* hist_;
//! 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_;
ClassDef(Lau2DHistDPPdf,0) // 2D Histogram utility class for DP analyses
};
#endif
diff --git a/inc/Lau2DHistPdf.hh b/inc/Lau2DHistPdf.hh
index 259793c..7610970 100644
--- a/inc/Lau2DHistPdf.hh
+++ b/inc/Lau2DHistPdf.hh
@@ -1,207 +1,208 @@
// Copyright University of Warwick 2008 - 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 Lau2DHistPdf.hh
\brief File containing declaration of Lau2DHistPdf class.
*/
/*! \class Lau2DHistPdf
\brief Class for defining a 2D histogram PDF.
Class for defining a 2D histogram PDF.
Employs linear interpolation to get the histogram value based on how far away a
point is to nearby bin centres. The returned values are normalised to
the total area.
*/
#ifndef LAU_2DHIST_PDF
#define LAU_2DHIST_PDF
#include "LauAbsPdf.hh"
class TH2;
class TH1;
class Lau1DHistPdf;
class Lau2DHistPdf : public LauAbsPdf {
public:
//! Constructor
/*!
\param [in] theVarNames the names of the abscissa variables
\param [in] hist the 2D histogram from which the PDF should be constructed
\param [in] minVals the minimum values of the abscissas
\param [in] maxVals the maximum values of the abscissas
\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] 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.
*/
Lau2DHistPdf(const std::vector<TString>& theVarNames, const TH2* hist,
const LauFitData& minVals, const LauFitData& maxVals,
Bool_t useInterpolation = kTRUE, Bool_t fluctuateBins = kFALSE);
//! Destructor
virtual ~Lau2DHistPdf();
//! Copy constructor
Lau2DHistPdf(const Lau2DHistPdf& other);
//! Get the minimum value of x-axis abscissa
Double_t getMinX() const {return minX_;}
//! Get the maximum value of x-axis abscissa
Double_t getMaxX() const {return maxX_;}
//! Get the range of x-axis abscissa values
Double_t getRangeX() const {return rangeX_;}
//! Get the minimum value of y-axis abscissa
Double_t getMinY() const {return minY_;}
//! Get the maximum value of y-axis abscissa
Double_t getMaxY() const {return maxY_;}
//! Get the range of y-axis abscissa values
Double_t getRangeY() const {return rangeY_;}
//! Cache information from data
/*!
\param [in] inputData the input data
*/
virtual void cacheInfo(const LauFitDataTree& inputData);
//! Calculate the likelihood (and intermediate info) for a given value of the abscissas
/*!
\param [in] abscissas the abscissa values
*/
virtual void calcLikelihoodInfo(const LauAbscissas& abscissas);
//! Calculate the likelihood (and intermediate info) for a given event number
/*!
\param [in] iEvt event number
*/
virtual void calcLikelihoodInfo(UInt_t iEvt);
//! Get likelihood only for the given variable
/*!
\param [in] theVarName the variable name (must be one of the two dimensions of this histogram!)
\return likelihood value
*/
virtual Double_t getLikelihood( const TString& theVarName ) const;
using LauAbsPdf::getLikelihood;
//! Check that PDF is positive
/*!
Check that PDF is positive.
Dealt with in getBinHistValue, so nothing to do here.
*/
virtual void checkPositiveness() {}; // Nothing to check here.
//! Calculate the normalisation
virtual void calcNorm();
//! Calculate the PDF height
/*!
\param [in] kinematics the current DP kinematics (not strictly required in this case since PDF has no DP dependence)
*/
virtual void calcPDFHeight( const LauKinematics* kinematics );
//! Generate an event from the PDF
/*!
\param [in] kinematics the current DP kinematics (used to determine the DP position by some PDFs that have dependence on it)
\return generated data
*/
virtual LauFitData generate(const LauKinematics* kinematics);
protected:
//! Fluctuate the histogram bin contents in accordance with their errors
void doBinFluctuation();
//! Check the normalisation calculation
void checkNormalisation();
//! Get the bin content from the histogram
/*!
\param [in] i the x-axis bin number
\param [in] j the y-axis bin number
\return the bin conent
*/
Double_t getBinHistValue(Int_t i, Int_t j) const;
//! Perform the interpolation (unnormalised)
/*!
\param [in] x the x-axis abscissa value
\param [in] y the y-aixs abscissa 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;
private:
//! The underlying histogram
TH2* hist_;
//! Projection of histogram x-axis
TH1* xProj_;
//! Projection of histogram y-axis
TH1* yProj_;
//! 1D PDF for x variable
Lau1DHistPdf* xVarPdf_;
//! 1D PDF for y variable
Lau1DHistPdf* yVarPdf_;
//! x variable name
TString xName_;
//! y variable name
TString yName_;
//! 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 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_;
//! Control boolean for using the linear interpolation
Bool_t useInterpolation_;
//! Control boolean for performing the fluctuation of the histogram bin contents
Bool_t fluctuateBins_;
ClassDef(Lau2DHistPdf,0) // 2D histogram pdf class
};
#endif
diff --git a/inc/Lau2DSplineDP.hh b/inc/Lau2DSplineDP.hh
index 9ecaba7..2e52e68 100644
--- a/inc/Lau2DSplineDP.hh
+++ b/inc/Lau2DSplineDP.hh
@@ -1,76 +1,78 @@
// 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)
- \param [in] avEff the desired average efficiency - see Lau2DSplineDP::raiseOrLowerBins
+ \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);
//! 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/Lau2DSplineDPPdf.hh b/inc/Lau2DSplineDPPdf.hh
index 1595154..fd51eb1 100644
--- a/inc/Lau2DSplineDPPdf.hh
+++ b/inc/Lau2DSplineDPPdf.hh
@@ -1,91 +1,92 @@
// 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] 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] 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 d1af760..169a568 100644
--- a/inc/LauBkgndDPModel.hh
+++ b/inc/LauBkgndDPModel.hh
@@ -1,150 +1,152 @@
// 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 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] 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] 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] 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] 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
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_;
//! Flag to track whether a warning has been issued for bin values less than zero
mutable Bool_t lowBinWarningIssued_;
ClassDef(LauBkgndDPModel,0) // DP background model
};
#endif
diff --git a/inc/LauEffModel.hh b/inc/LauEffModel.hh
index decf7f2..8ef388d 100644
--- a/inc/LauEffModel.hh
+++ b/inc/LauEffModel.hh
@@ -1,137 +1,141 @@
// 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.
- \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 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 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.
- \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 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);
//! 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/Lau1DHistPdf.cc b/src/Lau1DHistPdf.cc
index 36ca29a..7a1c8d2 100644
--- a/src/Lau1DHistPdf.cc
+++ b/src/Lau1DHistPdf.cc
@@ -1,260 +1,260 @@
// 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 Lau1DHistPdf.cc
\brief File containing implementation of Lau1DHistPdf class.
*/
#include <cstdlib>
#include <iostream>
#include <vector>
#include "TAxis.h"
#include "TH1.h"
#include "TRandom.h"
#include "TSystem.h"
#include "Lau1DHistPdf.hh"
#include "LauRandom.hh"
class LauParameter;
ClassImp(Lau1DHistPdf)
Lau1DHistPdf::Lau1DHistPdf(const TString& theVarName, const TH1* hist, Double_t minAbscissa, Double_t maxAbscissa,
Bool_t useInterpolation, Bool_t fluctuateBins) :
LauAbsPdf(theVarName, std::vector<LauParameter*>(), minAbscissa, maxAbscissa),
hist_(hist ? dynamic_cast<TH1*>(hist->Clone()) : 0),
useInterpolation_(useInterpolation),
fluctuateBins_(fluctuateBins),
nBins_(0),
axisMin_(0.0),
axisMax_(0.0),
axisRange_(0.0)
{
// Constructor
// Set the directory for the histogram
hist_->SetDirectory(0);
// Save various attributes of the histogram
nBins_ = hist_->GetNbinsX();
TAxis* xAxis = hist_->GetXaxis();
axisMin_ = xAxis->GetXmin();
axisMax_ = xAxis->GetXmax();
axisRange_ = axisMax_ - axisMin_;
// Check that axis range corresponds to range of abscissa
if (TMath::Abs(this->getMinAbscissa() - axisMin_)>1e-6) {
std::cerr << "ERROR in Lau1DHistPdf::Lau1DHistPdf : Histogram axis minimum: " << axisMin_ <<
" does not correspond to abscissa minimum: " << this->getMinAbscissa() << "." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
if (TMath::Abs(this->getMaxAbscissa() - axisMax_)>1e-6) {
std::cerr << "ERROR in Lau1DHistPdf::Lau1DHistPdf : Histogram axis maximum: " << axisMax_ <<
" does not correspond to abscissa maximum: " << this->getMaxAbscissa() << "." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
// If the bins are to be fluctuated then do so now before
// calculating anything that depends on the bin content.
if (fluctuateBins) {
this->doBinFluctuation();
}
// Calculate the PDF normalisation.
this->calcNorm();
// And check it is OK.
this->checkNormalisation();
}
Lau1DHistPdf::~Lau1DHistPdf()
{
// Destructor
delete hist_; hist_ = 0;
}
Lau1DHistPdf::Lau1DHistPdf(const Lau1DHistPdf& other) : LauAbsPdf(other.varName(), other.getParameters(), other.getMinAbscissa(), other.getMaxAbscissa())
{
hist_ = other.hist_ ? dynamic_cast<TH1*>(other.hist_->Clone()) : 0;
useInterpolation_ = other.useInterpolation_;
fluctuateBins_ = other.fluctuateBins_;
this->setRandomFun(other.getRandomFun());
this->calcNorm();
}
void Lau1DHistPdf::calcPDFHeight( const LauKinematics* /*kinematics*/ )
{
if (this->heightUpToDate()) {
return;
}
// Get the maximum height of the histogram
Int_t maxBin = hist_->GetMaximumBin();
Double_t height = hist_->GetBinContent(maxBin);
this->setMaxHeight(height);
}
void Lau1DHistPdf::calcNorm()
{
// Calculate the histogram normalisation.
// Loop over the range to get the total area.
// Just sum the contributions up using 1e-3 increments of the range.
// Multiply the end result by dx.
Double_t dx(1e-3*axisRange_);
Double_t area(0.0);
Double_t x(axisMin_ + dx/2.0);
while (x > axisMin_ && x < axisMax_) {
area += this->interpolate(x);
x += dx;
}
Double_t norm = area*dx;
this->setNorm(norm);
}
void Lau1DHistPdf::checkNormalisation()
{
Double_t dx(1e-3*axisRange_);
Double_t area(0.0);
Double_t areaNoNorm(0.0);
Double_t x(axisMin_ + dx/2.0);
while (x > axisMin_ && x < axisMax_) {
area += this->interpolateNorm(x);
useInterpolation_ = kFALSE;
areaNoNorm += this->interpolate(x);
useInterpolation_ = kTRUE;
x += dx;
}
Double_t norm = area*dx;
std::cout << "INFO in Lau1DHistPdf::checkNormalisation : Area = " << area << ", dx = " << dx << std::endl;
std::cout << " : Area with no norm = " << areaNoNorm << "*dx = " << areaNoNorm*dx << std::endl;
std::cout << " : The total area of the normalised histogram PDF is " << norm << std::endl;
}
Double_t Lau1DHistPdf::getBinHistValue(Int_t bin) const
{
// Check that bin is in range [1 , nBins_]
if ((bin < 1) || (bin > nBins_)) {
return 0.0;
}
Double_t value = static_cast<Double_t>(hist_->GetBinContent(bin));
// protect against negative values
if ( value < 0.0 ) {
std::cerr << "WARNING in Lau1DHistPdf::getBinHistValue : Negative bin content set to zero!" << std::endl;
value = 0.0;
}
return value;
}
Double_t Lau1DHistPdf::interpolateNorm(Double_t x) const
{
// Get the normalised interpolated value.
Double_t value = this->interpolate(x);
Double_t norm = this->getNorm();
return value/norm;
}
void Lau1DHistPdf::calcLikelihoodInfo(const LauAbscissas& abscissas)
{
// Check that the given abscissa is within the allowed range
if (!this->checkRange(abscissas)) {
gSystem->Exit(EXIT_FAILURE);
}
// Get our abscissa
Double_t abscissa = abscissas[0];
// Calculate the interpolated value
Double_t value = this->interpolate(abscissa);
this->setUnNormPDFVal(value);
}
Double_t Lau1DHistPdf::interpolate(Double_t x) const
{
// This function returns the interpolated value of the histogram function
// for the given value of x by finding the adjacent bins and extrapolating
// using weights based on the inverse distance of the point from the adajcent
// bin centres.
// Find the histogram bin
Int_t bin = hist_->FindFixBin(x);
// Ask whether we want to do the interpolation, since this method is
// not reliable for low statistics histograms.
if (useInterpolation_ == kFALSE) {
return this->getBinHistValue(bin);
}
// Find the bin centres (actual co-ordinate positions, not histogram indices)
Double_t cbinx = hist_->GetBinCenter(bin);
// Find the adjacent bins
Double_t deltax = x - cbinx;
Int_t bin_adj(0);
if (deltax > 0.0) {
bin_adj = bin + 1;
} else {
bin_adj = bin - 1;
}
Bool_t isBoundary(kFALSE);
if ( bin_adj > nBins_ || bin_adj < 1 ) {
isBoundary = kTRUE;
}
// At the edges, do no interpolation, use entry in bin.
if (isBoundary == kTRUE) {
return this->getBinHistValue(bin);
}
// Linear interpolation using inverse distance as weights.
// Find the adjacent bin centre
Double_t cbinx_adj = hist_->GetBinCenter(bin_adj);
Double_t deltax_adj = cbinx_adj - x;
Double_t dx0 = TMath::Abs(deltax);
Double_t dx1 = TMath::Abs(deltax_adj);
Double_t denom = dx0 + dx1;
Double_t value0 = this->getBinHistValue(bin);
Double_t value1 = this->getBinHistValue(bin_adj);
Double_t value = value0*dx1 + value1*dx0;
value /= denom;
return value;
}
void Lau1DHistPdf::doBinFluctuation()
{
- TRandom* random = LauRandom::zeroSeedRandom();
+ TRandom* random = LauRandom::randomFun();
for (Int_t bin(0); bin<nBins_; bin++) {
Double_t currentContent = hist_->GetBinContent(bin+1);
Double_t currentError = hist_->GetBinError(bin+1);
Double_t newContent = random->Gaus(currentContent,currentError);
if (newContent<0.0) {
hist_->SetBinContent(bin+1,0.0);
} else {
hist_->SetBinContent(bin+1,newContent);
}
}
}
diff --git a/src/Lau2DAbsHistDP.cc b/src/Lau2DAbsHistDP.cc
index a379fd1..93226a8 100644
--- a/src/Lau2DAbsHistDP.cc
+++ b/src/Lau2DAbsHistDP.cc
@@ -1,123 +1,123 @@
// 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::zeroSeedRandom();
+ 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::raiseOrLowerBins(TH2* hist, Double_t avEff, Double_t avEffError)
{
- TRandom* random = LauRandom::zeroSeedRandom();
+ 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/Lau2DAbsHistDPPdf.cc b/src/Lau2DAbsHistDPPdf.cc
index ac1f4e0..3ed0c8e 100644
--- a/src/Lau2DAbsHistDPPdf.cc
+++ b/src/Lau2DAbsHistDPPdf.cc
@@ -1,103 +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();
+ 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);
}
}
}
}
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/Lau2DHistPdf.cc b/src/Lau2DHistPdf.cc
index 92f1ea5..a68b095 100644
--- a/src/Lau2DHistPdf.cc
+++ b/src/Lau2DHistPdf.cc
@@ -1,501 +1,501 @@
// Copyright University of Warwick 2008 - 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 Lau2DHistPdf.cc
\brief File containing implementation of Lau2DHistPdf class.
*/
#include <cstdlib>
#include <iostream>
#include <vector>
#include "TAxis.h"
#include "TH2.h"
#include "TRandom.h"
#include "TSystem.h"
#include "Lau1DHistPdf.hh"
#include "Lau2DHistPdf.hh"
#include "LauRandom.hh"
class LauParameter;
ClassImp(Lau2DHistPdf)
Lau2DHistPdf::Lau2DHistPdf(const std::vector<TString>& theVarNames, const TH2* hist,
const LauFitData& minVals, const LauFitData& maxVals,
Bool_t useInterpolation, Bool_t fluctuateBins) :
LauAbsPdf(theVarNames, std::vector<LauParameter*>(), minVals, maxVals),
hist_(hist ? dynamic_cast<TH2*>(hist->Clone()) : 0),
xProj_(0),
yProj_(0),
xVarPdf_(0),
yVarPdf_(0),
xName_(""),
yName_(""),
nBinsX_(0),
nBinsY_(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),
invBinXWidth_(0.0),
invBinYWidth_(0.0),
useInterpolation_(useInterpolation),
fluctuateBins_(fluctuateBins)
{
// Constructor
// Check that we have two variables
if ( this->nInputVars() != 2 ) {
std::cerr << "ERROR in Lau2DHistPdf::Lau2DHistPdf : Have not been provided with exactly two variables." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
// Set the variable names from the abstract class
xName_ = this->varNames()[0];
yName_ = this->varNames()[1];
// Set the variable limits from the abstract class
minX_ = this->getMinAbscissa(xName_);
maxX_ = this->getMaxAbscissa(xName_);
minY_ = this->getMinAbscissa(yName_);
maxY_ = this->getMaxAbscissa(yName_);
rangeX_ = maxX_ - minX_;
rangeY_ = maxY_ - minY_;
// Have we got a valid histogram
if ( hist_ == 0 ) {
std::cerr << "ERROR in Lau2DHistPdf::Lau2DHistPdf : Histogram pointer is null." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
// Set the directory for the histogram
hist_->SetDirectory(0);
// Save various attributes of the histogram
nBinsX_ = hist_->GetNbinsX();
nBinsY_ = hist_->GetNbinsY();
TAxis* xAxis = hist_->GetXaxis();
Double_t xAxisMin = xAxis->GetXmin();
Double_t xAxisMax = xAxis->GetXmax();
TAxis* yAxis = hist_->GetYaxis();
Double_t yAxisMin = yAxis->GetXmin();
Double_t yAxisMax = yAxis->GetXmax();
// Check that axis ranges corresponds to ranges of abscissas
if (TMath::Abs(minX_ - xAxisMin)>1e-6) {
std::cerr << "ERROR in Lau2DHistPdf::Lau2DHistPdf : Histogram x-axis minimum: " << xAxisMin <<
" does not correspond to " << xName_ << " minimum: " << minX_ << "." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
if (TMath::Abs(maxX_ - xAxisMax)>1e-6) {
std::cerr << "ERROR in Lau2DHistPdf::Lau2DHistPdf : Histogram x-axis maximum: " << xAxisMax <<
" does not correspond to " << xName_ << " maximum: " << maxX_ << "." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
if (TMath::Abs(minY_ - yAxisMin)>1e-6) {
std::cerr << "ERROR in Lau2DHistPdf::Lau2DHistPdf : Histogram y-axis minimum: " << yAxisMin <<
" does not correspond to " << yName_ << " minimum: " << minY_ << "." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
if (TMath::Abs(maxY_ - yAxisMax)>1e-6) {
std::cerr << "ERROR in Lau2DHistPdf::Lau2DHistPdf : Histogram y-axis maximum: " << yAxisMax <<
" does not correspond to " << yName_ << " maximum: " << maxY_ << "." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
// Calculate the bin widths and inverse bin widths
binXWidth_ = static_cast<Double_t>(TMath::Abs(rangeX_)/(nBinsX_*1.0));
binYWidth_ = static_cast<Double_t>(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();
}
// Create the projections and 1D PDFs from those
xProj_ = hist_->ProjectionX();
yProj_ = hist_->ProjectionY();
if ( !xProj_ || !yProj_ ) {
std::cerr << "ERROR in Lau2DHistPdf::Lau2DHistPdf : Can't get X or Y projection from 2D histogram." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
xVarPdf_ = new Lau1DHistPdf(xName_, xProj_, minX_, maxX_, useInterpolation_, fluctuateBins_);
yVarPdf_ = new Lau1DHistPdf(yName_, yProj_, minY_, maxY_, useInterpolation_, fluctuateBins_);
// Calculate the PDF normalisation.
this->calcNorm();
// And check it is OK.
this->checkNormalisation();
}
Lau2DHistPdf::~Lau2DHistPdf()
{
// Destructor
delete hist_; hist_ = 0;
delete xProj_; xProj_ = 0;
delete yProj_; yProj_ = 0;
delete xVarPdf_; xVarPdf_ = 0;
delete yVarPdf_; yVarPdf_ = 0;
}
Lau2DHistPdf::Lau2DHistPdf(const Lau2DHistPdf& other) : LauAbsPdf(other.varName(), other.getParameters(), other.getMinAbscissa(), other.getMaxAbscissa())
{
hist_ = dynamic_cast<TH2*>(other.hist_->Clone());
xName_ = other.xName_;
yName_ = other.yName_;
nBinsX_ = other.nBinsX_;
nBinsY_ = other.nBinsY_;
minX_ = other.minX_;
maxX_ = other.maxX_;
minY_ = other.minY_;
maxY_ = other.maxY_;
rangeX_ = other.rangeX_;
rangeY_ = other.rangeY_;
binXWidth_ = other.binXWidth_;
binYWidth_ = other.binYWidth_;
invBinXWidth_ = other.invBinXWidth_;
invBinYWidth_ = other.invBinYWidth_;
useInterpolation_ = other.useInterpolation_;
fluctuateBins_ = other.fluctuateBins_;
xProj_ = hist_->ProjectionX();
yProj_ = hist_->ProjectionY();
xVarPdf_ = new Lau1DHistPdf(xName_, xProj_, minX_, maxX_, useInterpolation_, fluctuateBins_);
yVarPdf_ = new Lau1DHistPdf(yName_, yProj_, minY_, maxY_, useInterpolation_, fluctuateBins_);
this->setRandomFun(other.getRandomFun());
this->calcNorm();
}
void Lau2DHistPdf::calcPDFHeight( const LauKinematics* /*kinematics*/ )
{
if (this->heightUpToDate()) {
return;
}
// Get the maximum height of the histogram
Int_t maxBin = hist_->GetMaximumBin();
Double_t height = hist_->GetBinContent(maxBin);
this->setMaxHeight(height);
}
void Lau2DHistPdf::calcNorm()
{
// Calculate the histogram normalisation.
// 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.
Double_t dx(1e-3*rangeX_), dy(1e-3*rangeY_);
Double_t area(0.0);
Double_t x(minX_ + dx/2.0);
while (x > minX_ && x < maxX_) {
Double_t y(minY_ + dy/2.0);
while (y > minY_ && y < maxY_) {
area += this->interpolateXY(x,y);
y += dy;
} // y while loop
x += dx;
} // x while loop
Double_t norm = area*dx*dy;
this->setNorm(norm);
}
void Lau2DHistPdf::checkNormalisation()
{
Double_t dx(1e-3*rangeX_), dy(1e-3*rangeY_);
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(minX_ + dx/2.0);
while (x > minX_ && x < maxX_) {
Double_t y(minY_ + dy/2.0);
while (y > minY_ && y < maxY_) {
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 Lau2DHistPdf::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;
}
Double_t Lau2DHistPdf::getBinHistValue(Int_t i, Int_t j) const
{
// Check that x co-ord is in range [1 , nBinsX_]
if ((i < 1) || (i > nBinsX_)) {
return 0.0;
}
// Check that y co-ord is in range [1 , nBinsY_]
if ((j < 1) || (j > nBinsY_)) {
return 0.0;
}
Double_t value = hist_->GetBinContent(i,j);
// protect against negative values
if ( value < 0.0 ) {
std::cerr << "WARNING in Lau2DHistPdf::getBinHistValue : Negative bin content set to zero!" << std::endl;
value = 0.0;
}
return value;
}
Double_t Lau2DHistPdf::interpolateXYNorm(Double_t x, Double_t y) const
{
// Get the normalised interpolated value.
Double_t value = this->interpolateXY(x,y);
Double_t norm = this->getNorm();
return value/norm;
}
Double_t Lau2DHistPdf::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.
// Find the 2D histogram bin for x and y
Int_t i = static_cast<Int_t>((x - minX_)*invBinXWidth_)+1;
Int_t j = static_cast<Int_t>((y - minY_)*invBinYWidth_)+1;
if (i > nBinsX_) {i = nBinsX_;}
if (j > nBinsY_) {j = nBinsY_;}
// 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 = static_cast<Double_t>(i-0.5)*binXWidth_ + minX_;
Double_t cbiny = static_cast<Double_t>(j-0.5)*binYWidth_ + minY_;
// 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 < 1) {isYBoundary = kTRUE;}
if (j_adj > nBinsY_ || j_adj < 1) {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 = static_cast<Double_t>(i_adj-0.5)*binXWidth_ + minX_;
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 = static_cast<Double_t>(j_adj-0.5)*binYWidth_ + minY_;
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 = static_cast<Double_t>(i_adj-0.5)*binXWidth_ + minX_;
Double_t cbiny_adj = static_cast<Double_t>(j_adj-0.5)*binYWidth_ + minY_;
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 Lau2DHistPdf::calcLikelihoodInfo(const LauAbscissas& abscissas)
{
// Check that the given abscissa is within the allowed range
if (!this->checkRange(abscissas)) {
gSystem->Exit(EXIT_FAILURE);
}
// Calculate the interpolated value
Double_t x = abscissas[0];
Double_t y = abscissas[1];
Double_t value = this->interpolateXY(x,y);
this->setUnNormPDFVal(value);
// TODO - do we need this? I think probably not.
// Have the 1D PDFs calculate their likelihoods as well
//xVarPdf_->calcLikelihoodInfo(x);
//yVarPdf_->calcLikelihoodInfo(y);
}
void Lau2DHistPdf::calcLikelihoodInfo(UInt_t iEvt)
{
// Call the base class method
this->LauAbsPdf::calcLikelihoodInfo(iEvt);
// Have the 1D PDFs retrieve their values as well
xVarPdf_->calcLikelihoodInfo(iEvt);
yVarPdf_->calcLikelihoodInfo(iEvt);
}
Double_t Lau2DHistPdf::getLikelihood( const TString& theVarName ) const
{
if ( theVarName == xName_ ) {
return xVarPdf_->getLikelihood();
} else if ( theVarName == yName_ ) {
return yVarPdf_->getLikelihood();
} else {
std::cerr << "ERROR in Lau2DHistPdf::getLikelihood : Unrecognised variable name \"" << theVarName << "\", cannot determine likelihood." << std::endl;
return 0.0;
}
}
void Lau2DHistPdf::cacheInfo(const LauFitDataTree& inputData)
{
// Call the base class method
this->LauAbsPdf::cacheInfo(inputData);
// Have the 1D PDFs cache their info as well
xVarPdf_->cacheInfo(inputData);
yVarPdf_->cacheInfo(inputData);
}
void Lau2DHistPdf::doBinFluctuation()
{
- TRandom* random(LauRandom::zeroSeedRandom());
+ TRandom* random(LauRandom::randomFun());
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);
}
}
}
}
LauFitData Lau2DHistPdf::generate(const LauKinematics* kinematics)
{
this->withinGeneration(kTRUE);
// Check that the PDF height is up to date
this->calcPDFHeight( kinematics );
// Generate the value of the abscissa.
Bool_t gotAbscissa(kFALSE);
if (this->getRandomFun() == 0) {
std::cerr << "ERROR in Lau2DHistPdf::generate : Please set the random number generator for this PDF by using the setRandomFun(TRandom*) function." << std::endl;
this->withinGeneration(kFALSE);
gSystem->Exit(EXIT_FAILURE);
}
Double_t genPDFVal(0.0);
LauAbscissas genAbscissas(2);
Double_t PDFheight = this->getMaxHeight()*(1.0+1e-11);
while (!gotAbscissa) {
genAbscissas[0] = this->getRandomFun()->Rndm()*this->getRange(xName_) + this->getMinAbscissa(xName_);
genAbscissas[1] = this->getRandomFun()->Rndm()*this->getRange(yName_) + this->getMinAbscissa(yName_);
this->calcLikelihoodInfo(genAbscissas);
genPDFVal = this->getUnNormLikelihood();
if (this->getRandomFun()->Rndm() <= genPDFVal/PDFheight) {gotAbscissa = kTRUE;}
if (genPDFVal > PDFheight) {
std::cerr << "WARNING in LauAbsPdf::generate : genPDFVal = " << genPDFVal << " is larger than the specified PDF height " << this->getMaxHeight() << " for (x,y) = (" << genAbscissas[0] << "," << genAbscissas[1] << ")." << std::endl;
std::cerr << " : Need to reset height to be larger than " << genPDFVal << " by using the setMaxHeight(Double_t) function and re-run the Monte Carlo generation!" << std::endl;
}
}
LauFitData genData;
genData[ xName_ ] = genAbscissas[0];
genData[ yName_ ] = genAbscissas[1];
this->withinGeneration(kFALSE);
return genData;
}

File Metadata

Mime Type
text/x-diff
Expires
Sun, Feb 23, 2:32 PM (15 h, 42 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4486636
Default Alt Text
(75 KB)

Event Timeline