Page MenuHomeHEPForge
Files F10881363

No OneTemporary
Actions

View Options

This file is larger than 256 KB, so syntax highlighting was skipped.
diff --git a/doc/mainpage.dox.in b/doc/mainpage.dox.in
index 3099fbc..9185666 100644
--- a/doc/mainpage.dox.in
+++ b/doc/mainpage.dox.in
@@ -1,33 +1,33 @@
/*!
\copyright
-Copyright 2004-2019 University of Warwick<br>
+Copyright 2004-2023 University of Warwick<br>
Licensed under the Apache License, Version 2.0 (the "License");
you may not use these files except in compliance with the License.<br>
You may obtain a copy of the License at<br>
http://www.apache.org/licenses/LICENSE-2.0<br>
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.<br>
See the License for the specific language governing permissions and
limitations under the License.
*/
/*!
\authors
Laura++ package authors:<br>
John Back,
Paul Harrison,
Thomas Latham
*/
/*!
\mainpage
Welcome to the Laura++ documentation, please use the links at the top of the page to navigate to particular classes or namespaces.
The Laura++ project webpages can be found at @PROJECT_HOMEPAGE_URL@
See also the [README](@PROJECT_SOURCE_DIR@/README.md) and [release notes](@PROJECT_SOURCE_DIR@/doc/ReleaseNotes.md) for more information.
*/
diff --git a/inc/Lau1DCubicSpline.hh b/inc/Lau1DCubicSpline.hh
index 95d0241..10e8893 100644
--- a/inc/Lau1DCubicSpline.hh
+++ b/inc/Lau1DCubicSpline.hh
@@ -1,206 +1,207 @@
/*
Copyright 2015 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file Lau1DCubicSpline.hh
\brief File containing declaration of Lau1DCubicSpline class.
*/
/*! \class Lau1DCubicSpline
\brief Class for defining a 1D cubic spline based on a set of knots.
Class for defining a 1D cubic spline based on a set of knots.
Interpolation between the knots is performed either by one of
two types of cubic spline (standard or Akima) or by linear interpolation.
The splines are defined by a piecewise cubic function which, between knots i and i+1, has the form
f_i(x) = (1 - t)*y_i + t*y_i+1 + t*(1 - t)(a*(1 - t) + b*t)
where t = (x - x_i)/(x_i+1 - x_i),
a = k_i *(x_i+1 - x_i) - (y_i+1 - y_i),
b = -k_i+1*(x_i+1 - x_i) + (y_i+1 - y_i)
and k_i is (by construction) the first derivative at knot i.
f(x) and f'(x) are continuous at the internal knots by construction.
For the standard splines, f''(x) is required to be continuous
at all internal knots placing n-2 constraints on the n parameters, k_i.
The final two constraints are set by the boundary conditions.
At each boundary, the function may be:
(i) Clamped : f'(x) = C at the last knot
(ii) Natural : f''(x) = 0 at the last knot
(iii) Not a knot : f'''(x) continuous at the second last knot
The algorithms used in these splines can be found on:
http://en.wikipedia.org/wiki/Spline_interpolation#Algorithm_to_find_the_interpolating_cubic_spline
http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm
For the Akima splines, the values of k_i are determined from the slopes of the four nearest segments (a_i-1, a_i, a_i+1 and a_i+2) as
k_i = ( | a_i+2 - a_i+1 | a_i + | a_i - a_i-1 | a_i+1 ) / ( | a_i+2 - a_i+1 | + | a_i - a_i-1 | )
and as
k_i = ( a_i + a_i+1 ) / 2
in the special case a_i-1 == a_i != a_i+1 == a_i+2.
Boundary conditions are specified by the relations
a_2 - a_1 = a_1 - a_0 = a_0 - a_-1 and
a_n-1 - a_n-2 = a_n - a_n-1 = a_n+1 - a_n
The algorithms used in these splines can be found in:
J.ACM vol. 17 no. 4 pp 589-602
*/
#ifndef LAU_1DCUBICSPLINE
#define LAU_1DCUBICSPLINE
#include <vector>
#include "Rtypes.h"
class Lau1DCubicSpline {
public:
//! Define the allowed interpolation types
enum LauSplineType {
StandardSpline, /*!< standard cubic splines with f''(x) continuous at all internal knots */
AkimaSpline, /*!< Akima cubic splines with f'(x) at each knot defined locally by the positions of only five knots */
LinearInterpolation /*! Linear interpolation between each pair of knots */
};
//! Define the allowed boundary condition types
/*!
These are only supported by standard splines
*/
enum LauSplineBoundaryType {
Clamped, /*!< clamped boundary - f'(x) = C */
Natural, /*!< natural boundary - f''(x) = 0 */
NotAKnot /*!< 'not a knot' boundary - f'''(x) continuous at second last knot */
};
//! Constructor
/*!
/param [in] xs the x-values of the knots
/param [in] ys the y-values of the knots
+ /param [in] type the type of the spline (e.g. StandardSpline, AkimaSpline, LinearInterpolation)
/param [in] leftBound the left-hand boundary condition
/param [in] rightBound the right-hand boundary condition
/param [in] dydx0 the gradient at the left-hand boundary of a clamped spline
/param [in] dydxn the gradient at the right-hand boundary of a clamped spline
*/
Lau1DCubicSpline(const std::vector<Double_t>& xs, const std::vector<Double_t>& ys,
LauSplineType type = Lau1DCubicSpline::StandardSpline,
LauSplineBoundaryType leftBound = Lau1DCubicSpline::NotAKnot,
LauSplineBoundaryType rightBound = Lau1DCubicSpline::NotAKnot,
Double_t dydx0 = 0.0, Double_t dydxn = 0.0);
//! Destructor
virtual ~Lau1DCubicSpline();
//! Evaluate the function at given point
/*!
\param [in] x the x-coordinate
\return the value of the spline at x
*/
Double_t evaluate(Double_t x) const;
//! Update the y-values of the knots
/*!
\param [in] ys the y-values of the knots
*/
void updateYValues(const std::vector<Double_t>& ys);
//! Update the type of interpolation to perform
/*!
\param [in] type the type of interpolation
*/
void updateType(LauSplineType type);
//! Update the boundary conditions for the spline
/*!
/param [in] leftBound the left-hand boundary condition
/param [in] rightBound the right-hand boundary condition
/param [in] dydx0 the gradient at the left-hand boundary of a clamped spline
/param [in] dydxn the gradient at the right-hand boundary of a clamped spline
*/
void updateBoundaryConditions(LauSplineBoundaryType leftBound,
LauSplineBoundaryType rightBound,
Double_t dydx0 = 0.0,
Double_t dydxn = 0.0);
private:
//! Copy constructor - not implemented
Lau1DCubicSpline( const Lau1DCubicSpline& rhs );
//! Copy assignment operator - not implemented
Lau1DCubicSpline& operator=(const Lau1DCubicSpline& rhs);
//! Initialise the class
void init();
//! Calculate the first derivative at each knot
void calcDerivatives();
//! Calculate the first derivatives according to the standard method
void calcDerivativesStandard();
//! Calculate the first derivatives according to the Akima method
void calcDerivativesAkima();
//! The number of knots in the spline
const UInt_t nKnots_;
//! The x-value at each knot
std::vector<Double_t> x_;
//! The y-value at each knot
std::vector<Double_t> y_;
//! The first derivative at each knot
std::vector<Double_t> dydx_;
//! The 'a' coefficients used to determine the derivatives
std::vector<Double_t> a_;
//! The 'b' coefficients used to determine the derivatives
std::vector<Double_t> b_;
//! The 'c' coefficients used to determine the derivatives
std::vector<Double_t> c_;
//! The 'd' coefficients used to determine the derivatives
std::vector<Double_t> d_;
//! The type of interpolation to be performed
LauSplineType type_;
//! The left-hand boundary condition on the spline
LauSplineBoundaryType leftBound_;
//! The right-hand boundary condition on the spline
LauSplineBoundaryType rightBound_;
//! The gradient at the left boundary for a clamped spline
Double_t dydx0_;
//! The gradient at the right boundary for a clamped spline
Double_t dydxn_;
ClassDef(Lau1DCubicSpline, 0); // Class for defining a 1D cubic spline
};
#endif
diff --git a/inc/LauAbsResonance.hh b/inc/LauAbsResonance.hh
index 02a6b1a..adf337a 100644
--- a/inc/LauAbsResonance.hh
+++ b/inc/LauAbsResonance.hh
@@ -1,579 +1,581 @@
/*
Copyright 2004 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauAbsResonance.hh
\brief File containing declaration of LauAbsResonance class.
*/
/*! \class LauAbsResonance
\brief Abstract class for defining type for resonance amplitude models (Breit-Wigner, Flatte etc.)
Abstract Class for defining the type for all classes used to model
resonances in the Dalitz plot, such as Breit-Wigner functions.
In addition, some common functionality is implemented, including data such
as the mass and width of the desired state.
*/
#ifndef LAU_ABS_RESONANCE
#define LAU_ABS_RESONANCE
#include "TString.h"
#include "LauBlattWeisskopfFactor.hh"
#include "LauComplex.hh"
#include "LauParameter.hh"
class LauDaughters;
class LauKinematics;
class LauResonanceInfo;
class LauAbsResonance {
public:
//! Define the allowed resonance types
enum LauResonanceModel {
BW, /*!< simple Breit-Wigner */
RelBW, /*!< relativistic Breit-Wigner */
GS, /*!< a modified Breit-Wigner from Gounaris-Sakurai */
Flatte, /*!< Flatte or coupled-channel Breit-Wigner */
Sigma, /*!< special shape for the sigma or f_0(600) */
Kappa, /*!< special shape for the kappa, a low-mass Kpi scalar */
Dabba, /*!< special shape for the dabba, a low-mass Dpi scalar */
LASS, /*!< the LASS amplitude to describe the Kpi S-wave */
LASS_BW, /*!< the resonant part of the LASS amplitude */
LASS_NR, /*!< the nonresonant part of the LASS amplitude */
EFKLLM, /*!< a form-factor-based description of the Kpi S-wave */
KMatrix, /*!< description using K-matrix and P-vector */
FlatNR, /*!< a uniform nonresonant amplitude */
NRModel, /*!< a theoretical model nonresonant amplitude */
BelleNR, /*!< an empirical exponential nonresonant amplitude */
PowerLawNR, /*!< an empirical power law nonresonant amplitude */
BelleSymNR, /*!< an empirical exponential nonresonant amplitude for symmetrised DPs */
BelleSymNRNoInter, /*!< an empirical exponential nonresonant amplitude for symmetrised DPs without interference */
TaylorNR, /*!< an empirical Taylor expansion nonresonant amplitude for symmetrised DPs */
PolNR, /*!< an empirical polynomial nonresonant amplitude */
Pole, /*!< scalar Pole lineshape */
PolarFFNR, /*!< Polar Form Factor nonresonant amplitude */
PolarFFSymNR, /*!< Polar Form Factor nonresonant amplitude for symmetrised DPs */
PolarFFSymNRNoInter, /*!< Polar Form Factor nonresonant amplitude for symmetrised DPs without interference */
Rescattering, /*!< KK-PiPi inelastic scattering amplitude */
Rescattering2, /*!< KK-PiPi inelastic scattering amplitude */
RescatteringNoInter, /*!< KK-PiPi inelastic scattering amplitude */
MIPW_MagPhase, /*!< a model independent partial wave - magnitude and phase representation */
MIPW_RealImag, /*!< a model independent partial wave - real and imaginary part representation */
GaussIncoh, /*!< an incoherent Gaussian shape */
RhoOmegaMix_GS, /*!< mass mixing model using GS for res 1 and RBW for res 2 */
RhoOmegaMix_RBW, /*!< mass mixing model using two RBWs */
RhoOmegaMix_GS_1, /*!< mass mixing model using GS for res 1 and RBW for res 2, with denominator factor = 1 */
RhoOmegaMix_RBW_1 /*!< mass mixing model using two RBWs, with denominator factor = 1 */
};
//! Define the allowed spin formalisms
enum LauSpinType {
Zemach_P, /*!< Zemach tensor formalism, bachelor momentum in resonance rest frame */
Zemach_Pstar, /*!< Zemach tensor formalism, bachelor momentum in parent rest frame */
Covariant, /*!< Covariant tensor formalism */
Legendre /*!< Legendre polynomials only */
};
//! Is the resonance model incoherent?
/*!
\param [in] model the resonance model
\return true if the model is incoherent
*/
static bool isIncoherentModel(LauResonanceModel model);
//! Constructor (for use by standard resonances)
/*!
\param [in] resInfo the object containing information on the resonance name, mass, width, spin, charge, etc.
\param [in] resPairAmpInt the number of the daughter not produced by the resonance
\param [in] daughters the daughter particles
*/
LauAbsResonance(LauResonanceInfo* resInfo, const Int_t resPairAmpInt, const LauDaughters* daughters);
//! Constructor (for use by K-matrix components)
/*!
\param [in] resName the name of the component
\param [in] resPairAmpInt the number of the daughter not produced by the resonance
\param [in] daughters the daughter particles
\param [in] resSpin the spin of the final channel into which the K-matrix scatters
*/
LauAbsResonance(const TString& resName, const Int_t resPairAmpInt, const LauDaughters* daughters, const Int_t resSpin);
//! Destructor
virtual ~LauAbsResonance();
//! Initialise the model
virtual void initialise() = 0;
//! Calculate the complex amplitude
/*!
\param [in] kinematics the kinematic variables of the current event
\return the complex amplitude
*/
virtual LauComplex amplitude(const LauKinematics* kinematics);
//! Get the resonance model type
/*!
\return the resonance model type
*/
virtual LauResonanceModel getResonanceModel() const = 0;
//! Get the spin type
/*!
\return the spin formalism
*/
LauSpinType getSpinType() const {return spinType_;}
//! Get the name of the resonance
/*!
\return the resonance name
*/
const TString& getResonanceName() const {return resName_;}
//! Get the name of the resonance
/*!
\return the resonance name
*/
const TString& getSanitisedName() const {return sanitisedName_;}
//! Get the integer to identify which DP axis the resonance belongs to
/*!
\return the DP axis identification number, the ID of the bachelor
*/
Int_t getPairInt() const {return resPairAmpInt_;}
//! Get the spin of the resonance
/*!
\return the resonance spin
*/
Int_t getSpin() const {return resSpin_;}
//! Get the charge of the resonance
/*!
\return the resonance charge
*/
Int_t getCharge() const {return resCharge_;}
//! Get the mass of the resonance
/*!
\return the resonance mass
*/
Double_t getMass() const {return (resMass_!=0) ? resMass_->unblindValue() : -1.0;}
//! Get the width of the resonance
/*!
\return the resonance width
*/
Double_t getWidth() const {return (resWidth_!=0) ? resWidth_->unblindValue() : -1.0;}
//! Get the mass parameter of the resonance
/*!
\return the resonance mass parameter
*/
LauParameter* getMassPar() {return resMass_;}
//! Get the width parameter of the resonance
/*!
\return the resonance width parameter
*/
LauParameter* getWidthPar() {return resWidth_;}
//! Retrieve the resonance parameters, e.g. so that they can be loaded into a fit
/*!
\return floating parameters of the resonance
*/
virtual const std::vector<LauParameter*>& getFloatingParameters() { return this->getParameters(); };
//! Is the amplitude pre-symmetrised?
/*!
The default value is kFALSE, so pre-symmetrised lineshapes should override this.
\return whether the amplitude is already symmetrised
*/
virtual Bool_t preSymmetrised() const {return kFALSE;}
//! Get the helicity flip flag
/*!
\return the flip helicity flag
*/
Bool_t flipHelicity() const {return flipHelicity_;}
//! Set the helicity flip flag
/*!
\param [in] boolean the helicity flip status
*/
void flipHelicity(const Bool_t boolean) {flipHelicity_ = boolean;}
//! Get the ignore momenta flag
/*!
Whether to ignore the momentum factors in both the spin factor and the mass-dependent width
\return the ignore momenta flag
*/
Bool_t ignoreMomenta() const {return ignoreMomenta_;}
//! Set the ignore momenta flag
/*!
Whether to ignore the momentum factors in both the spin factor and the mass-dependent width
\param [in] boolean the ignore momenta status
*/
void ignoreMomenta(const Bool_t boolean) {ignoreMomenta_ = boolean;}
//! Get the ignore spin flag
/*!
Whether to set the spinTerm to unity always
\return the ignore spin flag
*/
Bool_t ignoreSpin() const {return ignoreSpin_;}
//! Set the ignore spin flag
/*!
Whether to set the spinTerm to unity always
\param [in] boolean the ignore spin status
*/
void ignoreSpin(const Bool_t boolean) {ignoreSpin_ = boolean;}
//! Get the ignore barrier factor scaling flag
/*!
Whether to ignore barrier factor scaling in the amplitude numerator, they are still used for the mass-dependent width
\return the ignore barrier amplitude scaling flag
*/
Bool_t ignoreBarrierScaling() const {return ignoreBarrierScaling_;}
//! Set the ignore barrier factor scaling flag
/*!
Whether to ignore barrier factor scaling in the amplitude numerator, they are still used for the mass-dependent width
\param [in] boolean the ignore barrier factor scaling status
*/
void ignoreBarrierScaling(const Bool_t boolean) {ignoreBarrierScaling_ = boolean;}
//! Allow the mass, width and spin of the resonance to be changed
/*!
Negative values wil be ignored, so if, for example, you
want to only change the spin you can provide negative
values for the mass and width
\param [in] newMass new value of the resonance mass
\param [in] newWidth new value of the resonance width
\param [in] newSpin new value of the resonance spin
*/
void changeResonance(const Double_t newMass, const Double_t newWidth, const Int_t newSpin);
//! Allow the Blatt-Weisskopf radius for the resonance and parent factors to be changed
/*!
Negative values wil be ignored, so if, for example, you
want to only change the parent radius you can provide a
negative value for the resonance radius
\param [in] resRadius new value of the resonance radius
\param [in] parRadius new value of the parent radius
*/
void changeBWBarrierRadii(const Double_t resRadius, const Double_t parRadius);
//! Set value of the various parameters
/*!
\param [in] name the name of the parameter to be changed
\param [in] value the new parameter value
*/
virtual void setResonanceParameter(const TString& name, const Double_t value);
//! Allow the various parameters to float in the fit
/*!
\param [in] name the name of the parameter to be floated
*/
virtual void floatResonanceParameter(const TString& name);
//! Access the given resonance parameter
/*!
\param [in] name the name of the parameter
\return the corresponding parameter
*/
virtual LauParameter* getResonanceParameter(const TString& name);
//! Fix or release the resonance mass
/*!
\param [in] parFixed new status of mass
*/
void fixMass(const Bool_t parFixed) { if (resMass_!=0) { resMass_->fixed(parFixed); } }
//! Fix or release the resonance width
/*!
\param [in] parFixed new status of width
*/
void fixWidth(const Bool_t parFixed) { if (resWidth_!=0) { resWidth_->fixed(parFixed); } }
//! Get the status of resonance mass (fixed or released)
/*!
\return the status of resonance mass (fixed or released)
*/
Bool_t fixMass() const { return (resMass_!=0) ? resMass_->fixed() : kTRUE; }
//! Get the status of resonance width (fixed or released)
/*!
\return the status of resonance width (fixed or released)
*/
Bool_t fixWidth() const { return (resWidth_!=0) ? resWidth_->fixed() : kTRUE; }
//! Set the spin formalism to be used
/*!
\param [in] spinType the spin formalism
*/
void setSpinType(const LauSpinType spinType) {spinType_ = spinType;}
//! Set the form factor model and parameters
/*!
\param [in] resFactor the barrier factor for the resonance decay
\param [in] parFactor the barrier factor for the parent decay
*/
void setBarrierRadii(LauBlattWeisskopfFactor* resFactor, LauBlattWeisskopfFactor* parFactor)
{
resBWFactor_ = resFactor;
parBWFactor_ = parFactor;
}
//! Fix or release the Blatt-Weisskopf barrier radii
void fixBarrierRadii(const Bool_t fixResRadius, const Bool_t fixParRadius);
//! Get the status of resonance barrier radius (fixed or released)
Bool_t fixResRadius() const;
//! Get the status of parent barrier radius (fixed or released)
Bool_t fixParRadius() const;
//! Get the radius of the resonance barrier factor
Double_t getResRadius() const;
//! Get the radius of the parent barrier factor
Double_t getParRadius() const;
protected:
//! Get the name of the parent particle
TString getNameParent() const;
//! Get the name of the first daughter of the resonance
TString getNameDaug1() const;
//! Get the name of the second daughter of the resonance
TString getNameDaug2() const;
//! Get the name of the daughter that does not originate form the resonance
TString getNameBachelor() const;
//! Get the parent particle mass
Double_t getMassParent() const;
//! Get the mass of daughter 1
Double_t getMassDaug1() const;
//! Get the mass of daughter 2
Double_t getMassDaug2() const;
//! Get the mass of the bachelor daughter
Double_t getMassBachelor() const;
//! Get the Charge of the parent particle
Int_t getChargeParent() const;
//! Get the charge of daughter 1
Int_t getChargeDaug1() const;
//! Get the charge of daughter 2
Int_t getChargeDaug2() const;
//! Get the charge of the bachelor daughter
Int_t getChargeBachelor() const;
//! Get the current value of the daughter momentum in the resonance rest frame
Double_t getQ() const {return q_;}
//! Get the current value of the bachelor momentum in the resonance rest frame
Double_t getP() const {return p_;}
//! Get the current value of the bachelor momentum in the parent rest frame
Double_t getPstar() const {return pstar_;}
//! Get the current value of the full spin-dependent covariant factor
Double_t getCovFactor() const {return covFactor_;}
//! Get the centrifugal barrier for the parent decay
LauBlattWeisskopfFactor* getParBWFactor() {return parBWFactor_;}
+ //! Get the centrifugal barrier for the parent decay
const LauBlattWeisskopfFactor* getParBWFactor() const {return parBWFactor_;}
//! Get the centrifugal barrier for the resonance decay
LauBlattWeisskopfFactor* getResBWFactor() {return resBWFactor_;}
+ //! Get the centrifugal barrier for the resonance decay
const LauBlattWeisskopfFactor* getResBWFactor() const {return resBWFactor_;}
//! Access the resonance info object
LauResonanceInfo* getResInfo() const {return resInfo_;}
//! Access the daughters object
const LauDaughters* getDaughters() const {return daughters_;}
//! Calculate the amplitude spin term using the Zemach tensor formalism
/*!
\param [in] pProd the momentum factor (either q * p or q * pstar)
*/
Double_t calcZemachSpinFactor( const Double_t pProd ) const;
//! Calculate the amplitude spin term using the covariant tensor formalism
/*!
\param [in] pProd the momentum factor (q * pstar)
*/
Double_t calcCovSpinFactor( const Double_t pProd );
//! Calculate the spin-dependent covariant factor
/*!
\param [in] erm E_ij in the parent rest-frame divided by m_ij (equivalent to sqrt(1 + p^2/mParent^2))
*/
void calcCovFactor( const Double_t erm );
//! Calculate the Legendre polynomial for the spin factor
/*!
Uses the current-event value of cosHel_
*/
Double_t calcLegendrePoly() const;
//! Calculate the Legendre polynomial for the spin factor (specifying the cosHel value)
/*!
\param [in] cosHel the cosine of the helicity angle
*/
Double_t calcLegendrePoly( const Double_t cosHel );
//! Complex resonant amplitude
/*!
\param [in] mass appropriate invariant mass for the resonance
\param [in] spinTerm spin term
*/
virtual LauComplex resAmp(Double_t mass, Double_t spinTerm) = 0;
//! Clear list of floating parameters
void clearFloatingParameters() { resParameters_.clear(); }
//! Add parameter to the list of floating parameters
/*!
\param [in] param the parameter to be added to the list
*/
void addFloatingParameter( LauParameter* param );
//! Access the list of floating parameters
std::vector<LauParameter*>& getParameters() { return resParameters_; }
private:
//! Copy constructor (not implemented)
LauAbsResonance(const LauAbsResonance& rhs);
//! Copy assignment operator (not implemented)
LauAbsResonance& operator=(const LauAbsResonance& rhs);
//! Information on the resonance
LauResonanceInfo* resInfo_{0};
//! Information on the particles
const LauDaughters* daughters_;
//! Parent name
TString nameParent_{""};
//! Daughter 1 name
TString nameDaug1_{""};
//! Daughter 2 name
TString nameDaug2_{""};
//! Bachelor name
TString nameBachelor_{""};
//! Parent charge
Int_t chargeParent_{0};
//! Daughter 1 charge
Int_t chargeDaug1_{0};
//! Daughter 2 charge
Int_t chargeDaug2_{0};
//! Bachelor charge
Int_t chargeBachelor_{0};
//! Parent mass
Double_t massParent_{0.0};
//! Daughter 1 mass
Double_t massDaug1_{0.0};
//! Daughter 2 mass
Double_t massDaug2_{0.0};
- // Bachelor mass
+ //! Bachelor mass
Double_t massBachelor_{0.0};
//! Resonance name
TString resName_;
//! Resonance name with illegal characters removed
TString sanitisedName_;
//! Resonance mass
LauParameter* resMass_{0};
//! Resonance width
LauParameter* resWidth_{0};
//! All parameters of the resonance
std::vector<LauParameter*> resParameters_;
//! Resonance spin
Int_t resSpin_;
//! Resonance charge
Int_t resCharge_{0};
//! DP axis identifier
Int_t resPairAmpInt_;
//! Blatt Weisskopf barrier for parent decay
LauBlattWeisskopfFactor* parBWFactor_{0};
//! Blatt Weisskopf barrier for resonance decay
LauBlattWeisskopfFactor* resBWFactor_{0};
//! Spin formalism
LauSpinType spinType_{Zemach_P};
//! Boolean to flip helicity
Bool_t flipHelicity_{kFALSE};
//! Boolean to ignore the momentum factors in both the spin factor and the mass-dependent width
Bool_t ignoreMomenta_{kFALSE};
//! Boolean to set the spinTerm to unity always
Bool_t ignoreSpin_{kFALSE};
//! Boolean to ignore barrier factor scaling in the amplitude numerator, they are still used for the mass-dependent width
Bool_t ignoreBarrierScaling_{kFALSE};
// Event kinematics information
//! Invariant mass
Double_t mass_{0.0};
//! Helicity angle cosine
Double_t cosHel_{0.0};
//! Daughter momentum in resonance rest frame
Double_t q_{0.0};
//! Bachelor momentum in resonance rest frame
Double_t p_{0.0};
//! Bachelor momentum in parent rest frame
Double_t pstar_{0.0};
//! Covariant factor
/*!
sqrt(1 + z*z), where z = p / mParent
Can also be expressed as E_ij in the parent rest-frame divided by m_ij - indeed this is how LauKinematics calculates it.
\see LauKinematics::getcov12
\see LauKinematics::getcov13
\see LauKinematics::getcov23
*/
Double_t erm_{1.0};
//! Covariant factor (full spin-dependent expression)
Double_t covFactor_{1.0};
ClassDef(LauAbsResonance,0) // Abstract resonance class
};
#endif
diff --git a/inc/LauCPFitModel.hh b/inc/LauCPFitModel.hh
index 357e123..e25761d 100644
--- a/inc/LauCPFitModel.hh
+++ b/inc/LauCPFitModel.hh
@@ -1,730 +1,732 @@
/*
Copyright 2004 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauCPFitModel.hh
\brief File containing declaration of LauCPFitModel class.
*/
/*! \class LauCPFitModel
\brief Class for defining a CP fit model.
LauCPFitModel is a class that allows the user to define a three-body Dalitz
plot according to the isobar model, i.e. defining a set of resonances that
have complex amplitudes that can interfere with each other.
It extends the LauSimpleFitModel in that it allows simultaneous fitting of
two parent flavours simultaneously. By default, it assumes perfect tagging
of those flavours but it can also be used in an untagged scenario.
*/
#ifndef LAU_CP_FIT_MODEL
#define LAU_CP_FIT_MODEL
#include <vector>
#include "TString.h"
#include "LauAbsFitModel.hh"
#include "LauComplex.hh"
#include "LauParameter.hh"
class TH2;
class LauAbsBkgndDPModel;
class LauAbsCoeffSet;
class LauIsobarDynamics;
class LauAbsPdf;
class LauEffModel;
class LauEmbeddedData;
class LauKinematics;
class LauScfMap;
class LauCPFitModel : public LauAbsFitModel {
public:
//! Constructor
/*!
\param [in] negModel DP model for the antiparticle
\param [in] posModel DP model for the particle
\param [in] tagged is the analysis tagged or untagged?
\param [in] tagVarName the variable name in the data tree that specifies the event tag
*/
LauCPFitModel(LauIsobarDynamics* negModel, LauIsobarDynamics* posModel, Bool_t tagged = kTRUE, const TString& tagVarName = "charge");
//! Destructor
virtual ~LauCPFitModel();
//! Set the signal event yield
/*!
\param [in] nSigEvents contains the signal yield and option to fix it
*/
virtual void setNSigEvents(LauParameter* nSigEvents);
//! Set the signal event yield if there is an asymmetry
/*!
\param [in] nSigEvents contains the signal yield and option to fix it
\param [in] sigAsym contains the signal asymmetry and option to fix it
\param [in] forceAsym the option to force there to be an asymmetry
*/
virtual void setNSigEvents(LauParameter* nSigEvents, LauParameter* sigAsym, Bool_t forceAsym = kFALSE);
//! Set the background event yield(s)
/*!
The name of the parameter must be that of the corresponding background category (so that it can be correctly assigned)
\param [in] nBkgndEvents contains the name, yield and option to fix the yield of the background
*/
virtual void setNBkgndEvents(LauAbsRValue* nBkgndEvents);
//! Set the background event yield(s)
/*!
The names of the parameters must be that of the corresponding background category (so that they can be correctly assigned)
\param [in] nBkgndEvents contains the name, yield and option to fix the yield of the background
\param [in] bkgndAsym contains the background asymmetry and option to fix it
*/
virtual void setNBkgndEvents(LauAbsRValue* nBkgndEvents, LauAbsRValue* bkgndAsym);
//! Set the background DP models
/*!
\param [in] bkgndClass the name of the background class
\param [in] negModel the DP model of the B- background
\param [in] posModel the DP model of the B+ background
*/
void setBkgndDPModels(const TString& bkgndClass, LauAbsBkgndDPModel* negModel, LauAbsBkgndDPModel* posModel);
//! Split the signal component into well-reconstructed and mis-reconstructed parts
/*!
The nomenclature used here is TM (truth-matched) and SCF (self cross feed)
In this option, the SCF fraction is DP-dependent
Can also optionally provide a smearing matrix to smear the SCF DP PDF
\param [in] dpHisto the DP histogram of the SCF fraction value
\param [in] upperHalf boolean flag to specify that the supplied histogram contains only the upper half of a symmetric DP (or lower half if using square DP coordinates)
\param [in] fluctuateBins whether the bins on the histogram should be varied in accordance with their uncertainties (for evaluation of systematic uncertainties)
\param [in] scfMap the (optional) smearing matrix
*/
void splitSignalComponent( const TH2* dpHisto, const Bool_t upperHalf = kFALSE, const Bool_t fluctuateBins = kFALSE, LauScfMap* scfMap = 0 );
//! Split the signal component into well reconstructed and mis-reconstructed parts
/*!
The nomenclature used here is TM (truth-matched) and SCF (self cross feed)
In this option, the SCF fraction is a single global number
\param [in] scfFrac the SCF fraction value
\param [in] fixed whether the SCF fraction is fixed or floated in the fit
*/
void splitSignalComponent( const Double_t scfFrac, const Bool_t fixed );
//! Determine whether we are splitting the signal into TM and SCF parts
Bool_t useSCF() const { return useSCF_; }
//! Determine whether the SCF fraction is DP-dependent
Bool_t useSCFHist() const { return useSCFHist_; }
//! Determine if we are smearing the SCF DP PDF
Bool_t smearSCFDP() const { return (scfMap_ != 0); }
// Set the DeltaE and mES models, i.e. give us the PDFs
//! Set the signal PDFs
/*!
\param [in] negPdf the PDF to be added to the B- signal model
\param [in] posPdf the PDF to be added to the B+ signal model
*/
void setSignalPdfs(LauAbsPdf* negPdf, LauAbsPdf* posPdf);
//! Set the SCF PDF for a given variable
/*!
\param [in] negPdf the PDF to be added to the B- signal model
\param [in] posPdf the PDF to be added to the B+ signal model
*/
void setSCFPdfs(LauAbsPdf* negPdf, LauAbsPdf* posPdf);
//! Set the background PDFs
/*!
\param [in] bkgndClass the name of the background class
\param [in] negPdf the PDF to be added to the B- background model
\param [in] posPdf the PDF to be added to the B+ background model
*/
void setBkgndPdfs(const TString& bkgndClass, LauAbsPdf* negPdf, LauAbsPdf* posPdf);
//! Embed full simulation events for the B- signal, rather than generating toy from the PDFs
/*!
\param [in] fileName the name of the file containing the events
\param [in] treeName the name of the tree
\param [in] reuseEventsWithinEnsemble sample with replacement but only replace events once each experiment has been generated
\param [in] reuseEventsWithinExperiment sample with immediate replacement
\param [in] useReweighting perform an accept/reject routine using the configured signal amplitude model based on the MC-truth DP coordinate
*/
void embedNegSignal(const TString& fileName, const TString& treeName,
Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment = kFALSE,
Bool_t useReweighting = kFALSE);
//! Embed full simulation events for the given background class, rather than generating toy from the PDFs
/*!
\param [in] bgClass the name of the background class
\param [in] fileName the name of the file containing the events
\param [in] treeName the name of the tree
\param [in] reuseEventsWithinEnsemble sample with replacement but only replace events once each experiment has been generated
\param [in] reuseEventsWithinExperiment sample with immediate replacement
*/
void embedNegBkgnd(const TString& bgClass, const TString& fileName, const TString& treeName,
Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment = kFALSE);
//! Embed full simulation events for the B+ signal, rather than generating toy from the PDFs
/*!
\param [in] fileName the name of the file containing the events
\param [in] treeName the name of the tree
\param [in] reuseEventsWithinEnsemble sample with replacement but only replace events once each experiment has been generated
\param [in] reuseEventsWithinExperiment sample with immediate replacement
\param [in] useReweighting perform an accept/reject routine using the configured signal amplitude model based on the MC-truth DP coordinate
*/
void embedPosSignal(const TString& fileName, const TString& treeName,
Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment = kFALSE,
Bool_t useReweighting = kFALSE);
//! Embed full simulation events for the given background class, rather than generating toy from the PDFs
/*!
\param [in] bgClass the name of the background class
\param [in] fileName the name of the file containing the events
\param [in] treeName the name of the tree
\param [in] reuseEventsWithinEnsemble sample with replacement but only replace events once each experiment has been generated
\param [in] reuseEventsWithinExperiment sample with immediate replacement
*/
void embedPosBkgnd(const TString& bgClass, const TString& fileName, const TString& treeName,
Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment = kFALSE);
//! Set the DP amplitude coefficients
/*!
The name of the coeffSet must match the name of one of the resonances in the DP model for the antiparticle (the name of the conjugate state in the model for the particle will be automatically determined).
The supplied order of coefficients will be rearranged to match the order in which the resonances are stored in the dynamics, see LauIsobarDynamics::addResonance.
\param [in] coeffSet the set of coefficients
*/
virtual void setAmpCoeffSet(LauAbsCoeffSet* coeffSet);
//! Calculate the DP amplitude(s) for a given DP position
/*!
If not already done, this function will initialise the DP models
\param [in] m13Sq the invariant mass squared of children 1 and 3
\param [in] m23Sq the invariant mass squared of children 2 and 3
\return a container of complex amplitudes, labelled to indicate to which component they belong
*/
virtual std::map<TString, LauComplex> getDPAmps( const Double_t m13Sq, const Double_t m23Sq );
//! Calculate the DP likelihood(s) for a given DP position
/*!
If not already done, this function will initialise the DP models
\param [in] m13Sq the invariant mass squared of children 1 and 3
\param [in] m23Sq the invariant mass squared of children 2 and 3
\return a container of likelihood values, labelled to indicate to which component they belong
*/
virtual std::map<TString, Double_t> getDPLikelihoods( const Double_t m13Sq, const Double_t m23Sq );
protected:
//! Define a map to be used to store a category name and numbers
typedef std::map< std::pair<TString,Int_t>, std::pair<Int_t,Double_t> > LauGenInfo;
//! Typedef for a vector of background DP models
typedef std::vector<LauAbsBkgndDPModel*> LauBkgndDPModelList;
//! Typedef for a vector of background PDFs
typedef std::vector<LauPdfList> LauBkgndPdfsList;
//! Typedef for a vector of background yields
typedef std::vector<LauAbsRValue*> LauBkgndYieldList;
//! Typedef for a vector of embedded data objects
typedef std::vector<LauEmbeddedData*> LauBkgndEmbDataList;
//! Typedef for a vector of booleans to flag if events are reused
typedef std::vector<Bool_t> LauBkgndReuseEventsList;
//! Weight events based on the DP model
/*!
\param [in] dataFileName the name of the data file
\param [in] dataTreeName the name of the data tree
*/
virtual void weightEvents( const TString& dataFileName, const TString& dataTreeName );
//! Initialise the fit
virtual void initialise();
//! Initialise the signal DP models
virtual void initialiseDPModels();
//! Recalculate Normalization the signal DP models
virtual void recalculateNormalisation();
//! Update the coefficients
virtual void updateCoeffs();
//! Toy MC generation and fitting overloaded functions
virtual Bool_t genExpt();
//! Calculate things that depend on the fit parameters after they have been updated by Minuit
virtual void propagateParUpdates();
//! Read in the input fit data variables, e.g. m13Sq and m23Sq
virtual void cacheInputFitVars();
//! Check the initial fit parameters
virtual void checkInitFitParams();
//! Get the fit results and store them
/*!
\param [in] tablePrefixName prefix for the name of the output file
*/
virtual void finaliseFitResults(const TString& tablePrefixName);
//! Print the fit fractions, total DP rate and mean efficiency
/*!
\param [out] output the stream to which to print
*/
virtual void printFitFractions(std::ostream& output);
//! Print the asymmetries
/*!
\param [out] output the stream to which to print
*/
virtual void printAsymmetries(std::ostream& output);
//! Write the fit results in latex table format
/*!
\param [in] outputFile the name of the output file
*/
virtual void writeOutTable(const TString& outputFile);
//! Save the pdf Plots for all the resonances
/*!
\param [in] label prefix for the file name to be saved
*/
virtual void savePDFPlots(const TString& label);
//! Save the pdf Plots for the sum of resonances of a given spin
/*!
\param [in] label prefix for the file name to be saved
\param [in] spin spin of the wave to be saved
*/
virtual void savePDFPlotsWave(const TString& label, const Int_t& spin);
//! Store the per event likelihood values
virtual void storePerEvtLlhds();
// Methods to do with calculating the likelihood functions
// and manipulating the fitting parameters.
//! Get the total likelihood for each event
/*!
\param [in] iEvt the event number
*/
virtual Double_t getTotEvtLikelihood(UInt_t iEvt);
//! Calculate the signal and background likelihoods for the DP for a given event
/*!
\param [in] iEvt the event number
*/
virtual void getEvtDPLikelihood(UInt_t iEvt);
//! Calculate the SCF likelihood for the DP for a given event
/*!
\param [in] iEvt the event number
*/
virtual Double_t getEvtSCFDPLikelihood(UInt_t iEvt);
//! Determine the signal and background likelihood for the extra variables for a given event
/*!
\param [in] iEvt the event number
*/
virtual void getEvtExtraLikelihoods(UInt_t iEvt);
//! Get the total number of events
/*!
\return the total number of events
*/
virtual Double_t getEventSum() const;
//! Set the fit parameters for the DP model
void setSignalDPParameters();
//! Set the fit parameters for the extra PDFs
void setExtraPdfParameters();
//! Set the initial yields
void setFitNEvents();
//! Set-up other parameters that are derived from the fit results, e.g. fit fractions
void setExtraNtupleVars();
//! Randomise the initial fit parameters
void randomiseInitFitPars();
//! Calculate the CP-conserving and CP-violating fit fractions
/*!
\param [in] initValues is this before or after the fit
*/
void calcExtraFractions(Bool_t initValues = kFALSE);
//! Calculate the CP asymmetries
/*!
\param [in] initValues is this before or after the fit
*/
void calcAsymmetries(Bool_t initValues = kFALSE);
//! Define the length of the background vectors
virtual void setupBkgndVectors();
//! Determine the number of events to generate for each hypothesis
std::pair<LauGenInfo,Bool_t> eventsToGenerate();
//! Generate signal event
Bool_t generateSignalEvent();
//! Generate background event
/*!
\param [in] bgID ID number of the background class
*/
Bool_t generateBkgndEvent(UInt_t bgID);
//! Setup the required ntuple branches
void setupGenNtupleBranches();
//! Store all of the DP information
void setDPBranchValues();
//! Generate from the extra PDFs
/*!
\param [in] extraPdfs the list of extra PDFs
\param [in] embeddedData the embedded data sample
*/
void generateExtraPdfValues(LauPdfList* extraPdfs, LauEmbeddedData* embeddedData);
//! Store the MC truth info on the TM/SCF nature of the embedded signal event
/*!
\param [in] embeddedData the full simulation information
*/
Bool_t storeSignalMCMatch(LauEmbeddedData* embeddedData);
//! Add sPlot branches for the extra PDFs
/*!
\param [in] extraPdfs the list of extra PDFs
\param [in] prefix the list of prefixes for the branch names
*/
void addSPlotNtupleBranches(const LauPdfList* extraPdfs, const TString& prefix);
//! Set the branches for the sPlot ntuple with extra PDFs
/*!
\param [in] extraPdfs the list of extra PDFs
\param [in] prefix the list of prefixes for the branch names
\param [in] iEvt the event number
*/
Double_t setSPlotNtupleBranchValues(LauPdfList* extraPdfs, const TString& prefix, UInt_t iEvt);
//! Update the signal events after Minuit sets background parameters
void updateSigEvents();
//! Add branches to store experiment number and the event number within the experiment
virtual void setupSPlotNtupleBranches();
//! Returns the names of all variables in the fit
virtual LauSPlot::NameSet variableNames() const;
//! Returns the names and yields of species that are free in the fit
virtual LauSPlot::NumbMap freeSpeciesNames() const;
//! Returns the names and yields of species that are fixed in the fit
virtual LauSPlot::NumbMap fixdSpeciesNames() const;
//! Returns the species and variables for all 2D PDFs in the fit
virtual LauSPlot::TwoDMap twodimPDFs() const;
//! Check if the signal is split into well-reconstructed and mis-reconstructed types
virtual Bool_t splitSignal() const {return this->useSCF();}
//! Check if the mis-reconstructed signal is to be smeared in the DP
virtual Bool_t scfDPSmear() const {return (scfMap_ != 0);}
//! Append fake data points to the inputData for each bin in the SCF smearing matrix
/*!
We'll be caching the DP amplitudes and efficiencies of the centres of the true bins.
To do so, we attach some fake points at the end of inputData, the number of the entry
minus the total number of events corresponding to the number of the histogram for that
given true bin in the LauScfMap object. (What this means is that when Laura is provided with
the LauScfMap object by the user, it's the latter who has to make sure that it contains the
right number of histograms and in exactly the right order!)
\param [in] inputData the fit data
*/
void appendBinCentres( LauFitDataTree* inputData );
+ //! Retrieve the signal DP model for the B-
LauIsobarDynamics* getNegSigModel() {return negSigModel_;}
+ //! Retrieve the signal DP model for the B+
LauIsobarDynamics* getPosSigModel() {return posSigModel_;}
//! Retrieve a named parameter from a TTree
/*!
\param [in] tree a reference to the tree from which to obtain the parameters
\param [in] name the name of the parameter to retrive
*/
Double_t getParamFromTree( TTree& tree, const TString& name );
//! Set a LauParameter to a given value
/*!
\param [in] param a pointer to the LauParameter to set
\param [in] val the value to set
\param [in] fix whether to fix the LauParameter or leave it floating
*/
void fixParam( LauParameter* param, const Double_t val, const Bool_t fix );
//! Set a vector of LauParameters according to the specified method
/*!
\param [in] params the vector of pointers to LauParameter to set values of
*/
void fixParams( std::vector<LauParameter*>& params );
private:
//! Copy constructor (not implemented)
LauCPFitModel(const LauCPFitModel& rhs);
//! Copy assignment operator (not implemented)
LauCPFitModel& operator=(const LauCPFitModel& rhs);
//! The B- signal Dalitz plot model
LauIsobarDynamics *negSigModel_;
//! The B+ signal Dalitz plot model
LauIsobarDynamics *posSigModel_;
//! The B- background Dalitz plot models
LauBkgndDPModelList negBkgndDPModels_;
//! The B+ background Dalitz plot models
LauBkgndDPModelList posBkgndDPModels_;
//! The B- Dalitz plot kinematics object
LauKinematics *negKinematics_;
//! The B+ Dalitz plot kinematics object
LauKinematics *posKinematics_;
//! The B- signal PDFs
LauPdfList negSignalPdfs_;
//! The B+ signal PDFs
LauPdfList posSignalPdfs_;
//! The B- SCF PDFs
LauPdfList negScfPdfs_;
//! The B+ SCF PDFs
LauPdfList posScfPdfs_;
//! The B- background PDFs
LauBkgndPdfsList negBkgndPdfs_;
//! The B+ background PDFs
LauBkgndPdfsList posBkgndPdfs_;
//! Background boolean
Bool_t usingBkgnd_;
//! Number of signal components
UInt_t nSigComp_;
//! Number of signal DP parameters
UInt_t nSigDPPar_;
//! Number of extra PDF parameters
UInt_t nExtraPdfPar_;
//! Number of normalisation parameters (yields, asymmetries)
UInt_t nNormPar_;
//! Magnitudes and Phases
std::vector<LauAbsCoeffSet*> coeffPars_;
//! The B- fit fractions
LauParArray negFitFrac_;
//! The B+ fit fractions
LauParArray posFitFrac_;
//! Fit B- fractions (uncorrected for the efficiency)
LauParArray negFitFracEffUnCorr_;
//! Fit B+ fractions (uncorrected for the efficiency)
LauParArray posFitFracEffUnCorr_;
//! The CP violating fit fraction
LauParArray CPVFitFrac_;
//! The CP conserving fit fraction
LauParArray CPCFitFrac_;
//! The fit fraction asymmetries
std::vector<LauParameter> fitFracAsymm_;
//! A_CP parameter
std::vector<LauParameter> acp_;
//! The mean efficiency for B- model
LauParameter negMeanEff_;
//! The mean efficiency for B+ model
LauParameter posMeanEff_;
//! The average DP rate for B-
LauParameter negDPRate_;
//! The average DP rate for B+
LauParameter posDPRate_;
//! Signal yield
LauParameter* signalEvents_;
//! Signal asymmetry
LauParameter* signalAsym_;
//! Option to force an asymmetry
Bool_t forceAsym_;
//! Background yield(s)
LauBkgndYieldList bkgndEvents_;
//! Background asymmetries(s)
LauBkgndYieldList bkgndAsym_;
//! IS the analysis tagged?
const Bool_t tagged_;
//! Event charge
const TString tagVarName_;
//! Current event charge
Int_t curEvtCharge_;
//! Vector to store event charges
std::vector<Int_t> evtCharges_;
//! Is the signal split into TM and SCF
Bool_t useSCF_;
//! Is the SCF fraction DP-dependent
Bool_t useSCFHist_;
//! The (global) SCF fraction parameter
LauParameter scfFrac_;
//! The histogram giving the DP-dependence of the SCF fraction
LauEffModel* scfFracHist_;
//! The smearing matrix for the SCF DP PDF
LauScfMap* scfMap_;
//! The cached values of the SCF fraction for each event
std::vector<Double_t> recoSCFFracs_;
//! The cached values of the SCF fraction for each bin centre
std::vector<Double_t> fakeSCFFracs_;
//! The cached values of the sqDP jacobians for each event
std::vector<Double_t> recoJacobians_;
//! The cached values of the sqDP jacobians for each true bin
std::vector<Double_t> fakeJacobians_;
//! Run choice variables
Bool_t compareFitData_;
//! Name of the parent particle
TString negParent_;
//! Name of the parent particle
TString posParent_;
//! The complex coefficients for B-
std::vector<LauComplex> negCoeffs_;
//! The complex coefficients for B+
std::vector<LauComplex> posCoeffs_;
// Embedding full simulation events
//! The B- signal event tree
LauEmbeddedData *negSignalTree_;
//! The B+ signal event tree
LauEmbeddedData *posSignalTree_;
//! The B- background event tree
LauBkgndEmbDataList negBkgndTree_;
//! The B+ background event tree
LauBkgndEmbDataList posBkgndTree_;
//! Boolean to reuse signal events
Bool_t reuseSignal_;
//! Boolean to use reweighting for B-
Bool_t useNegReweighting_;
//! Boolean to use reweighting for B+
Bool_t usePosReweighting_;
//! Vector of booleans to reuse background events
LauBkgndReuseEventsList reuseBkgnd_;
// Likelihood values
//! Signal DP likelihood value
Double_t sigDPLike_;
//! SCF DP likelihood value
Double_t scfDPLike_;
//! Background DP likelihood value(s)
std::vector<Double_t> bkgndDPLike_;
//! Signal likelihood from extra PDFs
Double_t sigExtraLike_;
//! SCF likelihood from extra PDFs
Double_t scfExtraLike_;
//! Background likelihood value(s) from extra PDFs
std::vector<Double_t> bkgndExtraLike_;
//! Total signal likelihood
Double_t sigTotalLike_;
//! Total SCF likelihood
Double_t scfTotalLike_;
//! Total background likelihood(s)
std::vector<Double_t> bkgndTotalLike_;
ClassDef(LauCPFitModel,0) // CP fit/ToyMC model
};
#endif
diff --git a/inc/LauCalcChiSq.hh b/inc/LauCalcChiSq.hh
index 7fec5cd..e6683c6 100644
--- a/inc/LauCalcChiSq.hh
+++ b/inc/LauCalcChiSq.hh
@@ -1,163 +1,163 @@
/*
Copyright 2008 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
#include "TH2Poly.h"
#include "TString.h"
#include <vector>
/*! \file LauCalcChiSq.hh
\brief File containing declaration of LauCalcChiSq class.
*/
/*! \class LauCalcChiSq
\brief Utility class to allow the calculation of the chisq of the fit to the Dalitz plot
A utility class to allow the calculation of the chisq of the fit to the Dalitz plot.
A text config file is provided that gives the datasets for the data and toy
MC generated from the fit results. These can be in the traditional DP or
the square DP.
A sample config file is provided. The fields are:
- <name of file containing data histogram> <name of data tree> <name of the x axis variable in data> <name of the y axis variable in data>
- <name of file containing toy MC histogram> <name of toy MC tree> <name of the x axis variable in toy MC> <name of the y axis variable in toy MC>
- <the minimum content for any bin> <the number of free parameter in the fit> <the scalefactor by which the toy MC bin content should be multiplied> <minimum x value> <maximum x value> <minimum y value> <maximum y value>
+ \<name of file containing data histogram\> \<name of data tree\> \<name of the x axis variable in data\> \<name of the y axis variable in data\><br>
+ \<name of file containing toy MC histogram\> \<name of toy MC tree\> \<name of the x axis variable in toy MC\> \<name of the y axis variable in toy MC\><br>
+ \<the minimum content for any bin\> \<the number of free parameter in the fit\> \<the scalefactor by which the toy MC bin content should be multiplied\> \<minimum x value\> \<maximum x value\> \<minimum y value\> \<maximum y value\>
*/
class LauCalcChiSq {
public:
//! Constructor
/*!
\param [in] inputFileName name of the config file
*/
LauCalcChiSq(const TString& inputFileName = "chiSqInput.txt");
//! Destructor
virtual ~LauCalcChiSq();
//! Toggle verbose printout
/*!
\param [in] flag true to enable verbose printout, false to disable
*/
inline void setVerbose(const Bool_t flag) {verbose_ = flag;}
//! Run the calculations
void run();
private:
//! Read the config file, read the data and create histograms
void initialiseHistos();
//! Choose the binning scheme
/*!
\param [in] xs x coordinates of low statistics sample
\param [in] ys y coordinates of low statistics sample
\param [in] nEntries number of entries in low statistics sample
\param [out] divisions resulting binning scheme
*/
void pickBinning(const Double_t* xs, const Double_t* ys, const Int_t nEntries, std::vector<Int_t>& divisions);
//! Create the template histogram based on the binning scheme
/*!
This function is called recursively to perform subdivisions
\param [in] xMin the minimum x coordinate of the region to be (sub)divided
\param [in] xMax the maximum x coordinate of the region to be (sub)divided
\param [in] yMin the minimum y coordinate of the region to be (sub)divided
\param [in] yMax the maximum y coordinate of the region to be (sub)divided
\param [in] xs x coordinates of low statistics sample
\param [in] ys y coordinates of low statistics sample
\param [in] nEntries number of entries in low statistics sample
\param [in] divisions the binning scheme
\param [in] iter indicates depth of the subdivisions
*/
void getHisto(const Double_t xMin, const Double_t xMax, const Double_t yMin, const Double_t yMax, const Double_t* xs, const Double_t* ys, const Int_t nEntries, const std::vector<Int_t>& divisions, const UInt_t iter=0);
//! Calculate the chisq from the data histograms
void calculateChiSq();
//! Create plots
void makePlots();
//! Name of the config file
TString inputFileName_;
//! Name of the low stats data file
TString fileName1_;
//! Name of the high stats data file
TString fileName2_;
//! Name of the low stats data tree
TString treeName1_;
//! Name of the high stats data tree
TString treeName2_;
//! Name of the x-coordinate branch in tree 1
TString xName1_;
//! Name of the x-coordinate branch in tree 2
TString xName2_;
//! Name of the y-coordinate branch in tree 1
TString yName1_;
//! Name of the y-coordinate branch in tree 2
TString yName2_;
//! The minimum bin content
Float_t minContent_;
//! Template histogram constructed from the binning scheme
TH2Poly* theHisto_;
//! Histogram (constructed from template) filled from tree 1
TH2Poly* histo1_;
//! Histogram (constructed from template) filled from tree 2
TH2Poly* histo2_;
//! Histogram (constructed from template) filled with pulls of tree1 vs tree2
TH2Poly* pullHisto_;
//! Histogram (constructed from template) filled with chisq of tree1 vs tree2
TH2Poly* chiSqHisto_;
//! Histogram (constructed from template) filled with signed chisq of tree1 vs tree2
TH2Poly* chiSqSignedHisto_;
//! Minimum x coordinate of histograms
Float_t xMin_;
//! Maximum x coordinate of histograms
Float_t xMax_;
//! Minimum y coordinate of histograms
Float_t yMin_;
//! Maximum y coordinate of histograms
Float_t yMax_;
//! Number of free parameters in fit (used for calculating the ndof)
Int_t nParams_;
//! Scalefactor between low and high stats data samples
Float_t scaleFactor_;
//! Verbose flag
Bool_t verbose_;
ClassDef(LauCalcChiSq,0)
};
diff --git a/inc/LauComplex.hh b/inc/LauComplex.hh
index 935a283..6b69e8b 100644
--- a/inc/LauComplex.hh
+++ b/inc/LauComplex.hh
@@ -1,357 +1,358 @@
/*
Copyright 2004 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauComplex.hh
\brief File containing declaration of LauComplex class.
*/
/*! \class LauComplex
\brief Class for defining a complex number
Class for complex number manipulation.
In the function descriptions, the form (a,b) is used to represent a complex
number. This is equivalent to the mathematical expression a + ib.
*/
/*****************************************************************************
* Class based on RooFit/RooComplex. *
* Original copyright given below. *
*****************************************************************************
* Authors: *
* WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
* DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
* *
* Copyright (c) 2000-2005, Regents of the University of California *
* and Stanford University. All rights reserved. *
* *
* Redistribution and use in source and binary forms, *
* with or without modification, are permitted according to the terms *
* listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
*****************************************************************************/
#ifndef LAU_COMPLEX
#define LAU_COMPLEX
#include <iosfwd>
#include "Rtypes.h"
#include "TMath.h"
class LauComplex {
public:
//! Default Constructor
inline LauComplex() : re_(0.0), im_(0.0) {}
//! Constructor
/*!
\param [in] a the value corresponding to the real part of the complex number
\param [in] b the value corresponding to the imaginary part of the complex number
*/
inline LauComplex(Double_t a, Double_t b) : re_(a), im_(b) {}
//! Destructor
virtual ~LauComplex() {}
//! Copy constructor
/*!
\param [in] other the complex number to be copied
*/
inline LauComplex(const LauComplex& other) : re_(other.re_), im_(other.im_) {}
//! Copy assignment operator
/*!
\param [in] other the complex number to be copied
\return the assigned complex number
*/
inline LauComplex& operator=(const LauComplex& other)
{
if ( &other != this ) {
re_ = other.re_;
im_ = other.im_;
}
return *this;
}
//! Unary minus operator
/*!
\return the negated complex number
*/
inline LauComplex operator-() const
{
return LauComplex(-re_,-im_);
}
//! Addition operator
/*!
\param [in] other the object to added to this one
\return the sum of the two complex numbers
*/
inline LauComplex operator+(const LauComplex& other) const
{
LauComplex tmpCmplx( *this );
tmpCmplx += other;
return tmpCmplx;
}
//! Subtraction operator
/*!
\param [in] other the object to be subtracted from this one
\return the difference of the two complex numbers
*/
inline LauComplex operator-(const LauComplex& other) const
{
LauComplex tmpCmplx( *this );
tmpCmplx -= other;
return tmpCmplx;
}
//! Multiplication operator
/*!
\param [in] other the object this one is to be multiplied by
\return the product of the two complex numbers
*/
inline LauComplex operator*(const LauComplex& other) const
{
LauComplex tmpCmplx( *this );
tmpCmplx *= other;
return tmpCmplx;
}
//! Division operator
/*!
\param [in] other the object this one is to be divided by
\return the ratio of the two complex numbers
*/
inline LauComplex operator/(const LauComplex& other) const
{
LauComplex tmpCmplx( *this );
tmpCmplx /= other;
return tmpCmplx;
}
//! Addition assignment operator
/*!
\param [in] other the object to be added to this one
\return the result of the addition
*/
inline LauComplex operator+=(const LauComplex& other)
{
this->re_ += other.re_;
this->im_ += other.im_;
return (*this);
}
//! Subtraction assignment operator
/*!
\param [in] other the object to be subtracted from this one
\return the result of the subtraction
*/
inline LauComplex operator-=(const LauComplex& other)
{
this->re_ -= other.re_;
this->im_ -= other.im_;
return (*this);
}
//! Multiplication assignment operator
/*!
\param [in] other the object this one is to be multiplied by
\return the result of the multiplication
*/
inline LauComplex operator*=(const LauComplex& other)
{
Double_t realPart = this->re_*other.re_ - this->im_*other.im_;
Double_t imagPart = this->re_*other.im_ + this->im_*other.re_;
this->setRealImagPart( realPart, imagPart );
return (*this);
}
//! Division assignment operator
/*!
\param [in] other the object this one is to be divided by
\return the results of the division
*/
inline LauComplex operator/=(const LauComplex& other)
{
Double_t x(other.abs2());
Double_t realPart = (this->re_*other.re_ + this->im_*other.im_)/x;
Double_t imagPart = (this->im_*other.re_ - this->re_*other.im_)/x;
this->setRealImagPart( realPart, imagPart );
return (*this);
}
//! Boolean comparison operator
/*!
\param [in] other the object to compared with this one
\return true/false for the comparison
*/
inline Bool_t operator==(const LauComplex& other) const
{
return (re_==other.re_ && im_==other.im_) ;
}
//! Get the real part
/*!
\return the real part of the complex number
*/
inline Double_t re() const
{
return re_;
}
//! Get the imaginary part
/*!
\return the imaginary part of the complex number
*/
inline Double_t im() const
{
return im_;
}
//! Obtain the absolute value of the complex number
/*!
\return the absolute value (magnitude)
*/
inline Double_t abs() const
{
return TMath::Sqrt( this->abs2() );
}
//! Obtain the square of the absolute value of the complex number
/*!
\return the square of the absolute value (magnitude^2)
*/
inline Double_t abs2() const
{
return re_*re_ + im_*im_;
}
//! Obtain the phase angle of the complex number
/*!
\return the phase angle of the complex number
*/
inline Double_t arg() const
{
return TMath::ATan2( im_, re_ );
}
//! Obtain the exponential of the complex number
/*!
\return the exponential of the complex number
*/
inline LauComplex exp() const
{
Double_t mag(TMath::Exp(re_));
return LauComplex(mag*TMath::Cos(im_),mag*TMath::Sin(im_));
}
//! Obtain the complex conjugate
/*!
\return the complex conjugate
*/
inline LauComplex conj() const
{
return LauComplex(re_,-im_);
}
//! Transform this to its complex conjugate
inline void makeConj()
{
im_ = -im_;
}
//! Obtain the complex number scaled by some factor
/*!
\param [in] scaleVal the value used to scale the complex number
\return the complex number scaled by the scaleVal
*/
inline LauComplex scale(Double_t scaleVal) const
{
return LauComplex(scaleVal*re_, scaleVal*im_);
}
//! Scale this by a factor
/*!
\param [in] scaleVal the value used to scale the complex number
*/
inline void rescale(Double_t scaleVal)
{
re_ *= scaleVal;
im_ *= scaleVal;
}
//! Set the real part
/*!
\param [in] realpart the value to be set as the real part
*/
inline void setRealPart(Double_t realpart)
{
re_ = realpart;
}
//! Set the imaginary part
/*!
\param [in] imagpart the value to be set as the imaginary part
*/
inline void setImagPart(Double_t imagpart)
{
im_ = imagpart;
}
//! Set both real and imaginary part
/*!
\param [in] realpart the value to be set as the real part
\param [in] imagpart the value to be set as the imaginary part
*/
inline void setRealImagPart(Double_t realpart, Double_t imagpart)
{
this->setRealPart( realpart );
this->setImagPart( imagpart );
}
//! Set both real and imaginary part to zero
inline void zero()
{
this->setRealImagPart(0.0,0.0);
}
//! Print the complex number
void print() const;
private:
//! The real part
Double_t re_;
//! The imaginary part
Double_t im_;
ClassDef(LauComplex,0) // a non-persistent bare-bones complex class
};
-//! input/output operator formatting of a complex number
+//! input operator formatting of a complex number
std::istream& operator>>(std::istream& os, LauComplex& z);
+//! output operator formatting of a complex number
std::ostream& operator<<(std::ostream& os, const LauComplex& z);
#endif
diff --git a/inc/LauCruijffPdf.hh b/inc/LauCruijffPdf.hh
index f54998b..28b6d6b 100644
--- a/inc/LauCruijffPdf.hh
+++ b/inc/LauCruijffPdf.hh
@@ -1,116 +1,114 @@
/*
Copyright 2008 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauCruijffPdf.hh
\brief File containing declaration of LauCruijffPdf class.
*/
/*! \class LauCruijffPdf
\brief Class for defining a Cruijff PDF.
Class that allows the user to define a Cruijff PDF, a bifurcated Gaussian with asymmetric tails.
The guts of the implementation have been copied from Wouter Hulsbergen's RooFit class.
*/
/*****************************************************************************
* Class based on RooFit/RooCruijff. *
* Original copyright given below. *
*****************************************************************************
* Authors: *
* WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
* DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
* *
* Copyright (c) 2000-2005, Regents of the University of California *
* and Stanford University. All rights reserved. *
* *
* Redistribution and use in source and binary forms, *
* with or without modification, are permitted according to the terms *
* listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
*****************************************************************************/
#ifndef LAU_CRUIJFF_PDF
#define LAU_CRUIJFF_PDF
#include "TString.h"
#include "TRandom.h"
#include "LauAbsPdf.hh"
#include "LauParameter.hh"
#include <vector>
-using std::vector;
-
class LauCruijffPdf : public LauAbsPdf {
public:
//! Constructor
/*!
\param [in] theVarName the name of the abscissa variable
\param [in] params the PDF parameters - mean, sigmaR, sigmaL, alphaR and alphaL
\param [in] minAbscissa the minimum value of the abscissa
\param [in] maxAbscissa the maximum value of the abscissa
*/
- LauCruijffPdf(const TString& theVarName, const vector<LauAbsRValue*>& params, Double_t minAbscissa, Double_t maxAbscissa);
+ LauCruijffPdf(const TString& theVarName, const std::vector<LauAbsRValue*>& params, Double_t minAbscissa, Double_t maxAbscissa);
//! Destructor
virtual ~LauCruijffPdf();
//! 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;
//! Calculate the PDF height
/*!
\param [in] kinematics the current DP kinematics
*/
virtual void calcPDFHeight( const LauKinematics* kinematics );
protected:
private:
//! Copy constructor (not implemented)
LauCruijffPdf(const LauCruijffPdf& other);
//! Copy assignment operator (not implemented)
LauCruijffPdf& operator=(const LauCruijffPdf& other);
//! Gaussian mean
LauAbsRValue* mean_;
//! Sigma of left Gaussian
LauAbsRValue* sigmaL_;
//! Sigma of right Gaussian
LauAbsRValue* sigmaR_;
//! Alpha of left Gaussian
LauAbsRValue* alphaL_;
//! Alpha of right Gaussian
LauAbsRValue* alphaR_;
ClassDef(LauCruijffPdf,0) // Define the Cruijff PDF
};
#endif
diff --git a/inc/LauDPPartialIntegralInfo.hh b/inc/LauDPPartialIntegralInfo.hh
index b50fb8f..2faf60d 100644
--- a/inc/LauDPPartialIntegralInfo.hh
+++ b/inc/LauDPPartialIntegralInfo.hh
@@ -1,275 +1,276 @@
/*
Copyright 2014 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauDPPartialIntegralInfo.hh
\brief File containing declaration of LauDPPartialIntegralInfo class.
*/
/*! \class LauDPPartialIntegralInfo
\brief Class for defining (a section of) the Dalitz plot integration binning scheme
Defines the range and bin size of the integration grid.
Stores the weights and Jacobian for each grid point.
Also stores the amplitude values for each model component at each grid point.
*/
#ifndef LAU_DPPARTIALINTEGRAL_INFO
#define LAU_DPPARTIALINTEGRAL_INFO
#include <iosfwd>
#include "TString.h"
#include "LauComplex.hh"
class LauKinematics;
class LauDPPartialIntegralInfo {
public:
//! Constructor
/*!
\param [in] minm13 the minimum of the m13 range
\param [in] maxm13 the maximum of the m13 range
\param [in] minm23 the minimum of the m23 range
\param [in] maxm23 the maximum of the m23 range
\param [in] m13BinWidth the m13 bin width
\param [in] m23BinWidth the m23 bin width
\param [in] precision the precision required for the Gauss-Legendre weights
\param [in] nAmp the number of coherent amplitude components
\param [in] nIncohAmp the number of incoherent amplitude components
\param [in] squareDP whether or not to use the square DP for the integration - if so, m13 is actually mPrime and m23 is actually thetaPrime
\param [in] kinematics the kinematics object to use to calculate the Jacobians (only relevant if squareDP is true)
*/
LauDPPartialIntegralInfo(const Double_t minm13, const Double_t maxm13,
const Double_t minm23, const Double_t maxm23,
const Double_t m13BinWidth, const Double_t m23BinWidth,
const Double_t precision,
const UInt_t nAmp,
const UInt_t nIncohAmp,
const Bool_t squareDP = kFALSE,
const LauKinematics* kinematics = 0);
//! Destructor
virtual ~LauDPPartialIntegralInfo();
//! Retrieve the minm13 of DP
/*!
\return the the minm13 of DP
*/
inline Double_t getMinm13() const {return minm13_;}
//! Retrieve the maxm13 of DP
/*!
\return the the maxm13 of DP
*/
inline Double_t getMaxm13() const {return maxm13_;}
//! Retrieve the minm23 of DP
/*!
\return the the minm23 of DP
*/
inline Double_t getMinm23() const {return minm23_;}
//! Retrieve the maxm23 of DP
/*!
\return the the maxm23 of DP
*/
inline Double_t getMaxm23() const {return maxm23_;}
//! Retrieve the m13BinWidth of DP
/*!
\return the the m13BinWidth of DP
*/
inline Double_t getM13BinWidth() const {return m13BinWidth_;}
//! Retrieve the m23BinWidth of DP
/*!
\return the the m23BinWidth of DP
*/
inline Double_t getM23BinWidth() const {return m23BinWidth_;}
//! Retrieve the number of bins in m13
/*!
\return the number of bins in m13
*/
inline UInt_t getnm13Points() const {return nm13Points_;}
//! Retrieve the number of bins in m23
/*!
\return the number of bins in m23
*/
inline UInt_t getnm23Points() const {return nm23Points_;}
//! Retrieve the square DP flag
/*!
\return whether or not the integration is performed in the square DP
*/
inline Bool_t getSquareDP() const {return squareDP_;}
//! Retrieve the weight for the given grid point
/*!
\param [in] m13Point the grid index in m13
\param [in] m23Point the grid index in m23
\return the value of the weight
*/
inline Double_t getWeight(const UInt_t m13Point, const UInt_t m23Point) const {return weights_[m13Point][m23Point];}
//! Retrieve the m13 value at the given grid point
/*!
\param [in] m13Point the grid index in m13
\return the m13 value
*/
inline Double_t getM13Value(const UInt_t m13Point) const {return m13Points_[m13Point];}
//! Retrieve the m23 value at the given grid point
/*!
\param [in] m23Point the grid index in m23
\return the m23 value
*/
inline Double_t getM23Value(const UInt_t m23Point) const {return m23Points_[m23Point];}
//! Retrieve the efficiency for the given grid point
/*!
\param [in] m13Point the grid index in m13
\param [in] m23Point the grid index in m23
\return the efficiency value
*/
inline Double_t getEfficiency(const UInt_t m13Point, const UInt_t m23Point) const { return efficiencies_[m13Point][m23Point]; }
//! Store the efficiency for the given grid point
/*!
\param [in] m13Point the grid index in m13
\param [in] m23Point the grid index in m23
\param [in] efficiency the new efficiency value
*/
inline void storeEfficiency(const UInt_t m13Point, const UInt_t m23Point, const Double_t efficiency) { efficiencies_[m13Point][m23Point] = efficiency; }
//! Retrieve the amplitude for the given grid point and amplitude index
/*!
\param [in] m13Point the grid index in m13
\param [in] m23Point the grid index in m23
\param [in] iAmp the amplitude index
\return the amplitude value
*/
inline const LauComplex& getAmplitude(const UInt_t m13Point, const UInt_t m23Point, const UInt_t iAmp) const { return amplitudes_[m13Point][m23Point][iAmp]; }
//! Store the amplitude for the given grid point and amplitude index
/*!
\param [in] m13Point the grid index in m13
\param [in] m23Point the grid index in m23
\param [in] iAmp the amplitude index
\param [in] amplitude the new amplitude value
*/
inline void storeAmplitude(const UInt_t m13Point, const UInt_t m23Point, const UInt_t iAmp, const LauComplex& amplitude) { amplitudes_[m13Point][m23Point][iAmp] = amplitude; }
//! Retrieve the intensity for the given grid point and intensity index
/*!
\param [in] m13Point the grid index in m13
\param [in] m23Point the grid index in m23
\param [in] iAmp the intensity index
\return the intensity value
*/
inline Double_t getIntensity(const UInt_t m13Point, const UInt_t m23Point, const UInt_t iAmp) const { return incohIntensities_[m13Point][m23Point][iAmp]; }
//! Store the intensity for the given grid point and intensity index
/*!
\param [in] m13Point the grid index in m13
\param [in] m23Point the grid index in m23
\param [in] iAmp the intensity index
\param [in] intensity the new intensity value
*/
inline void storeIntensity(const UInt_t m13Point, const UInt_t m23Point, const UInt_t iAmp, const Double_t intensity) { incohIntensities_[m13Point][m23Point][iAmp] = intensity; }
private:
//! Copy constructor (not implemented)
LauDPPartialIntegralInfo( const LauDPPartialIntegralInfo& other );
//! Copy assignment operator (not implemented)
LauDPPartialIntegralInfo& operator=( const LauDPPartialIntegralInfo& other );
//! The minimum of the m13 range
const Double_t minm13_;
//! The maximum of the m13 range
const Double_t maxm13_;
//! The minimum of the m23 range
const Double_t minm23_;
//! The maximum of the m23 range
const Double_t maxm23_;
//! The bin width for m13
const Double_t m13BinWidth_;
//! The bin width for m23
const Double_t m23BinWidth_;
//! The number of bins in m13
const UInt_t nm13Points_;
//! The number of bins in m23
const UInt_t nm23Points_;
//! The number of amplitude components
const UInt_t nAmp_;
//! The number of amplitude components
const UInt_t nIncohAmp_;
//! Flag whether or not we're using the square DP for the integration
const Bool_t squareDP_;
//! The m13 positions of the grid points
std::vector<Double_t> m13Points_;
//! The m23 positions of the grid points
std::vector<Double_t> m23Points_;
//! The Gauss-Legendre weights of the m13 grid points
std::vector<Double_t> m13Weights_;
//! The Gauss-Legendre weights of the m23 grid points
std::vector<Double_t> m23Weights_;
//! The combined weights at each 2D grid point
std::vector< std::vector<Double_t> > weights_;
//! The efficiency at each 2D grid point
std::vector< std::vector<Double_t> > efficiencies_;
//! The amplitude values at each 2D grid point
std::vector< std::vector< std::vector<LauComplex> > > amplitudes_;
//! The incoherent intensity values at each 2D grid point
std::vector< std::vector< std::vector<Double_t> > > incohIntensities_;
ClassDef(LauDPPartialIntegralInfo, 0)
};
+//! output operator formatting
std::ostream& operator<<( std::ostream& stream, const LauDPPartialIntegralInfo& infoRecord );
#endif
diff --git a/inc/LauKinematics.hh b/inc/LauKinematics.hh
index 1ab595b..aad5d56 100644
--- a/inc/LauKinematics.hh
+++ b/inc/LauKinematics.hh
@@ -1,675 +1,673 @@
/*
Copyright 2004 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauKinematics.hh
\brief File containing declaration of LauKinematics class.
*/
/*! \class LauKinematics
\brief Class for calculating 3-body kinematic quantities.
Class for defining the many routines related to the three body kinematics.
Given the two DP co-ordinates, all other useful quantities e.g. the helicity angles are calculated.
Optionally, the so-called ``square Dalitz plot'' quantities can also be calculated.
*/
#ifndef LAU_KINEMATICS
#define LAU_KINEMATICS
#include <vector>
#include "Rtypes.h"
#include "TMath.h"
class LauKinematics {
public:
//! Constructor
/*!
\param [in] m1 the first daughter mass
\param [in] m2 the second daughter mass
\param [in] m3 the third daughter mass
\param [in] mParent the parent particle mass
\param [in] calcSquareDPCoords boolean flag to enable/disable calculation of the square Dalitz plot co-ordinates
\param [in] symmetricalDP boolean flag to indicate whether the DP is symmetric (i.e. two identical particle in final state)
\param [in] fullySymmetricDP boolean flag to indicate whether the DP is fully symmetric (i.e. all three final-state particles are identical)
*/
LauKinematics(const Double_t m1, const Double_t m2, const Double_t m3, const Double_t mParent,
const Bool_t calcSquareDPCoords = kFALSE, const Bool_t symmetricalDP = kFALSE, const Bool_t fullySymmetricDP = kFALSE);
//! Destructor
virtual ~LauKinematics();
//! Enable/disable the calculation of square Dalitz plot co-ordinates
/*!
\param [in] calcSquareDPCoords kTRUE/kFALSE to enable/disable calculation of the square DP co-ordinates
*/
inline void squareDP( const Bool_t calcSquareDPCoords ) { squareDP_ = calcSquareDPCoords; }
//! Are the square Dalitz plot co-ordinates being calculated?
/*!
\return kTRUE if the square Dalitz plot co-ordinates are being calculated, kFALSE otherwise
*/
inline Bool_t squareDP() const { return squareDP_; }
//! Is the DP symmetric?
/*!
\return kTRUE if the DP is symmetric (i.e. daughters 1 and 2 are identical), kFALSE otherwise
*/
inline Bool_t gotSymmetricalDP() const { return (symmetricalDP_ && ! fullySymmetricDP_); }
//! Is the DP fully symmetric?
/*!
\return kTRUE if the DP is fully symmetric (i.e. daughters 1, 2 and 3 are identical), kFALSE otherwise
*/
inline Bool_t gotFullySymmetricDP() const { return fullySymmetricDP_; }
//! Enable/disable warning messages
inline void warningMessages(const Bool_t boolean) { warnings_ = boolean; }
//! Update all kinematic quantities based on the DP co-ordinates m13Sq and m23Sq
/*!
It can be useful to first check that the point is within the kinematic boundary (using LauKinematics::withinDPLimits) before calling this method.
\param [in] m13Sq the invariant mass squared of daughters 1 and 3
\param [in] m23Sq the invariant mass squared of daughters 2 and 3
*/
void updateKinematics(const Double_t m13Sq, const Double_t m23Sq);
//! Update all kinematic quantities based on the square DP co-ordinates m' and theta'
/*!
It can be useful to first check that the point is within the kinematic boundary (using LauKinematics::withinSqDPLimits) before calling this method.
\param [in] mPrime the m' co-ordinate
\param [in] thetaPrime the theta' co-ordinate
*/
void updateSqDPKinematics(const Double_t mPrime, const Double_t thetaPrime);
//! Update all kinematic quantities based on m23 and cos(theta23)
/*!
\param [in] m23 the invariant mass of daughters 2 and 3
\param [in] c23 the cosine of the helicity angle theta23, \see getc23
*/
void updateKinematicsFrom23(const Double_t m23, const Double_t c23);
//! Update all kinematic quantities based on m13 and cos(theta13)
/*!
\param [in] m13 the invariant mass of daughters 1 and 3
\param [in] c13 the cosine of the helicity angle theta13, \see getc13
*/
void updateKinematicsFrom13(const Double_t m13, const Double_t c13);
//! Update all kinematic quantities based on m12 and cos(theta12)
/*!
\param [in] m12 the invariant mass of daughters 1 and 2
\param [in] c12 the cosine of the helicity angle theta12, \see getc12
*/
void updateKinematicsFrom12(const Double_t m12, const Double_t c12);
//! Calculate the Jacobian for the transformation m23^2, m13^2 -> m', theta' (square DP) at the given point in the square DP
/*!
\param [in] mPrime the m' co-ordinate
\param [in] thPrime the theta' co-ordinate
\return the jacobian of the transformation
*/
Double_t calcSqDPJacobian(const Double_t mPrime, const Double_t thPrime) const;
//! Calculate the Jacobian for the transformation m23^2, m13^2 -> m', theta' (square DP) at the currently stored point in the square DP
/*!
\return the jacobian of the transformation
*/
Double_t calcSqDPJacobian() const;
//! Routine to generate events flat in phase-space
/*!
\param [out] m13Sq the invariant mass squared of daughters 1 and 3
\param [out] m23Sq the invariant mass squared of daughters 2 and 3
*/
void genFlatPhaseSpace(Double_t& m13Sq, Double_t& m23Sq) const;
//! Routine to generate events flat in the square Dalitz plot
/*!
\param [out] mPrime the m' variable
\param [out] thetaPrime the theta' variable
*/
void genFlatSqDP(Double_t& mPrime, Double_t& thetaPrime) const;
//! Check whether a given (m13Sq,m23Sq) point is within the kinematic limits of the Dalitz plot
/*!
This method first checks that m13Sq is within its absolute
min and max and then for the given m13Sq calculates the
local min and max of m23Sq and checks whether the given
m23Sq satisfies these bounds.
\param [in] m13Sq the invariant mass squared of daughters 1 and 3
\param [in] m23Sq the invariant mass squared of daughters 2 and 3
\return kTRUE if the event is inside the kinematic limit, kFALSE otherwise
*/
Bool_t withinDPLimits(const Double_t m13Sq, const Double_t m23Sq) const;
//! Check whether a given (m13Sq,m23Sq) point is within the kinematic limits of the Dalitz plot (alternate method)
/*!
This method first checks that m23Sq is within its absolute
min and max and then for the given m23Sq calculates the
local min and max of m13Sq and checks whether the given
m13Sq satisfies these bounds.
\param [in] m13Sq the m13 invariant mass pair squared
\param [in] m23Sq the m23 invariant mass pair squared
\return kTRUE if the event is inside the kinematic limit, kFALSE otherwise
*/
Bool_t withinDPLimits2(const Double_t m13Sq, const Double_t m23Sq) const;
//! Check whether a given (m',theta') point is within the kinematic limits of the Dalitz plot
/*!
\param [in] mPrime the m' co-ordinate
\param [in] thetaPrime the theta' co-ordinate
\return kTRUE if the event is inside the kinematic limit, kFALSE otherwise
*/
Bool_t withinSqDPLimits(const Double_t mPrime, const Double_t thetaPrime) const;
//! Calculate the third invariant mass square from the two provided (e.g. mjkSq from mijSq and mikSq)
/*!
\param [in] firstMassSq the first invariant mass squared
\param [in] secondMassSq the second invariant mass squared
\return the third invariant mass square
*/
Double_t calcThirdMassSq(const Double_t firstMassSq, const Double_t secondMassSq) const;
//! Calculate the distance from the currently set (m13Sq, m23Sq) point to the centre of the Dalitz plot (which is defined as the point where m12=m13=m23)
/*!
\return the distance to the DP centre
*/
Double_t distanceFromDPCentre() const;
//! Calculate the distance from a given (m13Sq, m23Sq) point to the centre of the Dalitz plot (which is defined as the point where m12=m13=m23)
/*!
\return the distance to the DP centre
*/
Double_t distanceFromDPCentre(const Double_t m13Sq, const Double_t m23Sq) const;
//! Get the m12 invariant mass
/*!
\return the m12 invariant mass
*/
inline Double_t getm12() const {return m12_;}
//! Get the m23 invariant mass
/*!
\return the m23 invariant mass
*/
inline Double_t getm23() const {return m23_;}
//! Get the m13 invariant mass
/*!
\return the m13 invariant mass
*/
inline Double_t getm13() const {return m13_;}
//! Get the m12 invariant mass square
/*!
\return the m12 invariant mass square
*/
inline Double_t getm12Sq() const {return m12Sq_;}
//! Get the m23 invariant mass square
/*!
\return the m23 invariant mass square
*/
inline Double_t getm23Sq() const {return m23Sq_;}
//! Get the m13 invariant mass square
/*!
\return the m13 invariant mass square
*/
inline Double_t getm13Sq() const {return m13Sq_;}
//! Get the cosine of the helicity angle theta12
/*!
theta12 is defined as the angle between 1&3 in the rest frame of 1&2
\return the cosine of the helicity angle theta12
*/
inline Double_t getc12() const {return c12_;}
//! Get the cosine of the helicity angle theta23
/*!
theta23 is defined as the angle between 3&1 in the rest frame of 2&3
\return the cosine of the helicity angle theta23
*/
inline Double_t getc23() const {return c23_;}
//! Get the cosine of the helicity angle theta13
/*!
theta13 is defined as the angle between 3&2 in the rest frame of 1&3
\return the cosine of the helicity angle theta13
*/
inline Double_t getc13() const {return c13_;}
//! Get m' value
/*!
\return m' value
*/
inline Double_t getmPrime() const {return mPrime_;}
//! Get theta' value
/*!
\return theta' value
*/
inline Double_t getThetaPrime() const {return thetaPrime_;}
//! Get parent mass
/*!
\return parent mass
*/
inline Double_t getmParent() const {return mParent_;}
//! Get parent mass squared
/*!
\return parent mass squared
*/
inline Double_t getmParentSq() const {return mParentSq_;}
//! Get the box area defined from the kinematic bounds
/*!
The box area is defined as:\n
[(M-m3)^2 - (m1+m2)^2]*[(M-m2)^2 - (m1+m3)^2] .:. (m13SqMax - m13SqMin)*(m23SqMax - m23SqMin)
\return the Dalitz plot box area
*/
inline Double_t getDPBoxArea() const {return (mSqMax_[1] - mSqMin_[1])*(mSqMax_[0] - mSqMin_[0]);}
//! Flips the DP variables m13^2 <-> m23^2 and recalculates all kinematic quantities
/*!
Useful in the case of symmetrical Dalitz plots, i.e. when two final state particles are identical
*/
void flipAndUpdateKinematics();
//! Cyclically rotates the DP variables (m12 -> m23, m23 -> m13, m13 -> m12) and recalculates all kinematic quantities
/*!
Useful in the case of a fully symmetric Dalitz plot, i.e. when all three final state particles are identical
*/
void rotateAndUpdateKinematics();
//! Get the m1 mass
/*!
\return the m1 mass
*/
inline Double_t getm1() const {return m1_;}
//! Get the m2 mass
/*!
\return the m2 mass
*/
inline Double_t getm2() const {return m2_;}
//! Get the m3 mass
/*!
\return the m3 mass
*/
inline Double_t getm3() const {return m3_;}
//! Get the m23 minimum defined as (m2 + m3)
/*!
\return the minimum value for m23
*/
inline Double_t getm23Min() const {return TMath::Sqrt(mSqMin_[0]);}
//! Get the m13 minimum defined as (m1 + m3)
/*!
\return the minimum value for m13
*/
inline Double_t getm13Min() const {return TMath::Sqrt(mSqMin_[1]);}
//! Get the m12 minimum defined as (m1 + m2)
/*!
\return the minimum value for m12
*/
inline Double_t getm12Min() const {return TMath::Sqrt(mSqMin_[2]);}
//! Get the m23 maximum defined as (mParent - m1)
/*!
\return the maximum value for m23
*/
inline Double_t getm23Max() const {return TMath::Sqrt(mSqMax_[0]);}
//! Get the m13 maximum defined as (mParent - m2)
/*!
\return the maximum value for m13
*/
inline Double_t getm13Max() const {return TMath::Sqrt(mSqMax_[1]);}
//! Get the m12 maximum defined as (mParent - m3)
/*!
\return the maximum value for m12
*/
inline Double_t getm12Max() const {return TMath::Sqrt(mSqMax_[2]);}
//! Get the m23Sq minimum, (m2 + m3)^2
/*!
\return the minimum value for m23Sq
*/
inline Double_t getm23SqMin() const {return mSqMin_[0];}
//! Get the m13Sq minimum, (m1 + m3)^2
/*!
\return the minimum value for m13Sq
*/
inline Double_t getm13SqMin() const {return mSqMin_[1];}
//! Get the m12Sq minimum, (m1 + m2)^2
/*!
\return the minimum value for m12Sq
*/
inline Double_t getm12SqMin() const {return mSqMin_[2];}
//! Get the m23Sq maximum, (mParent - m1)^2
/*!
\return the maximum value for m23Sq
*/
inline Double_t getm23SqMax() const {return mSqMax_[0];}
//! Get the m13Sq maximum, (mParent - m2)^2
/*!
\return the maximum value for m13Sq
*/
inline Double_t getm13SqMax() const {return mSqMax_[1];}
//! Get the m12Sq maximum, (mParent - m3)^2
/*!
\return the maximum value for m12Sq
*/
inline Double_t getm12SqMax() const {return mSqMax_[2];}
//! Get the momentum of the track 1 in 12 rest frame
/*!
\return the momentum of track 1 in 12 rest frame
*/
inline Double_t getp1_12() const {return p1_12_;}
//! Get the momentum of the track 3 in 12 rest frame
/*!
\return the momentum of track 3 in 12 rest frame
*/
inline Double_t getp3_12() const {return p3_12_;}
//! Get the momentum of the track 2 in 23 rest frame
/*!
\return the momentum of track 2 in 23 rest frame
*/
inline Double_t getp2_23() const {return p2_23_;}
//! Get the momentum of the track 1 in 23 rest frame
/*!
\return the momentum of track 1 in 23 rest frame
*/
inline Double_t getp1_23() const {return p1_23_;}
//! Get the momentum of the track 1 in 13 rest frame
/*!
\return the momentum of track 1 in 13 rest frame
*/
inline Double_t getp1_13() const {return p1_13_;}
//! Get the momentum of the track 2 in 13 rest frame
/*!
\return the momentum of track 2 in 13 rest frame
*/
inline Double_t getp2_13() const {return p2_13_;}
//! Get the momentum of the track 1 in parent rest frame
/*!
\return the momentum of track 1 in parent rest frame
*/
inline Double_t getp1_Parent() const {return p1_Parent_;}
//! Get the momentum of the track 2 in parent rest frame
/*!
\return the momentum of track 2 in parent rest frame
*/
inline Double_t getp2_Parent() const {return p2_Parent_;}
//! Get the momentum of the track 3 in parent rest frame
/*!
\return the momentum of track 3 in parent rest frame
*/
inline Double_t getp3_Parent() const {return p3_Parent_;}
//! Method to draw the Dalitz plot contours on the top of the histo previously drawn
/*!
\param [in] orientation orientation used for the draw, with default set to 1323 that corresponds x = m13, y = m23
\param [in] nbins number of bins for the contour to be sampled over (default = 100)
*/
void drawDPContour(const Int_t orientation = 1323, const Int_t nbins = 100);
//! Get covariant factor in 12 axis
/*!
\return covariant factor in 12 axis
*/
inline Double_t getcov12() const {return (mParentSq_ + m12Sq_ - m3Sq_)/(2.*mParent_*m12_);}
//! Get covariant factor in 13 axis
/*!
\return covariant factor in 13 axis
*/
inline Double_t getcov13() const {return (mParentSq_ + m13Sq_ - m2Sq_)/(2.*mParent_*m13_);}
//! Get covariant factor in 23 axis
/*!
\return covariant factor in 23 axis
*/
inline Double_t getcov23() const {return (mParentSq_ + m23Sq_ - m1Sq_)/(2.*mParent_*m23_);}
protected:
//! Update the variables m23Sq_ and m13Sq_ given the invariant mass m12 and the cosine of the helicity angle c12
/*!
\param [in] m12 the invariant mass m12
\param [in] c12 the cosine of the helicity angle c12
*/
void updateMassSq_m12(const Double_t m12, const Double_t c12);
//! Update the variables m12Sq_ and m13Sq_ given the invariant mass m23 and the cosine of the helicity angle c23
/*!
\param [in] m23 the invariant mass m12
\param [in] c23 the cosine of the helicity angle c23
*/
void updateMassSq_m23(const Double_t m23, const Double_t c23);
//! Update the variables m12Sq_ and m23Sq_ given the invariant mass m13 and the cosine of the helicity angle c13
/*!
\param [in] m13 the invariant mass m13
\param [in] c13 the cosine of the helicity angle c13
*/
void updateMassSq_m13(const Double_t m13, const Double_t c13);
//! Update some kinematic quantities based on the DP co-ordinates m13Sq and m23Sq
/*!
Only the three invariant masses and their squares, plus the parent rest-frame momenta are updated.
\param [in] m13Sq the invariant mass squared of daughters 1 and 3
\param [in] m23Sq the invariant mass squared of daughters 2 and 3
*/
void updateMassSquares(const Double_t m13Sq, const Double_t m23Sq);
//! Update some kinematic quantities based on the square DP co-ordinates m' and theta'
/*!
Only m', theta', the three invariant masses and their squares, plus the parent rest-frame momenta are updated.
\param [in] mPrime the m' co-ordinate
\param [in] thetaPrime the theta' co-ordinate
*/
void updateSqDPMassSquares(const Double_t mPrime, const Double_t thetaPrime);
//! General method to calculate the cos(helicity) variables from the masses of the particles
/*!
\param [in] mijSq the mij invariant mass square
\param [in] mikSq the mik invariant mass square
\param [in] mij the mij invariant mass
\param [in] i the index for the first track
\param [in] j the index for the second track
\param [in] k the index for the third track
\return helicity angle in the ij rest frame
*/
Double_t cFromM(const Double_t mijSq, const Double_t mikSq, const Double_t mij, const Int_t i, const Int_t j, const Int_t k) const;
//! General method to calculate mikSq given mijSq and cosHel_ij
/*!
\param [in] mijSq the mij invariant mass square
\param [in] cij the helicity angle for the pair which is made from tracks i and j
\param [in] mij the mij invariant mass
\param [in] i the index for the first track
\param [in] j the index for the second track
\param [in] k the index for the third track
\return the invariant mass square mikSq
*/
Double_t mFromC(const Double_t mijSq, const Double_t cij, const Double_t mij, const Int_t i, const Int_t j, const Int_t k) const;
//! General method to calculate the momentum of a particle, given its energy and invariant mass squared.
/*!
\param [in] energy the energy of the particle
\param [in] massSq the invariant mass squared of the particle
\return the momentum of the particle
*/
Double_t pCalc(const Double_t energy, const Double_t massSq) const;
//! Randomly generate the invariant mass squared m13Sq
/*!
\return the invariant mass squared m13Sq
*/
Double_t genm13Sq() const;
//! Randomly generate the invariant mass squared m23Sq
/*!
\return the invariant mass squared m23Sq
*/
Double_t genm23Sq() const;
//! Randomly generate the invariant mass squared m12Sq
/*!
\return the invariant mass squared m12Sq
*/
Double_t genm12Sq() const;
//! Calculate m12Sq from m13Sq and m23Sq
void calcm12Sq();
//! Calculate cosines of the helicity angles, momenta of daughters and bachelor in various ij rest frames
void calcHelicities();
//! Calculate the m' and theta' variables for the square Dalitz plot
void calcSqDPVars();
//! Calculate the momenta of each daughter in the parent rest frame
void calcParentFrameMomenta();
private:
//! Copy constructor (not implemented)
LauKinematics(const LauKinematics& rhs);
//! Copy assignment operator (not implemented)
LauKinematics& operator=(const LauKinematics& rhs);
//! Symmetrical DP
const Bool_t symmetricalDP_;
//! Fully-symmetrical DP
const Bool_t fullySymmetricDP_;
//! Mass of particle 1
const Double_t m1_;
//! Mass of particle 2
const Double_t m2_;
//! Mass of particle 3
const Double_t m3_;
//! Mass of parent particle
const Double_t mParent_;
//! Mass of particle 1 squared
const Double_t m1Sq_;
//! Mass of particle 2 squared
const Double_t m2Sq_;
//! Mass of particle 3 squared
const Double_t m3Sq_;
//! Mass of parent particle squared
const Double_t mParentSq_;
//! Vector of daughter particles masses
std::vector<Double_t> mass_;
//! Vector of the minimum mij values
std::vector<Double_t> mMin_;
//! Vector of the maximum mij values
std::vector<Double_t> mMax_;
//! Vector of the difference between the mMax and mMin
std::vector<Double_t> mDiff_;
//! Vector of daughter particles masses squared
std::vector<Double_t> mSq_;
//! Vector of the minimum mijSq values
std::vector<Double_t> mSqMin_;
//! Vector of the maximum mijSq values
std::vector<Double_t> mSqMax_;
//! Vector of the difference between the mSqMax and mSqMin
std::vector<Double_t> mSqDiff_;
//! Sum of the daughter masses
const Double_t mDTot_;
//! Mass difference between the parent particle and the sum of the daughter particles
const Double_t massInt_;
//! Sum of the squares of the daughter masses
const Double_t mSqDTot_;
//! Invariant mass m12
Double_t m12_;
//! Invariant mass m23
Double_t m23_;
//! Invariant mass m13
Double_t m13_;
//! Invariant mass m12 squared
Double_t m12Sq_;
//! Invariant mass m23 squared
Double_t m23Sq_;
//! Invariant mass m13 squared
Double_t m13Sq_;
//! Cosine of the helicity angle theta12, which is defined as the angle between 1&3 in the rest frame of 1&2
Double_t c12_;
//! Cosine of the helicity angle theta23, which is defined as the angle between 1&2 in the rest frame of 2&3
Double_t c23_;
//! Cosine of the helicity angle theta13, which is defined as the angle between 1&2 in the rest frame of 1&3
Double_t c13_;
//! m' co-ordinate
Double_t mPrime_;
//! theta' co-ordinate
Double_t thetaPrime_;
//! Momentum q of particle i
mutable Double_t qi_;
//! Momentum q of particle k
mutable Double_t qk_;
//! Momentum of track 1 in 1-2 rest frame
Double_t p1_12_;
//! Momentum of track 3 in 1-2 rest frame
Double_t p3_12_;
//! Momentum of track 2 in 2-3 rest frame
Double_t p2_23_;
//! Momentum of track 1 in 2-3 rest frame
Double_t p1_23_;
//! Momentum of track 1 in 1-3 rest frame
Double_t p1_13_;
//! Momentum of track 2 in 1-3 rest frame
Double_t p2_13_;
//! Momentum of track 1 in parent rest frame
Double_t p1_Parent_;
//! Momentum of track 2 in parent rest frame
Double_t p2_Parent_;
//! Momentum of track 3 in parent rest frame
Double_t p3_Parent_;
//! Should we calculate the square DP co-ordinates or not?
Bool_t squareDP_;
//! Enable/disable warning messages
Bool_t warnings_;
ClassDef(LauKinematics,0)
};
-Double_t dal(Double_t* x, Double_t* par);
-
#endif
diff --git a/inc/LauLinearPdf.hh b/inc/LauLinearPdf.hh
index 242bc49..5f3a1e6 100644
--- a/inc/LauLinearPdf.hh
+++ b/inc/LauLinearPdf.hh
@@ -1,94 +1,95 @@
/*
Copyright 2006 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauLinearPdf.hh
\brief File containing declaration of LauLinearPdf class.
*/
/*! \class LauLinearPdf
\brief Class for defining a straight line PDF.
Class that allows the user to define a straight line PDF.
*/
#ifndef LAU_LINEAR_PDF
#define LAU_LINEAR_PDF
#include <vector>
#include "TString.h"
#include "LauAbsPdf.hh"
class LauLinearPdf : public LauAbsPdf {
public:
//! Constructor
/*!
\param [in] theVarName the name of the abscissa variable
\param [in] params the PDF parameter - slope
\param [in] minAbscissa the minimum value of the abscissa
\param [in] maxAbscissa the maximum value of the abscissa
*/
LauLinearPdf(const TString& theVarName, const std::vector<LauAbsRValue*>& params, Double_t minAbscissa, Double_t maxAbscissa);
//! Destructor
virtual ~LauLinearPdf();
//! 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;
//! Calculate the normalisation
virtual void calcNorm();
//! Calculate the PDF height
/*!
\param [in] kinematics the current DP kinematics
*/
virtual void calcPDFHeight( const LauKinematics* kinematics );
protected:
private:
//! Copy constructor (not implemented)
LauLinearPdf(const LauLinearPdf& rhs);
//! Copy assignment operator (not implemented)
LauLinearPdf& operator=(const LauLinearPdf& rhs);
//! Line slope
LauAbsRValue* slope_;
+ //! Flag to control printing of warnings about the PDF going negative
Bool_t posflag_;
ClassDef(LauLinearPdf,0) // Define the Linear PDF
};
#endif
diff --git a/inc/LauMergeDataFiles.hh b/inc/LauMergeDataFiles.hh
index 940d6e0..6c49f20 100644
--- a/inc/LauMergeDataFiles.hh
+++ b/inc/LauMergeDataFiles.hh
@@ -1,127 +1,127 @@
/*
Copyright 2008 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
#include <map>
#include <cstdlib>
#include "TFile.h"
#include "TString.h"
#include "TTree.h"
/*! \file LauMergeDataFiles.hh
\brief File containing declaration of LauMergeDataFiles class.
*/
-/*! \class LauResultsExtractor
+/*! \class LauMergeDataFiles
\brief Utility class to allow the merging of data files on a expt-by-expt basis
The files are merged such that events for expt 0 from tree 1 will be followed
by events for expt 0 from tree 2, then expt 1 from tree1, expt 1 from tree 2, etc.
*/
class LauMergeDataFiles
{
public:
//! Constructor
/*!
\param [in] fileName1 name of first file to be merged
\param [in] fileName2 name of second file to be merged
\param [in] treeName name of the tree to read from the input files
*/
LauMergeDataFiles(const TString& fileName1, const TString& fileName2, const TString& treeName);
//! Destructor
virtual ~LauMergeDataFiles();
//! Do the merge
/*!
\param [in] fileName name of file to which to write the merged tree
*/
void process(const TString& fileName);
protected:
//! Type to relate leaf names with their double-precision value
typedef std::map<TString,Double_t> LeafDoubleMap;
//! Type to relate leaf names with their integer value
typedef std::map<TString,Int_t> LeafIntegerMap;
//! Type to hold for each experiment the first and last entry numbers in a tree
typedef std::map< Int_t,std::pair<Int_t,Int_t> > ExptsMap;
//! Open the specified input files and check that the trees can be read
void openInputFiles();
//! Read the structure of the input trees, create appropriate storage and set the branch addresses
void setupInputTrees();
//! Create the structure of the output tree
void setupOutputTree();
//! Determine the experiments stored a given tree
void findExperiments(TTree* tree, ExptsMap& exptsMap);
//! Check that the experiments in each tree match
Bool_t checkExperimentMaps() const;
//! Read the entries for a given experiment from the given tree and store in the output tree
void readExperiment(TTree* tree, const ExptsMap::const_iterator& exptsMap, Int_t offset);
//! Write the output file
void writeFile();
private:
//! Name of file 1
TString fileName1_;
//! Name of file 2
TString fileName2_;
//! Name of the tree
TString treeName_;
//! Input file 1
TFile * inputFile1_;
//! Input file 2
TFile * inputFile2_;
//! Input tree 1
TTree * inputTree1_;
//! Input tree 2
TTree * inputTree2_;
//! Output file
TFile * outputFile_;
//! Output tree
TTree * outputTree_;
// Tree variables
//! Storage for the experiment index variable
Int_t iExpt_;
//! Storage for the event-within-experiment index variable
Int_t iEvtWithinExpt_;
//! Storage for double-precision leaves
LeafDoubleMap doubleVars_;
//! Storage for integer leaves
LeafIntegerMap integerVars_;
//! Experiment -> first and last tree entry for tree 1
ExptsMap tree1Expts_;
//! Experiment -> first and last tree entry for tree 2
ExptsMap tree2Expts_;
ClassDef(LauMergeDataFiles,0)
};
diff --git a/inc/LauNovosibirskPdf.hh b/inc/LauNovosibirskPdf.hh
index 0e88b81..6ed8649 100644
--- a/inc/LauNovosibirskPdf.hh
+++ b/inc/LauNovosibirskPdf.hh
@@ -1,95 +1,93 @@
/*
Copyright 2008 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauNovosibirskPdf.hh
\brief File containing declaration of LauNovosibirskPdf class.
*/
/*! \class LauNovosibirskPdf
\brief Class for defining a Novosibirsk function PDF.
Class that allows the user to define a Novosibirsk function PDF.
*/
#ifndef LAU_NOVOSIBIRSK_PDF
#define LAU_NOVOSIBIRSK_PDF
#include "TString.h"
#include "TRandom.h"
#include "LauAbsPdf.hh"
#include "LauParameter.hh"
#include <vector>
-using std::vector;
-
class LauNovosibirskPdf : public LauAbsPdf {
public:
//! Constructor
/*!
\param [in] theVarName the name of the abscissa variable
\param [in] params the PDF parameters - mean, sigma and tail
\param [in] minAbscissa the minimum value of the abscissa
\param [in] maxAbscissa the maximum value of the abscissa
*/
- LauNovosibirskPdf(const TString& theVarName, const vector<LauAbsRValue*>& params, Double_t minAbscissa, Double_t maxAbscissa);
+ LauNovosibirskPdf(const TString& theVarName, const std::vector<LauAbsRValue*>& params, Double_t minAbscissa, Double_t maxAbscissa);
//! Destructor
virtual ~LauNovosibirskPdf();
//! 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;
//! Calculate the PDF height
/*!
\param [in] kinematics the current DP kinematics
*/
virtual void calcPDFHeight( const LauKinematics* kinematics );
protected:
private:
//! Copy constructor (not implemented)
LauNovosibirskPdf(const LauNovosibirskPdf& rhs);
//! Copy assignment operator (not implemented)
LauNovosibirskPdf& operator=(const LauNovosibirskPdf& rhs);
//! Gaussian mean
LauAbsRValue* mean_;
//! Gaussian sigma
LauAbsRValue* sigma_;
//! Gaussian tail
LauAbsRValue* tail_;
ClassDef(LauNovosibirskPdf,0) // Define the Novosibirsk PDF
};
#endif
diff --git a/inc/LauPolarFormFactorNR.hh b/inc/LauPolarFormFactorNR.hh
index a4d248f..568da95 100644
--- a/inc/LauPolarFormFactorNR.hh
+++ b/inc/LauPolarFormFactorNR.hh
@@ -1,141 +1,142 @@
/*
Copyright 2018 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauPolarFormFactorNR.hh
\brief File containing declaration of LauPolarFormFactorNR class.
*/
/*! \class LauPolarFormFactorNR
\brief Class for defining a nonresonant form factor model
Defines the nonresonant form factor model from:
Nogueira, Bediaga, Cavalcante, Frederico, Lorenco: Phys. Rev. D92 (2015) 054010, arXiv:1506.08332 [hep-ph]
Pelaez, Yndurain: Phys. Rev. D71 (2005) 074016, arXiv:hep-ph/0411334
*/
#ifndef LAU_POLAR_FORM_FACTOR_NR
#define LAU_POLAR_FORM_FACTOR_NR
#include "TString.h"
#include "LauComplex.hh"
#include "LauAbsResonance.hh"
class LauKinematics;
class LauParameter;
class LauPolarFormFactorNR : public LauAbsResonance {
public:
//! Constructor
/*!
\param [in] resInfo the object containing information on the resonance name, mass, width, spin, charge, etc.
\param [in] resType the model of the resonance
\param [in] resPairAmpInt the number of the daughter not produced by the resonance
\param [in] daughters the daughter particles
*/
LauPolarFormFactorNR(LauResonanceInfo* resInfo, const LauAbsResonance::LauResonanceModel resType,
const Int_t resPairAmpInt, const LauDaughters* daughters);
//! Destructor
virtual ~LauPolarFormFactorNR();
//! Initialise the model
virtual void initialise();
//! Get the resonance model type
/*!
\return the resonance model type
*/
virtual LauAbsResonance::LauResonanceModel getResonanceModel() const {return model_;}
//! Set value of the various parameters
/*!
\param [in] name the name of the parameter to be changed
\param [in] value the new parameter value
*/
virtual void setResonanceParameter(const TString& name, const Double_t value);
//! Allow the various parameters to float in the fit
/*!
\param [in] name the name of the parameter to be floated
*/
virtual void floatResonanceParameter(const TString& name);
//! Access the given resonance parameter
/*!
\param [in] name the name of the parameter
\return the corresponding parameter
*/
virtual LauParameter* getResonanceParameter(const TString& name);
//! Retrieve the resonance parameters, e.g. so that they can be loaded into a fit
/*!
\return floating parameters of the resonance
*/
virtual const std::vector<LauParameter*>& getFloatingParameters();
protected:
//! Set the parameter lambda, the NR shape parameter
/*!
\param [in] lambda the NR shape parameter
*/
void setLambda(const Double_t lambda);
//! Get the parameter lambda, the NR shape parameter
/*!
\return lambda, the NR shape parameter
*/
Double_t getLambda() const {return (lambda_!=0) ? lambda_->value() : 0.0;}
//! See if the lambda parameter is fixed or floating
/*!
\return kTRUE if the lambda parameter is fixed, kFALSE otherwise
*/
Bool_t fixLambda() const {return (lambda_!=0) ? lambda_->fixed() : kTRUE;}
//! Complex resonant amplitude
/*!
\param [in] mass appropriate invariant mass for the resonance
\param [in] spinTerm Zemach spin term
*/
virtual LauComplex resAmp(Double_t mass, Double_t spinTerm);
private:
//! Copy constructor (not implemented)
LauPolarFormFactorNR(const LauPolarFormFactorNR& rhs);
//! Copy assignment operator (not implemented)
LauPolarFormFactorNR& operator=(const LauPolarFormFactorNR& rhs);
+ //! The NR shape parameter
LauParameter* lambda_;
//! The model to use
LauAbsResonance::LauResonanceModel model_;
ClassDef(LauPolarFormFactorNR,0)
};
#endif
diff --git a/inc/LauPolarFormFactorSymNR.hh b/inc/LauPolarFormFactorSymNR.hh
index f79de32..111b3a2 100644
--- a/inc/LauPolarFormFactorSymNR.hh
+++ b/inc/LauPolarFormFactorSymNR.hh
@@ -1,149 +1,150 @@
/*
Copyright 2018 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauPolarFormFactorSymNR.hh
\brief File containing declaration of LauPolarFormFactorSymNR class.
*/
/*! \class LauPolarFormFactorSymNR
\brief Class for defining a nonresonant form factor model
Defines the nonresonant form factor model from:
Nogueira, Bediaga, Cavalcante, Frederico, Lorenco: Phys. Rev. D92 (2015) 054010, arXiv:1506.08332 [hep-ph]
Pelaez, Yndurain: Phys. Rev. D71 (2005) 074016, arXiv:hep-ph/0411334
modified for symmetric DPs.
*/
#ifndef LAU_POLAR_FORM_FACTOR_SYM_NR
#define LAU_POLAR_FORM_FACTOR_SYM_NR
#include "TString.h"
#include "LauComplex.hh"
#include "LauAbsResonance.hh"
class LauKinematics;
class LauParameter;
class LauPolarFormFactorSymNR : public LauAbsResonance {
public:
//! Constructor
/*!
\param [in] resInfo the object containing information on the resonance name, mass, width, spin, charge, etc.
\param [in] resType the model of the resonance
\param [in] resPairAmpInt the number of the daughter not produced by the resonance
\param [in] daughters the daughter particles
*/
LauPolarFormFactorSymNR(LauResonanceInfo* resInfo, const LauAbsResonance::LauResonanceModel resType,
const Int_t resPairAmpInt, const LauDaughters* daughters);
//! Destructor
virtual ~LauPolarFormFactorSymNR();
//! Initialise the model
virtual void initialise();
//! Get the complex dynamical amplitude
/*!
\param [in] kinematics the kinematic variables of the current event
\return the complex amplitude
*/
virtual LauComplex amplitude(const LauKinematics* kinematics);
//! Get the resonance model type
/*!
\return the resonance model type
*/
virtual LauAbsResonance::LauResonanceModel getResonanceModel() const {return model_;}
//! Set value of the various parameters
/*!
\param [in] name the name of the parameter to be changed
\param [in] value the new parameter value
*/
virtual void setResonanceParameter(const TString& name, const Double_t value);
//! Allow the various parameters to float in the fit
/*!
\param [in] name the name of the parameter to be floated
*/
virtual void floatResonanceParameter(const TString& name);
//! Access the given resonance parameter
/*!
\param [in] name the name of the parameter
\return the corresponding parameter
*/
virtual LauParameter* getResonanceParameter(const TString& name);
//! Retrieve the resonance parameters, e.g. so that they can be loaded into a fit
/*!
\return floating parameters of the resonance
*/
virtual const std::vector<LauParameter*>& getFloatingParameters();
protected:
//! Set the parameter lambda, the NR shape parameter
/*!
\param [in] lambda the NR shape parameter
*/
void setLambda(const Double_t lambda);
//! Get the parameter lambda, the NR shape parameter
/*!
\return lambda the NR shape parameter
*/
Double_t getLambda() const {return (lambda_!=0) ? lambda_->value() : 0.0;}
//! See if the lambda parameter is fixed or floating
/*!
\return kTRUE if the lambda parameter is fixed, kFALSE otherwise
*/
Bool_t fixLambda() const {return (lambda_!=0) ? lambda_->fixed() : kTRUE;}
//! Complex resonant amplitude
/*!
\param [in] mass appropriate invariant mass for the resonance
\param [in] spinTerm Zemach spin term
*/
virtual LauComplex resAmp(Double_t mass, Double_t spinTerm);
private:
//! Copy constructor (not implemented)
LauPolarFormFactorSymNR(const LauPolarFormFactorSymNR& rhs);
//! Copy assignment operator (not implemented)
LauPolarFormFactorSymNR& operator=(const LauPolarFormFactorSymNR& rhs);
+ //! The NR shape parameter
LauParameter* lambda_;
//! The model to use
LauAbsResonance::LauResonanceModel model_;
ClassDef(LauPolarFormFactorSymNR,0)
};
#endif
diff --git a/inc/LauResonanceInfo.hh b/inc/LauResonanceInfo.hh
index cb7d701..3df334d 100644
--- a/inc/LauResonanceInfo.hh
+++ b/inc/LauResonanceInfo.hh
@@ -1,204 +1,205 @@
/*
Copyright 2006 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauResonanceInfo.hh
\brief File containing declaration of LauResonanceInfo class.
*/
/*! \class LauResonanceInfo
\brief Class for defining the properties of a resonant particle.
*/
#ifndef LAU_RESONANCE_INFO
#define LAU_RESONANCE_INFO
#include <iosfwd>
#include <set>
#include "TString.h"
#include "LauBlattWeisskopfFactor.hh"
class LauParameter;
class LauResonanceInfo {
public:
//! Constructor
/*!
\param [in] name the name of the resonant particle
\param [in] mass the mass of the resonant particle
\param [in] width the width of the resonant particle
\param [in] spin the spin of the resonant particle
\param [in] charge the charge of the resonant particle
\param [in] bwCategory the Blatt-Weisskopf barrier factor category
\param [in] bwRadius the Blatt-Weisskopf radius of the resonant particle
*/
LauResonanceInfo(const TString& name, const Double_t mass, const Double_t width, const Int_t spin, const Int_t charge, const LauBlattWeisskopfFactor::BlattWeisskopfCategory bwCategory, const Double_t bwRadius = 4.0);
//! Destructor
virtual ~LauResonanceInfo();
//! Retrieve the name of the resonant particle
/*!
\return the name of the resonant particle
*/
TString getName() const {return name_;}
//! Retrieve the sanitised name of the resonant particle
/*!
Removes/replaces characters from the name that cause
problems when used in TBranch names
\return the sanitised name of the resonant particle
*/
TString getSanitisedName() const {return sanitisedName_;}
//! Retrieve the mass of the resonant particle
/*!
\return the mass of the resonant particle
*/
LauParameter* getMass() const {return mass_;}
//! Retrieve the width of the resonant particle
/*!
\return the width of the resonant particle
*/
LauParameter* getWidth() const {return width_;}
//! Retrieve the spin of the resonant particle
/*!
\return the spin of the resonant particle
*/
UInt_t getSpin() const {return spin_;}
//! Retrieve the charge of the resonant particle
/*!
\return the charge of the resonant particle
*/
Int_t getCharge() const {return charge_;}
//! Retrieve the BW category of the resonant particle
/*!
\return the BW category of the resonant particle
*/
LauBlattWeisskopfFactor::BlattWeisskopfCategory getBWCategory() const {return bwCategory_;}
//! Retrieve the BW radius of the resonant particle
/*!
\return the BW radius of the resonant particle
*/
Double_t getBWRadius() const {return bwRadius_;}
//! Create the charge conjugate particle info record
/*!
The mass and width parameters are cloned
*/
LauResonanceInfo* createChargeConjugate();
//! Create another record that will share parameters with this one
/*!
The name needs to be specified.
The spin and charge are assumed to be the same.
The mass, width and other parameters will be cloned.
\param [in] name the name of the resonant particle
*/
LauResonanceInfo* createSharedParameterRecord( const TString& name );
//! Retrieve an extra parameter of the resonance
/*!
\return the extra parameter (or null pointer if not found)
*/
LauParameter* getExtraParameter( const TString& parName );
//! Add an extra parameter of the resonance
/*!
\param [in] param the extra parameter to be added
\param [in] independentPar governs whether any info record that usually shares parameters with this one should also share this one (the default) or make its own independent version
*/
void addExtraParameter( LauParameter* param, const Bool_t independentPar = kFALSE );
protected:
//! Add a clone of an extra parameter of the resonance
/*!
\param [in] param the extra parameter to be added
\param [in] copyNotClone should we create an unlinked copy instead of cloning - default is to clone
*/
void addCloneOfExtraParameter( LauParameter* param, const Bool_t copyNotClone = kFALSE );
private:
//! Copy constructor (not implemented)
LauResonanceInfo( const LauResonanceInfo& other );
//! Copy constructor (with new name and charge)
LauResonanceInfo( const LauResonanceInfo& other, const TString& newName, const Int_t newCharge );
//! Copy assignment operator (not implemented)
LauResonanceInfo& operator=( const LauResonanceInfo& other );
//! Create the sanitised name by removing characters that are illegal in TBranch names
void sanitiseName();
//! The name of the resonant particle
TString name_;
//! The name of the resonant particle with illegal characters removed
TString sanitisedName_;
//! The mass of the resonant particle
LauParameter* mass_;
//! The width of the resonant particle
LauParameter* width_;
//! The spin of the resonant particle
UInt_t spin_;
//! The charge of the resonant particle
Int_t charge_;
//! The Blatt-Weisskopf barrier factor category
LauBlattWeisskopfFactor::BlattWeisskopfCategory bwCategory_;
//! The Blatt-Weisskopf radius of the resonant particle
Double_t bwRadius_;
//! The conjugate info object
LauResonanceInfo* conjugate_;
//! Other info objects that share parameters with this one
std::vector<LauResonanceInfo*> sharedParRecords_;
//! Extra parameters
std::set<LauParameter*> extraPars_;
ClassDef(LauResonanceInfo, 0) // Specify each allowed resonance
};
+//! output operator formatting of an info record
std::ostream& operator<<( std::ostream& stream, const LauResonanceInfo& infoRecord );
#endif
diff --git a/src/LauCruijffPdf.cc b/src/LauCruijffPdf.cc
index 6fa8b41..ce3cc69 100644
--- a/src/LauCruijffPdf.cc
+++ b/src/LauCruijffPdf.cc
@@ -1,168 +1,167 @@
/*
Copyright 2008 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauCruijffPdf.cc
\brief File containing implementation of LauCruijffPdf class.
*/
/*****************************************************************************
* Class based on RooFit/RooCruijff. *
* Original copyright given below. *
*****************************************************************************
* Authors: *
* WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
* DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
* *
* Copyright (c) 2000-2005, Regents of the University of California *
* and Stanford University. All rights reserved. *
* *
* Redistribution and use in source and binary forms, *
* with or without modification, are permitted according to the terms *
* listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
*****************************************************************************/
#include <iostream>
#include <vector>
using std::cout;
using std::cerr;
using std::endl;
-using std::vector;
#include "TMath.h"
#include "TSystem.h"
#include "LauCruijffPdf.hh"
#include "LauConstants.hh"
ClassImp(LauCruijffPdf)
-LauCruijffPdf::LauCruijffPdf(const TString& theVarName, const vector<LauAbsRValue*>& params, Double_t minAbscissa, Double_t maxAbscissa) :
+LauCruijffPdf::LauCruijffPdf(const TString& theVarName, const std::vector<LauAbsRValue*>& params, Double_t minAbscissa, Double_t maxAbscissa) :
LauAbsPdf(theVarName, params, minAbscissa, maxAbscissa),
mean_(0),
sigmaL_(0),
sigmaR_(0),
alphaL_(0),
alphaR_(0)
{
// Constructor for the Cruijff PDF.
//
// The parameters in params are the mean, sigmaR, sigmaL, alphaR
// and alphaL.
// The last two arguments specify the range in which the PDF is defined, and the PDF
// will be normalised w.r.t. these limits.
mean_ = this->findParameter("mean");
sigmaR_ = this->findParameter("sigmaR");
sigmaL_ = this->findParameter("sigmaL");
alphaR_ = this->findParameter("alphaR");
alphaL_ = this->findParameter("alphaL");
if ((this->nParameters() != 5) || (mean_ == 0) || (sigmaR_ == 0) || (sigmaL_ == 0) || (alphaL_ == 0) || (alphaR_ == 0)) {
cerr<<"ERROR in LauCruijffPdf constructor: LauCruijffPdf requires 5 parameters: \"mean\", \"sigmaL\", \"sigmaR\", \"alphaR\" and \"alphaL\"."<<endl;
gSystem->Exit(EXIT_FAILURE);
}
// Cache the normalisation factor
this->calcNorm();
}
LauCruijffPdf::~LauCruijffPdf()
{
// Destructor
}
void LauCruijffPdf::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];
// Get the up to date parameter values
Double_t mean = mean_->unblindValue();
Double_t sigmaL = sigmaL_->unblindValue();
Double_t sigmaR = sigmaR_->unblindValue();
Double_t alphaL = alphaL_->unblindValue();
Double_t alphaR = alphaR_->unblindValue();
// Evaluate the LauCruijff PDF value
Double_t arg = abscissa - mean;
Double_t coef(0.0);
Double_t value(0.0);
if (arg < 0.0){
if (TMath::Abs(sigmaL) > 1e-30) {
coef = -1.0/(2.0*sigmaL*sigmaL + alphaL*arg*arg);
}
} else {
if (TMath::Abs(sigmaR) > 1e-30) {
coef = -1.0/(2.0*sigmaR*sigmaR+ alphaR*arg*arg);
}
}
value = TMath::Exp(coef*arg*arg);
// if the parameters are floating then we
// need to recalculate the normalisation
if (!this->cachePDF() && !this->withinNormCalc() && !this->withinGeneration()) {
this->calcNorm();
}
this->setUnNormPDFVal(value);
}
void LauCruijffPdf::calcPDFHeight( const LauKinematics* /*kinematics*/ )
{
if (this->heightUpToDate()) {
return;
}
// Get the up to date parameter values
Double_t mean = mean_->unblindValue();
LauAbscissas maxPoint(1);
maxPoint[0] = mean;
// Calculate the PDF height
if (mean < this->getMinAbscissa()) {
maxPoint[0] = this->getMinAbscissa();
} else if (mean > this->getMaxAbscissa()) {
maxPoint[0] = this->getMaxAbscissa();
}
this->calcLikelihoodInfo(maxPoint);
Double_t height = this->getUnNormLikelihood();
// Multiply by a small factor to avoid problems from rounding errors
height *= (1.0 + 1e-1);
this->setMaxHeight(height);
}
diff --git a/src/LauCrystalBallPdf.cc b/src/LauCrystalBallPdf.cc
index 9c7da10..b3d8e38 100644
--- a/src/LauCrystalBallPdf.cc
+++ b/src/LauCrystalBallPdf.cc
@@ -1,230 +1,229 @@
/*
Copyright 2006 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauCrystalBallPdf.cc
\brief File containing implementation of LauCrystalBallPdf class.
*/
/*****************************************************************************
* Class based on RooFit/RooCBShape. *
* Original copyright given below. *
*****************************************************************************
* Authors: *
* WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
* DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
* *
* Copyright (c) 2000-2005, Regents of the University of California *
* and Stanford University. All rights reserved. *
* *
* Redistribution and use in source and binary forms, *
* with or without modification, are permitted according to the terms *
* listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
*****************************************************************************/
#include <iostream>
#include <vector>
using std::cout;
using std::cerr;
using std::endl;
-using std::vector;
#include "TMath.h"
#include "TSystem.h"
#include "LauConstants.hh"
#include "LauCrystalBallPdf.hh"
ClassImp(LauCrystalBallPdf)
-LauCrystalBallPdf::LauCrystalBallPdf(const TString& theVarName, const vector<LauAbsRValue*>& params, Double_t minAbscissa, Double_t maxAbscissa) :
+LauCrystalBallPdf::LauCrystalBallPdf(const TString& theVarName, const std::vector<LauAbsRValue*>& params, Double_t minAbscissa, Double_t maxAbscissa) :
LauAbsPdf(theVarName, params, minAbscissa, maxAbscissa),
mean_(0),
sigma_(0),
alpha_(0),
n_(0)
{
// Constructor for the Crystal Ball PDF, which is a gaussian and a decaying tail
// smoothly matched up. The tail goes as a 1/x^n
//
// The parameters in params are the mean and the sigma (half the width) of the gaussian,
// the distance from the mean in which the gaussian and the tail are matched up (which
// can be negative or positive), and the power "n" for the tail.
// The last two arguments specify the range in which the PDF is defined, and the PDF
// will be normalised w.r.t. these limits.
mean_ = this->findParameter("mean");
sigma_ = this->findParameter("sigma");
alpha_ = this->findParameter("alpha");
n_ = this->findParameter("order");
if ((this->nParameters() != 4) || (mean_ == 0) || (sigma_ == 0) || (alpha_ == 0) || (n_ == 0)) {
cerr<<"ERROR in LauCrystalBallPdf constructor: LauCrystalBallPdf requires 4 parameters: \"mean\", \"sigma\", \"alpha\" and \"order\"."<<endl;
gSystem->Exit(EXIT_FAILURE);
}
// Cache the normalisation factor
this->calcNorm();
}
LauCrystalBallPdf::~LauCrystalBallPdf()
{
// Destructor
}
void LauCrystalBallPdf::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];
// Get the up to date parameter values
Double_t mean = mean_->unblindValue();
Double_t sigma = sigma_->unblindValue();
Double_t alpha = alpha_->unblindValue();
Double_t n = n_->unblindValue();
Double_t result(0.0);
Double_t t = (abscissa - mean)/sigma;
if (alpha < 0.0) {
t = -t;
}
Double_t absAlpha = TMath::Abs(alpha);
if (t >= -absAlpha) {
result = TMath::Exp(-0.5*t*t);
} else {
Double_t a = TMath::Power(n/absAlpha,n)*TMath::Exp(-0.5*absAlpha*absAlpha);
Double_t b = n/absAlpha - absAlpha;
result = a/TMath::Power(b - t, n);
}
this->setUnNormPDFVal(result);
// if the parameters are floating then we
// need to recalculate the normalisation
if (!this->cachePDF() && !this->withinNormCalc() && !this->withinGeneration()) {
this->calcNorm();
}
}
void LauCrystalBallPdf::calcNorm()
{
// Get the up to date parameter values
Double_t mean = mean_->unblindValue();
Double_t sigma = sigma_->unblindValue();
Double_t alpha = alpha_->unblindValue();
Double_t n = n_->unblindValue();
Double_t result = 0.0;
Bool_t useLog = kFALSE;
if ( TMath::Abs(n-1.0) < 1.0e-05 ) {
useLog = kTRUE;
}
Double_t sig = TMath::Abs(sigma);
Double_t tmin = (this->getMinAbscissa() - mean)/sig;
Double_t tmax = (this->getMaxAbscissa() - mean)/sig;
if (alpha < 0) {
Double_t tmp = tmin;
tmin = -tmax;
tmax = -tmp;
}
Double_t absAlpha = TMath::Abs(alpha);
if ( tmin >= -absAlpha ) {
result += sig*LauConstants::rootPiBy2*( this->approxErf( tmax / LauConstants::root2 )
- approxErf( tmin / LauConstants::root2 ) );
} else if ( tmax <= -absAlpha ) {
Double_t a = TMath::Power(n/absAlpha, n)*TMath::Exp( -0.5*absAlpha*absAlpha);
Double_t b = n/absAlpha - absAlpha;
if ( useLog == kTRUE ) {
result += a*sig*( TMath::Log(b - tmin) - TMath::Log(b - tmax) );
} else {
result += a*sig/(1.0 - n)*( 1.0/(TMath::Power( b - tmin, n - 1.0))
- 1.0/(TMath::Power( b - tmax, n - 1.0)) );
}
} else {
Double_t a = TMath::Power(n/absAlpha, n)*TMath::Exp( -0.5*absAlpha*absAlpha );
Double_t b = n/absAlpha - absAlpha;
Double_t term1 = 0.0;
if ( useLog == kTRUE )
term1 = a*sig*( TMath::Log(b - tmin) - TMath::Log(n / absAlpha));
else
term1 = a*sig/(1.0 - n)*( 1.0/(TMath::Power( b - tmin, n - 1.0))
- 1.0/(TMath::Power( n/absAlpha, n - 1.0)) );
Double_t term2 = sig*LauConstants::rootPiBy2*( this->approxErf( tmax / LauConstants::root2 )
- this->approxErf( -absAlpha / LauConstants::root2 ) );
result += term1 + term2;
}
this->setNorm(result);
}
void LauCrystalBallPdf::calcPDFHeight( const LauKinematics* /*kinematics*/ )
{
if (this->heightUpToDate()) {
return;
}
// Get the up to date parameter values
Double_t mean = mean_->unblindValue();
// The Crystall Ball function is a Gaussian with an exponentially decaying tail
// Therefore, calculate the PDF height for the Gaussian function.
LauAbscissas abscissa(1);
abscissa[0] = mean;
this->calcLikelihoodInfo(abscissa);
Double_t height = this->getUnNormLikelihood();
this->setMaxHeight(height);
}
Double_t LauCrystalBallPdf::approxErf(Double_t arg) const
{
static const Double_t erflim = 5.0;
if ( arg > erflim ) {
return 1.0;
}
if ( arg < -erflim ) {
return -1.0;
}
return TMath::Erf(arg);
}
diff --git a/src/LauDPDepBifurGaussPdf.cc b/src/LauDPDepBifurGaussPdf.cc
index 9974d6b..3444785 100644
--- a/src/LauDPDepBifurGaussPdf.cc
+++ b/src/LauDPDepBifurGaussPdf.cc
@@ -1,283 +1,282 @@
/*
Copyright 2007 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauDPDepBifurGaussPdf.cc
\brief File containing implementation of LauDPDepBifurGaussPdf class.
*/
#include <iostream>
#include <vector>
using std::cout;
using std::cerr;
using std::endl;
-using std::vector;
#include "TMath.h"
#include "TRandom.h"
#include "TSystem.h"
#include "Lau1DHistPdf.hh"
#include "LauConstants.hh"
#include "LauComplex.hh"
#include "LauDaughters.hh"
#include "LauFitDataTree.hh"
#include "LauKinematics.hh"
#include "LauDPDepBifurGaussPdf.hh"
ClassImp(LauDPDepBifurGaussPdf)
-LauDPDepBifurGaussPdf::LauDPDepBifurGaussPdf(const TString& theVarName, const vector<LauAbsRValue*>& params,
+LauDPDepBifurGaussPdf::LauDPDepBifurGaussPdf(const TString& theVarName, const std::vector<LauAbsRValue*>& params,
Double_t minAbscissa, Double_t maxAbscissa,
const LauDaughters* daughters,
- const vector<Double_t>& meanCoeffs,
- const vector<Double_t>& sigmaLCoeffs,
- const vector<Double_t>& sigmaRCoeffs,
+ const std::vector<Double_t>& meanCoeffs,
+ const std::vector<Double_t>& sigmaLCoeffs,
+ const std::vector<Double_t>& sigmaRCoeffs,
DPAxis dpAxis) :
LauAbsPdf(theVarName, params, minAbscissa, maxAbscissa),
kinematics_(daughters ? daughters->getKinematics() : 0),
mean_(0),
sigmaL_(0),
sigmaR_(0),
meanVal_(0.0),
sigmaLVal_(0.0),
sigmaRVal_(0.0),
meanCoeffs_(meanCoeffs),
sigmaLCoeffs_(sigmaLCoeffs),
sigmaRCoeffs_(sigmaRCoeffs),
dpAxis_(dpAxis),
scaleMethod_(poly)
{
// Check we have a valid kinematics object
if ( ! kinematics_ ) {
cerr<<"ERROR in LauDPDepBifurGaussPdf::LauDPDepBifurGaussPdf : Have not been provided with a valid DP kinematics object."<<endl;
gSystem->Exit(EXIT_FAILURE);
}
// The parameters are:
// - the mean, the sigmaL and the sigmaR of the Bifurcated Gaussian
//
// The next two arguments specify the range in which the PDF is defined,
// and the PDF will be normalised w.r.t. these limits.
//
// The final three argument define whether the BF Gaussian parameters should be scaled or not
mean_ = this->findParameter("mean");
sigmaL_ = this->findParameter("sigmaL");
sigmaR_ = this->findParameter("sigmaR");
if ((this->nParameters() != 3) || (mean_ == 0) || (sigmaL_ == 0) || (sigmaR_ == 0) ) {
cerr<<"ERROR in LauDPDepBifurGaussPdf constructor: LauDPDepBifurGaussPdf requires 3 parameters:"
<<" \"mean\", \"sigmaL\" and \"sigmaR\"."<<endl;
gSystem->Exit(EXIT_FAILURE);
}
// Cache the normalisation factor
this->calcNorm();
}
LauDPDepBifurGaussPdf::~LauDPDepBifurGaussPdf()
{
// Destructor
}
void LauDPDepBifurGaussPdf::calcLikelihoodInfo(const LauAbscissas& abscissas)
{
// Check that the given abscissas are within the allowed range
if (!this->checkRange(abscissas)) {
gSystem->Exit(EXIT_FAILURE);
}
// Get our abscissa
Double_t abscissa = abscissas[0];
// Get the up to date parameter values
meanVal_ = mean_->unblindValue();
sigmaLVal_ = sigmaL_->unblindValue();
sigmaRVal_ = sigmaR_->unblindValue();
// Find out the DP position
Double_t dpPos(0.0);
UInt_t nVars = this->nInputVars();
if ( abscissas.size() == nVars+2 ) {
Double_t m13Sq = abscissas[nVars];
Double_t m23Sq = abscissas[nVars+1];
if ( dpAxis_ == M12 ) {
dpPos = kinematics_->calcThirdMassSq(m13Sq,m23Sq);
} else if ( dpAxis_ == M13 ) {
dpPos = m13Sq;
} else if ( dpAxis_ == M23 ) {
dpPos = m23Sq;
} else if ( dpAxis_ == MMIN ) {
dpPos = TMath::Min( m13Sq, m23Sq );
} else if ( dpAxis_ == MMAX ) {
dpPos = TMath::Max( m13Sq, m23Sq );
} else {
dpPos = kinematics_->distanceFromDPCentre(m13Sq,m23Sq);
}
}
// Scale the gaussian parameters by the dpPos (if appropriate)
ScaleMethod scale = this->scaleMethod();
if ( scale==poly ) {
this->scalePars_poly(dpPos);
} else if ( scale==polyNegPower ) {
this->scalePars_polyNegPower(dpPos);
} else {
cerr<<"Scaling method unknown! Methods available: <poly> and <polyNegPower>."<<endl;
gSystem->Exit(EXIT_FAILURE);
}
// Evaluate the Birfucated Gaussian PDF value
Double_t arg = abscissa - meanVal_;
Double_t coef(0.0);
Double_t value(0.0);
if (arg < 0.0) {
if (TMath::Abs(sigmaLVal_) > 1e-30) {
coef = -0.5/(sigmaLVal_*sigmaLVal_);
}
} else {
if (TMath::Abs(sigmaRVal_) > 1e-30) {
coef = -0.5/(sigmaRVal_*sigmaRVal_);
}
}
value = TMath::Exp(coef*arg*arg);
// Calculate the norm
Double_t xscaleL = LauConstants::root2*sigmaLVal_;
Double_t xscaleR = LauConstants::root2*sigmaRVal_;
Double_t integral(0.0);
Double_t norm(0.0);
Double_t result(0.0);
if ( this->getMaxAbscissa() < meanVal_ ) {
integral = sigmaLVal_ * ( TMath::Erf((this->getMaxAbscissa() - meanVal_)/xscaleL) - TMath::Erf((this->getMinAbscissa() - meanVal_)/xscaleL));
} else if ( this->getMinAbscissa() > meanVal_ ) {
integral = sigmaRVal_ * (TMath::Erf((this->getMaxAbscissa() - meanVal_)/xscaleR) - TMath::Erf((this->getMinAbscissa() - meanVal_)/xscaleR));
} else {
integral = sigmaRVal_*TMath::Erf((this->getMaxAbscissa() -meanVal_)/xscaleR) - sigmaLVal_*TMath::Erf((this->getMinAbscissa() - meanVal_)/xscaleL);
}
norm = LauConstants::rootPiBy2*integral;
// the result
result = value/norm;
this->setUnNormPDFVal(result);
}
void LauDPDepBifurGaussPdf::scalePars_poly(Double_t dpPos)
{
Int_t power = 1;
- for (vector<Double_t>::const_iterator iter = meanCoeffs_.begin(); iter != meanCoeffs_.end(); ++iter) {
+ for (std::vector<Double_t>::const_iterator iter = meanCoeffs_.begin(); iter != meanCoeffs_.end(); ++iter) {
Double_t coeff = (*iter);
meanVal_ += coeff * TMath::Power(dpPos,power);
++power;
}
power = 1;
- for (vector<Double_t>::const_iterator iter = sigmaLCoeffs_.begin(); iter != sigmaLCoeffs_.end(); ++iter) {
+ for (std::vector<Double_t>::const_iterator iter = sigmaLCoeffs_.begin(); iter != sigmaLCoeffs_.end(); ++iter) {
Double_t coeff = (*iter);
sigmaLVal_ += coeff * TMath::Power(dpPos,power);
++power;
}
power = 1;
- for (vector<Double_t>::const_iterator iter = sigmaRCoeffs_.begin(); iter != sigmaRCoeffs_.end(); ++iter) {
+ for (std::vector<Double_t>::const_iterator iter = sigmaRCoeffs_.begin(); iter != sigmaRCoeffs_.end(); ++iter) {
Double_t coeff = (*iter);
sigmaRVal_ += coeff * TMath::Power(dpPos,power);
++power;
}
}
void LauDPDepBifurGaussPdf::scalePars_polyNegPower(Double_t dpPos)
{
Int_t power = -1;
- for (vector<Double_t>::const_iterator iter = meanCoeffs_.begin(); iter != meanCoeffs_.end(); ++iter) {
+ for (std::vector<Double_t>::const_iterator iter = meanCoeffs_.begin(); iter != meanCoeffs_.end(); ++iter) {
Double_t coeff = (*iter);
meanVal_ += coeff * TMath::Power(dpPos,power);
--power;
}
power = -1;
- for (vector<Double_t>::const_iterator iter = sigmaLCoeffs_.begin(); iter != sigmaLCoeffs_.end(); ++iter) {
+ for (std::vector<Double_t>::const_iterator iter = sigmaLCoeffs_.begin(); iter != sigmaLCoeffs_.end(); ++iter) {
Double_t coeff = (*iter);
sigmaLVal_ += coeff * TMath::Power(dpPos,power);
--power;
}
power = -1;
- for (vector<Double_t>::const_iterator iter = sigmaRCoeffs_.begin(); iter != sigmaRCoeffs_.end(); ++iter) {
+ for (std::vector<Double_t>::const_iterator iter = sigmaRCoeffs_.begin(); iter != sigmaRCoeffs_.end(); ++iter) {
Double_t coeff = (*iter);
sigmaRVal_ += coeff * TMath::Power(dpPos,power);
--power;
}
}
void LauDPDepBifurGaussPdf::calcNorm()
{
this->setNorm(1.0);
}
void LauDPDepBifurGaussPdf::calcPDFHeight( const LauKinematics* kinematics )
{
// Get the up to date parameter values
meanVal_ = mean_->unblindValue();
sigmaLVal_ = sigmaL_->unblindValue();
sigmaRVal_ = sigmaR_->unblindValue();
// Scale the gaussian parameters by the dpCentreDist (if appropriate)
Double_t dpCentreDist = kinematics->distanceFromDPCentre();
ScaleMethod scale = this->scaleMethod();
if ( scale==poly ) {
this->scalePars_poly(dpCentreDist);
} else if ( scale==polyNegPower ) {
this->scalePars_polyNegPower(dpCentreDist);
} else {
cerr<<"Scaling method unknown! Methods available: <poly> and <polyNegPower>."<<endl;
gSystem->Exit(EXIT_FAILURE);
}
// Calculate the PDF height for the Bifurcated Gaussian function.
LauAbscissas maxPoint(3);
maxPoint[0] = meanVal_;
maxPoint[1] = kinematics->getm13Sq();
maxPoint[2] = kinematics->getm23Sq();
if ( meanVal_ > this->getMaxAbscissa() ) {
maxPoint[0] = this->getMaxAbscissa();
} else if ( meanVal_ < this->getMinAbscissa() ) {
maxPoint[0] = this->getMinAbscissa();
}
this->calcLikelihoodInfo(maxPoint);
Double_t height = this->getUnNormLikelihood();
// Multiply by a small factor to avoid problems from rounding errors
height *= 1.001;
this->setMaxHeight(height);
}
diff --git a/src/LauDPDepSumPdf.cc b/src/LauDPDepSumPdf.cc
index 47f0815..0695409 100644
--- a/src/LauDPDepSumPdf.cc
+++ b/src/LauDPDepSumPdf.cc
@@ -1,432 +1,431 @@
/*
Copyright 2009 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauDPDepSumPdf.cc
\brief File containing implementation of LauDPDepSumPdf class.
*/
#include <iostream>
#include <vector>
using std::cout;
using std::cerr;
using std::endl;
-using std::vector;
#include "TMath.h"
#include "TSystem.h"
#include "LauConstants.hh"
#include "LauDaughters.hh"
#include "LauEffModel.hh"
#include "LauParameter.hh"
#include "LauDPDepSumPdf.hh"
ClassImp(LauDPDepSumPdf)
LauDPDepSumPdf::LauDPDepSumPdf(LauAbsPdf* pdf1, LauAbsPdf* pdf2,
const LauDaughters* daughters,
const TH2* dpHisto, Bool_t upperHalf, Bool_t useSpline) :
- LauAbsPdf(pdf1 ? pdf1->varNames() : vector<TString>(), vector<LauAbsRValue*>(), pdf1 ? pdf1->getMinAbscissas() : LauFitData(), pdf1 ? pdf1->getMaxAbscissas() : LauFitData()),
+ LauAbsPdf(pdf1 ? pdf1->varNames() : std::vector<TString>(), std::vector<LauAbsRValue*>(), pdf1 ? pdf1->getMinAbscissas() : LauFitData(), pdf1 ? pdf1->getMaxAbscissas() : LauFitData()),
daughters_( new LauDaughters(*daughters) ),
pdf1_(pdf1),
pdf2_(pdf2),
frac_(0),
fracVal_(0.5),
dpDependence_( new LauEffModel(daughters, 0) ),
dpAxis_( CentreDist )
{
// Constructor for the sum PDF.
// We are defining the sum as:
// f x (PDF1/S(PDF1)) + (1-f) x (PDF2/S(PDF2))
// where f is the fraction, x is multiplication, PDFi is the i'th PDF,
// and S(PDFi) is the integral of the i'th PDF.
// The value of the fraction is read from the DP histogram.
if(useSpline) {
dpDependence_->setEffSpline( dpHisto, kFALSE, 0.0, 0.0, upperHalf, daughters->squareDP());
} else {
dpDependence_->setEffHisto( dpHisto, kTRUE, kFALSE, 0.0, 0.0, upperHalf, daughters->squareDP() );
}
// So the first thing we have to do is check the pointers are all valid.
if (!pdf1 || !pdf2) {
cerr<<"ERROR in LauDPDepSumPdf constructor: one of the 2 PDF pointers is null."<<endl;
gSystem->Exit(EXIT_FAILURE);
}
// Next check that the abscissa ranges are the same for each PDF
if (pdf1->getMinAbscissa() != pdf2->getMinAbscissa()) {
cerr<<"ERROR in LauDPDepSumPdf constructor: minimum abscissa values not the same for the two PDFs."<<endl;
gSystem->Exit(EXIT_FAILURE);
}
if (pdf1->getMaxAbscissa() != pdf2->getMaxAbscissa()) {
cerr<<"ERROR in LauDPDepSumPdf constructor: maximum abscissa values not the same for the two PDFs."<<endl;
gSystem->Exit(EXIT_FAILURE);
}
// Also check that both PDFs are expecting the same number of input variables
if (pdf1->nInputVars() != pdf2->nInputVars()) {
cerr<<"ERROR in LauDPDepSumPdf constructor: number of input variables not the same for the two PDFs."<<endl;
gSystem->Exit(EXIT_FAILURE);
}
// Also check that both PDFs are expecting the same variable name(s)
if (pdf1->varNames() != pdf2->varNames()) {
cerr<<"ERROR in LauDPDepSumPdf constructor: variable name(s) not the same for the two PDFs."<<endl;
gSystem->Exit(EXIT_FAILURE);
}
// Then we need to grab all the parameters and pass them to the base class.
// This is so that when we are asked for them they can be put into the fit.
// The number of parameters is the number in PDF1 + the number in PDF2.
UInt_t nPar(pdf1->nParameters()+pdf2->nParameters());
- vector<LauAbsRValue*> params; params.reserve(nPar);
- vector<LauAbsRValue*>& pdf1pars = pdf1->getParameters();
- vector<LauAbsRValue*>& pdf2pars = pdf2->getParameters();
- for (vector<LauAbsRValue*>::iterator iter = pdf1pars.begin(); iter != pdf1pars.end(); ++iter) {
+ std::vector<LauAbsRValue*> params; params.reserve(nPar);
+ std::vector<LauAbsRValue*>& pdf1pars = pdf1->getParameters();
+ std::vector<LauAbsRValue*>& pdf2pars = pdf2->getParameters();
+ for (std::vector<LauAbsRValue*>::iterator iter = pdf1pars.begin(); iter != pdf1pars.end(); ++iter) {
params.push_back(*iter);
}
- for (vector<LauAbsRValue*>::iterator iter = pdf2pars.begin(); iter != pdf2pars.end(); ++iter) {
+ for (std::vector<LauAbsRValue*>::iterator iter = pdf2pars.begin(); iter != pdf2pars.end(); ++iter) {
params.push_back(*iter);
}
this->addParameters(params);
// Cache the normalisation factor
this->calcNorm();
}
LauDPDepSumPdf::LauDPDepSumPdf(LauAbsPdf* pdf1, LauAbsPdf* pdf2,
LauParameter* frac,
const LauDaughters* daughters,
const std::vector<Double_t>& fracCoeffs,
DPAxis dpAxis) :
- LauAbsPdf(pdf1 ? pdf1->varNames() : vector<TString>(), vector<LauAbsRValue*>(), pdf1 ? pdf1->getMinAbscissas() : LauFitData(), pdf1 ? pdf1->getMaxAbscissas() : LauFitData()),
+ LauAbsPdf(pdf1 ? pdf1->varNames() : std::vector<TString>(), std::vector<LauAbsRValue*>(), pdf1 ? pdf1->getMinAbscissas() : LauFitData(), pdf1 ? pdf1->getMaxAbscissas() : LauFitData()),
daughters_( new LauDaughters(*daughters) ),
pdf1_(pdf1),
pdf2_(pdf2),
frac_(frac),
fracVal_( frac ? frac->unblindValue() : 0.0 ),
dpDependence_( 0 ),
fracCoeffs_( fracCoeffs ),
dpAxis_( dpAxis )
{
// Constructor for the sum PDF.
// We are defining the sum as:
// f x (PDF1/S(PDF1)) + (1-f) x (PDF2/S(PDF2))
// where f is the fraction, x is multiplication, PDFi is the i'th PDF,
// and S(PDFi) is the integral of the i'th PDF.
// The value of the fraction has a polynomial dependence on one of
// the DP axes.
// So the first thing we have to do is check the pointers are all valid.
if (!pdf1 || !pdf2) {
cerr<<"ERROR in LauDPDepSumPdf constructor: one of the 2 PDF pointers is null."<<endl;
gSystem->Exit(EXIT_FAILURE);
}
if ( !frac ) {
cerr<<"ERROR in LauDPDepSumPdf constructor: the fraction parameter pointer is null."<<endl;
gSystem->Exit(EXIT_FAILURE);
}
// Next check that the abscissa ranges are the same for each PDF
if (pdf1->getMinAbscissa() != pdf2->getMinAbscissa()) {
cerr<<"ERROR in LauDPDepSumPdf constructor: minimum abscissa values not the same for the two PDFs."<<endl;
gSystem->Exit(EXIT_FAILURE);
}
if (pdf1->getMaxAbscissa() != pdf2->getMaxAbscissa()) {
cerr<<"ERROR in LauDPDepSumPdf constructor: maximum abscissa values not the same for the two PDFs."<<endl;
gSystem->Exit(EXIT_FAILURE);
}
// Also check that both PDFs are expecting the same number of input variables
if (pdf1->nInputVars() != pdf2->nInputVars()) {
cerr<<"ERROR in LauDPDepSumPdf constructor: number of input variables not the same for the two PDFs."<<endl;
gSystem->Exit(EXIT_FAILURE);
}
// Also check that both PDFs are expecting the same variable name(s)
if (pdf1->varNames() != pdf2->varNames()) {
cerr<<"ERROR in LauDPDepSumPdf constructor: variable name(s) not the same for the two PDFs."<<endl;
gSystem->Exit(EXIT_FAILURE);
}
// Then we need to grab all the parameters and pass them to the base class.
// This is so that when we are asked for them they can be put into the fit.
// The number of parameters is the number in PDF1 + the number in PDF2.
UInt_t nPar( pdf1->nParameters() + pdf2->nParameters() + 1 );
- vector<LauAbsRValue*> params; params.reserve(nPar);
+ std::vector<LauAbsRValue*> params; params.reserve(nPar);
params.push_back(frac);
- vector<LauAbsRValue*>& pdf1pars = pdf1->getParameters();
- vector<LauAbsRValue*>& pdf2pars = pdf2->getParameters();
- for (vector<LauAbsRValue*>::iterator iter = pdf1pars.begin(); iter != pdf1pars.end(); ++iter) {
+ std::vector<LauAbsRValue*>& pdf1pars = pdf1->getParameters();
+ std::vector<LauAbsRValue*>& pdf2pars = pdf2->getParameters();
+ for (std::vector<LauAbsRValue*>::iterator iter = pdf1pars.begin(); iter != pdf1pars.end(); ++iter) {
params.push_back(*iter);
}
- for (vector<LauAbsRValue*>::iterator iter = pdf2pars.begin(); iter != pdf2pars.end(); ++iter) {
+ for (std::vector<LauAbsRValue*>::iterator iter = pdf2pars.begin(); iter != pdf2pars.end(); ++iter) {
params.push_back(*iter);
}
this->addParameters(params);
// Now check that we can find the fraction parameter ok
frac_ = this->findParameter("frac");
if (frac_ == 0) {
cerr<<"ERROR in LauDPDepSumPdf constructor: parameter \"frac\" not found."<<endl;
gSystem->Exit(EXIT_FAILURE);
}
// Cache the normalisation factor
this->calcNorm();
}
LauDPDepSumPdf::~LauDPDepSumPdf()
{
// Destructor
delete daughters_; daughters_ = 0;
delete dpDependence_; dpDependence_ = 0;
}
void LauDPDepSumPdf::calcLikelihoodInfo(const LauAbscissas& abscissas)
{
// Check that the given abscissa is within the allowed range
if (!this->checkRange(abscissas)) {
gSystem->Exit(EXIT_FAILURE);
}
LauAbscissas noDPVars(1);
noDPVars[0] = abscissas[0];
// Evaluate the normalised PDF values
if ( pdf1_->isDPDependent() ) {
pdf1_->calcLikelihoodInfo(abscissas);
} else {
pdf1_->calcLikelihoodInfo(noDPVars);
}
if ( pdf2_->isDPDependent() ) {
pdf2_->calcLikelihoodInfo(abscissas);
} else {
pdf2_->calcLikelihoodInfo(noDPVars);
}
Double_t result1 = pdf1_->getLikelihood();
Double_t result2 = pdf2_->getLikelihood();
// Determine the fraction
// The DP variables will be abscissas[nInputVars] and
// abscissas[nInputVars+1] (if present).
UInt_t nVars = this->nInputVars();
if ( abscissas.size() == nVars+2 ) {
Double_t m13Sq = abscissas[nVars];
Double_t m23Sq = abscissas[nVars+1];
LauKinematics* kinematics = daughters_->getKinematics();
if ( dpDependence_ ) {
kinematics->updateKinematics( m13Sq, m23Sq );
fracVal_ = dpDependence_->calcEfficiency( kinematics );
} else {
fracVal_ = frac_->unblindValue();
Double_t dpPos(0.0);
if ( dpAxis_ == M12 ) {
dpPos = kinematics->calcThirdMassSq(m13Sq,m23Sq);
} else if ( dpAxis_ == M13 ) {
dpPos = m13Sq;
} else if ( dpAxis_ == M23 ) {
dpPos = m23Sq;
} else if ( dpAxis_ == MMIN ) {
dpPos = TMath::Min( m13Sq, m23Sq );
} else if ( dpAxis_ == MMAX ) {
dpPos = TMath::Max( m13Sq, m23Sq );
} else {
dpPos = kinematics->distanceFromDPCentre(m13Sq,m23Sq);
}
this->scaleFrac( dpPos );
}
}
// Add them together
Double_t result = fracVal_ * result1 + (1.0-fracVal_) * result2;
this->setUnNormPDFVal(result);
}
void LauDPDepSumPdf::scaleFrac( Double_t dpPos )
{
Int_t power = 1;
for (std::vector<Double_t>::const_iterator iter = fracCoeffs_.begin(); iter != fracCoeffs_.end(); ++iter) {
Double_t coeff = (*iter);
fracVal_ += coeff * TMath::Power(dpPos,power);
++power;
}
}
void LauDPDepSumPdf::calcNorm()
{
// Nothing to do here, since it is already normalized
this->setNorm(1.0);
}
void LauDPDepSumPdf::calcPDFHeight( const LauKinematics* kinematics )
{
// This method gives you the maximum possible height of the PDF.
// It combines the maximum heights of the two individual PDFs.
// So it would give the true maximum if the two individual maxima coincided.
// It is guaranteed to always be >= the true maximum.
// Update the heights of the individual PDFs
pdf1_->calcPDFHeight( kinematics );
pdf2_->calcPDFHeight( kinematics );
// Get the up to date parameter values
if ( dpDependence_ ) {
fracVal_ = dpDependence_->calcEfficiency( kinematics );
} else {
fracVal_ = frac_->unblindValue();
Double_t dpPos(0.0);
if ( dpAxis_ == M12 ) {
dpPos = kinematics->getm12Sq();
} else if ( dpAxis_ == M13 ) {
dpPos = kinematics->getm13Sq();
} else if ( dpAxis_ == M23 ) {
dpPos = kinematics->getm23Sq();
} else if ( dpAxis_ == MMIN ) {
Double_t m13Sq = kinematics->getm13Sq();
Double_t m23Sq = kinematics->getm23Sq();
dpPos = TMath::Min( m13Sq, m23Sq );
} else if ( dpAxis_ == MMAX ) {
Double_t m13Sq = kinematics->getm13Sq();
Double_t m23Sq = kinematics->getm23Sq();
dpPos = TMath::Max( m13Sq, m23Sq );
} else {
dpPos = kinematics->distanceFromDPCentre();
}
this->scaleFrac( dpPos );
}
// Find the (un-normalised) individual PDF maxima
Double_t height1 = pdf1_->getMaxHeight();
Double_t height2 = pdf2_->getMaxHeight();
// Get the individual PDF normalisation factors
Double_t norm1 = pdf1_->getNorm();
Double_t norm2 = pdf2_->getNorm();
// Calculate the normalised individual PDF maxima
height1 /= norm1;
height2 /= norm2;
// Combine these heights together
Double_t height = fracVal_ * height1 + (1-fracVal_) * height2;
this->setMaxHeight(height);
}
// Override the base class methods for cacheInfo and calcLikelihoodInfo(UInt_t iEvt).
// For both of these we delegate to the two constituent PDFs.
void LauDPDepSumPdf::cacheInfo(const LauFitDataTree& inputData)
{
// delegate to the two sub-PDFs to cache their information
pdf1_->cacheInfo(inputData);
pdf2_->cacheInfo(inputData);
// now check that the DP variables are included in the data
Bool_t hasBranch = inputData.haveBranch( "m13Sq" );
hasBranch &= inputData.haveBranch( "m23Sq" );
if (!hasBranch) {
cerr<<"ERROR in LauDPDepSumPdf::cacheInfo : Input data does not contain Dalitz plot variables m13Sq and/or m23Sq."<<endl;
return;
}
// determine whether we are caching our PDF value
Bool_t doCaching( this->nFixedParameters() == this->nParameters() );
this->cachePDF( doCaching );
if ( !doCaching ) {
// in this case we seem to be doing a fit where the parameters are floating
// so need to mark that the PDF height is no longer up to date
this->heightUpToDate(kFALSE);
}
// clear the vector and reserve enough space
UInt_t nEvents = inputData.nEvents();
fractions_.clear();
fractions_.reserve(nEvents);
// loop through the events, determine the fraction and store
for (UInt_t iEvt = 0; iEvt < nEvents; iEvt++) {
const LauFitData& dataValues = inputData.getData(iEvt);
Double_t m13Sq = dataValues.find("m13Sq")->second;
Double_t m23Sq = dataValues.find("m23Sq")->second;
LauKinematics* kinematics = daughters_->getKinematics();
if ( dpDependence_ ) {
// if we're using the histogram then just
// determine the fraction and store
kinematics->updateKinematics( m13Sq, m23Sq );
fracVal_ = dpDependence_->calcEfficiency( kinematics );
} else {
// if we're scaling the fraction parameter then we
// just store the scaling info since the parameter
// might be floating
fracVal_ = frac_->unblindValue();
Double_t dpPos(0.0);
if ( dpAxis_ == M12 ) {
dpPos = kinematics->calcThirdMassSq(m13Sq,m23Sq);
} else if ( dpAxis_ == M13 ) {
dpPos = m13Sq;
} else if ( dpAxis_ == M23 ) {
dpPos = m23Sq;
} else if ( dpAxis_ == MMIN ) {
dpPos = TMath::Min( m13Sq, m23Sq );
} else if ( dpAxis_ == MMAX ) {
dpPos = TMath::Max( m13Sq, m23Sq );
} else {
dpPos = kinematics->distanceFromDPCentre(m13Sq,m23Sq);
}
this->scaleFrac( dpPos );
fracVal_ -= frac_->unblindValue();
}
fractions_.push_back( fracVal_ );
}
}
void LauDPDepSumPdf::calcLikelihoodInfo(UInt_t iEvt)
{
// Get the fraction value for this event
fracVal_ = fractions_[iEvt];
if ( frac_ ) {
// if we're scaling the parameter then need to add the
// current value of the parameter
fracVal_ += frac_->unblindValue();
}
// Evaluate the normalised PDF values
pdf1_->calcLikelihoodInfo(iEvt);
pdf2_->calcLikelihoodInfo(iEvt);
Double_t result1 = pdf1_->getLikelihood();
Double_t result2 = pdf2_->getLikelihood();
// Add them together
Double_t result = fracVal_ * result1 + (1-fracVal_) * result2;
this->setUnNormPDFVal(result);
}
diff --git a/src/LauExponentialPdf.cc b/src/LauExponentialPdf.cc
index 3f07b92..e89aaef 100644
--- a/src/LauExponentialPdf.cc
+++ b/src/LauExponentialPdf.cc
@@ -1,140 +1,139 @@
/*
Copyright 2011 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauExponentialPdf.cc
\brief File containing implementation of LauExponentialPdf class.
*/
#include <iostream>
#include <vector>
using std::cout;
using std::cerr;
using std::endl;
-using std::vector;
#include "TMath.h"
#include "TSystem.h"
#include "LauConstants.hh"
#include "LauExponentialPdf.hh"
ClassImp(LauExponentialPdf)
-LauExponentialPdf::LauExponentialPdf(const TString& theVarName, const vector<LauAbsRValue*>& params, Double_t minAbscissa, Double_t maxAbscissa) :
+LauExponentialPdf::LauExponentialPdf(const TString& theVarName, const std::vector<LauAbsRValue*>& params, Double_t minAbscissa, Double_t maxAbscissa) :
LauAbsPdf(theVarName, params, minAbscissa, maxAbscissa),
slope_(0)
{
// Constructor for the Exponential PDF.
//
// The parameter in param is the slope of the exponential ie exp(slope*x).
// The last two arguments specify the range in which the PDF is defined, and the PDF
// will be normalised w.r.t. these limits.
slope_ = this->findParameter("slope");
if ((this->nParameters() != 1) || (slope_ == 0)) {
cerr<<"ERROR in LauExponentialPdf constructor: LauExponentialPdf requires 1 parameter: \"slope\"."<<endl;
gSystem->Exit(EXIT_FAILURE);
}
// Cache the normalisation factor
this->calcNorm();
}
LauExponentialPdf::~LauExponentialPdf()
{
// Destructor
}
void LauExponentialPdf::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];
// Get the up to date parameter values
Double_t slope = slope_->unblindValue();
// Calculate the value of the Exponential for the given value of the abscissa.
Double_t exponent(0.0);
if(slope != 0){
exponent = slope*abscissa;
}
Double_t value = TMath::Exp(exponent);
this->setUnNormPDFVal(value);
// if the parameters are floating then we
// need to recalculate the normalisation
if (!this->cachePDF() && !this->withinNormCalc() && !this->withinGeneration()) {
this->calcNorm();
}
}
void LauExponentialPdf::calcNorm()
{
// Get the up to date parameter values
Double_t slope = slope_->unblindValue();
// Calculate the normalisation of the exponential and cache it.
Double_t norm(0.0);
if(slope != 0){
norm = (1/slope)*(TMath::Exp(slope*this->getMaxAbscissa()) - TMath::Exp(slope*this->getMinAbscissa())) ;
}
this->setNorm(norm);
}
void LauExponentialPdf::calcPDFHeight(const LauKinematics* /*kinematics*/)
{
if (this->heightUpToDate()) {
return;
}
// Get the up to date parameter values
Double_t slope = slope_->unblindValue();
LauAbscissas maxPoint(1);
maxPoint[0] = 0;
// Calculate the PDF height for the Exponential function.
if (slope > 0) {
maxPoint[0] = this->getMaxAbscissa();
} else if (slope < 0) {
maxPoint[0] = this->getMinAbscissa();
}
this->calcLikelihoodInfo(maxPoint);
Double_t height = this->getUnNormLikelihood();
this->setMaxHeight(height);
}
diff --git a/src/LauGaussPdf.cc b/src/LauGaussPdf.cc
index 0c519bc..7f48922 100644
--- a/src/LauGaussPdf.cc
+++ b/src/LauGaussPdf.cc
@@ -1,146 +1,145 @@
/*
Copyright 2006 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauGaussPdf.cc
\brief File containing implementation of LauGaussPdf class.
*/
#include <iostream>
#include <vector>
using std::cout;
using std::cerr;
using std::endl;
-using std::vector;
#include "TMath.h"
#include "TSystem.h"
#include "LauConstants.hh"
#include "LauGaussPdf.hh"
ClassImp(LauGaussPdf)
-LauGaussPdf::LauGaussPdf(const TString& theVarName, const vector<LauAbsRValue*>& params, Double_t minAbscissa, Double_t maxAbscissa) :
+LauGaussPdf::LauGaussPdf(const TString& theVarName, const std::vector<LauAbsRValue*>& params, Double_t minAbscissa, Double_t maxAbscissa) :
LauAbsPdf(theVarName, params, minAbscissa, maxAbscissa),
mean_(0),
sigma_(0)
{
// Constructor for the Gaussian PDF.
//
// The parameters in params are the mean and the sigma (half the width) of the gaussian.
// The last two arguments specify the range in which the PDF is defined, and the PDF
// will be normalised w.r.t. these limits.
mean_ = this->findParameter("mean");
sigma_ = this->findParameter("sigma");
if ((this->nParameters() != 2) || (mean_ == 0) || (sigma_ == 0)) {
cerr<<"ERROR in LauGaussPdf constructor: LauGaussPdf requires 2 parameters: \"mean\" and \"sigma\"."<<endl;
gSystem->Exit(EXIT_FAILURE);
}
// Cache the normalisation factor
this->calcNorm();
}
LauGaussPdf::~LauGaussPdf()
{
// Destructor
}
void LauGaussPdf::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];
// Get the up to date parameter values
Double_t mean = mean_->unblindValue();
Double_t sigma = sigma_->unblindValue();
// Calculate the value of the Gaussian for the given value of the abscissa.
Double_t arg = abscissa - mean;
Double_t exponent(0.0);
if (TMath::Abs(sigma) > 1e-10) {
exponent = -0.5*arg*arg/(sigma*sigma);
}
Double_t value = TMath::Exp(exponent);
this->setUnNormPDFVal(value);
// if the parameters are floating then we
// need to recalculate the normalisation
if (!this->cachePDF() && !this->withinNormCalc() && !this->withinGeneration()) {
this->calcNorm();
}
}
void LauGaussPdf::calcNorm()
{
// Get the up to date parameter values
Double_t mean = mean_->unblindValue();
Double_t sigma = sigma_->unblindValue();
// Calculate the normalisation of the gaussian and cache it.
Double_t scale = LauConstants::root2*sigma;
Double_t norm(0.0);
if (TMath::Abs(sigma) > 1e-10) {
norm = LauConstants::rootPiBy2*sigma*(TMath::Erf((this->getMaxAbscissa() - mean)/scale) - TMath::Erf((this->getMinAbscissa() - mean)/scale));
}
this->setNorm(norm);
}
void LauGaussPdf::calcPDFHeight(const LauKinematics* /*kinematics*/)
{
if (this->heightUpToDate()) {
return;
}
// Get the up to date parameter values
Double_t mean = mean_->unblindValue();
LauAbscissas maxPoint(1);
maxPoint[0] = mean;
// Calculate the PDF height for the Gaussian function.
if (mean>this->getMaxAbscissa()) {
maxPoint[0] = this->getMaxAbscissa();
} else if (mean<this->getMinAbscissa()) {
maxPoint[0] = this->getMinAbscissa();
}
this->calcLikelihoodInfo(maxPoint);
Double_t height = this->getUnNormLikelihood();
this->setMaxHeight(height);
}
diff --git a/src/LauKinematics.cc b/src/LauKinematics.cc
index eb7d9f4..a263816 100644
--- a/src/LauKinematics.cc
+++ b/src/LauKinematics.cc
@@ -1,744 +1,743 @@
/*
Copyright 2004 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauKinematics.cc
\brief File containing implementation of LauKinematics class.
*/
#include <iostream>
#include "TF2.h"
#include "TMath.h"
#include "TRandom.h"
#include "LauConstants.hh"
#include "LauKinematics.hh"
#include "LauRandom.hh"
ClassImp(LauKinematics)
LauKinematics::LauKinematics(const Double_t m1, const Double_t m2, const Double_t m3, const Double_t mParent, const Bool_t calcSquareDPCoords, const Bool_t symmetricalDP, const Bool_t fullySymmetricDP) :
symmetricalDP_(symmetricalDP),
fullySymmetricDP_(fullySymmetricDP),
m1_(m1), m2_(m2), m3_(m3), mParent_(mParent),
m1Sq_(m1*m1), m2Sq_(m2*m2), m3Sq_(m3*m3), mParentSq_(mParent*mParent),
mDTot_(m1 + m2 + m3),
massInt_(mParent - (m1+m2+m3)),
mSqDTot_(m1*m1 + m2*m2 + m3*m3),
m12_(0.0), m23_(0.0), m13_(0.0),
m12Sq_(0.0), m23Sq_(0.0), m13Sq_(0.0),
c12_(0.0), c23_(0.0), c13_(0.0),
mPrime_(0.0), thetaPrime_(0.0),
qi_(0.0), qk_(0.0),
p1_12_(0.0), p3_12_(0.0),
p2_23_(0.0), p1_23_(0.0),
p1_13_(0.0), p2_13_(0.0),
p1_Parent_(0.0), p2_Parent_(0.0), p3_Parent_(0.0),
squareDP_(calcSquareDPCoords), warnings_(kTRUE)
{
// Constructor
mass_.clear(); mMin_.clear(); mMax_.clear(); mSqMin_.clear(); mSqMax_.clear();
mSq_.clear();
mSqDiff_.clear();
mDiff_.clear();
mass_.push_back(m1_);
mass_.push_back(m2_);
mass_.push_back(m3_);
mSq_.push_back(m1Sq_);
mSq_.push_back(m2Sq_);
mSq_.push_back(m3Sq_);
// DP max and min kinematic boundary for circumscribed box
// (see figure in PDG book).
for (Int_t i = 0; i < 3; i++) {
mMin_.push_back(mDTot_ - mass_[i]);
mMax_.push_back(mParent_ - mass_[i]);
mSqMin_.push_back(mMin_[i]*mMin_[i]);
mSqMax_.push_back(mMax_[i]*mMax_[i]);
mSqDiff_.push_back(mSqMax_[i] - mSqMin_[i]);
mDiff_.push_back(mMax_[i] - mMin_[i]);
}
if (this->squareDP()) {
std::cout<<"INFO in LauKinematics::LauKinematics : squareDP = kTRUE"<<std::endl;
} else {
std::cout<<"INFO in LauKinematics::LauKinematics : squareDP = kFALSE"<<std::endl;
}
// add covariant factor calculation
}
LauKinematics::~LauKinematics()
{
// Destructor
}
void LauKinematics::updateKinematics(const Double_t m13Sq, const Double_t m23Sq)
{
// Sets the internal private data members
// m13Sq and m23Sq, as well as m13 and m23.
// Also calculate m12Sq and m12, given mParent defined in the constructor.
// Lastly, calculate the helicity cosines.
// Update the 3 mass-squares
this->updateMassSquares(m13Sq, m23Sq);
// Now update the helicity cosines
this->calcHelicities();
// Also calculate the m', theta' variables
if (squareDP_) {this->calcSqDPVars();}
}
void LauKinematics::updateSqDPKinematics(const Double_t mPrime, const Double_t thetaPrime)
{
// From square DP co-ordinates, calculate remaining kinematic variables
this->updateSqDPMassSquares(mPrime, thetaPrime);
// Finally calculate the helicities and track momenta
this->calcHelicities();
}
void LauKinematics::calcSqDPVars()
{
// For given m_12 and cos(theta_12) values, calculate m' and theta' for the square Dalitz plot
Double_t value = (2.0*(m12_ - mMin_[2])/mDiff_[2]) - 1.0;
mPrime_ = LauConstants::invPi*TMath::ACos(value);
thetaPrime_ = LauConstants::invPi*TMath::ACos(c12_);
// Sometimes events are assigned exactly thetaPrime = 0.0 or 1.0
// which gives problems with efficiency and other histograms
if (thetaPrime_ == 0.0)
{
thetaPrime_ += 1.0e-10;
} else if (thetaPrime_ == 1.0)
{
thetaPrime_ -= 1.0e-10;
}
}
Double_t LauKinematics::calcSqDPJacobian() const
{
return this->calcSqDPJacobian(mPrime_,thetaPrime_);
}
Double_t LauKinematics::calcSqDPJacobian(const Double_t mPrime, const Double_t thPrime) const
{
// Calculate the Jacobian for the transformation
// m23^2, m13^2 -> m', theta' (square DP)
const Double_t m12 = 0.5*mDiff_[2]*(1.0 + TMath::Cos(LauConstants::pi*mPrime)) + mMin_[2];
const Double_t m12Sq = m12*m12;
const Double_t e1Cms12 = (m12Sq - m2Sq_ + m1Sq_)/(2.0*m12);
const Double_t e3Cms12 = (mParentSq_ - m12Sq - m3Sq_)/(2.0*m12);
const Double_t p1Cms12 = this->pCalc(e1Cms12, m1Sq_);
const Double_t p3Cms12 = this->pCalc(e3Cms12, m3Sq_);
const Double_t deriv1 = LauConstants::piBy2*mDiff_[2]*TMath::Sin(LauConstants::pi*mPrime);
const Double_t deriv2 = LauConstants::pi*TMath::Sin(LauConstants::pi*thPrime);
const Double_t jacobian = 4.0*p1Cms12*p3Cms12*m12*deriv1*deriv2;
return jacobian;
}
void LauKinematics::updateMassSquares(const Double_t m13Sq, const Double_t m23Sq)
{
m13Sq_ = m13Sq;
if (m13Sq_ > 0.0) {
m13_ = TMath::Sqrt(m13Sq_);
} else {
m13_ = 0.0;
}
m23Sq_ = m23Sq;
if (m23Sq_ > 0.0) {
m23_ = TMath::Sqrt(m23Sq_);
} else {
m23_ = 0.0;
}
// Now calculate m12Sq and m12.
this->calcm12Sq();
// Calculate momenta of tracks in parent (B, D etc.) rest-frame
this->calcParentFrameMomenta();
}
void LauKinematics::updateSqDPMassSquares(const Double_t mPrime, const Double_t thetaPrime)
{
// From square DP co-ordinates, calculate only the mass-squares
// First set the square-DP variables
mPrime_ = mPrime; thetaPrime_ = thetaPrime;
// Next calculate m12 and c12
Double_t m12 = 0.5*mDiff_[2]*(1.0 + TMath::Cos(LauConstants::pi*mPrime)) + mMin_[2];
Double_t c12 = TMath::Cos(LauConstants::pi*thetaPrime);
// From these calculate m13Sq and m23Sq
this->updateMassSq_m12(m12, c12);
// Calculate momenta of tracks in parent (B, D etc.) rest-frame
this->calcParentFrameMomenta();
}
void LauKinematics::calcm12Sq()
{
// Calculate m12Sq from m13Sq and m23Sq.
m12Sq_ = mParentSq_ + mSqDTot_ - m13Sq_ - m23Sq_;
// If m12Sq is too low, return lower limit,
// and modify m13Sq accordingly.
if (m12Sq_ < mSqMin_[2]) {
m12Sq_ = mSqMin_[2] + 1.0e-3;
m13Sq_ = mParentSq_ + mSqDTot_ - m12Sq_ - m23Sq_;
}
if (m12Sq_ > 0.0) {
m12_ = TMath::Sqrt(m12Sq_);
} else {
m12_ = 0.0;
}
}
void LauKinematics::calcParentFrameMomenta()
{
Double_t e1 = (mParentSq_ + m1Sq_ - m23Sq_) / (2.0*mParent_);
Double_t e2 = (mParentSq_ + m2Sq_ - m13Sq_) / (2.0*mParent_);
Double_t e3 = (mParentSq_ + m3Sq_ - m12Sq_) / (2.0*mParent_);
p1_Parent_ = this->pCalc(e1, m1Sq_);
p2_Parent_ = this->pCalc(e2, m2Sq_);
p3_Parent_ = this->pCalc(e3, m3Sq_);
}
void LauKinematics::calcHelicities()
{
// Calculate helicity angle cosines, given m12Sq, m13Sq and m23Sq.
// cij_ is the cosine of the helicity angle in the rest frame of the
// system of particles i and j.
// It is (but note the caveat below) the angle between tracks i and k
// in the ij rest frame with indices permuted cyclically.
// However, it is important to note that it is not exactly a cyclic
// permutation (there is a special treatment for c23 inside the cFromM
// function) for reasons of preserving the symmetry about m13=m23 for
// symmetric final states.
// The precise definitions are:
// theta12 is defined as the angle between 1&3 in the rest frame of 1&2
// theta23 is defined as the angle between 3&1 in the rest frame of 2&3
// theta13 is defined as the angle between 3&2 in the rest frame of 1&3
//
// It is a prerequisite that all mij_ and mijSq_ variables be correctly set.
Int_t zero(0), one(1), two(2);
c12_ = cFromM(m12Sq_, m13Sq_, m12_, zero, one, two);
p1_12_ = qi_; p3_12_ = qk_; // i, j = 12 (rest frame), k = 3
c23_ = cFromM(m23Sq_, m12Sq_, m23_, one, two, zero);
p2_23_ = qi_; p1_23_ = qk_; // i, j = 23 (rest frame), k = 1
c13_ = cFromM(m13Sq_, m23Sq_, m13_, two, zero, one);
p1_13_ = qi_; p2_13_ = qk_; // i, j = 31 (rest frame), k = 2
}
Double_t LauKinematics::cFromM(const Double_t mijSq, const Double_t mikSq, const Double_t mij, const Int_t i, const Int_t j, const Int_t k) const
{
// Routine to calculate the cos(helicity) variables from the masses of the particles.
// (See comment in LauKinematics::calcHelicities for futher information.)
Double_t EiCmsij = (mijSq - mSq_[j] + mSq_[i])/(2.0*mij);
Double_t EkCmsij = (mParentSq_ - mijSq - mSq_[k])/(2.0*mij);
if (EiCmsij < mass_[i]) {
if (warnings_) {
std::cerr<<"WARNING in LauKinematics::cFromM : EiCmsij = "<<EiCmsij<<" too small < mass_["<<i<<"] = "<<mass_[i]<<" in cFromM, i, j, k = "<<i<<", "<<j<<", "<<k<<std::endl;
std::cerr<<" : mijSq = "<<mijSq<<"; mij = "<<mij<<"; mSq_["<<j<<"] = "<<mSq_[j]<<"; mSq_["<<i<<"] = "<<mSq_[i]<<std::endl;
}
return 0.0;
}
if (EkCmsij < mass_[k]) {
if (warnings_) {
std::cerr<<"WARNING in LauKinematics::cFromM : EkCmsij = "<<EkCmsij<<" too small < mass_["<<k<<"] = "<<mass_[k]<<" in cFromM, i, j, k = "<<i<<", "<<j<<", "<<k<<std::endl;
std::cerr<<" : mijSq = "<<mijSq<<"; mij = "<<mij<<"; mSq_["<<j<<"] = "<<mSq_[j]<<"; mSq_["<<i<<"] = "<<mSq_[i]<<std::endl;
}
return 0.0;
}
// Find track i and k momenta in ij rest frame
qi_ = this->pCalc(EiCmsij, mSq_[i]);
qk_ = this->pCalc(EkCmsij, mSq_[k]);
// Find ij helicity angle
Double_t cosHel = -(mikSq - mSq_[i] - mSq_[k] - 2.0*EiCmsij*EkCmsij)/(2.0*qi_*qk_);
if (cosHel > 1.0) {
cosHel = 1.0;
} else if (cosHel < -1.0) {
cosHel = -1.0;
}
if (i == 1) {cosHel *= -1.0;}
return cosHel;
}
Double_t LauKinematics::mFromC(const Double_t mijSq, const Double_t cij, const Double_t mij, const Int_t i, const Int_t j, const Int_t k) const
{
// Returns the mass-squared for a pair of particles, i,j, given their
// invariant mass (squared) and the helicity angle.
// cij is helicity angle for the pair which is made from tracks i and j.
// It is the angle between tracks i and k in the ij rest frame.
Double_t cosHel( cij );
if (i == 1) {cosHel *= -1.0;}
Double_t EiCmsij = (mijSq - mSq_[j] + mSq_[i])/(2.0*mij);
Double_t EkCmsij = (mParentSq_ - mijSq - mSq_[k])/(2.0*mij);
if (TMath::Abs(EiCmsij - mass_[i]) > 1e-6 && EiCmsij < mass_[i]) {
if (warnings_) {
std::cerr<<"WARNING in LauKinematics::mFromC : EiCmsij = "<<EiCmsij<<" < "<<mass_[i]<<" in mFromC, i, j, k = "<<i<<", "<<j<<", "<<k<<std::endl;
}
return 0.0;
}
if (TMath::Abs(EkCmsij - mass_[k]) > 1e-6 && EkCmsij < mass_[k]) {
if (warnings_) {
std::cerr<<"WARNING in LauKinematics::mFromC : EkCmsij = "<<EkCmsij<<" < "<<mass_[k]<<" in mFromC, i, j, k = "<<i<<", "<<j<<", "<<k<<std::endl;
}
return 0.0;
}
// Find track i and k momenta in ij rest fram
Double_t qi = this->pCalc(EiCmsij, mSq_[i]);
Double_t qk = this->pCalc(EkCmsij, mSq_[k]);
// Find mikSq
Double_t massSq = mSq_[i] + mSq_[k] + 2.0*EiCmsij*EkCmsij - 2.0*qi*qk*cosHel;
if (massSq < mSqMin_[j]) {
if (warnings_) {
std::cerr<<"WARNING in LauKinematics::mFromC : mFromC below bound: i, j, k = "<<i<<", "<<j<<", "<<k<<std::endl;
}
massSq = mSqMin_[j];
}
return massSq;
}
void LauKinematics::genFlatPhaseSpace(Double_t& m13Sq, Double_t& m23Sq) const
{
// Routine to generate flat phase-space events.
// DP max kinematic boundaries in circumscribed box
// See DP figure in PDG book.
// m13max=mbrec-mass(2)
// m13sqmax=m13max*m13max
// m23max=mbrec-mass(1)
// m23sqmax=m23max*m23max
// Generate m13Sq and m23Sq flat within DP overall rectangular box
// Loop if the generated point is not within the DP kinematic boundary
do {
m13Sq = mSqMin_[1] + LauRandom::randomFun()->Rndm()*mSqDiff_[1];
m23Sq = mSqMin_[0] + LauRandom::randomFun()->Rndm()*mSqDiff_[0];
} while ( ! this->withinDPLimits( m13Sq, m23Sq ) );
}
void LauKinematics::genFlatSqDP(Double_t& mPrime, Double_t& thetaPrime) const
{
// Generate random event in the square Dalitz plot
mPrime = LauRandom::randomFun()->Rndm();
thetaPrime = LauRandom::randomFun()->Rndm();
}
Bool_t LauKinematics::withinDPLimits(const Double_t m13Sq, const Double_t m23Sq) const
{
// Find out whether the point (m13Sq,m23Sq) is within the limits of the
// Dalitz plot. The limits are specified by the invariant masses
// of the parent (e.g. B) and its three daughters that were
// defined in the constructor of this class. Here
// m_13Sq = square of invariant mass of daughters 1 and 3
// m_23Sq = square of invariant mass of daughters 2 and 3.
Bool_t withinDP = kFALSE;
// First check that m13Sq is within its absolute min and max
if (!((m13Sq > mSqMin_[1]) && (m13Sq < mSqMax_[1]))) {
return kFALSE;
}
// Now for the given value of m13Sq calculate the local min and max of m23Sq
Double_t m13 = TMath::Sqrt(m13Sq);
Double_t e3Cms13 = (m13Sq - m1Sq_ + m3Sq_)/(2.0*m13);
Double_t p3Cms13 = this->pCalc(e3Cms13, m3Sq_);
Double_t e2Cms13 = (mParentSq_ - m13Sq - m2Sq_)/(2.0*m13);
Double_t p2Cms13 = this->pCalc(e2Cms13, m2Sq_);
Double_t term = 2.0*e2Cms13*e3Cms13 + m2Sq_ + m3Sq_;
Double_t m23SqLocMin = term - 2.0*p2Cms13*p3Cms13;
Double_t m23SqLocMax = term + 2.0*p2Cms13*p3Cms13;
// Check whether the given value of m23Sq satisfies these bounds
if (m23Sq > m23SqLocMin && m23Sq < m23SqLocMax) {
withinDP = kTRUE;
}
return withinDP;
}
Bool_t LauKinematics::withinDPLimits2(const Double_t m13Sq, const Double_t m23Sq) const
{
// Same as withinDPLimits, but this time testing whether the m13Sq
// variable is within the kinematic boundary for the given m23Sq value
Bool_t withinDP = kFALSE;
// First check that m13Sq is within its absolute min and max
if (!((m23Sq > mSqMin_[0]) && (m23Sq < mSqMax_[0]))) {
return kFALSE;
}
// Now for the given value of m13Sq calculate the local min and max of m23Sq
Double_t m23 = TMath::Sqrt(m23Sq);
Double_t e3Cms23 = (m23Sq - m2Sq_ + m3Sq_)/(2.0*m23);
Double_t p3Cms23 = this->pCalc(e3Cms23, m3Sq_);
Double_t e1Cms23 = (mParentSq_ - m23Sq - m1Sq_)/(2.0*m23);
Double_t p1Cms23 = this->pCalc(e1Cms23, m1Sq_);
Double_t term = 2.0*e3Cms23*e1Cms23 + m1Sq_ + m3Sq_;
Double_t m13SqLocMin = term - 2.0*p3Cms23*p1Cms23;
Double_t m13SqLocMax = term + 2.0*p3Cms23*p1Cms23;
// Check whether the given value of m23Sq satisfies these bounds
if (m13Sq > m13SqLocMin && m13Sq < m13SqLocMax) {
withinDP = kTRUE;
}
return withinDP;
}
Bool_t LauKinematics::withinSqDPLimits(const Double_t mPrime, const Double_t thetaPrime) const
{
// Check whether m' and theta' are between 0 and 1
Bool_t withinDP(kFALSE);
if (mPrime > 0.0 && mPrime < 1.0 &&
thetaPrime > 0.0 && thetaPrime < 1.0) {
withinDP = kTRUE;
}
return withinDP;
}
Double_t LauKinematics::calcThirdMassSq(const Double_t firstMassSq, const Double_t secondMassSq) const
{
// Calculate one massSq from the other two
return mParentSq_ + mSqDTot_ - firstMassSq - secondMassSq;
}
Double_t LauKinematics::distanceFromDPCentre() const
{
return this->distanceFromDPCentre(m13Sq_,m23Sq_);
}
Double_t LauKinematics::distanceFromDPCentre(const Double_t m13Sq, const Double_t m23Sq) const
{
// DP centre is defined as the point where m12=m13=m23
// First find the m^2_ij value for the centre
Double_t centreMijSq = (mParentSq_ + m1Sq_ + m2Sq_ + m3Sq_)/3.0;
// Then find the difference between this and the two provided co-ords
Double_t diff13 = m13Sq - centreMijSq;
Double_t diff23 = m23Sq - centreMijSq;
// Calculate the total distance
Double_t distance = TMath::Sqrt(diff13*diff13 + diff23*diff23);
return distance;
}
Double_t LauKinematics::pCalc(const Double_t energy, const Double_t massSq) const
{
// Routine to calculate the momentum of a particle, given its energy and
// invariant mass (squared).
Double_t arg = energy*energy - massSq;
if (arg < 0.0) {
//if (warnings_) {
//std::cerr<<"WARNING in LauKinematics::pCalc : arg < 0.0: "<<arg<<" for e = "<<energy<<", mSq = "<<massSq<<std::endl;
//}
arg = 0.0;
}
Double_t pCalcVal = TMath::Sqrt(arg);
return pCalcVal;
}
void LauKinematics::flipAndUpdateKinematics()
{
// Flips the DP variables m13^2 <-> m23^2.
// Used in the case of symmetrical Dalitz plots (i.e. when two of the
// daughter tracks are the same type) within the
// LauIsobarDynamics::resAmp function.
this->updateKinematics(m23Sq_, m13Sq_);
}
void LauKinematics::rotateAndUpdateKinematics()
{
// Cyclically rotates the DP variables (m12 -> m23, m23 -> m13, m13 -> m12)
// Used in the case of fully symmetric Dalitz plots (i.e. when all
// three of the daughter tracks are the same type) within the
// LauIsobarDynamics::resAmp function.
this->updateKinematics(m12Sq_, m13Sq_);
}
void LauKinematics::updateMassSq_m23(const Double_t m23, const Double_t c23)
{
// Update the variables m12Sq_ and m13Sq_ given m23 and c23.
m23_ = m23; m23Sq_ = m23*m23; c23_ = c23;
const Int_t zero(0), one(1), two(2);
m12Sq_ = this->mFromC(m23Sq_, c23_, m23_, one, two, zero);
m13Sq_ = mParentSq_ - m23Sq_ - m12Sq_ + mSqDTot_;
m12_ = TMath::Sqrt(m12Sq_);
m13_ = TMath::Sqrt(m13Sq_);
}
void LauKinematics::updateMassSq_m13(const Double_t m13, const Double_t c13)
{
// Update the variables m12Sq_ and m23Sq_ given m13 and c13.
m13_ = m13; m13Sq_ = m13*m13; c13_ = c13;
const Int_t zero(0), one(1), two(2);
m23Sq_ = this->mFromC(m13Sq_, c13_, m13_, two, zero, one);
m12Sq_ = mParentSq_ - m23Sq_ - m13Sq_ + mSqDTot_;
m23_ = TMath::Sqrt(m23Sq_);
m12_ = TMath::Sqrt(m12Sq_);
}
void LauKinematics::updateMassSq_m12(const Double_t m12, const Double_t c12)
{
// Update the variables m23Sq_ and m13Sq_ given m12 and c12.
m12_ = m12; m12Sq_ = m12*m12; c12_ = c12;
const Int_t zero(0), one(1), two(2);
m13Sq_ = this->mFromC(m12Sq_, c12_, m12_, zero, one, two);
m23Sq_ = mParentSq_ - m12Sq_ - m13Sq_ + mSqDTot_;
m13_ = TMath::Sqrt(m13Sq_);
m23_ = TMath::Sqrt(m23Sq_);
}
void LauKinematics::updateKinematicsFrom23(const Double_t m23, const Double_t c23)
{
// Calculate the other mass squares
this->updateMassSq_m23(m23,c23);
// Calculate the remaining helicity angles
this->calcHelicities();
// Calculate momenta of tracks in parent (B, D etc.) rest-frame
this->calcParentFrameMomenta();
// Also calculate the m', theta' variables
if (squareDP_) {this->calcSqDPVars();}
}
void LauKinematics::updateKinematicsFrom13(const Double_t m13, const Double_t c13)
{
// Calculate the other mass squares
this->updateMassSq_m13(m13,c13);
// Calculate the remaining helicity angles
this->calcHelicities();
// Calculate momenta of tracks in parent (B, D etc.) rest-frame
this->calcParentFrameMomenta();
// Also calculate the m', theta' variables
if (squareDP_) {this->calcSqDPVars();}
}
void LauKinematics::updateKinematicsFrom12(const Double_t m12, const Double_t c12)
{
// Calculate the other mass squares
this->updateMassSq_m12(m12,c12);
// Calculate the remaining helicity angles
this->calcHelicities();
// Calculate momenta of tracks in parent (B, D etc.) rest-frame
this->calcParentFrameMomenta();
// Also calculate the m', theta' variables
if (squareDP_) {this->calcSqDPVars();}
}
Double_t LauKinematics::genm13Sq() const
{
Double_t m13Sq = mSqMin_[1] + LauRandom::randomFun()->Rndm()*mSqDiff_[1];
return m13Sq;
}
Double_t LauKinematics::genm23Sq() const
{
Double_t m23Sq = mSqMin_[0] + LauRandom::randomFun()->Rndm()*mSqDiff_[0];
return m23Sq;
}
Double_t LauKinematics::genm12Sq() const
{
Double_t m12Sq = mSqMin_[2] + LauRandom::randomFun()->Rndm()*mSqDiff_[2];
return m12Sq;
}
+Double_t dal(Double_t* x, Double_t* par)
+{
+ Double_t mParent = par[0];
+ Double_t mi = par[1];
+ Double_t mj = par[2];
+ Double_t mk = par[3];
+
+ Double_t mikSq=x[0];
+ Double_t mjkSq=x[1];
+ Double_t mik = TMath::Sqrt(mikSq);
+ Double_t mjk = TMath::Sqrt(mjkSq);
+
+ Double_t ejcmsik = (mParent*mParent-mj*mj-mik*mik)/(2.0*mik);
+ Double_t ekcmsik = (mik*mik+mk*mk-mi*mi)/(2.0*mik);
+ if (ekcmsik<mk || ejcmsik<mj) return 2.0;
+
+ Double_t pj = TMath::Sqrt(ejcmsik*ejcmsik-mj*mj);
+ Double_t pk = TMath::Sqrt(ekcmsik*ekcmsik-mk*mk);
+ Double_t coshelik = (mjk*mjk - mk*mk - mj*mj - 2.0*ejcmsik*ekcmsik)/(2.0*pj*pk);
+
+ Double_t coshelikSq = coshelik*coshelik;
+ return coshelikSq;
+}
void LauKinematics::drawDPContour(Int_t orientation, Int_t nbins)
{
// orientation -
// 1323 : x = m13, y = m23
// etc.
Double_t m1 = this->getm1();
Double_t m2 = this->getm2();
Double_t m3 = this->getm3();
Double_t mParent = this->getmParent();
Double_t m13SqMin = this->getm13SqMin();
Double_t m23SqMin = this->getm23SqMin();
Double_t m12SqMin = this->getm12SqMin();
Double_t m13SqMax = this->getm13SqMax();
Double_t m23SqMax = this->getm23SqMax();
Double_t m12SqMax = this->getm12SqMax();
Double_t xMin(0.0);
Double_t xMax(mParent*mParent);
Double_t yMin(0.0);
Double_t yMax(mParent*mParent);
if (orientation == 1323) {
xMin = m13SqMin-1.0; xMax = m13SqMax+1.0;
yMin = m23SqMin-1.0; yMax = m23SqMax+1.0;
} else if (orientation == 2313) {
xMin = m23SqMin-1.0; xMax = m23SqMax+1.0;
yMin = m13SqMin-1.0; yMax = m13SqMax+1.0;
} else if (orientation == 1213) {
xMin = m12SqMin-1.0; xMax = m12SqMax+1.0;
yMin = m13SqMin-1.0; yMax = m13SqMax+1.0;
} else if (orientation == 1312) {
xMin = m13SqMin-1.0; xMax = m13SqMax+1.0;
yMin = m12SqMin-1.0; yMax = m12SqMax+1.0;
} else if (orientation == 1223) {
xMin = m12SqMin-1.0; xMax = m12SqMax+1.0;
yMin = m23SqMin-1.0; yMax = m23SqMax+1.0;
} else if (orientation == 2312) {
xMin = m23SqMin-1.0; xMax = m23SqMax+1.0;
yMin = m12SqMin-1.0; yMax = m12SqMax+1.0;
} else {
std::cerr<<"ERROR in LauKinematics::drawDPContour : Unrecognised orientation, not drawing contour."<<std::endl;
return;
}
Int_t npar = 4;
TF2 * f2 = new TF2("contour", &dal, xMin, xMax, yMin, yMax, npar);
// Set the parameters
f2->SetParameter(0,mParent);
if (orientation == 1323) {
f2->SetParameter(1,m1);
f2->SetParameter(2,m2);
f2->SetParameter(3,m3);
} else if (orientation == 2313) {
f2->SetParameter(1,m2);
f2->SetParameter(2,m1);
f2->SetParameter(3,m3);
} else if (orientation == 1213) {
f2->SetParameter(1,m2);
f2->SetParameter(2,m3);
f2->SetParameter(3,m1);
} else if (orientation == 1312) {
f2->SetParameter(1,m3);
f2->SetParameter(2,m2);
f2->SetParameter(3,m1);
} else if (orientation == 1223) {
f2->SetParameter(1,m1);
f2->SetParameter(2,m3);
f2->SetParameter(3,m2);
} else if (orientation == 2312) {
f2->SetParameter(1,m3);
f2->SetParameter(2,m1);
f2->SetParameter(3,m2);
}
// Set up the contour to be drawn when the value of the function == 1.0
Double_t b[]={1.0};
f2->SetContour(1,b);
// Set the number of bins for the contour to be sampled over
f2->SetNpx(nbins);
f2->SetNpy(nbins);
// and the line style
f2->SetLineWidth(3);
f2->SetLineStyle(kSolid);
// Draw the contour on top of the histo that should already have been drawn
f2->DrawCopy("same");
delete f2;
}
-Double_t dal(Double_t* x, Double_t* par)
-{
- Double_t mParent = par[0];
- Double_t mi = par[1];
- Double_t mj = par[2];
- Double_t mk = par[3];
-
- Double_t mikSq=x[0];
- Double_t mjkSq=x[1];
- Double_t mik = TMath::Sqrt(mikSq);
- Double_t mjk = TMath::Sqrt(mjkSq);
-
- Double_t ejcmsik = (mParent*mParent-mj*mj-mik*mik)/(2.0*mik);
- Double_t ekcmsik = (mik*mik+mk*mk-mi*mi)/(2.0*mik);
- if (ekcmsik<mk || ejcmsik<mj) return 2.0;
-
- Double_t pj = TMath::Sqrt(ejcmsik*ejcmsik-mj*mj);
- Double_t pk = TMath::Sqrt(ekcmsik*ekcmsik-mk*mk);
- Double_t coshelik = (mjk*mjk - mk*mk - mj*mj - 2.0*ejcmsik*ekcmsik)/(2.0*pj*pk);
-
- Double_t coshelikSq = coshelik*coshelik;
- return coshelikSq;
-}
-
diff --git a/src/LauLinearPdf.cc b/src/LauLinearPdf.cc
index 5ae16bc..b0904c0 100644
--- a/src/LauLinearPdf.cc
+++ b/src/LauLinearPdf.cc
@@ -1,129 +1,128 @@
/*
Copyright 2006 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauLinearPdf.cc
\brief File containing implementation of LauLinearPdf class.
*/
#include <iostream>
#include <vector>
using std::cout;
using std::cerr;
using std::endl;
-using std::vector;
#include "TMath.h"
#include "TSystem.h"
#include "LauConstants.hh"
#include "LauLinearPdf.hh"
ClassImp(LauLinearPdf)
- LauLinearPdf::LauLinearPdf(const TString& theVarName, const vector<LauAbsRValue*>& params, Double_t minAbscissa, Double_t maxAbscissa) :
+ LauLinearPdf::LauLinearPdf(const TString& theVarName, const std::vector<LauAbsRValue*>& params, Double_t minAbscissa, Double_t maxAbscissa) :
LauAbsPdf(theVarName, params, minAbscissa, maxAbscissa),
slope_(0),posflag_(kTRUE)
{
// Constructor for the linear PDF.
//
// The parameters in params are the slope and y-intercept of the line.
// The last two arguments specify the range in which the PDF is defined, and the PDF
// will be normalised w.r.t. these limits.
slope_ = this->findParameter("slope");
if ((this->nParameters() != 1) || (slope_ == 0)) {
cerr<<"Warning. LauLinearPdf requires 1 parameter: \"slope\"."<<endl;
gSystem->Exit(EXIT_FAILURE);
}
// Cache the normalisation factor.
this->calcNorm();
}
LauLinearPdf::~LauLinearPdf()
{
// Destructor
}
void LauLinearPdf::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];
// Get the up to date parameter values
Double_t slope = slope_->unblindValue();
// Calculate the value of the straight line for the given value of the abscissa.
Double_t constVal = 1.0/(this->getMaxAbscissa() - this->getMinAbscissa());
constVal -= slope * (this->getMaxAbscissa() + this->getMinAbscissa()) * 0.5;
Double_t value = slope*abscissa + constVal;
// Make sure the PDF doesn't go negative
if ( value < 0.0 ) {
if ( posflag_ ) {
std::cerr << "WARNING in LauLinearPdf::calcLikelihoodInfo : The PDF is negative, setting to zero" << std::endl;
+ posflag_= kFALSE;
}
value = 0.0;
- posflag_= kFALSE;
}
this->setUnNormPDFVal(value);
}
void LauLinearPdf::calcNorm()
{
// Nothing to calculate here since the PDF is already normalised to 1
this->setNorm(1.0);
}
void LauLinearPdf::calcPDFHeight( const LauKinematics* /*kinematics*/ )
{
if (this->heightUpToDate()) {
return;
}
// Get the up to date parameter values
Double_t slope = slope_->unblindValue();
// Calculate the PDF height for the straight line
LauAbscissas maxPoint(1);
if (slope>0.0) {
maxPoint[0] = this->getMaxAbscissa();
} else {
maxPoint[0] = this->getMinAbscissa();
}
this->calcLikelihoodInfo(maxPoint);
Double_t height = this->getUnNormLikelihood();
this->setMaxHeight(height);
}
diff --git a/src/LauNovosibirskPdf.cc b/src/LauNovosibirskPdf.cc
index 3e363e4..ab191eb 100644
--- a/src/LauNovosibirskPdf.cc
+++ b/src/LauNovosibirskPdf.cc
@@ -1,152 +1,152 @@
/*
Copyright 2008 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauNovosibirskPdf.cc
\brief File containing implementation of LauNovosibirskPdf class.
*/
#include <iostream>
using std::cout;
using std::cerr;
using std::endl;
#include "TMath.h"
#include "TSystem.h"
#include "LauNovosibirskPdf.hh"
#include "LauConstants.hh"
ClassImp(LauNovosibirskPdf)
-LauNovosibirskPdf::LauNovosibirskPdf(const TString& theVarName, const vector<LauAbsRValue*>& params, Double_t minAbscissa, Double_t maxAbscissa) :
+LauNovosibirskPdf::LauNovosibirskPdf(const TString& theVarName, const std::vector<LauAbsRValue*>& params, Double_t minAbscissa, Double_t maxAbscissa) :
LauAbsPdf(theVarName, params, minAbscissa, maxAbscissa),
mean_(0),
sigma_(0),
tail_(0)
{
// Constructor for the Novosibirsk PDF.
//
// The parameters in params are the mean, sigma and tail.
// The last two arguments specify the range in which the PDF is defined, and the PDF
// will be normalised w.r.t. these limits.
mean_ = this->findParameter("mean");
sigma_ = this->findParameter("sigma");
tail_ = this->findParameter("tail");
if ((this->nParameters() != 3) || (mean_ == 0) || (sigma_ == 0) || (tail_ == 0)) {
cerr<<"ERROR in LauNovosibirskPdf constructor: LauNovosibirskPdf requires 3 parameters: \"mean\", \"sigma\"and \"tail\" "<<endl;
gSystem->Exit(EXIT_FAILURE);
}
// Cache the normalisation factor
this->calcNorm();
}
LauNovosibirskPdf::~LauNovosibirskPdf()
{
// Destructor
}
void LauNovosibirskPdf::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];
// Get the up to date parameter values
Double_t mean = mean_->unblindValue();
Double_t sigma = sigma_->unblindValue();
Double_t tail = tail_->unblindValue();
// Evaluate the Novosibirsk PDF value
Double_t qa(0.0), qb(0.0), qx(0.0), qy(0.0);
Double_t arg(0.0);
Double_t value(0.0);
if(TMath::Abs(tail) < 1.e-7)
arg = 0.5*((abscissa - mean)/sigma)*((abscissa - mean)/sigma);
else {
qa = tail*TMath::Sqrt(LauConstants::log4);
qb = TMath::SinH(qa)/qa;
qx = (abscissa - mean)/sigma*qb;
qy = 1.0+ tail*qx;
//---- Cutting curve from right side
if( qy > 1.E-7) {
arg = 0.5*( (log(qy)/tail)*(log(qy)/tail) + tail*tail);
}else{
arg = 15.0;
}
}
value = TMath::Exp(-arg);
// if the parameters are floating then we
// need to recalculate the normalisation
if (!this->cachePDF() && !this->withinNormCalc() && !this->withinGeneration()) {
this->calcNorm();
}
this->setUnNormPDFVal(value);
}
void LauNovosibirskPdf::calcPDFHeight( const LauKinematics* /*kinematics*/ )
{
if (this->heightUpToDate()) {
return;
}
// Get the up to date parameter values
Double_t mean = mean_->unblindValue();
LauAbscissas maxPoint(1);
maxPoint[0] = mean;
// Calculate the PDF height for the Bifurcated Gaussian function.
if (mean < this->getMinAbscissa()) {
maxPoint[0] = this->getMinAbscissa();
} else if (mean > this->getMaxAbscissa()) {
maxPoint[0] = this->getMaxAbscissa();
}
this->calcLikelihoodInfo(maxPoint);
Double_t height = this->getUnNormLikelihood();
// Multiply by a small factor to avoid problems from rounding errors
height *= (1.0 + 1e-1);
this->setMaxHeight(height);
}
diff --git a/src/LauParameter.cc b/src/LauParameter.cc
index 535e017..ab7b48a 100644
--- a/src/LauParameter.cc
+++ b/src/LauParameter.cc
@@ -1,762 +1,761 @@
/*
Copyright 2006 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauParameter.cc
\brief File containing implementation of LauParameter class.
*/
#include <iostream>
#include <map>
using std::cout;
using std::cerr;
using std::endl;
-using std::map;
#include "TRandom.h"
#include "LauParameter.hh"
#include "LauRandom.hh"
ClassImp(LauParameter)
LauParameter::LauParameter() :
name_(""),
value_(0.0),
error_(0.0),
negError_(0.0),
posError_(0.0),
genValue_(0.0),
initValue_(0.0),
minValue_(0.0),
maxValue_(0.0),
fixed_(kTRUE),
secondStage_(kFALSE),
gaussConstraint_(kFALSE),
constraintTrueMean_(0.0),
constraintMean_(0.0),
constraintWidth_(0.0),
gcc_(0.0),
bias_(0.0),
pull_(0.0),
clone_(kFALSE),
parent_(0),
blinder_(0)
{
}
LauParameter::LauParameter(const TString& parName) :
name_(parName),
value_(0.0),
error_(0.0),
negError_(0.0),
posError_(0.0),
genValue_(0.0),
initValue_(0.0),
minValue_(0.0),
maxValue_(0.0),
fixed_(kTRUE),
secondStage_(kFALSE),
gaussConstraint_(kFALSE),
constraintTrueMean_(0.0),
constraintMean_(0.0),
constraintWidth_(0.0),
gcc_(0.0),
bias_(0.0),
pull_(0.0),
clone_(kFALSE),
parent_(0),
blinder_(0)
{
}
LauParameter::LauParameter(Double_t parValue) :
name_(""),
value_(parValue),
error_(0.0),
negError_(0.0),
posError_(0.0),
genValue_(parValue),
initValue_(parValue),
minValue_(parValue-1e-6),
maxValue_(parValue+1e-6),
fixed_(kTRUE),
secondStage_(kFALSE),
gaussConstraint_(kFALSE),
constraintTrueMean_(0.0),
constraintMean_(0.0),
constraintWidth_(0.0),
gcc_(0.0),
bias_(0.0),
pull_(0.0),
clone_(kFALSE),
parent_(0),
blinder_(0)
{
}
LauParameter::LauParameter(const TString& parName, Double_t parValue) :
name_(parName),
value_(parValue),
error_(0.0),
negError_(0.0),
posError_(0.0),
genValue_(parValue),
initValue_(parValue),
minValue_(parValue-1e-6),
maxValue_(parValue+1e-6),
fixed_(kTRUE),
secondStage_(kFALSE),
gaussConstraint_(kFALSE),
constraintTrueMean_(0.0),
constraintMean_(0.0),
constraintWidth_(0.0),
gcc_(0.0),
bias_(0.0),
pull_(0.0),
clone_(kFALSE),
parent_(0),
blinder_(0)
{
}
LauParameter::LauParameter(Double_t parValue, Double_t min, Double_t max) :
name_(""),
value_(parValue),
error_(0.0),
negError_(0.0),
posError_(0.0),
genValue_(parValue),
initValue_(parValue),
minValue_(min),
maxValue_(max),
fixed_(kTRUE),
secondStage_(kFALSE),
gaussConstraint_(kFALSE),
constraintTrueMean_(0.0),
constraintMean_(0.0),
constraintWidth_(0.0),
gcc_(0.0),
bias_(0.0),
pull_(0.0),
clone_(kFALSE),
parent_(0),
blinder_(0)
{
this->checkRange();
}
LauParameter::LauParameter(Double_t parValue, Double_t parError, Double_t min, Double_t max) :
name_(""),
value_(parValue),
error_(parError),
negError_(0.0),
posError_(0.0),
genValue_(parValue),
initValue_(parValue),
minValue_(min),
maxValue_(max),
fixed_(kTRUE),
secondStage_(kFALSE),
gaussConstraint_(kFALSE),
constraintTrueMean_(0.0),
constraintMean_(0.0),
constraintWidth_(0.0),
gcc_(0.0),
bias_(0.0),
pull_(0.0),
clone_(kFALSE),
parent_(0),
blinder_(0)
{
this->checkRange();
}
LauParameter::LauParameter(Double_t parValue, Double_t min, Double_t max, Bool_t parFixed) :
name_(""),
value_(parValue),
error_(0.0),
negError_(0.0),
posError_(0.0),
genValue_(parValue),
initValue_(parValue),
minValue_(min),
maxValue_(max),
fixed_(parFixed),
secondStage_(kFALSE),
gaussConstraint_(kFALSE),
constraintTrueMean_(0.0),
constraintMean_(0.0),
constraintWidth_(0.0),
gcc_(0.0),
bias_(0.0),
pull_(0.0),
clone_(kFALSE),
parent_(0),
blinder_(0)
{
this->checkRange();
}
LauParameter::LauParameter(const TString& parName, Double_t parValue, Double_t min, Double_t max) :
name_(parName),
value_(parValue),
error_(0.0),
negError_(0.0),
posError_(0.0),
genValue_(parValue),
initValue_(parValue),
minValue_(min),
maxValue_(max),
fixed_(kTRUE),
secondStage_(kFALSE),
gaussConstraint_(kFALSE),
constraintTrueMean_(0.0),
constraintMean_(0.0),
constraintWidth_(0.0),
gcc_(0.0),
bias_(0.0),
pull_(0.0),
clone_(kFALSE),
parent_(0),
blinder_(0)
{
this->checkRange();
}
LauParameter::LauParameter(const TString& parName, Double_t parValue, Double_t min, Double_t max, Bool_t parFixed) :
name_(parName),
value_(parValue),
error_(0.0),
negError_(0.0),
posError_(0.0),
genValue_(parValue),
initValue_(parValue),
minValue_(min),
maxValue_(max),
fixed_(parFixed),
secondStage_(kFALSE),
gaussConstraint_(kFALSE),
constraintTrueMean_(0.0),
constraintMean_(0.0),
constraintWidth_(0.0),
gcc_(0.0),
bias_(0.0),
pull_(0.0),
clone_(kFALSE),
parent_(0),
blinder_(0)
{
this->checkRange();
}
LauParameter::LauParameter(const TString& parName, Double_t parValue, Double_t parError, Double_t min, Double_t max) :
name_(parName),
value_(parValue),
error_(parError),
negError_(0.0),
posError_(0.0),
genValue_(parValue),
initValue_(parValue),
minValue_(min),
maxValue_(max),
fixed_(kTRUE),
secondStage_(kFALSE),
gaussConstraint_(kFALSE),
constraintTrueMean_(0.0),
constraintMean_(0.0),
constraintWidth_(0.0),
gcc_(0.0),
bias_(0.0),
pull_(0.0),
clone_(kFALSE),
parent_(0),
blinder_(0)
{
this->checkRange();
}
LauParameter::LauParameter(const LauParameter& rhs) : TObject(rhs), LauAbsRValue(rhs),
name_(rhs.name_),
value_(rhs.value_),
error_(rhs.error_),
negError_(rhs.negError_),
posError_(rhs.posError_),
genValue_(rhs.genValue_),
initValue_(rhs.initValue_),
minValue_(rhs.minValue_),
maxValue_(rhs.maxValue_),
fixed_(rhs.fixed_),
secondStage_(rhs.secondStage_),
gaussConstraint_(rhs.gaussConstraint_),
constraintTrueMean_(rhs.constraintTrueMean_),
constraintMean_(rhs.constraintMean_),
constraintWidth_(rhs.constraintWidth_),
gcc_(rhs.gcc_),
bias_(rhs.bias_),
pull_(rhs.pull_),
clone_(rhs.clone_),
parent_(rhs.parent_),
clones_(rhs.clones_),
blinder_((rhs.blinder_==0) ? 0 : new LauBlind(*(rhs.blinder_)))
{
}
LauParameter& LauParameter::operator=(const LauParameter& rhs)
{
if (&rhs != this) {
TObject::operator=(rhs);
LauAbsRValue::operator=(rhs);
name_ = rhs.name_;
value_ = rhs.value_;
error_ = rhs.error_;
negError_ = rhs.negError_;
posError_ = rhs.posError_;
genValue_ = rhs.genValue_;
initValue_ = rhs.initValue_;
minValue_ = rhs.minValue_;
maxValue_ = rhs.maxValue_;
fixed_ = rhs.fixed_;
secondStage_ = rhs.secondStage_;
gaussConstraint_ = rhs.gaussConstraint_;
constraintTrueMean_ = rhs.constraintTrueMean_;
constraintMean_ = rhs.constraintMean_;
constraintWidth_ = rhs.constraintWidth_;
gcc_ = rhs.gcc_;
bias_ = rhs.bias_;
pull_ = rhs.pull_;
clone_ = rhs.clone_;
parent_ = rhs.parent_;
clones_ = rhs.clones_;
delete blinder_;
blinder_ = (rhs.blinder_==0) ? 0 : new LauBlind(*(rhs.blinder_));
}
return *this;
}
LauParameter::~LauParameter()
{
delete blinder_;
}
std::vector<LauParameter*> LauParameter::getPars()
{
std::vector<LauParameter*> list;
list.push_back(this);
return list;
}
void LauParameter::value(Double_t newValue)
{
if (this->clone()) {
parent_->value(newValue);
} else {
this->checkRange(newValue,this->minValue(),this->maxValue());
this->updateClones(kTRUE);
}
}
void LauParameter::error(Double_t newError)
{
if (this->clone()) {
parent_->error(newError);
} else {
error_ = TMath::Abs(newError);
this->updateClones(kFALSE);
}
}
void LauParameter::negError(Double_t newNegError)
{
if (this->clone()) {
parent_->negError(newNegError);
} else {
negError_ = TMath::Abs(newNegError);
this->updateClones(kFALSE);
}
}
void LauParameter::posError(Double_t newPosError)
{
if (this->clone()) {
parent_->posError(newPosError);
} else {
posError_ = TMath::Abs(newPosError);
this->updateClones(kFALSE);
}
}
void LauParameter::errors(Double_t newError, Double_t newNegError, Double_t newPosError)
{
if (this->clone()) {
parent_->errors(newError,newNegError,newPosError);
} else {
error_ = TMath::Abs(newError);
negError_ = TMath::Abs(newNegError);
posError_ = TMath::Abs(newPosError);
this->updateClones(kFALSE);
}
}
void LauParameter::valueAndErrors(Double_t newValue, Double_t newError, Double_t newNegError, Double_t newPosError)
{
if (this->clone()) {
parent_->valueAndErrors(newValue,newError,newNegError,newPosError);
} else {
this->checkRange(newValue,this->minValue(),this->maxValue());
error_ = TMath::Abs(newError);
negError_ = TMath::Abs(newNegError);
posError_ = TMath::Abs(newPosError);
this->updateClones(kFALSE);
}
}
void LauParameter::globalCorrelationCoeff(Double_t newGCCValue)
{
if (this->clone()) {
parent_->globalCorrelationCoeff(newGCCValue);
} else {
gcc_ = newGCCValue;
this->updateClones(kFALSE);
}
}
void LauParameter::genValue(Double_t newGenValue)
{
if (this->clone()) {
parent_->genValue(newGenValue);
} else {
genValue_ = newGenValue;
this->updateClones(kFALSE);
}
}
void LauParameter::initValue(Double_t newInitValue)
{
if (this->clone()) {
parent_->initValue(newInitValue);
} else {
initValue_ = newInitValue;
this->updateClones(kFALSE);
}
}
void LauParameter::minValue(Double_t newMinValue)
{
if (this->clone()) {
parent_->minValue(newMinValue);
} else {
this->checkRange(this->value(),newMinValue,this->maxValue());
this->updateClones(kFALSE);
}
}
void LauParameter::maxValue(Double_t newMaxValue)
{
if (this->clone()) {
parent_->maxValue(newMaxValue);
} else {
this->checkRange(this->value(),this->minValue(),newMaxValue);
this->updateClones(kFALSE);
}
}
void LauParameter::range(Double_t newMinValue, Double_t newMaxValue)
{
if (this->clone()) {
parent_->range(newMinValue,newMaxValue);
} else {
this->checkRange(this->value(),newMinValue,newMaxValue);
this->updateClones(kFALSE);
}
}
void LauParameter::valueAndRange(Double_t newValue, Double_t newMinValue, Double_t newMaxValue)
{
if (this->clone()) {
parent_->valueAndRange(newValue,newMinValue,newMaxValue);
} else {
this->checkRange(newValue,newMinValue,newMaxValue);
this->updateClones(kFALSE);
}
}
void LauParameter::name(const TString& newName)
{
// no need to update clones here
// clones are allowed to have different names
name_ = newName;
}
void LauParameter::fixed(Bool_t parFixed)
{
if (this->clone()) {
parent_->fixed(parFixed);
} else {
fixed_ = parFixed;
this->updateClones(kFALSE);
}
}
void LauParameter::secondStage(Bool_t secondStagePar)
{
if (this->clone()) {
parent_->secondStage(secondStagePar);
} else {
secondStage_ = secondStagePar;
this->updateClones(kFALSE);
}
}
void LauParameter::addGaussianConstraint(Double_t newGaussMean, Double_t newGaussWidth)
{
if (this->clone()) {
parent_->addGaussianConstraint(newGaussMean,newGaussWidth);
} else {
gaussConstraint_ = kTRUE;
constraintTrueMean_ = newGaussMean;
constraintMean_ = newGaussMean;
constraintWidth_ = newGaussWidth;
this->updateClones(kFALSE);
}
}
void LauParameter::removeGaussianConstraint()
{
if (this->clone()) {
parent_->removeGaussianConstraint();
} else {
gaussConstraint_ = kFALSE;
this->updateClones(kFALSE);
}
}
void LauParameter::generateConstraintMean()
{
if (this->clone()) {
parent_->generateConstraintMean();
} else {
constraintMean_ = LauRandom::randomFun()->Gaus( constraintTrueMean_, constraintWidth_ );
this->updateClones(kFALSE);
}
}
Double_t LauParameter::constraintPenalty() const
{
const Double_t val { this->unblindValue() };
const Double_t diff { val - constraintMean_ };
const Double_t term { diff * diff };
return term / ( 2.0 * constraintWidth_ * constraintWidth_ );
}
void LauParameter::blindParameter(const TString& blindingString, const Double_t width)
{
if (this->clone()) {
parent_->blindParameter(blindingString,width);
return;
}
if ( blinder_ != 0 ) {
std::cerr << "WARNING in LauParameter::blindParameter : blinding has already been set up for this parameter" << std::endl;
return;
}
blinder_ = new LauBlind(blindingString,width);
- for (map<LauParameter*,Double_t>::iterator iter = clones_.begin(); iter != clones_.end(); ++iter) {
+ for (std::map<LauParameter*,Double_t>::iterator iter = clones_.begin(); iter != clones_.end(); ++iter) {
LauParameter* clonePar = iter->first;
if ( clonePar->blinder_ != 0 ) {
std::cerr << "WARNING in LauParameter::blindParameter : blinding has already been set up for a clone of this parameter - it will be replaced!" << std::endl;
delete clonePar->blinder_;
clonePar->blinder_ = 0;
}
clonePar->blinder_ = new LauBlind(*blinder_);
}
}
void LauParameter::updatePull()
{
if (this->clone()) {
parent_->updatePull();
return;
}
// calculate the bias
bias_ = value_ - genValue_;
// if we have errors calculated then calculate
// the pull using the best error available
if ((bias_ > 0.0) && (negError_ > 1e-10)) {
pull_ = bias_ / negError_;
} else if ((bias_ < 0.0) && (posError_ > 1e-10)) {
pull_ = bias_ / posError_;
} else if (error_ > 1e-10) {
pull_ = bias_ / error_;
} else {
pull_ = 0.0;
}
this->updateClones(kFALSE);
}
void LauParameter::checkRange(Double_t val, Double_t minVal, Double_t maxVal)
{
// first check that min is less than max (or they are the same - this is allowed)
if (minVal > maxVal) {
cerr<<"ERROR in LauParameter::checkRange : minValue: "<<minVal<<" greater than maxValue: "<<maxVal<<endl;
if (minValue_ > maxValue_) {
minValue_ = maxValue_;
cerr<<" : Setting both to "<<maxValue_<<"."<<endl;
} else {
cerr<<" : Not changing anything."<<endl;
}
return;
}
minValue_ = minVal;
maxValue_ = maxVal;
// now check that the value is still within the range
if ((val < minVal) || (val > maxVal)) {
if (name_ != "") {
cerr<<"ERROR in LauParameter::checkRange : value: "<<val<<" not in allowed range ["<<minVal<<" , "<<maxVal<<"] for parameter \""<<name_<<"\"."<<endl;
} else {
cerr<<"ERROR in LauParameter::checkRange : value: "<<val<<" not in allowed range ["<<minVal<<" , "<<maxVal<<"]."<<endl;
}
if (val < minVal) {
cerr<<" : Setting value to minValue: "<<minVal<<endl;
val = minVal;
} else {
cerr<<" : Setting value to maxValue: "<<maxVal<<endl;
val = maxVal;
}
}
value_ = val;
}
LauParameter* LauParameter::createClone(Double_t constFactor)
{
// if we're a clone we mustn't be cloned ourselves
// but instead return another clone of our parent
// (this is so that the parent knows of all its clones)
if (this->clone()) {
LauParameter* clonePar = parent_->createClone(constFactor);
clonePar->name(this->name());
return clonePar;
}
// clone ourselves using the copy-constructor
LauParameter* clonePar = new LauParameter(*this);
Double_t newValue = clonePar->value() * constFactor;
clonePar->value( newValue );
clonePar->wipeClones();
clonePar->clone(this);
clones_.insert( std::make_pair( clonePar, constFactor ) );
return clonePar;
}
LauParameter* LauParameter::createClone(const TString& newName, Double_t constFactor)
{
// self message to create the clone
LauParameter* clonePar = this->createClone(constFactor);
// set the new name
clonePar->name(newName);
// and return
return clonePar;
}
void LauParameter::updateClones(Bool_t justValue)
{
// if we don't have any clones then there's nothing to do
if ( clones_.empty() ) {
return;
}
// we have to set the values directly rather than using member functions because otherwise we'd get into an infinite loop
if (justValue) {
for ( auto& [ clonePar, constFactor ] : clones_ ) {
clonePar->value_ = constFactor * value_;
}
} else {
for ( auto& [ clonePar, constFactor ] : clones_ ) {
clonePar->value_ = constFactor * value_;
clonePar->error_ = constFactor * error_;
clonePar->negError_ = constFactor * negError_;
clonePar->posError_ = constFactor * posError_;
clonePar->genValue_ = constFactor * genValue_;
clonePar->initValue_ = constFactor * initValue_;
clonePar->minValue_ = constFactor * minValue_;
clonePar->maxValue_ = constFactor * maxValue_;
clonePar->fixed_ = fixed_;
clonePar->secondStage_ = secondStage_;
clonePar->gaussConstraint_ = gaussConstraint_;
clonePar->constraintTrueMean_ = constraintTrueMean_;
clonePar->constraintMean_ = constraintMean_;
clonePar->constraintWidth_ = constraintWidth_;
clonePar->gcc_ = gcc_;
clonePar->bias_ = bias_;
clonePar->pull_ = pull_;
}
}
}
void LauParameter::randomiseValue()
{
this->randomiseValue(this->minValue(), this->maxValue());
}
void LauParameter::randomiseValue(Double_t minVal, Double_t maxVal)
{
// if we're fixed then do nothing
if (this->fixed()) {
return;
}
// check supplied values are sensible
if (maxVal < minVal) {
cerr<<"ERROR in LauParameter::randomiseValue : Supplied maximum value smaller than minimum value."<<endl;
return;
}
if (maxVal > this->maxValue()) {
maxVal = this->maxValue();
}
if (minVal < this->minValue()) {
minVal = this->minValue();
}
// use the zero-seed random number generator to get values that are
// uniformly distributed over the given range
Double_t randNo = LauRandom::zeroSeedRandom()->Rndm();
Double_t val = randNo*(maxVal - minVal) + minVal;
this->initValue(val);
}
// ostream operator
std::ostream& operator << (std::ostream& stream, const LauParameter& par)
{
stream << par.value();
return stream;
}
diff --git a/src/LauParametricStepFuncPdf.cc b/src/LauParametricStepFuncPdf.cc
index 3c9e1ed..19752c6 100644
--- a/src/LauParametricStepFuncPdf.cc
+++ b/src/LauParametricStepFuncPdf.cc
@@ -1,231 +1,230 @@
/*
Copyright 2008 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauParametricStepFuncPdf.cc
\brief File containing implementation of LauParametricStepFuncPdf class.
*/
/*****************************************************************************
* Class based on RooFit/RooParametricStepFunction. *
* Original copyright given below. *
*****************************************************************************
* Authors: *
* WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
* DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
* *
* Copyright (c) 2000-2005, Regents of the University of California *
* and Stanford University. All rights reserved. *
* *
* Redistribution and use in source and binary forms, *
* with or without modification, are permitted according to the terms *
* listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
*****************************************************************************/
#include <iostream>
#include <vector>
using std::cout;
using std::cerr;
using std::endl;
-using std::vector;
#include "TMath.h"
#include "TSystem.h"
#include "LauConstants.hh"
#include "LauParametricStepFuncPdf.hh"
ClassImp(LauParametricStepFuncPdf)
-LauParametricStepFuncPdf::LauParametricStepFuncPdf(const TString& theVarName, const vector<LauAbsRValue*>& params, const vector<Double_t>& limits, NormBin normalisationBin) :
+LauParametricStepFuncPdf::LauParametricStepFuncPdf(const TString& theVarName, const std::vector<LauAbsRValue*>& params, const std::vector<Double_t>& limits, NormBin normalisationBin) :
LauAbsPdf(theVarName, params, limits.front(), limits.back()),
normBin_(normalisationBin),
limits_(limits)
{
// Constructor for the PSF PDF.
//
// The parameters in params are the bin contents of all but the
// normalisation bin, so has N_bins-1 entries.
// The last argument specifies the limits of the bins and the range
// as a whole, so has N_bins+1 entries.
if (this->nParameters() != this->nBins()-1) {
cerr<<"ERROR in LauParametricStepFuncPdf constructor: LauParametricStepFuncPdf requires N-1 parameters, where N is the number of bins."<<endl;
gSystem->Exit(EXIT_FAILURE);
}
// Cache the normalisation factor
this->calcNorm();
}
LauParametricStepFuncPdf::~LauParametricStepFuncPdf()
{
// Destructor
}
void LauParametricStepFuncPdf::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];
// Get the parameters
- const vector<LauAbsRValue*>& pars = this->getParameters();
+ const std::vector<LauAbsRValue*>& pars = this->getParameters();
// Calculate value
Double_t value(0.0);
const UInt_t numBins = this->nBins();
if ( this->normBin() == Last ) {
// the last bin is our normalisation bin
for ( UInt_t i(1); i<=numBins; ++i ) {
if ( abscissa < limits_[i] ) {
// in bin i-1 (starting with bin 0)
if ( i < numBins ) {
// not in last bin
value = pars[i-1]->unblindValue();
break;
} else {
// in last bin
Double_t sum(0.0);
Double_t binSize(0.0);
for ( UInt_t j(1); j<numBins; ++j ) {
binSize = limits_[j] - limits_[j-1];
sum += ( pars[j-1]->unblindValue() * binSize );
}
binSize = limits_[numBins] - limits_[numBins-1];
value = ( 1.0 - sum ) / binSize;
break;
}
}
}
} else {
// the first bin is our normalisation bin
for ( UInt_t i(1); i<=numBins; ++i ) {
if ( abscissa < limits_[i] ) {
// in bin i-1 (starting with bin 0)
if ( i > 1 ) {
// not in first bin
value = pars[i-2]->unblindValue();
break;
} else {
// in first bin
Double_t sum(0.0);
Double_t binSize(0.0);
for ( UInt_t j(2); j<=numBins; ++j ) {
binSize = limits_[j] - limits_[j-1];
sum += ( pars[j-2]->unblindValue() * binSize );
}
binSize = limits_[1] - limits_[0];
value = ( 1.0 - sum ) / binSize;
break;
}
}
}
}
this->setUnNormPDFVal(value);
// if the parameters are floating then we
// need to recalculate the normalisation
if (!this->cachePDF() && !this->withinNormCalc() && !this->withinGeneration()) {
this->calcNorm();
}
}
void LauParametricStepFuncPdf::calcNorm()
{
this->setNorm(1.0);
}
void LauParametricStepFuncPdf::calcPDFHeight( const LauKinematics* /*kinematics*/ )
{
if (this->heightUpToDate()) {
return;
}
// Get the parameters
- const vector<LauAbsRValue*>& pars = this->getParameters();
+ const std::vector<LauAbsRValue*>& pars = this->getParameters();
// Find the PDF height
Double_t height(0.0);
Double_t value(0.0);
const UInt_t numBins = this->nBins();
if ( this->normBin() == Last ) {
// the last bin is our normalisation bin
// Check through all the parameterised bins
for ( UInt_t i(0); i<numBins-1; ++i ) {
value = pars[i]->unblindValue();
if ( height < value ) {
height = value;
}
}
// Check the last bin
Double_t sum(0.0);
Double_t binSize(0.0);
for ( UInt_t j(1); j<numBins; ++j ) {
binSize = limits_[j] - limits_[j-1];
sum += ( pars[j-1]->unblindValue() * binSize );
}
binSize = limits_[numBins] - limits_[numBins-1];
value = ( 1.0 - sum ) / binSize;
if ( height < value ) {
height = value;
}
} else {
// the first bin is our normalisation bin
// Check through all the parameterised bins
for ( UInt_t i(1); i<numBins; ++i ) {
value = pars[i-1]->unblindValue();
if ( height < value ) {
height = value;
}
}
// Check the first bin
Double_t sum(0.0);
Double_t binSize(0.0);
for ( UInt_t j(2); j<=numBins; ++j ) {
binSize = limits_[j] - limits_[j-1];
sum += ( pars[j-2]->unblindValue() * binSize );
}
binSize = limits_[1] - limits_[0];
value = ( 1.0 - sum ) / binSize;
if ( height < value ) {
height = value;
}
}
this->setMaxHeight(height);
}
diff --git a/src/LauSPlot.cc b/src/LauSPlot.cc
index 993a543..a7a8f09 100644
--- a/src/LauSPlot.cc
+++ b/src/LauSPlot.cc
@@ -1,1303 +1,1302 @@
/*
Copyright 2006 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauSPlot.cc
\brief File containing implementation of LauSPlot class.
Class for defining the SPlot technique based on TSplot from ROOT by the following authors:
Muriel Pivk, Anna Kreshuk (10/2005).
(Original copyright notice below)
*/
/**********************************************************************
* *
* Copyright (c) 2005 ROOT Foundation, CERN/PH-SFT *
* *
**********************************************************************/
#include <cfloat>
#include <iostream>
#include <vector>
using std::cout;
using std::cerr;
using std::endl;
-using std::vector;
#include "TEventList.h"
#include "TFile.h"
#include "TLeaf.h"
#include "TMath.h"
#include "TObjArray.h"
#include "TSystem.h"
#include "TTree.h"
#include "TVirtualFitter.h"
#include "LauSPlot.hh"
extern void Yields(Int_t &, Double_t *, Double_t &f, Double_t *x, Int_t iflag);
ClassImp(LauSPlot)
LauSPlot::LauSPlot(const TString& fileName, const TString& treeName,
Int_t firstExpt,
Int_t nExpt,
const NameSet& variableNames,
const NumbMap& freeSpecies,
const NumbMap& fixdSpecies,
const TwoDMap& twodimPDFs,
Bool_t sigSplit,
Bool_t scfDPSmeared) :
fileName_(fileName),
inputTreeName_(treeName),
cnTreeName_(""),
sweightTreeName_(""),
file_(0),
inputTree_(0),
cnTree_(0),
sweightTree_(0),
eventList_(0),
variableNames_(variableNames),
freeSpecies_(freeSpecies),
fixdSpecies_(fixdSpecies),
origFreeSpecies_(freeSpecies),
origFixdSpecies_(fixdSpecies),
twodimPDFs_(twodimPDFs),
signalSplit_(sigSplit),
scfDPSmear_(scfDPSmeared),
readInput_(kFALSE),
definedCNBranches_(kFALSE),
definedSWeightBranches_(kFALSE),
firstExpt_(firstExpt),
nExpt_(nExpt),
iExpt_(0),
nEvents_(0),
nDiscVars_(variableNames.size()),
nFreeSpecies_(freeSpecies.size()),
nFixdSpecies_(fixdSpecies.size()),
nSpecies_(freeSpecies.size()+fixdSpecies.size())
{
this->openInputFileAndTree();
this->readInputInfo();
this->createSWeightTree();
if (nFixdSpecies_>0) {
this->createCNTree();
}
}
LauSPlot::~LauSPlot()
{
// seems that closing the file deletes the tree
// so only delete if the file is still open
if (file_ && file_->IsOpen()) {
delete inputTree_; inputTree_ = 0;
delete sweightTree_; sweightTree_ = 0;
delete cnTree_; cnTree_ = 0;
}
delete file_; file_ = 0;
}
void LauSPlot::openInputFileAndTree()
{
// first check whether we've already opened up the file or not
if (!file_) {
// if not, first check the filename and if all ok create the file
if (fileName_ == "") {
cerr<<"ERROR in LauSPlot::createFileAndTree : Bad filename supplied, not creating file or tree."<<endl;
return;
}
file_ = TFile::Open(fileName_, "update");
if (!file_ || file_->IsZombie() || !file_->IsWritable()) {
cerr<<"ERROR in LauSPlot::createFileAndTree : Problem opening file \""<<fileName_<<"\" for updating, can't do anything."<<endl;
return;
}
}
// next open the input tree for reading
if (!inputTree_) {
file_->cd();
inputTree_ = dynamic_cast<TTree*>(file_->Get(inputTreeName_));
inputTree_->SetDirectory(file_);
}
}
void LauSPlot::readInputInfo()
{
// Read the tree and then setup the maps so we know which leaves to
// read from the tree to get the various PDF values
Bool_t inputOK = this->readInputLeaves();
if (!inputOK) {
this->readInput(inputOK);
return;
}
inputOK = this->checkLeaves();
this->readInput(inputOK);
return;
}
Bool_t LauSPlot::readInputLeaves()
{
// Read all the leaves in the tree and store them in the leaves map
if (!inputTree_) {
cerr<<"ERROR in LauSPlot::readInputInfo : Invalid pointer to data tree."<<endl;
return kFALSE;
}
Int_t nBranches = inputTree_ ? static_cast<Int_t>(inputTree_->GetNbranches()) : 0;
TObjArray* pLeaves = inputTree_->GetListOfLeaves();
if (!pLeaves) {
cerr<<"ERROR in LauSPlot::readInputInfo : Problem retrieving leaves from the tree."<<endl;
return kFALSE;
}
TObjArray& leaves = *pLeaves;
if (nBranches > leaves.GetSize()) {
cerr<<"ERROR in LauSPlot::readInputInfo : List of leaves is smaller than number of branches - this is strange!"<<endl;
return kFALSE;
}
for (Int_t iLeaf = 0; iLeaf < nBranches; ++iLeaf) {
TLeaf * leaf = dynamic_cast<TLeaf*>(leaves[iLeaf]);
// we can't deal with arrays
Int_t size = leaf->GetNdata();
if (size != 1) {
cerr<<"ERROR in LauSPlot::readInputInfo : Tree has array branches, can't deal with those."<<endl;
return kFALSE;
}
// find the name of the leaf
TString name = leaf->GetName();
// initialise an entry in the maps to hold the value
leaves_[name] = leaf;
}
return kTRUE;
}
Bool_t LauSPlot::checkLeaves() const
{
// Analyse the names of the leaves to check that we have a leaf for
// all bits of information we require, i.e. a likelihood value for
// each combination of variable and species
// If we have 2D PDFs then we have to look for some extra leaves
for ( TwoDMap::const_iterator twodim_iter = twodimPDFs_.begin(); twodim_iter != twodimPDFs_.end(); ++twodim_iter ) {
const TString& specName = twodim_iter->first;
const TString& firstVarName = twodim_iter->second.first;
const TString& secondVarName = twodim_iter->second.second;
TString expectedName(specName);
expectedName += firstVarName;
expectedName += secondVarName;
expectedName += "Like";
Bool_t found(kFALSE);
for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
const TString& leafName = leaf_iter->first;
const TLeaf* leaf = leaf_iter->second;
if ((leafName == expectedName) && (leaf != 0)) {
found = kTRUE;
break;
}
}
if (!found) {
cerr<<"ERROR in LauSPlot::checkLeaves : Could not find expected leaf \""<<expectedName<<"\"."<<endl;
return kFALSE;
}
}
// Search for all other "standard" leaves, i.e. <species><var>Like
for (NumbMap::const_iterator fixd_iter = fixdSpecies_.begin(); fixd_iter != fixdSpecies_.end(); ++fixd_iter) {
const TString& specName = fixd_iter->first;
// if the signal is split we need to do a dedicated search
// for sigTM and sigSCF, so skip the signal here
if ( specName == "sig" && this->signalSplit() ) {
continue;
}
for (NameSet::const_iterator vars_iter = variableNames_.begin(); vars_iter != variableNames_.end(); ++vars_iter) {
const TString& varName = (*vars_iter);
TString expectedName(specName);
expectedName += varName;
expectedName += "Like";
Bool_t found(kFALSE);
for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
const TString& leafName = leaf_iter->first;
const TLeaf* leaf = leaf_iter->second;
if ((leafName == expectedName) && (leaf != 0)) {
found = kTRUE;
break;
}
}
if (!found) {
cerr<<"ERROR in LauSPlot::checkLeaves : Could not find expected leaf \""<<expectedName<<"\"."<<endl;
return kFALSE;
}
}
}
for (NumbMap::const_iterator free_iter = freeSpecies_.begin(); free_iter != freeSpecies_.end(); ++free_iter) {
const TString& specName = free_iter->first;
// if the signal is split we need to do a dedicated search
// for sigTM and sigSCF, so skip the signal here
if ( specName == "sig" && this->signalSplit() ) {
continue;
}
for (NameSet::const_iterator vars_iter = variableNames_.begin(); vars_iter != variableNames_.end(); ++vars_iter) {
const TString& varName = (*vars_iter);
TString expectedName(specName);
expectedName += varName;
expectedName += "Like";
Bool_t found(kFALSE);
for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
const TString& leafName = leaf_iter->first;
const TLeaf* leaf = leaf_iter->second;
if ((leafName == expectedName) && (leaf != 0)) {
found = kTRUE;
break;
}
}
if (!found) {
cerr<<"ERROR in LauSPlot::checkLeaves : Could not find expected leaf \""<<expectedName<<"\"."<<endl;
return kFALSE;
}
}
}
if ( this->signalSplit() ) {
// now need to search for the sigTM and sigSCF leaves
for (NameSet::const_iterator vars_iter = variableNames_.begin(); vars_iter != variableNames_.end(); ++vars_iter) {
const TString& varName = (*vars_iter);
TString expectedName("sigTM");
expectedName += varName;
expectedName += "Like";
Bool_t found(kFALSE);
for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
const TString& leafName = leaf_iter->first;
const TLeaf* leaf = leaf_iter->second;
if ((leafName == expectedName) && (leaf != 0)) {
found = kTRUE;
break;
}
}
if (!found) {
cerr<<"ERROR in LauSPlot::checkLeaves : Could not find expected leaf \""<<expectedName<<"\"."<<endl;
return kFALSE;
}
expectedName = "sigSCF";
expectedName += varName;
expectedName += "Like";
found = kFALSE;
for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
const TString& leafName = leaf_iter->first;
const TLeaf* leaf = leaf_iter->second;
if ((leafName == expectedName) && (leaf != 0)) {
found = kTRUE;
break;
}
}
if (!found) {
cerr<<"ERROR in LauSPlot::checkLeaves : Could not find expected leaf \""<<expectedName<<"\"."<<endl;
return kFALSE;
}
expectedName = "sigSCFFrac";
found = kFALSE;
for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
const TString& leafName = leaf_iter->first;
const TLeaf* leaf = leaf_iter->second;
if ((leafName == expectedName) && (leaf != 0)) {
found = kTRUE;
break;
}
}
if (!found) {
cerr<<"ERROR in LauSPlot::checkLeaves : Could not find expected leaf \""<<expectedName<<"\"."<<endl;
return kFALSE;
}
}
}
return kTRUE;
}
void LauSPlot::createCNTree()
{
// check whether we've already created the tree
if (!cnTree_) {
// if not change to the file's directory and create the tree
cnTreeName_ = inputTreeName_;
cnTreeName_ += "_cNCoeffs";
file_->cd();
cnTree_ = new TTree(cnTreeName_, cnTreeName_);
cnTree_->SetDirectory(file_);
this->definedCNBranches(kFALSE);
}
}
void LauSPlot::createSWeightTree()
{
// check whether we've already created the tree
if (!sweightTree_) {
// if not change to the file's directory and create the tree
sweightTreeName_ = inputTreeName_;
sweightTreeName_ += "_sWeights";
file_->cd();
sweightTree_ = new TTree(sweightTreeName_, sweightTreeName_);
sweightTree_->SetDirectory(file_);
this->definedSWeightBranches(kFALSE);
}
}
void LauSPlot::defineCNBranches()
{
if (this->definedCNBranches()) {
cerr<<"ERROR in LauSPlot::defineCNBranches : Already defined branches, not doing it again."<<endl;
return;
}
if (cN_.empty()) {
cerr<<"ERROR in LauSPlot::defineCNBranches : No entries in the cN container, can't define branches."<<endl;
return;
}
if (!cnTree_) {
cerr<<"ERROR in LauSPlot::defineCNBranches : Invalid pointer to the tree, can't define branches."<<endl;
return;
}
// In the cN tree there is one entry per experiment, so need to know the experiment number.
cnTree_->Branch("iExpt", &iExpt_, "iExpt/I");
for (std::map<TString,NumbMap>::iterator excl_iter = cN_.begin(); excl_iter != cN_.end(); ++excl_iter) {
const TString& exclName = excl_iter->first;
NumbMap& species = excl_iter->second;
for (NumbMap::iterator spec_iter = species.begin(); spec_iter != species.end(); ++spec_iter) {
const TString& specName = spec_iter->first;
Double_t * pointer = &(spec_iter->second);
TString name(specName); name += "_cN";
if (exclName == "none") {
name += "_all";
} else {
name += "_no";
name += exclName;
}
TString thirdPart(name); thirdPart += "/D";
cnTree_->Branch(name, pointer, thirdPart);
}
}
this->definedCNBranches(kTRUE);
}
void LauSPlot::defineSWeightBranches()
{
if (this->definedSWeightBranches()) {
cerr<<"ERROR in LauSPlot::defineSWeightBranches : Already defined branches, not doing it again."<<endl;
return;
}
if (sWeights_.empty()) {
cerr<<"ERROR in LauSPlot::defineSWeightBranches : No entries in the sWeights container, can't define branches."<<endl;
return;
}
if (!sweightTree_) {
cerr<<"ERROR in LauSPlot::defineSWeightBranches : Invalid pointer to the tree, can't define branches."<<endl;
return;
}
for (std::map<TString,NumbMap>::iterator excl_iter = sWeightsCurrent_.begin(); excl_iter != sWeightsCurrent_.end(); ++excl_iter) {
const TString& exclName = excl_iter->first;
NumbMap& species = excl_iter->second;
for (NumbMap::iterator spec_iter = species.begin(); spec_iter != species.end(); ++spec_iter) {
const TString& specName = spec_iter->first;
Double_t * pointer = &(spec_iter->second);
TString name(specName); name += "_sWeight";
if (exclName == "none") {
name += "_all";
} else {
name += "_no";
name += exclName;
}
TString thirdPart(name); thirdPart += "/D";
sweightTree_->Branch(name, pointer, thirdPart);
}
}
this->definedSWeightBranches(kTRUE);
}
void LauSPlot::setExperiment(Int_t iExpt)
{
if (!inputTree_) {
cerr<<"ERROR in LauSPlot::setExperiment : Invalid pointer to data tree."<<endl;
return;
}
// create the event list if we haven't already done so
if (!eventList_) {
eventList_ = new TEventList("splotelist","splotelist");
eventList_->SetDirectory(file_);
}
// fill the event list with this experiment's events
TString listName(eventList_->GetName());
listName.Prepend(">>");
TString selection("iExpt==");
selection += iExpt;
cout<<"LauSPlot::setExperiment : Setting tree to experiment number "<<iExpt<<"."<<endl;
inputTree_->Draw(listName,selection);
// find how many events there are
nEvents_ = eventList_->GetN();
cout<<" Found "<<nEvents_<<" events."<<endl;
// make sure we have enough space in the per-event value vectors
pdfTot_.clear(); pdfTot_.resize(nEvents_);
discPdf_.clear(); discPdf_.resize(nEvents_);
scfFrac_.clear(); scfFrac_.resize(nEvents_);
sWeights_.clear(); sWeights_.resize(nEvents_);
// read the info for this experiment
this->readExpt();
}
void LauSPlot::readExpt()
{
for (Int_t iEvt(0); iEvt < nEvents_; ++iEvt) {
// Find which entry from the full tree contains the requested event
Long64_t iEntry = eventList_ ? eventList_->GetEntry(iEvt) : iEvt;
if (iEntry<0) { // this shouldn't happen, but just in case...
cerr<<"ERROR in LauSPlot::readExpt : Problem retrieving event."<<endl;
gSystem->Exit(EXIT_FAILURE);
}
// Then retrieve that entry from the tree
inputTree_->GetEntry(iEntry);
// If needed retrieve the SCF fraction values
if ( this->signalSplit() ) {
for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
const TString& leafName = leaf_iter->first;
const TLeaf* leaf = leaf_iter->second;
if ((leafName == "sigSCFFrac") && (leaf != 0)) {
scfFrac_[iEvt] = leaf->GetValue();
break;
}
}
}
// Copy the leaf values into discPdf_
for ( TwoDMap::const_iterator twodim_iter = twodimPDFs_.begin(); twodim_iter != twodimPDFs_.end(); ++twodim_iter ) {
const TString& specName = twodim_iter->first;
const TString& firstVarName = twodim_iter->second.first;
const TString& secondVarName = twodim_iter->second.second;
TString varName = firstVarName + secondVarName;
TString expectedName(specName);
expectedName += varName;
expectedName += "Like";
for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
const TString& leafName = leaf_iter->first;
const TLeaf* leaf = leaf_iter->second;
if ((leafName == expectedName) && (leaf != 0)) {
discPdf_[iEvt][specName][varName] = leaf->GetValue();
break;
}
}
}
Bool_t needSignalSearch(kFALSE);
for (NumbMap::const_iterator fixd_iter = fixdSpecies_.begin(); fixd_iter != fixdSpecies_.end(); ++fixd_iter) {
const TString& specName = fixd_iter->first;
// if the signal is split we need to do a dedicated search
// for sigTM and sigSCF, so skip the signal here
if ( specName == "sig" && this->signalSplit() ) {
needSignalSearch = kTRUE;
continue;
}
for (NameSet::const_iterator vars_iter = variableNames_.begin(); vars_iter != variableNames_.end(); ++vars_iter) {
const TString& varName = (*vars_iter);
TString expectedName(specName);
expectedName += varName;
expectedName += "Like";
for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
const TString& leafName = leaf_iter->first;
const TLeaf* leaf = leaf_iter->second;
if ((leafName == expectedName) && (leaf != 0)) {
discPdf_[iEvt][specName][varName] = leaf->GetValue();
break;
}
}
}
}
for (NumbMap::const_iterator free_iter = freeSpecies_.begin(); free_iter != freeSpecies_.end(); ++free_iter) {
const TString& specName = free_iter->first;
// if the signal is split we need to do a dedicated search
// for sigTM and sigSCF, so skip the signal here
if ( specName == "sig" && this->signalSplit() ) {
needSignalSearch = kTRUE;
continue;
}
for (NameSet::const_iterator vars_iter = variableNames_.begin(); vars_iter != variableNames_.end(); ++vars_iter) {
const TString& varName = (*vars_iter);
TString expectedName(specName);
expectedName += varName;
expectedName += "Like";
for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
const TString& leafName = leaf_iter->first;
const TLeaf* leaf = leaf_iter->second;
if ((leafName == expectedName) && (leaf != 0)) {
discPdf_[iEvt][specName][varName] = leaf->GetValue();
break;
}
}
}
}
if ( needSignalSearch ) {
for (NameSet::const_iterator vars_iter = variableNames_.begin(); vars_iter != variableNames_.end(); ++vars_iter) {
const TString& varName = (*vars_iter);
TString specName("sigTM");
TString expectedName(specName);
expectedName += varName;
expectedName += "Like";
for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
const TString& leafName = leaf_iter->first;
const TLeaf* leaf = leaf_iter->second;
if ((leafName == expectedName) && (leaf != 0)) {
discPdf_[iEvt][specName][varName] = leaf->GetValue();
break;
}
}
specName = "sigSCF";
expectedName = specName;
expectedName += varName;
expectedName += "Like";
for (LeafMap::const_iterator leaf_iter = leaves_.begin(); leaf_iter != leaves_.end(); ++leaf_iter) {
const TString& leafName = leaf_iter->first;
const TLeaf* leaf = leaf_iter->second;
if ((leafName == expectedName) && (leaf != 0)) {
discPdf_[iEvt][specName][varName] = leaf->GetValue();
break;
}
}
}
}
}
}
void LauSPlot::runCalculations(const TString& option)
{
// Calculates the sWeights and cN coeffs
// The option controls the print level:
// "Q" - no print out (default)
// "V" - prints the estimated # of events in species
// "VV" - as "V" + the MINUIT printing + sums of weights for control
if (!this->readInput()) {
cerr<<"ERROR in LauSPlot::runCalculations : The input ntuple has not been successfully read, can't calculate anything."<<endl;
return;
}
if (freeSpecies_.empty()) {
cerr<<"ERROR in LauSPlot::runCalculations : Numbers of events in species have not been set."<<endl;
return;
}
TString opt(option);
opt.ToUpper();
opt.ReplaceAll("VV", "W");
// Make sure that global fitter is MINUIT
this->checkFitter();
// Loop over experiments
for (iExpt_ = firstExpt_; iExpt_ < (firstExpt_+nExpt_); ++iExpt_) {
this->setExperiment(iExpt_);
if (nEvents_ < 1) {
cerr<<"ERROR in LauSPlot::runCalculations : Zero events in experiment "<<iExpt_<<", skipping..."<<endl;
continue;
}
// Now loop over the PDFs to exclude, including the case where none are to be excluded.
NameSet excludePdf;
if (variableNames_.size()<2) {
excludePdf.insert("none");
} else {
excludePdf = variableNames_;
excludePdf.insert("none");
}
for (NameSet::const_iterator excl_iter = excludePdf.begin(); excl_iter != excludePdf.end(); ++excl_iter) {
const TString& exclName = (*excl_iter);
cout<<"LauSPlot::runCalculations : Calculating sWeights, excluding PDF: "<<exclName<<"."<<endl;
// Calculate the per-event total PDF values for each species.
this->calcTotPDFValues(exclName);
// Reset the fitter
this->initialiseFitter(opt);
// Set the parameters to their initial values
this->setFitParameters();
// Run the fit
this->runFit();
// Get the fitted parameter values back
this->retrieveFittedParameters(opt);
// Get the covariance matrix
Bool_t covMatOK = this->calcCovMatrix();
Double_t * covmat(0);
if (!covMatOK) {
TVirtualFitter * fitter = TVirtualFitter::GetFitter();
covmat = fitter->GetCovarianceMatrix();
}
if (opt.Contains("W")) {
this->printCovMatrixElements(covmat);
}
// calculate the cN and sWeights from the covariance matrix
if (nFixdSpecies_ > 0) {
this->calcCNCoeffs(exclName, covmat);
}
this->calcSWeights(exclName, covmat);
// print verbose info if required
if (opt.Contains("W")) {
this->printSumOfWeights(exclName);
}
}
// Finally fill all the branches for this experiment
if (nFixdSpecies_ > 0) {
this->fillCNBranches();
}
this->fillSWeightBranches();
}
}
void LauSPlot::checkFitter() const
{
TString minuitName("TFitter");
TVirtualFitter * fitter = TVirtualFitter::GetFitter();
if (fitter) {
TString fitterName(fitter->IsA()->GetName());
if (fitterName != minuitName) {
delete fitter; fitter = 0;
}
}
fitter = TVirtualFitter::Fitter(0);
}
void LauSPlot::initialiseFitter(const TString& opt)
{
TVirtualFitter * fitter = TVirtualFitter::GetFitter();
fitter->Clear();
fitter->SetFCN(Yields);
fitter->SetObjectFit(this);
// Set the print level
Double_t arglist[10];
if (opt.Contains("Q")) {
arglist[0]=-1;
}
if (opt.Contains("V")) {
arglist[0]=0;
}
if (opt.Contains("W")) {
arglist[0]=1;
}
fitter->ExecuteCommand("SET PRINT", arglist, 1);
// Need to set the "SET ERR" command to +0.5 for +/-1 sigma errors
// for maximum likelihood fit.
arglist[0] = 0.5;
fitter->ExecuteCommand("SET ERR", arglist, 1);
}
void LauSPlot::setFitParameters() const
{
TVirtualFitter * fitter = TVirtualFitter::GetFitter();
// must add the parameters in the same order as they are stored in pdfTot_
Int_t ispecies(0);
const NumbMap& species = pdfTot_.front();
for (NumbMap::const_iterator spec_iter = species.begin(); spec_iter != species.end(); ++spec_iter) {
const TString& name(spec_iter->first);
// starting parameters should be the original values,
// not those that came out of the last fit
NumbMap::const_iterator free_iter = origFreeSpecies_.find(name);
NumbMap::const_iterator fixd_iter = origFixdSpecies_.find(name);
Bool_t fixed = fixd_iter != origFixdSpecies_.end();
Double_t value(0.0);
if (fixed) {
value = fixd_iter->second;
} else {
value = free_iter->second;
}
fitter->SetParameter(ispecies, name, value, 1, 0, 0);
if (fixed) {
fitter->FixParameter(ispecies);
}
++ispecies;
}
}
void LauSPlot::runFit()
{
TVirtualFitter * fitter = TVirtualFitter::GetFitter();
Double_t arglist[10];
arglist[0] = 1000*nFreeSpecies_; // maximum iterations
arglist[1] = 0.05; // tolerance : min EDM = 0.001*tolerance
// Execute MIGRAD
Int_t fitStatus = fitter->ExecuteCommand("MIGRAD", arglist, 2);
if (fitStatus != 0) {
cerr<<"ERROR in LauSPlot::runFit : Error during MIGRAD minimisation."<<endl;
} else {
// Execute HESSE
fitStatus = fitter->ExecuteCommand("HESSE", arglist, 1);
if (fitStatus != 0) {
cerr<<"ERROR in LauSPlot::runFit : Error during HESSE error calculation."<<endl;
}
}
}
void LauSPlot::printCovMatrixElements(const Double_t * covmat) const
{
Int_t a(0);
cout<<endl;
for (NumbMap::const_iterator iter_a = freeSpecies_.begin(); iter_a != freeSpecies_.end(); ++iter_a) {
Int_t b(0);
for (NumbMap::const_iterator iter_b = freeSpecies_.begin(); iter_b != freeSpecies_.end(); ++iter_b) {
if (covmat) {
cout<<"CovMat Elem: ["<<iter_a->first<<","<<iter_b->first<<"] = "<<covmat[a*nSpecies_+b]<<endl;
} else {
cout<<"CovMat Elem: ["<<iter_a->first<<","<<iter_b->first<<"] = "<<covMat_(a,b)<<endl;
}
++b;
}
++a;
}
cout<<endl;
if (!covmat) {
covMat_.Print();
}
}
void LauSPlot::retrieveFittedParameters(const TString& opt)
{
TVirtualFitter * fitter = TVirtualFitter::GetFitter();
// remember fit parameters are in same order as in pdfTot_
Int_t ispecies(0);
const NumbMap& species = pdfTot_.front();
for (NumbMap::const_iterator spec_iter = species.begin(); spec_iter != species.end(); ++spec_iter) {
const TString& name(spec_iter->first);
NumbMap::iterator free_iter = freeSpecies_.find(name);
if (free_iter != freeSpecies_.end()) {
free_iter->second = fitter->GetParameter(ispecies);
if (!opt.Contains("Q")) {
cout<<"Estimated # of events in species "<<name<<" = "<<free_iter->second<<endl;
}
}
++ispecies;
}
if (!opt.Contains("Q")) {cout<<endl;}
}
void LauSPlot::printSumOfWeights(const TString& exclName) const
{
for (NumbMap::const_iterator spec_iter = freeSpecies_.begin(); spec_iter != freeSpecies_.end(); ++spec_iter) {
const TString& specName = spec_iter->first;
Double_t sumweight(0.0);
for (Int_t iEvt(0); iEvt < nEvents_; ++iEvt) {
const NumbMap& specWeightMap = sWeights_[iEvt].find(exclName)->second;
Double_t weight = specWeightMap.find(specName)->second;
sumweight += weight;
}
cout<<"Sum of sWeights for species \""<<specName<<"\" = "<<sumweight<<endl;
}
cout<<endl;
}
Bool_t LauSPlot::calcCovMatrix()
{
// Calculate our inverse covariance matrix from the various PDFs
TMatrixD invMatrix(nFreeSpecies_,nFreeSpecies_);
// First calculate the denominator, which is common to all elements
- vector<Double_t> denominator(nEvents_);
+ std::vector<Double_t> denominator(nEvents_);
for (Int_t iEvt(0); iEvt<nEvents_; ++iEvt) {
denominator[iEvt] = 0.0;
for (NumbMap::const_iterator spec_iter = freeSpecies_.begin(); spec_iter != freeSpecies_.end(); ++spec_iter) {
const TString& specName = spec_iter->first;
Double_t specNumEvents = spec_iter->second;
denominator[iEvt] += specNumEvents * pdfTot_[iEvt][specName];
}
for (NumbMap::const_iterator spec_iter = fixdSpecies_.begin(); spec_iter != fixdSpecies_.end(); ++spec_iter) {
const TString& specName = spec_iter->first;
Double_t specNumEvents = spec_iter->second;
denominator[iEvt] += specNumEvents * pdfTot_[iEvt][specName];
}
// Square to get the final denominator
denominator[iEvt] *= denominator[iEvt];
}
// Then calculate each element
Int_t i(0);
for (NumbMap::const_iterator spec_iter_i = freeSpecies_.begin(); spec_iter_i != freeSpecies_.end(); ++spec_iter_i) {
const TString& specName_i = spec_iter_i->first;
Int_t j(0);
for (NumbMap::const_iterator spec_iter_j = freeSpecies_.begin(); spec_iter_j != freeSpecies_.end(); ++spec_iter_j) {
const TString& specName_j = spec_iter_j->first;
invMatrix(i,j) = 0.0;
for (Int_t iEvt(0); iEvt<nEvents_; ++iEvt) {
Double_t numerator = pdfTot_[iEvt][specName_i] * pdfTot_[iEvt][specName_j];
invMatrix(i,j) += numerator/denominator[iEvt];
}
++j;
}
++i;
}
// Check for a singular matrix
if (invMatrix.Determinant() < 1e-15) {
cerr<<"ERROR in LauSPlot::calcCovMatrix : Calculated inverse covariance matrix is singular, can't invert it."<<endl;
return kFALSE;
}
// Invert and store in the covariance matrix
covMat_.ResizeTo(nFreeSpecies_,nFreeSpecies_);
covMat_ = TMatrixD(TMatrixD::kInverted, invMatrix);
return kTRUE;
}
void LauSPlot::calcTotPDFValues(const TString& exclName)
{
for (Int_t iEvt(0); iEvt<nEvents_; iEvt++) {
Bool_t needSignalSearch(kFALSE);
for (NumbMap::const_iterator spec_iter = fixdSpecies_.begin(); spec_iter != fixdSpecies_.end(); ++spec_iter) {
const TString& specName = spec_iter->first;
// if the signal is split we need to do a dedicated search
// for sigTM and sigSCF, so skip the signal here
if ( specName == "sig" && this->signalSplit() ) {
needSignalSearch = kTRUE;
continue;
}
pdfTot_[iEvt][specName] = 1.0;
// loop through the 2D histo list
NameSet skipList;
for ( TwoDMap::const_iterator twodim_iter = twodimPDFs_.begin(); twodim_iter != twodimPDFs_.end(); ++twodim_iter ) {
// if the entry doesn't refer to this
// species then skip on
if ( specName != twodim_iter->first ) {
continue;
}
// retrieve the two variable names
const TString& firstVarName = twodim_iter->second.first;
const TString& secondVarName = twodim_iter->second.second;
if ( firstVarName != exclName && secondVarName != exclName ) {
// if neither name is the one being excluded then...
// add them both to the skip list
skipList.insert( firstVarName );
skipList.insert( secondVarName );
// and multiply the total by the combined PDF value
TString varName = firstVarName + secondVarName;
pdfTot_[iEvt][specName] *= discPdf_[iEvt][specName][varName];
}
}
// loop through all the variables
for (NameSet::const_iterator var_iter = variableNames_.begin(); var_iter != variableNames_.end(); ++var_iter) {
const TString& varName = (*var_iter);
// if the variable isn't the one being excluded
if (exclName != varName) {
// and it's not involved in a 2D PDF
if ( skipList.find(varName) == skipList.end() ) {
// multiply the total by its PDF value
pdfTot_[iEvt][specName] *= discPdf_[iEvt][specName][varName];
}
}
}
}
for (NumbMap::const_iterator spec_iter = freeSpecies_.begin(); spec_iter != freeSpecies_.end(); ++spec_iter) {
const TString& specName = spec_iter->first;
// if the signal is split we need to do a dedicated search
// for sigTM and sigSCF, so skip the signal here
if ( specName == "sig" && this->signalSplit() ) {
needSignalSearch = kTRUE;
continue;
}
pdfTot_[iEvt][specName] = 1.0;
// loop through the 2D histo list
NameSet skipList;
for ( TwoDMap::const_iterator twodim_iter = twodimPDFs_.begin(); twodim_iter != twodimPDFs_.end(); ++twodim_iter ) {
// if the entry doesn't refer to this
// species then skip on
if ( specName != twodim_iter->first ) {
continue;
}
// retrieve the two variable names
const TString& firstVarName = twodim_iter->second.first;
const TString& secondVarName = twodim_iter->second.second;
if ( firstVarName != exclName && secondVarName != exclName ) {
// if neither name is the one being excluded then...
// add them both to the skip list
skipList.insert( firstVarName );
skipList.insert( secondVarName );
// and multiply the total by the combined PDF value
TString varName = firstVarName + secondVarName;
pdfTot_[iEvt][specName] *= discPdf_[iEvt][specName][varName];
}
}
// loop through all the variables
for (NameSet::const_iterator var_iter = variableNames_.begin(); var_iter != variableNames_.end(); ++var_iter) {
const TString& varName = (*var_iter);
// if the variable isn't the one being excluded
if (exclName != varName) {
// and it's not involved in a 2D PDF
if ( skipList.find(varName) == skipList.end() ) {
// multiply the total by its PDF value
pdfTot_[iEvt][specName] *= discPdf_[iEvt][specName][varName];
}
}
}
}
if ( needSignalSearch ) {
// loop through the 2D histo list
TString specName("sigTM");
Double_t tmPDFVal(1.0);
NameSet skipList;
for ( TwoDMap::const_iterator twodim_iter = twodimPDFs_.begin(); twodim_iter != twodimPDFs_.end(); ++twodim_iter ) {
// if the entry doesn't refer to this
// species then skip on
if ( specName != twodim_iter->first ) {
continue;
}
// retrieve the two variable names
const TString& firstVarName = twodim_iter->second.first;
const TString& secondVarName = twodim_iter->second.second;
if ( firstVarName != exclName && secondVarName != exclName ) {
// if neither name is the one being excluded then...
// add them both to the skip list
skipList.insert( firstVarName );
skipList.insert( secondVarName );
// and multiply the total by the combined PDF value
TString varName = firstVarName + secondVarName;
tmPDFVal *= discPdf_[iEvt][specName][varName];
}
}
// loop through all the variables
for (NameSet::const_iterator var_iter = variableNames_.begin(); var_iter != variableNames_.end(); ++var_iter) {
const TString& varName = (*var_iter);
// if the variable isn't the one being excluded
if (exclName != varName) {
// and it's not involved in a 2D PDF
if ( skipList.find(varName) == skipList.end() ) {
// multiply the total by its PDF value
tmPDFVal *= discPdf_[iEvt][specName][varName];
}
}
}
tmPDFVal *= (1.0 - scfFrac_[iEvt]);
// loop through the 2D histo list
specName = "sigSCF";
Double_t scfPDFVal(1.0);
skipList.clear();
for ( TwoDMap::const_iterator twodim_iter = twodimPDFs_.begin(); twodim_iter != twodimPDFs_.end(); ++twodim_iter ) {
// if the entry doesn't refer to this
// species then skip on
if ( specName != twodim_iter->first ) {
continue;
}
// retrieve the two variable names
const TString& firstVarName = twodim_iter->second.first;
const TString& secondVarName = twodim_iter->second.second;
if ( firstVarName != exclName && secondVarName != exclName ) {
// if neither name is the one being excluded then...
// add them both to the skip list
skipList.insert( firstVarName );
skipList.insert( secondVarName );
// and multiply the total by the combined PDF value
TString varName = firstVarName + secondVarName;
scfPDFVal *= discPdf_[iEvt][specName][varName];
}
}
// loop through all the variables
for (NameSet::const_iterator var_iter = variableNames_.begin(); var_iter != variableNames_.end(); ++var_iter) {
const TString& varName = (*var_iter);
// if the variable isn't the one being excluded
if (exclName != varName) {
// and it's not involved in a 2D PDF
if ( skipList.find(varName) == skipList.end() ) {
// multiply the total by its PDF value
scfPDFVal *= discPdf_[iEvt][specName][varName];
}
}
}
if ( exclName == "DP" || !this->scfDPSmear() ) {
scfPDFVal *= scfFrac_[iEvt];
}
pdfTot_[iEvt]["sig"] = tmPDFVal + scfPDFVal;
}
}
}
void LauSPlot::calcCNCoeffs(const TString& exclName, const Double_t *covmat)
{
// Computes the cN for the extended sPlots from the covariance matrix
if (nFixdSpecies_ <= 0) {
return;
}
Int_t species_n(0);
for (NumbMap::const_iterator iter_n = freeSpecies_.begin(); iter_n != freeSpecies_.end(); ++iter_n) {
Int_t species_j(0);
const TString& specName = iter_n->first;
Double_t value = iter_n->second;
cN_[exclName][specName] = value;
for (NumbMap::const_iterator iter_j = freeSpecies_.begin(); iter_j != freeSpecies_.end(); ++iter_j) {
if (covmat) {
cN_[exclName][specName] -= covmat[species_n*nSpecies_+species_j];
} else {
cN_[exclName][specName] -= covMat_(species_n,species_j);
}
++species_j;
}
++species_n;
}
}
void LauSPlot::calcSWeights(const TString& exclName, Double_t *covmat)
{
// Computes the sWeights from the covariance matrix
// NB for the extended sPlot the sum in the denominator is still over all species,
// while that in the numerator is only over the free species.
// Similarly the sWeights can only be calculated for the free species.
Double_t numerator(0.0), denominator(0.0);
for (Int_t iEvent = 0; iEvent < nEvents_; ++iEvent) {
denominator = 0.0;
for (NumbMap::const_iterator free_iter = freeSpecies_.begin(); free_iter != freeSpecies_.end(); ++free_iter) {
denominator += free_iter->second * pdfTot_[iEvent][free_iter->first];
}
for (NumbMap::const_iterator fixd_iter = fixdSpecies_.begin(); fixd_iter != fixdSpecies_.end(); ++fixd_iter) {
denominator += fixd_iter->second * pdfTot_[iEvent][fixd_iter->first];
}
Int_t species_n(0);
for (NumbMap::const_iterator iter_n = freeSpecies_.begin(); iter_n != freeSpecies_.end(); ++iter_n) {
numerator = 0.0;
Int_t species_j(0);
for (NumbMap::const_iterator iter_j = freeSpecies_.begin(); iter_j != freeSpecies_.end(); ++iter_j) {
if (covmat) {
numerator += covmat[species_n*nSpecies_+species_j] * pdfTot_[iEvent][iter_j->first];
} else {
numerator += covMat_(species_n,species_j) * pdfTot_[iEvent][iter_j->first];
}
++species_j;
}
sWeights_[iEvent][exclName][iter_n->first] = numerator/denominator;
++species_n;
}
}
}
void LauSPlot::fillCNBranches()
{
if (!cnTree_) {
cerr<<"ERROR in LauSPlot::fillCNBranches : Tree not created, cannot fill branches."<<endl;
return;
} else if (!this->definedCNBranches()) {
this->defineCNBranches();
}
cnTree_->Fill();
}
void LauSPlot::fillSWeightBranches()
{
if (!sweightTree_) {
cerr<<"ERROR in LauSPlot::fillSWeightBranches : Tree not created, cannot fill branches."<<endl;
return;
} else if (sWeights_.empty()) {
cerr<<"ERROR in LauSPlot::fillSWeightBranches : No sWeights calculated, can't fill branches."<<endl;
return;
} else if (!this->definedSWeightBranches()) {
this->copyEventWeights(0);
this->defineSWeightBranches();
}
for (Int_t iEvent = 0; iEvent < nEvents_; ++iEvent) {
this->copyEventWeights(iEvent);
sweightTree_->Fill();
}
}
void LauSPlot::copyEventWeights(Int_t iEvent)
{
for (std::map<TString,NumbMap>::const_iterator excl_iter = sWeights_[iEvent].begin(); excl_iter != sWeights_[iEvent].end(); ++excl_iter) {
const TString& exclName = excl_iter->first;
const NumbMap& species = excl_iter->second;
for (NumbMap::const_iterator spec_iter = species.begin(); spec_iter != species.end(); ++spec_iter) {
const TString& specName = spec_iter->first;
sWeightsCurrent_[exclName][specName] = spec_iter->second;
}
}
}
void LauSPlot::writeOutResults()
{
// write out the results
// remove the transient objects that we don't want (re-)written to the file
if (eventList_) {
delete eventList_; eventList_ = 0;
}
if (inputTree_) {
delete inputTree_; inputTree_ = 0;
}
// first add the input tree as a friend of the output tree
this->addFriendTree();
// then write everything to the file and clean up
file_->cd();
file_->Write();
file_->Close();
delete file_; file_ = 0;
}
void LauSPlot::addFriendTree()
{
if (!sweightTree_) {
cerr<<"ERROR in LauSPlot::addFriendTree : Tree not created, cannot add friend."<<endl;
return;
}
sweightTree_->AddFriend(inputTreeName_,fileName_);
}
void Yields(Int_t &, Double_t *, Double_t &f, Double_t *x, Int_t /*iflag*/)
{
// FCN-function for Minuit
f = 0.0;
TVirtualFitter *fitter = TVirtualFitter::GetFitter();
LauSPlot* theModel = dynamic_cast<LauSPlot*>(fitter->GetObjectFit());
const std::vector<LauSPlot::NumbMap>& pdfTot = theModel->totalPdf();
Double_t ntot(0.0);
for (std::vector<LauSPlot::NumbMap>::const_iterator evt_iter = pdfTot.begin(); evt_iter != pdfTot.end(); ++evt_iter) { // loop over events
const LauSPlot::NumbMap& species = (*evt_iter);
Int_t ispecies(0);
Double_t lik(0.0);
ntot = 0.0;
for (LauSPlot::NumbMap::const_iterator spec_iter = species.begin(); spec_iter != species.end(); ++spec_iter) { // loop over species
lik += x[ispecies] * spec_iter->second;
ntot += x[ispecies];
++ispecies;
}
if (lik < 0.0) {
// make f the worst possible value to force MINUIT
// out of this region of parameter space
f = -DBL_MAX;
break;
} else {
f += TMath::Log(lik);
}
}
// extended likelihood
f = (ntot-f);
}
diff --git a/src/LauSigmoidPdf.cc b/src/LauSigmoidPdf.cc
index 886f474..0da1e05 100644
--- a/src/LauSigmoidPdf.cc
+++ b/src/LauSigmoidPdf.cc
@@ -1,159 +1,158 @@
/*
Copyright 2012 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/
/*! \file LauSigmoidPdf.cc
\brief File containing implementation of LauSigmoidPdf class.
*/
#include <iostream>
#include <vector>
using std::cout;
using std::cerr;
using std::endl;
-using std::vector;
#include "TMath.h"
#include "TSystem.h"
#include "LauConstants.hh"
#include "LauSigmoidPdf.hh"
ClassImp(LauSigmoidPdf)
-LauSigmoidPdf::LauSigmoidPdf(const TString& theVarName, const vector<LauAbsRValue*>& params, Double_t minAbscissa, Double_t maxAbscissa) :
+LauSigmoidPdf::LauSigmoidPdf(const TString& theVarName, const std::vector<LauAbsRValue*>& params, Double_t minAbscissa, Double_t maxAbscissa) :
LauAbsPdf(theVarName, params, minAbscissa, maxAbscissa),
a_(0),
b_(0)
{
// Constructor for the general form of a Sigmoid PDF.
//
// The parameters in params "a" and "b" define the steepness of the
// slope and shift of the distribution (negative parameter "a" flips
// the distribution around the y-axis).
// The last two arguments specify the range in which the PDF is
// defined, and the PDF will be normalised w.r.t. these limits.
a_ = this->findParameter("a");
b_ = this->findParameter("b");
if ((this->nParameters() != 2) || (a_ == 0) || (b_ == 0)) {
cerr<<"ERROR in LauSigmoidPdf constructor: LauSigmoidPdf requires 2 parameters: \"a\" and \"b\"."<<endl;
gSystem->Exit(EXIT_FAILURE);
}
// Cache the normalisation factor
this->calcNorm();
}
LauSigmoidPdf::~LauSigmoidPdf()
{
// Destructor
}
LauSigmoidPdf::LauSigmoidPdf(const LauSigmoidPdf& other) : LauAbsPdf(other.varName(), other.getParameters(), other.getMinAbscissa(), other.getMaxAbscissa())
{
// Copy constructor
this->setRandomFun(other.getRandomFun());
this->calcNorm();
}
void LauSigmoidPdf::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];
// Get the up to date parameter values
Double_t a = a_->unblindValue();
Double_t b = b_->unblindValue();
// Calculate the value of the exponent for the given value of the abscissa.
Double_t exponent = -a*abscissa + b;
Double_t value = 1.0/(1.0+ TMath::Exp(exponent));
this->setUnNormPDFVal(value);
// if the parameters are floating then we
// need to recalculate the normalisation
if (!this->cachePDF() && !this->withinNormCalc() && !this->withinGeneration()) {
this->calcNorm();
}
}
void LauSigmoidPdf::calcNorm()
{
// Get the up to date parameter values
Double_t a = a_->unblindValue();
Double_t b = b_->unblindValue();
// Calculate the normalisation of the sigmoid and cache it.
Double_t norm(0.0);
if (TMath::Abs(a) > 1e-10) {
Double_t expb = TMath::Exp(b);
norm = TMath::Log( TMath::Exp(a * this->getMaxAbscissa()) + expb ) -
TMath::Log( TMath::Exp(a * this->getMinAbscissa()) + expb );
norm /= a;
}
this->setNorm(norm);
}
void LauSigmoidPdf::calcPDFHeight(const LauKinematics* /*kinematics*/)
{
if (this->heightUpToDate()) {
return;
}
// Get the up to date parameter values
Double_t a = a_->unblindValue();
LauAbscissas maxPoint(1);
maxPoint[0] = 0;
// Calculate the PDF height for the Sigmoid function.
if (a > 0.0) {
maxPoint[0] = this->getMaxAbscissa();
} else {
maxPoint[0] = this->getMinAbscissa();
}
this->calcLikelihoodInfo(maxPoint);
Double_t height = this->getUnNormLikelihood();
this->setMaxHeight(height);
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, May 3, 6:15 AM (1 d, 13 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4982963
Default Alt Text
(281 KB)

Event Timeline