Page MenuHomeHEPForge
Files F7877267

No OneTemporary
Actions

View Options

diff --git a/Makefile b/Makefile
index d42bac3..d053b91 100644
--- a/Makefile
+++ b/Makefile
@@ -1,147 +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 #
# #
# ------------------------------- #
# Standalone Makefile for Laura++ #
# ------------------------------- #
# #
# Instructions #
# - Review 'external configuration' section below #
# to match systems compilers setup #
# #
# - Make sure the ROOTSYS environment variable is set and points #
# to your ROOT release or the root-config script is in your PATH #
# #
# - run 'make <target>' #
# #
# Build targets #
# lib - make libLaura++.a #
# shlib - make libLaura++.so (default) #
# clean - delete all intermediate and final build objects #
# #
##########################################################################
# --- External configuration ----------------------------------
# first check that ROOTSYS is defined
ifndef ROOTSYS
ROOTSYS := $(shell root-config --prefix)
ROOTBINDIR := $(shell root-config --bindir)
ifeq ($(ROOTSYS), )
$(error running of root-config failed or reported null value)
endif
else
ROOTBINDIR := $(ROOTSYS)/bin
endif
ROOTCONFIG := $(ROOTBINDIR)/root-config
ARCH := $(shell $(ROOTCONFIG) --arch)
PLATFORM := $(shell $(ROOTCONFIG) --platform)
INCLUDES =
SRCDIR = src
INCDIR = inc
LIBDIR = lib
WORKDIR = tmp
ifeq ($(findstring linux, $(ARCH)),linux)
# This set here should work for Linux.
CXX = g++
LD = g++
CXXFLAGS = -g -O2 -Wall -Wextra -Wshadow -Woverloaded-virtual -fPIC
MFLAGS = -MM
SOFLAGS = -shared
endif
ifeq ($(ARCH),macosx64)
# For Mac OS X you may need to put -m64 in CXXFLAGS and SOFLAGS.
CXX = g++
LD = g++
CXXFLAGS = -g -O3 -Wall -Wextra -Wshadow -Woverloaded-virtual -fPIC -m64
MFLAGS = -MM
SOFLAGS = -m64 -dynamiclib -single_module -undefined dynamic_lookup
endif
# --- Internal configuration ----------------------------------
PACKAGE=Laura++
DEPDIR=$(WORKDIR)/dependencies
OBJDIR=$(WORKDIR)/objects
INCLUDES += -I$(INCDIR) -I$(shell $(ROOTBINDIR)/root-config --incdir)
CXXFLAGS += $(INCLUDES)
SKIPLIST = test.cc
CINTFILE = $(WORKDIR)/$(PACKAGE)Cint.cc
CINTOBJ = $(OBJDIR)/$(PACKAGE)Cint.o
LIBFILE = $(LIBDIR)/lib$(PACKAGE).a
SHLIBFILE = $(LIBDIR)/lib$(PACKAGE).so
+ROOTMAPFILE := $(patsubst %.so,%.rootmap,$(SHLIBFILE))
default: shlib
# List of all header files
HHLIST:=$(wildcard $(INCDIR)/*.hh)
# List of all source files to build
CCLIST:=$(filter-out $(SKIPLIST),$(wildcard $(SRCDIR)/*.cc))
# List of all object files to build
OLIST:=$(patsubst %.cc,%.o,$(addprefix $(OBJDIR)/,$(notdir $(CCLIST))))
# List of all dependency files to make
DLIST:=$(patsubst %.cc,%.d,$(addprefix $(DEPDIR)/,$(notdir $(CCLIST))))
# Implicit rule making all dependency Makefiles included at the end of this makefile
$(DEPDIR)/%.d: $(SRCDIR)/%.cc
@echo "Making $@"
@mkdir -p $(DEPDIR)
@set -e; $(CXX) $(MFLAGS) $(CXXFLAGS) $< \
| sed 's#\($(notdir $*)\)\.o[ :]*#$(OBJDIR)/\1.o $@ : #g' > $@; \
[ -s $@ ] || rm -f $@
# Implicit rule to compile all classes
$(OBJDIR)/%.o : $(SRCDIR)/%.cc
@echo "Compiling $<"
@mkdir -p $(OBJDIR)
@$(CXX) $(CXXFLAGS) -c $< -o $@
# Rule to make ROOTCINT output file
$(CINTOBJ): $(HHLIST) $(INCDIR)/$(PACKAGE)_LinkDef.h
@mkdir -p $(OBJDIR)
@echo "Running rootcint"
@$(ROOTBINDIR)/rootcint -f $(CINTFILE) -c $(INCLUDES) $(notdir $(HHLIST)) $(INCDIR)/$(PACKAGE)_LinkDef.h
@echo "Compiling $(CINTFILE)"
@$(CXX) $(CXXFLAGS) -c $(CINTFILE) -o $(CINTOBJ)
# Rule to combine objects into a library
$(LIBFILE): $(OLIST) $(CINTOBJ)
@echo "Making $(LIBFILE)"
@mkdir -p $(LIBDIR)
@rm -f $(LIBFILE)
@ar rcs $(LIBFILE) $(OLIST) $(CINTOBJ)
# Rule to combine objects into a shared library
$(SHLIBFILE): $(OLIST) $(CINTOBJ)
@echo "Making $(SHLIBFILE)"
@mkdir -p $(LIBDIR)
@rm -f $(SHLIBFILE)
@$(CXX) $(OLIST) $(CINTOBJ) $(SOFLAGS) -o $(SHLIBFILE)
+ @rlibmap -f -o $(ROOTMAPFILE) -l $(SHLIBFILE) -d libMatrix libHist libTree -c $(INCDIR)/$(PACKAGE)_LinkDef.h
# Useful build targets
lib: $(LIBFILE)
shlib: $(SHLIBFILE)
clean:
rm -rf $(WORKDIR)
rm -f $(LIBFILE)
rm -f $(SHLIBFILE)
+ rm -f $(ROOTMAPFILE)
.PHONY : shlib lib default clean
-include $(DLIST)
diff --git a/doc/release.notes b/doc/release.notes
index c9d387b..4bf21ba 100644
--- a/doc/release.notes
+++ b/doc/release.notes
@@ -1,126 +1,145 @@
///////////////////////////////////////////////////////////////
/// ///
/// This is the History file for the Laura++ package. ///
/// ///
///////////////////////////////////////////////////////////////
+=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
+ 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
- - LauBkgndModel, modified to use Lau2DAbsDPPdf instead of Lau2DHistDPPdf
+ - 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
29th November 2013 Thomas Latham
* Fixed bug in ResultsExtractor where the output file was not written
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Laura++ v1r1p1
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
22nd November 2013 Thomas Latham
* Fixed bug in LauCPFitModel where values of q = -1 extra PDFs
were used regardless of the event charge.
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Laura++ v1r1
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
20th November 2013 Mark Whitehead
* Changed convention for the barrier factors, swapping from p* to p.
This seems to give more physically reasonable results.
The files updated are:
- src/LauGounarisSakuraiRes.cc
- src/LauRelBreitWignerRes.cc
18th October 2013 Thomas Latham
* Fix dependency problem in Makefile
8th October 2013 Thomas Latham
* Some fixes to yield implementation
* Minor bug fix in DP background histogram class
7th October 2013 Mark Whitehead
* Update to accept the yields and yield asymmetries as LauParameters.
All examples have been updated to reflect this change.
This updated the following files:
- inc/LauCPFitModel.hh
- inc/LauSimpleFitModel.hh
- inc/LauAbsFitModel.hh
- src/LauCPFitModel.cc
- src/LauSimpleFitModel.cc
* Addition of the following particles to src/LauResonanceMaker.cc
Ds*+-, Ds0*(2317)+-, Ds2*(2573)+-, Ds1*(2700)+- and Bs*0
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Laura++ v1r0
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
13th September 2013 Thomas Latham
* Initial import of the package into HEPforge
diff --git a/inc/Lau2DAbsHistDPPdf.hh b/inc/Lau2DAbsHistDPPdf.hh
index d647567..4937752 100644
--- a/inc/Lau2DAbsHistDPPdf.hh
+++ b/inc/Lau2DAbsHistDPPdf.hh
@@ -1,149 +1,149 @@
// Copyright University of Warwick 2004 - 2013.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
// Authors:
// Thomas Latham
// John Back
// Paul Harrison
/*! \file Lau2DAbsHistDPPdf.hh
\brief File containing declaration of Lau2DAbsHistDPPdf class.
*/
/*! \class Lau2DAbsHistDPPdf
\brief Abstract base class for defining a variation across a 2D DP based on a histogram.
Abstract base class for defining a normalised variation across a 2D DP based on a histogram.
The returned values are normalised to the total area.
The histogram can be defined in the conventional DP (m13Sq vs m23Sq) or in the square DP and
the upper half can be used to describe symmetric DPs.
*/
#ifndef LAU_2DABSHIST_DP_PDF
#define LAU_2DABSHIST_DP_PDF
#include "Lau2DAbsDPPdf.hh"
class TH2;
class LauKinematics;
class LauVetoes;
class Lau2DAbsHistDPPdf : public Lau2DAbsDPPdf {
public:
//! Constructor
/*!
\param [in] kinematics the current DP kinematics
\param [in] vetoes the vetoes within the DP
\param [in] useUpperHalfOnly boolean flag to specify that the supplied histogram contains only the upper half of a symmetric DP
\param [in] squareDP boolean flag to specify whether the supplied histogram is in square DP coordinates
*/
Lau2DAbsHistDPPdf(LauKinematics* kinematics, const LauVetoes* vetoes,
Bool_t useUpperHalfOnly = kFALSE, Bool_t squareDP = kFALSE);
//! Destructor
virtual ~Lau2DAbsHistDPPdf();
//! Retrieve maximum height of the PDF
/*!
\return the maximum height
*/
Double_t getMaxHeight() const {return maxHeight_;}
//! Retrieve PDF normalisation
/*!
\return the normalisation factor
*/
virtual Double_t getHistNorm() const=0;
//! Perform the interpolation (unnormalised)
/*!
\param [in] x the x-axis value
\param [in] y the y-aixs value
\return the unnormalised PDF value
*/
virtual Double_t interpolateXY(Double_t x, Double_t y) const=0;
//! Perform the interpolation and divide by the normalisation
/*!
\param [in] x the x-axis abscissa value
\param [in] y the y-axis abscissa value
\return the normalised PDF value
*/
virtual Double_t interpolateXYNorm(Double_t x, Double_t y) const=0;
protected:
//! Get the kinematics object
/*!
\return the kinematics
*/
const LauKinematics * getKinematics() const {return kinematics_;}
//! Get the vetoes object
/*!
\return the vetoes
*/
const LauVetoes * getVetoes() const {return vetoes_;}
//! Calculate maximum height
/*!
\param [in,out] hist the histogram
*/
void calcMaxHeight(TH2* hist);
//! Fluctuate the histogram bin contents in accordance with their errors
/*!
\param [in,out] hist the histogram
*/
void doBinFluctuation(TH2* hist);
- //! Check whether the given co-oridinates are within the kinematic boundary
+ //! 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-oridinates are within the kinematic boundary, otherwise false
+ \return true if the co-ordinates are within the kinematic boundary, otherwise false
*/
Bool_t withinDPBoundaries(Double_t x, Double_t y) const;
//! Update the current co-ordinates in the kinematic space
/*!
\param [in] x the x co-ordinate
\param [in] y the y co-ordinate
*/
void updateKinematics(Double_t x, Double_t y) const;
//! If only using the upper half of the (symmetric) DP then transform into the correct half
/*!
\param [in,out] x the x co-ordinate
\param [in,out] y the y co-ordinate
*/
void getUpperHalf(Double_t& x, Double_t& y) const;
private:
//! Copy constructor - not implemented
Lau2DAbsHistDPPdf( const Lau2DAbsHistDPPdf& rhs );
//! Copy assignment operator - not implemented
Lau2DAbsHistDPPdf& operator=(const Lau2DAbsHistDPPdf& rhs);
//! DP kinematics
LauKinematics* kinematics_;
//! Vetos within DP
const LauVetoes* vetoes_;
//! The maximum height of 2D histogram
Double_t maxHeight_;
//! Boolean for using the upper half of DP
Bool_t upperHalf_;
//! Boolean for using square DP variables
Bool_t squareDP_;
ClassDef(Lau2DAbsHistDPPdf,0) // Abstract base class for 2D DP variations based on a histogram
};
#endif
diff --git a/inc/LauBkgndDPModel.hh b/inc/LauBkgndDPModel.hh
index cd947b7..d1af760 100644
--- a/inc/LauBkgndDPModel.hh
+++ b/inc/LauBkgndDPModel.hh
@@ -1,147 +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 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] useUpperHalfOnly boolean flag to determine whether, in the case of a symmetric DP, the histogram supplied only includes the upper half of the DP
\param [in] squareDP boolean flag to determine whether the supplied histogram is given in square DP coordinates
*/
void setBkgndHisto(const TH2* histo, Bool_t useInterpolation,
Bool_t fluctuateBins, Bool_t useUpperHalfOnly,
Bool_t squareDP = kFALSE);
//! Set the background histogram and generate a spline
/*!
\param [in] histo the 2D histogram describing the DP distribution
\param [in] fluctuateBins boolean flag to determine whether the bin contents should be fluctuated in accordance with their errors
\param [in] useUpperHalfOnly boolean flag to determine whether, in the case of a symmetric DP, the histogram supplied only includes the upper half of the DP
\param [in] squareDP boolean flag to determine whether the supplied histogram is given in square DP coordinates
*/
void setBkgndSpline(const TH2* histo, Bool_t fluctuateBins,
Bool_t useUpperHalfOnly, Bool_t squareDP);
//! Initialise the model
virtual void initialise();
//! Generate a toy MC event from the model
/*!
\return success/failure flag
*/
virtual Bool_t generate();
//! Get unnormalised likelihood for a given event
/*!
\param [in] iEvt the event number
\return the unnormalised likelihood value
*/
virtual Double_t getUnNormValue(UInt_t iEvt);
//! Get PDF normalisation constant
/*!
\return the PDF normalisation constant
*/
virtual Double_t getPdfNorm() const {return pdfNorm_;}
//! Get likelihood for a given event
/*!
\param [in] iEvt the event number
\return the likelihood value
*/
virtual Double_t getLikelihood(UInt_t iEvt);
//! Cache the input data and (if appropriate) the per-event likelihood values
/*!
\param [in] fitDataTree the input data
*/
virtual void fillDataTree(const LauFitDataTree& fitDataTree);
protected:
//! Calulate histogram value at a given point
/*!
\param [in] xVal the x-value
\param [in] yVal the y-value
\return the histogram value
*/
Double_t calcHistValue(Double_t xVal, Double_t yVal);
//! Set data event number
/*!
\param [in] iEvt the event number
*/
virtual void setDataEventNo(UInt_t iEvt);
private:
//! Flags whether the DP is symmetrical or not
Bool_t symmetricalDP_;
//! Flags whether or not to work in square DP coordinates
Bool_t squareDP_;
//! PDF of Dalitz plot background, from a 2D histogram
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/src/Lau2DSplineDPPdf.cc b/src/Lau2DSplineDPPdf.cc
index 524db7c..5b71111 100644
--- a/src/Lau2DSplineDPPdf.cc
+++ b/src/Lau2DSplineDPPdf.cc
@@ -1,101 +1,98 @@
// Copyright University of Warwick 2004 - 2013.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
// Authors:
// Thomas Latham
// John Back
// Paul Harrison
/*! \file Lau2DSplineDPPdf.cc
\brief File containing implementation of Lau2DSplineDPPdf class.
*/
#include <iostream>
#include "TAxis.h"
#include "TH2.h"
#include "TRandom.h"
#include "TSystem.h"
#include "Lau2DSplineDPPdf.hh"
#include "Lau2DCubicSpline.hh"
#include "LauDaughters.hh"
#include "LauKinematics.hh"
#include "LauRandom.hh"
ClassImp(Lau2DSplineDPPdf)
Lau2DSplineDPPdf::Lau2DSplineDPPdf(const TH2* hist, LauKinematics* kinematics, const LauVetoes* vetoes,
Bool_t fluctuateBins, Bool_t useUpperHalfOnly, Bool_t squareDP) :
Lau2DAbsHistDPPdf(kinematics,vetoes,useUpperHalfOnly,squareDP),
spline_(0)
{
//We may need to modify the histogram so clone it
TH2* tempHist(hist ? dynamic_cast<TH2*>(hist->Clone()) : 0);
if ( ! tempHist ) {
std::cerr << "ERROR in Lau2DSplineDPPdf constructor : the histogram pointer is null." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
if (fluctuateBins) {
this->doBinFluctuation(tempHist);
}
spline_ = new Lau2DCubicSpline(*tempHist);
// Calculate the PDF normalisation.
- //this->calcHistNorm();
-
- // And check it is OK.
- //this->checkNormalisation();
+ this->calcHistNorm();
// Also obtain the maximum height
this->calcMaxHeight(tempHist);
delete tempHist;
}
Lau2DSplineDPPdf::~Lau2DSplineDPPdf()
{
delete spline_;
spline_ = 0;
}
Double_t Lau2DSplineDPPdf::interpolateXYNorm(Double_t x, Double_t y) const
{
// Get the normalised interpolated value.
Double_t value = this->interpolateXY(x,y);
return value/norm_;
}
Double_t Lau2DSplineDPPdf::interpolateXY(Double_t x, Double_t y) const
{
// This function returns the interpolated value of the histogram function
// for the given values of x and y by finding the adjacent bins and extrapolating
// using weights based on the inverse distance of the point from the adajcent
// bin centres.
// Here, x = m13^2, y = m23^2, or m', theta' for square DP co-ordinates
// If we're only using one half then flip co-ordinates
// appropriately for conventional or square DP
getUpperHalf(x,y);
// First ask whether the point is inside the kinematic region.
if (withinDPBoundaries(x,y) == kFALSE) {
std::cerr << "WARNING in Lau2DSplineDPPdf::interpolateXY : Given position is outside the DP boundary, returning 0.0." << std::endl;
return 0.0;
}
return spline_->evaluate(x,y);
}
void Lau2DSplineDPPdf::calcHistNorm()
{
norm_ = spline_->analyticalIntegral();
}
diff --git a/src/LauBkgndDPModel.cc b/src/LauBkgndDPModel.cc
index 633eeef..fba3b84 100644
--- a/src/LauBkgndDPModel.cc
+++ b/src/LauBkgndDPModel.cc
@@ -1,272 +1,288 @@
// Copyright University of Warwick 2004 - 2013.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
// Authors:
// Thomas Latham
// John Back
// Paul Harrison
/*! \file LauBkgndDPModel.cc
\brief File containing implementation of LauBkgndDPModel class.
*/
#include <cstdlib>
#include <iostream>
using std::cout;
using std::cerr;
using std::endl;
#include "TRandom.h"
#include "TSystem.h"
#include "Lau2DHistDPPdf.hh"
#include "Lau2DSplineDPPdf.hh"
#include "LauBkgndDPModel.hh"
#include "LauDaughters.hh"
#include "LauFitDataTree.hh"
#include "LauKinematics.hh"
#include "LauRandom.hh"
#include "LauVetoes.hh"
ClassImp(LauBkgndDPModel)
LauBkgndDPModel::LauBkgndDPModel(LauDaughters* daughters, LauVetoes* vetoes) :
LauAbsBkgndDPModel(daughters, vetoes),
symmetricalDP_(kFALSE),
squareDP_(kFALSE),
bgHistDPPdf_(0),
curEvtHistVal_(0.0),
maxPdfHeight_(1.0),
pdfNorm_(1.0),
- doneGenWarning_(kFALSE)
+ doneGenWarning_(kFALSE),
+ lowBinWarningIssued_(kFALSE)
{
if (daughters != 0) {
symmetricalDP_ = daughters->gotSymmetricalDP();
}
}
LauBkgndDPModel::~LauBkgndDPModel()
{
if (bgHistDPPdf_ != 0) {
delete bgHistDPPdf_; bgHistDPPdf_ = 0;
}
}
void LauBkgndDPModel::setBkgndHisto(const TH2* histo, Bool_t useInterpolation, Bool_t fluctuateBins, Bool_t useUpperHalfOnly, Bool_t squareDP)
{
Bool_t upperHalf = kFALSE;
if (symmetricalDP_ == kTRUE && useUpperHalfOnly == kTRUE) {upperHalf = kTRUE;}
cout<<"Bg histogram has upperHalf = "<<static_cast<Int_t>(upperHalf)<<endl;
squareDP_ = squareDP;
LauKinematics* kinematics = this->getKinematics();
const LauVetoes* vetoes = this->getVetoes();
bgHistDPPdf_ = new Lau2DHistDPPdf(histo, kinematics, vetoes, useInterpolation, fluctuateBins, upperHalf, squareDP);
maxPdfHeight_ = bgHistDPPdf_->getMaxHeight();
pdfNorm_ = bgHistDPPdf_->getHistNorm();
}
void LauBkgndDPModel::setBkgndSpline(const TH2* histo, Bool_t fluctuateBins, Bool_t useUpperHalfOnly, Bool_t squareDP)
{
Bool_t upperHalf = kFALSE;
if (symmetricalDP_ == kTRUE && useUpperHalfOnly == kTRUE) {upperHalf = kTRUE;}
cout<<"Bg histogram has upperHalf = "<<static_cast<Int_t>(upperHalf)<<endl;
squareDP_ = squareDP;
LauKinematics* kinematics = this->getKinematics();
const LauVetoes* vetoes = this->getVetoes();
bgHistDPPdf_ = new Lau2DSplineDPPdf(histo, kinematics, vetoes, fluctuateBins, upperHalf, squareDP);
maxPdfHeight_ = bgHistDPPdf_->getMaxHeight();
pdfNorm_ = bgHistDPPdf_->getHistNorm();
}
Double_t LauBkgndDPModel::calcHistValue(Double_t xVal, Double_t yVal)
{
// Get the likelihood value of the background in the Dalitz plot.
// Check that we have a valid histogram PDF
if (bgHistDPPdf_ == 0) {
cerr << "WARNING in LauBkgndDPModel::calcHistValue : We don't have a histogram so assuming the likelihood is flat in the Dalitz plot." << endl;
this->setBkgndHisto( 0, kFALSE, kFALSE, kFALSE, kFALSE );
}
// Find out the un-normalised PDF value
Double_t value = bgHistDPPdf_->interpolateXY(xVal, yVal);
+ // Check that the value is greater than zero
+ // If we're using a spline then negative values can be caused by adjacent bins that all contain a value of zero.
+ // The spline requires the value, 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. Such dips are unavoidable but are correctly removed here.
+ if ( value < 0.0 ) {
+ if(!lowBinWarningIssued_) {
+ std::cerr << "WARNING in LauBkgndDPModel::calcHistValue : Value " << value << " 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;
+ }
+ return 0.0;
+ }
+
LauKinematics* kinematics = this->getKinematics();
// For square DP co-ordinates, need to divide by Jacobian
if (squareDP_ == kTRUE) {
// Make sure that square DP kinematics are correct, then get Jacobian
kinematics->updateSqDPKinematics(xVal, yVal);
Double_t jacobian = kinematics->calcSqDPJacobian();
value /= jacobian;
}
return value;
}
void LauBkgndDPModel::initialise()
{
}
Bool_t LauBkgndDPModel::generate()
{
// Routine to generate the background, using data provided by an
// already defined histogram.
LauKinematics* kinematics = this->getKinematics();
Bool_t gotBG(kFALSE);
while (gotBG == kFALSE) {
if (squareDP_ == kTRUE) {
// Generate a point in m', theta' space. By construction, this point
// is already within the DP region.
Double_t mPrime(0.0), thetaPrime(0.0);
kinematics->genFlatSqDP(mPrime, thetaPrime);
// If we're in a symmetrical DP then we should only generate events in one half
if ( symmetricalDP_ && thetaPrime > 0.5 ) {
thetaPrime = 1.0 - thetaPrime;
}
// Calculate histogram height for DP point and
// compare with the maximum height
if ( bgHistDPPdf_ != 0 ) {
Double_t bgContDP = bgHistDPPdf_->interpolateXY(mPrime, thetaPrime)/maxPdfHeight_;
if (LauRandom::randomFun()->Rndm() < bgContDP) {
kinematics->updateSqDPKinematics(mPrime, thetaPrime);
gotBG = kTRUE;
}
} else {
if ( !doneGenWarning_ ) {
cerr << "WARNING in LauBkgndDPModel::generate : We don't have a histogram so generating flat in the square DP, which won't be flat in the conventional DP!" << endl;
cerr << "WARNING in LauBkgndDPModel::generate : This should never happen!! What have you done?!" << endl;
doneGenWarning_ = kTRUE;
}
kinematics->updateSqDPKinematics(mPrime, thetaPrime);
gotBG = kTRUE;
}
} else {
// Generate a point in the Dalitz plot (phase-space).
Double_t m13Sq(0.0), m23Sq(0.0);
kinematics->genFlatPhaseSpace(m13Sq, m23Sq);
// If we're in a symmetrical DP then we should only generate events in one half
if ( symmetricalDP_ && m13Sq > m23Sq ) {
Double_t tmpSq = m13Sq;
m13Sq = m23Sq;
m23Sq = tmpSq;
}
// Calculate histogram height for DP point and
// compare with the maximum height
if ( bgHistDPPdf_ != 0 ) {
Double_t bgContDP = bgHistDPPdf_->interpolateXY(m13Sq, m23Sq)/maxPdfHeight_;
if (LauRandom::randomFun()->Rndm() < bgContDP) {
kinematics->updateKinematics(m13Sq, m23Sq);
gotBG = kTRUE;
}
} else {
if ( !doneGenWarning_ ) {
cerr << "WARNING in LauBkgndDPModel::generate : We don't have a histogram so generating flat in the DP." << endl;
doneGenWarning_ = kTRUE;
}
kinematics->updateKinematics(m13Sq, m23Sq);
gotBG = kTRUE;
}
}
}
// Implement veto
Bool_t vetoOK(kTRUE);
const LauVetoes* vetoes = this->getVetoes();
if (vetoes) {
vetoOK = vetoes->passVeto(kinematics);
}
// Call this function recusively until we pass the veto.
if (vetoOK == kFALSE) {
this->generate();
}
return kTRUE;
}
void LauBkgndDPModel::fillDataTree(const LauFitDataTree& inputFitTree)
{
// In LauFitDataTree, the first two variables should always be
// m13^2 and m23^2. Other variables follow thus: charge.
Int_t nBranches = inputFitTree.nBranches();
if (nBranches < 2) {
cout<<"Error in LauBkgndDPModel::fillDataTree. Expecting at least 2 variables "
<<"in input data tree, but there are "<<nBranches<<"! Make sure you have "
<<"the right number of variables in your input data file!"<<endl;
return;
}
Double_t m13Sq(0.0), m23Sq(0.0), mPrime(0.0), thetaPrime(0.0);
UInt_t nEvents = inputFitTree.nEvents();
LauKinematics* kinematics = this->getKinematics();
// clear and resize the data vector
bgData_.clear();
bgData_.resize(nEvents);
for (UInt_t iEvt = 0; iEvt < nEvents; iEvt++) {
const LauFitData& dataValues = inputFitTree.getData(iEvt);
LauFitData::const_iterator iter = dataValues.find("m13Sq");
m13Sq = iter->second;
iter = dataValues.find("m23Sq");
m23Sq = iter->second;
// Update the kinematics. This will also update m' and theta' if squareDP = true
kinematics->updateKinematics(m13Sq, m23Sq);
if (squareDP_ == kTRUE) {
mPrime = kinematics->getmPrime();
thetaPrime = kinematics->getThetaPrime();
curEvtHistVal_ = this->calcHistValue(mPrime, thetaPrime);
} else {
curEvtHistVal_ = this->calcHistValue(m13Sq, m23Sq);
}
bgData_[iEvt] = curEvtHistVal_;
}
}
Double_t LauBkgndDPModel::getUnNormValue(UInt_t iEvt)
{
// Retrieve the likelihood for the given event
this->setDataEventNo(iEvt);
return curEvtHistVal_;
}
Double_t LauBkgndDPModel::getLikelihood(UInt_t iEvt)
{
// Retrieve the likelihood for the given event
this->setDataEventNo(iEvt);
Double_t llhd = curEvtHistVal_ / this->getPdfNorm();
return llhd;
}
void LauBkgndDPModel::setDataEventNo(UInt_t iEvt)
{
// Retrieve the data for event iEvt
if (bgData_.size() > iEvt) {
curEvtHistVal_ = bgData_[iEvt];
} else {
cerr<<"ERROR in LauBkgndDPModel::setDataEventNo : Event index too large: "<<iEvt<<" >= "<<bgData_.size()<<"."<<endl;
}
}
diff --git a/src/LauEffModel.cc b/src/LauEffModel.cc
index cff25dc..464d101 100644
--- a/src/LauEffModel.cc
+++ b/src/LauEffModel.cc
@@ -1,191 +1,191 @@
// 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::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;
}
}
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;
+ << " : 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;
+ << " : 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, 3:00 PM (1 d, 14 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3799080
Default Alt Text
(39 KB)

Event Timeline