Page MenuHomeHEPForge
Files F7878661

No OneTemporary
Actions

View Options

diff --git a/inc/LauAbsResonance.hh b/inc/LauAbsResonance.hh
index 4a30e65..d536de4 100644
--- a/inc/LauAbsResonance.hh
+++ b/inc/LauAbsResonance.hh
@@ -1,481 +1,502 @@
// Copyright University of Warwick 2004 - 2014.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
// Authors:
// Thomas Latham
// John Back
// Paul Harrison
/*! \file 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, /*!< S-wave 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 */
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 */
};
//! 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
*/
LauAbsResonance(const TString& resName, const Int_t resPairAmpInt, const LauDaughters* daughters);
//! 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 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
/*!
\return the ignore momenta flag
*/
Bool_t ignoreMomenta() const {return ignoreMomenta_;}
//! Set the ignore p_ and q_ flag
/*!
\param [in] boolean the ignore momenta status
*/
void ignoreMomenta(const Bool_t boolean) {ignoreMomenta_ = boolean;}
//! Get the ignore spin flag
/*!
\return the ignore spin flag
*/
Bool_t ignoreSpin() const {return ignoreSpin_;}
//! Set the ignore p_ and q_ flag
/*!
\param [in] boolean the ignore spin status
*/
void ignoreSpin(const Bool_t boolean) {ignoreSpin_ = boolean;}
//! Get the ignore barrier factor scaling flag
/*!
\return the ignore barrier amplitude scaling flag
*/
Bool_t ignoreBarrierScaling() const {return ignoreBarrierScaling_;}
//! Set the ignore barrier factor scaling flag
/*!
\param [in] boolean the ignore barrier factor scaling status
*/
void ignoreBarrierScaling(const Bool_t boolean) {ignoreBarrierScaling_ = boolean;}
+ //! Get the ignore covariant flag
+ /*!
+ \return the ignore covariant flag
+ */
+ Bool_t ignoreCovariant() const {return ignoreCovariant_;}
+
+
+ //! Set the ignore covariant flag
+ /*!
+ \param [in] boolean the ignore covariant status
+ */
+ void ignoreCovariant(const Bool_t boolean) {ignoreCovariant_ = 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 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 centrifugal barrier for the parent decay
LauBlattWeisskopfFactor* getParBWFactor() {return parBWFactor_;}
const LauBlattWeisskopfFactor* getParBWFactor() const {return parBWFactor_;}
//! Get the centrifugal barrier for the resonance decay
LauBlattWeisskopfFactor* getResBWFactor() {return resBWFactor_;}
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
/*!
\param [in] cosHel the cosine of the helicity angle
\param [in] pProd the momentum factor (typically q * p)
*/
Double_t calcSpinTerm( const Double_t cosHel, const Double_t pProd ) const;
//! 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) = 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_; }
+ //! calculate covariant factor
+ Double_t calcCovFactor( const Double_t Erm) const;
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_;
//! 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_;
//! Daughter 1 charge
Int_t chargeDaug1_;
//! Daughter 2 charge
Int_t chargeDaug2_;
//! Bachelor charge
Int_t chargeBachelor_;
//! Parent mass
Double_t massParent_;
//! Daughter 1 mass
Double_t massDaug1_;
//! Daughter 2 mass
Double_t massDaug2_;
// Bachelor mass
Double_t massBachelor_;
//! Resonance name
TString resName_;
//! Resonance name with illegal characters removed
TString sanitisedName_;
//! Resonance mass
LauParameter* resMass_;
//! Resonance width
LauParameter* resWidth_;
//! All parameters of the resonance
std::vector<LauParameter*> resParameters_;
//! Resonance spin
Int_t resSpin_;
//! Resonance charge
Int_t resCharge_;
//! DP axis identifier
Int_t resPairAmpInt_;
//! Blatt Weisskopf barrier for parent decay
LauBlattWeisskopfFactor* parBWFactor_;
//! Blatt Weisskopf barrier for resonance decay
LauBlattWeisskopfFactor* resBWFactor_;
//! Boolean to flip helicity
Bool_t flipHelicity_;
//! Boolean to ignore q_ and p_ in spinTerm, as well as ignoring momentum dependent widths
Bool_t ignoreMomenta_;
//! Boolean to set the spinTerm to unity always
Bool_t ignoreSpin_;
//! Boolean to ignore barrier factor amplitude scaling; they are still used for the width
Bool_t ignoreBarrierScaling_;
+ //! Boolean to ignore covariant factor
+ Bool_t ignoreCovariant_;
+
//! Daughter momentum in resonance rest frame
Double_t q_;
//! Bachelor momentum in resonance rest frame
Double_t p_;
//! Bachelor momentum in parent rest frame
Double_t pstar_;
+ //! Covariant factor
+ Double_t erm_;
+
ClassDef(LauAbsResonance,0) // Abstract resonance class
};
#endif
diff --git a/inc/LauKinematics.hh b/inc/LauKinematics.hh
index 75865cb..e892fa4 100644
--- a/inc/LauKinematics.hh
+++ b/inc/LauKinematics.hh
@@ -1,593 +1,611 @@
// Copyright University of Warwick 2004 - 2013.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
// Authors:
// Thomas Latham
// John Back
// Paul Harrison
/*! \file 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 "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
*/
LauKinematics(Double_t m1, Double_t m2, Double_t m3, Double_t mParent, Bool_t calcSquareDPCoords = kFALSE);
//! Destructor
virtual ~LauKinematics();
//! 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_; }
//! Enable/disable warning messages
inline void warningMessages(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(Double_t m13Sq, 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(Double_t mPrime, Double_t thetaPrime);
//! 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(Double_t m13Sq, 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(Double_t m13Sq, 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(Double_t mPrime, 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(Double_t firstMassSq, 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(Double_t m13Sq, 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(Int_t orientation = 1323, 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(Double_t m12, 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(Double_t m23, 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(Double_t m13, 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(Double_t m13Sq, 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(Double_t mPrime, 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(Double_t mijSq, Double_t mikSq, Double_t mij, Int_t i, Int_t j, 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(Double_t mijSq, Double_t cij, Double_t mij, Int_t i, Int_t j, 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(Double_t energy, 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);
//! 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/src/LauAbsResonance.cc b/src/LauAbsResonance.cc
index 1c325b4..b2a77dd 100644
--- a/src/LauAbsResonance.cc
+++ b/src/LauAbsResonance.cc
@@ -1,599 +1,633 @@
// Copyright University of Warwick 2004 - 2014.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
// Authors:
// Thomas Latham
// John Back
// Paul Harrison
/*! \file LauAbsResonance.cc
\brief File containing implementation of LauAbsResonance class.
*/
#include <iostream>
#include "TSystem.h"
#include "LauAbsResonance.hh"
#include "LauConstants.hh"
#include "LauDaughters.hh"
#include "LauKinematics.hh"
#include "LauParameter.hh"
#include "LauResonanceInfo.hh"
ClassImp(LauAbsResonance)
bool LauAbsResonance::isIncoherentModel(LauResonanceModel model) {
switch(model) {
case BW:
case RelBW:
case GS:
case Flatte:
case Sigma:
case Kappa:
case Dabba:
case LASS:
case LASS_BW:
case LASS_NR:
case EFKLLM:
case KMatrix:
case FlatNR:
case NRModel:
case BelleNR:
case PowerLawNR:
case BelleSymNR:
case BelleSymNRNoInter:
case TaylorNR:
case PolNR:
case MIPW_MagPhase:
case MIPW_RealImag:
case RhoOmegaMix_GS:
case RhoOmegaMix_RBW:
case RhoOmegaMix_GS_1:
case RhoOmegaMix_RBW_1:
break;
case GaussIncoh:
return true;
}
return false;
}
// Constructor
LauAbsResonance::LauAbsResonance(LauResonanceInfo* resInfo, const Int_t resPairAmpInt, const LauDaughters* daughters) :
resInfo_(resInfo),
daughters_(daughters),
nameParent_(""), nameDaug1_(""), nameDaug2_(""), nameBachelor_(""),
chargeParent_(0), chargeDaug1_(0), chargeDaug2_(0), chargeBachelor_(0),
massParent_(0.0), massDaug1_(0.0), massDaug2_(0.0), massBachelor_(0.0),
resName_( (resInfo!=0) ? resInfo->getName() : "" ),
sanitisedName_( (resInfo!=0) ? resInfo->getSanitisedName() : "" ),
resMass_( (resInfo!=0) ? resInfo->getMass() : 0 ),
resWidth_( (resInfo!=0) ? resInfo->getWidth() : 0 ),
resSpin_( (resInfo!=0) ? resInfo->getSpin() : 0 ),
resCharge_( (resInfo!=0) ? resInfo->getCharge() : 0 ),
resPairAmpInt_(resPairAmpInt),
parBWFactor_(0),
resBWFactor_(0),
flipHelicity_(kFALSE),
ignoreMomenta_(kFALSE),
ignoreSpin_(kFALSE),
ignoreBarrierScaling_(kFALSE),
+ ignoreCovariant_(kFALSE),
q_(0.0),
p_(0.0),
- pstar_(0.0)
+ pstar_(0.0),
+ erm_(1.0)
{
if ( resInfo == 0 ) {
std::cerr << "ERROR in LauAbsResonance constructor : null LauResonanceInfo object provided" << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
if ( daughters_ == 0 ) {
std::cerr << "ERROR in LauAbsResonance constructor : null LauDaughters object provided" << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
nameParent_ = this->getNameParent();
nameDaug1_ = this->getNameDaug1();
nameDaug2_ = this->getNameDaug2();
nameBachelor_ = this->getNameBachelor();
massParent_ = this->getMassParent();
massDaug1_ = this->getMassDaug1();
massDaug2_ = this->getMassDaug2();
massBachelor_ = this->getMassBachelor();
chargeParent_ = this->getChargeParent();
chargeDaug1_ = this->getChargeDaug1();
chargeDaug2_ = this->getChargeDaug2();
chargeBachelor_ = this->getChargeBachelor();
// check that the total charge adds up to that of the resonance
Int_t totalCharge = chargeDaug1_ + chargeDaug2_;
if ( (totalCharge != resCharge_) && (resPairAmpInt_ != 0) ) {
std::cerr << "ERROR in LauAbsResonance : Total charge of daughters = " << totalCharge << ". Resonance charge = " << resCharge_ << "." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
}
// Constructor
LauAbsResonance::LauAbsResonance(const TString& resName, const Int_t resPairAmpInt, const LauDaughters* daughters) :
resInfo_(0),
daughters_(daughters),
nameParent_(""), nameDaug1_(""), nameDaug2_(""), nameBachelor_(""),
chargeParent_(0), chargeDaug1_(0), chargeDaug2_(0), chargeBachelor_(0),
massParent_(0.0), massDaug1_(0.0), massDaug2_(0.0), massBachelor_(0.0),
resName_(resName),
sanitisedName_(resName),
resMass_(0),
resWidth_(0),
resSpin_(0),
resCharge_(0),
resPairAmpInt_(resPairAmpInt),
parBWFactor_(0),
resBWFactor_(0),
flipHelicity_(kFALSE),
ignoreMomenta_(kFALSE),
ignoreSpin_(kFALSE),
ignoreBarrierScaling_(kFALSE),
+ ignoreCovariant_(kFALSE),
q_(0.0),
p_(0.0),
- pstar_(0.0)
+ pstar_(0.0),
+ erm_(1.0)
{
if ( daughters_ == 0 ) {
std::cerr << "ERROR in LauAbsResonance constructor : null LauDaughters object provided" << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
nameParent_ = this->getNameParent();
nameDaug1_ = this->getNameDaug1();
nameDaug2_ = this->getNameDaug2();
nameBachelor_ = this->getNameBachelor();
massParent_ = this->getMassParent();
massDaug1_ = this->getMassDaug1();
massDaug2_ = this->getMassDaug2();
massBachelor_ = this->getMassBachelor();
chargeParent_ = this->getChargeParent();
chargeDaug1_ = this->getChargeDaug1();
chargeDaug2_ = this->getChargeDaug2();
chargeBachelor_ = this->getChargeBachelor();
// check that the total charge adds up to that of the resonance
Int_t totalCharge = chargeDaug1_ + chargeDaug2_;
if ( (totalCharge != resCharge_) && (resPairAmpInt_ != 0) ) {
std::cerr << "ERROR in LauAbsResonance : Total charge of daughters = " << totalCharge << ". Resonance charge = " << resCharge_ << "." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
}
// Destructor
LauAbsResonance::~LauAbsResonance()
{
}
LauComplex LauAbsResonance::amplitude(const LauKinematics* kinematics)
{
// Use LauKinematics interface for amplitude
Double_t mass(0.0), cosHel(0.0);
// For resonance made from tracks i, j, we need the momenta
// of tracks i and k in the i-j rest frame for spin helicity calculations
// in the Zemach tensor formalism.
// Also need the momentum of track k in the parent rest-frame for
// calculation of the Blatt-Weisskopf factors.
q_ = 0.0; p_ = 0.0; pstar_ = 0.0;
+ erm_ = 1.0;
if (resPairAmpInt_ == 1) {
mass = kinematics->getm23();
cosHel = kinematics->getc23();
q_ = kinematics->getp2_23();
p_ = kinematics->getp1_23();
pstar_ = kinematics->getp1_Parent();
+ erm_ = kinematics->getcov23();
} else if (resPairAmpInt_ == 2) {
mass = kinematics->getm13();
cosHel = kinematics->getc13();
q_ = kinematics->getp1_13();
p_ = kinematics->getp2_13();
pstar_ = kinematics->getp2_Parent();
+ erm_ = kinematics->getcov13();
} else if (resPairAmpInt_ == 3) {
mass = kinematics->getm12();
cosHel = kinematics->getc12();
q_ = kinematics->getp1_12();
p_ = kinematics->getp3_12();
pstar_ = kinematics->getp3_Parent();
+ erm_ = kinematics->getcov12();
} else {
std::cerr << "ERROR in LauAbsResonance::amplitude : Nonsense setup of resPairAmp array." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
if (this->flipHelicity()) {
cosHel *= -1.0;
}
if (this->ignoreMomenta()) {
q_ = 1.0;
p_ = 1.0;
}
// Calculate the spin factors
Double_t spinTerm(1.0);
+ Double_t covFactor(1.0);
if (!this->ignoreSpin()) {
Double_t pProd = q_*p_;
spinTerm = this->calcSpinTerm( cosHel, pProd );
+ covFactor = this->calcCovFactor( erm_ );
}
+
+ if (!this->ignoreCovariant()) {
+ covFactor = this->calcCovFactor( erm_ );
+ }
+
// Calculate the full amplitude
- LauComplex resAmplitude = this->resAmp(mass, spinTerm);
+ LauComplex resAmplitude = this->resAmp(mass, spinTerm*covFactor);
return resAmplitude;
}
Double_t LauAbsResonance::calcSpinTerm( const Double_t cosHel, const Double_t pProd ) const
{
// Calculate the spin factors
//
// These are calculated as follows
//
// -2^j * (q*p)^j * cj * Pj(cosHel)
//
// where Pj(coshHel) is the jth order Legendre polynomial and
//
// cj = j! / (2j-1)!!
Double_t spinTerm = 1.0;
if (resSpin_ == 1) {
// Calculate vector resonance Zemach helicity factor
spinTerm = -2.0*pProd*cosHel;
} else if (resSpin_ == 2) {
// Calculate tensor resonance Zemach helicity factor
spinTerm = 4.0*(pProd*pProd)*(3.0*cosHel*cosHel - 1.0)/3.0;
} else if (resSpin_ == 3) {
// Calculate spin 3 resonance Zemach helicity factor
spinTerm = -8.0*(pProd*pProd*pProd)*(5.0*cosHel*cosHel*cosHel - 3.0*cosHel)/5.0;
} else if (resSpin_ == 4) {
// Calculate spin 4 resonance Zemach helicity factor
spinTerm = 16.0*(pProd*pProd*pProd*pProd)*(35.0*cosHel*cosHel*cosHel*cosHel - 30.0*cosHel*cosHel + 3.0)/35.0;
} else if (resSpin_ == 5) {
// Calculate spin 5 resonance Zemach helicity factor
spinTerm = -32.0*(pProd*pProd*pProd*pProd*pProd)*(63.0*cosHel*cosHel*cosHel*cosHel*cosHel - 70.0*cosHel*cosHel*cosHel + 15.0*cosHel)/63.0;
}
return spinTerm;
}
void LauAbsResonance::changeResonance(const Double_t newMass, const Double_t newWidth, const Int_t newSpin)
{
if (newMass > 0.0) {
resMass_->valueAndRange(newMass,0.0,3.0*newMass);
resMass_->initValue(newMass);
resMass_->genValue(newMass);
std::cout << "INFO in LauAbsResonance::changeResonance : Setting mass to " << resMass_->value() << std::endl;
}
if (newWidth > 0.0) {
resWidth_->valueAndRange(newWidth,0.0,3.0*newWidth);
resWidth_->initValue(newWidth);
resWidth_->genValue(newWidth);
std::cout << "INFO in LauAbsResonance::changeResonance : Setting width to " << resWidth_->value() << std::endl;
}
if (newSpin > -1) {
resSpin_ = newSpin;
std::cout << "INFO in LauAbsResonance::changeResonance : Setting spin to " << resSpin_ << std::endl;
}
}
void LauAbsResonance::changeBWBarrierRadii(const Double_t resRadius, const Double_t parRadius)
{
if ( resRadius >= 0.0 && resBWFactor_ != 0 ) {
LauParameter* resBWRadius = resBWFactor_->getRadiusParameter();
resBWRadius->value(resRadius);
resBWRadius->initValue(resRadius);
resBWRadius->genValue(resRadius);
std::cout << "INFO in LauAbsResonance::changeBWBarrierRadii : Setting resonance factor radius to " << resBWRadius->value() << std::endl;
}
if ( parRadius >= 0.0 && parBWFactor_ != 0 ) {
LauParameter* parBWRadius = parBWFactor_->getRadiusParameter();
parBWRadius->value(parRadius);
parBWRadius->initValue(parRadius);
parBWRadius->genValue(parRadius);
std::cout << "INFO in LauAbsResonance::changeBWBarrierRadii : Setting parent factor radius to " << parBWRadius->value() << std::endl;
}
}
void LauAbsResonance::setResonanceParameter(const TString& name, const Double_t value)
{
//This function should always be overwritten if needed in classes inheriting from LauAbsResonance.
std::cerr << "WARNING in LauAbsResonance::setResonanceParameter : Unable to set parameter \"" << name << "\" to value: " << value << "." << std::endl;
}
void LauAbsResonance::floatResonanceParameter(const TString& name)
{
//This function should always be overwritten if needed in classes inheriting from LauAbsResonance.
std::cerr << "WARNING in LauAbsResonance::floatResonanceParameter : Unable to release parameter \"" << name << "\"." << std::endl;
}
LauParameter* LauAbsResonance::getResonanceParameter(const TString& name)
{
//This function should always be overwritten if needed in classes inheriting from LauAbsResonance.
std::cerr << "WARNING in LauAbsResonance::getResonanceParameter : Unable to get parameter \"" << name << "\"." << std::endl;
return 0;
}
void LauAbsResonance::addFloatingParameter( LauParameter* param )
{
if ( param == 0 ) {
return;
}
if ( param->clone() ) {
resParameters_.push_back( param->parent() );
} else {
resParameters_.push_back( param );
}
}
void LauAbsResonance::fixBarrierRadii(const Bool_t fixResRad, const Bool_t fixParRad)
{
if ( resBWFactor_ == 0 ) {
std::cerr << "WARNING in LauAbsResonance::fixBarrierRadii : resonance barrier factor not present, cannot fix/float it" << std::endl;
return;
}
if ( parBWFactor_ == 0 ) {
std::cerr << "WARNING in LauAbsResonance::fixBarrierRadii : parent barrier factor not present, cannot fix/float it" << std::endl;
return;
}
LauParameter* resBWRadius = resBWFactor_->getRadiusParameter();
resBWRadius->fixed(fixResRad);
LauParameter* parBWRadius = parBWFactor_->getRadiusParameter();
parBWRadius->fixed(fixParRad);
}
Bool_t LauAbsResonance::fixResRadius() const
{
if ( resBWFactor_ == 0 ) {
std::cerr << "WARNING in LauAbsResonance::fixResRadius : resonance barrier factor not present" << std::endl;
return kTRUE;
}
LauParameter* bwRadius = resBWFactor_->getRadiusParameter();
return bwRadius->fixed();
}
Bool_t LauAbsResonance::fixParRadius() const
{
if ( parBWFactor_ == 0 ) {
std::cerr << "WARNING in LauAbsResonance::fixParRadius : parent barrier factor not present" << std::endl;
return kTRUE;
}
LauParameter* bwRadius = parBWFactor_->getRadiusParameter();
return bwRadius->fixed();
}
Double_t LauAbsResonance::getResRadius() const
{
if ( resBWFactor_ == 0 ) {
std::cerr << "WARNING in LauAbsResonance::getResRadius : resonance barrier factor not present" << std::endl;
return -1.0;
}
LauParameter* bwRadius = resBWFactor_->getRadiusParameter();
return bwRadius->unblindValue();
}
Double_t LauAbsResonance::getParRadius() const
{
if ( parBWFactor_ == 0 ) {
std::cerr << "WARNING in LauAbsResonance::getParRadius : parent barrier factor not present" << std::endl;
return -1.0;
}
LauParameter* bwRadius = parBWFactor_->getRadiusParameter();
return bwRadius->unblindValue();
}
Double_t LauAbsResonance::getMassParent() const
{
// Get the parent mass
Double_t mass(LauConstants::mB);
if (daughters_) {
mass = daughters_->getMassParent();
}
return mass;
}
Double_t LauAbsResonance::getMassDaug1() const
{
// Get the daughter mass
Double_t mass(LauConstants::mPi);
if (daughters_) {
if (resPairAmpInt_ == 1) {
mass = daughters_->getMassDaug2();
} else if (resPairAmpInt_ == 2) {
mass = daughters_->getMassDaug1();
} else if (resPairAmpInt_ == 3) {
mass = daughters_->getMassDaug1();
}
}
return mass;
}
Double_t LauAbsResonance::getMassDaug2() const
{
// Get the daughter mass
Double_t mass(LauConstants::mPi);
if (daughters_) {
if (resPairAmpInt_ == 1) {
mass = daughters_->getMassDaug3();
} else if (resPairAmpInt_ == 2) {
mass = daughters_->getMassDaug3();
} else if (resPairAmpInt_ == 3) {
mass = daughters_->getMassDaug2();
}
}
return mass;
}
Double_t LauAbsResonance::getMassBachelor() const
{
// Get the bachelor mass
Double_t mass(LauConstants::mPi);
if (daughters_) {
if (resPairAmpInt_ == 1) {
mass = daughters_->getMassDaug1();
} else if (resPairAmpInt_ == 2) {
mass = daughters_->getMassDaug2();
} else if (resPairAmpInt_ == 3) {
mass = daughters_->getMassDaug3();
}
}
return mass;
}
Int_t LauAbsResonance::getChargeParent() const
{
// Get the parent charge
Int_t charge(0);
if (daughters_) {
charge = daughters_->getChargeParent();
}
return charge;
}
Int_t LauAbsResonance::getChargeDaug1() const
{
// Get the daughter charge
Int_t charge(0);
if (daughters_) {
if (resPairAmpInt_ == 1) {
charge = daughters_->getChargeDaug2();
} else if (resPairAmpInt_ == 2) {
charge = daughters_->getChargeDaug1();
} else if (resPairAmpInt_ == 3) {
charge = daughters_->getChargeDaug1();
}
}
return charge;
}
Int_t LauAbsResonance::getChargeDaug2() const
{
// Get the daughter charge
Int_t charge(0);
if (daughters_) {
if (resPairAmpInt_ == 1) {
charge = daughters_->getChargeDaug3();
} else if (resPairAmpInt_ == 2) {
charge = daughters_->getChargeDaug3();
} else if (resPairAmpInt_ == 3) {
charge = daughters_->getChargeDaug2();
}
}
return charge;
}
Int_t LauAbsResonance::getChargeBachelor() const
{
// Get the bachelor charge
Int_t charge(0);
if (daughters_) {
if (resPairAmpInt_ == 1) {
charge = daughters_->getChargeDaug1();
} else if (resPairAmpInt_ == 2) {
charge = daughters_->getChargeDaug2();
} else if (resPairAmpInt_ == 3) {
charge = daughters_->getChargeDaug3();
}
}
return charge;
}
TString LauAbsResonance::getNameParent() const
{
// Get the parent name
TString name("");
if (daughters_) {
name = daughters_->getNameParent();
}
return name;
}
TString LauAbsResonance::getNameDaug1() const
{
// Get the daughter name
TString name("");
if (daughters_) {
if (resPairAmpInt_ == 1) {
name = daughters_->getNameDaug2();
} else if (resPairAmpInt_ == 2) {
name = daughters_->getNameDaug1();
} else if (resPairAmpInt_ == 3) {
name = daughters_->getNameDaug1();
}
}
return name;
}
TString LauAbsResonance::getNameDaug2() const
{
// Get the daughter name
TString name("");
if (daughters_) {
if (resPairAmpInt_ == 1) {
name = daughters_->getNameDaug3();
} else if (resPairAmpInt_ == 2) {
name = daughters_->getNameDaug3();
} else if (resPairAmpInt_ == 3) {
name = daughters_->getNameDaug2();
}
}
return name;
}
TString LauAbsResonance::getNameBachelor() const
{
// Get the bachelor name
TString name("");
if (daughters_) {
if (resPairAmpInt_ == 1) {
name = daughters_->getNameDaug1();
} else if (resPairAmpInt_ == 2) {
name = daughters_->getNameDaug2();
} else if (resPairAmpInt_ == 3) {
name = daughters_->getNameDaug3();
}
}
return name;
}
+Double_t LauAbsResonance::calcCovFactor( const Double_t Erm) const
+{
+ Double_t covFactor = 1.0;
+
+ if (resSpin_ == 1) {
+ covFactor = Erm;
+ } else if (resSpin_ == 2) {
+ covFactor = Erm*Erm + 0.5;
+ } else if (resSpin_ == 3) {
+ covFactor = Erm*(Erm*Erm + 1.5);
+ } else if (resSpin_ == 4) {
+ covFactor = (8.*Erm*Erm*Erm*Erm + 24.*Erm*Erm + 3.)/35.;
+ } else {
+ std::cerr << "cov factor not calculated for spin 5 "<<std::endl;
+ }
+
+ return covFactor;
+}
+
diff --git a/src/LauKinematics.cc b/src/LauKinematics.cc
index b93bc41..c14884f 100644
--- a/src/LauKinematics.cc
+++ b/src/LauKinematics.cc
@@ -1,675 +1,677 @@
// Copyright University of Warwick 2004 - 2013.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
// Authors:
// Thomas Latham
// John Back
// Paul Harrison
/*! \file 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(Double_t m1, Double_t m2, Double_t m3, Double_t mParent, Bool_t calcSquareDPCoords) :
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(Double_t m13Sq, 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(Double_t mPrime, 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(Double_t m13Sq, 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(Double_t mPrime, 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.
// The notation is confusing for the helicity angle cosines:
// 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
// Make sure that setMassSqVars is run first.
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(Double_t mijSq, Double_t mikSq, Double_t mij, Int_t i, Int_t j, Int_t k) const
{
// Routine to calculate the cos(helicity) variables from the masses of the particles.
// 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 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(Double_t mijSq, Double_t cij, Double_t mij, Int_t i, Int_t j, 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 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*cij;
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(Double_t m13Sq, 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(Double_t m13Sq, 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(Double_t mPrime, 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(Double_t firstMassSq, 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(Double_t m13Sq, 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(Double_t energy, 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(Double_t m23, Double_t c23)
{
// Update the variables m12Sq_ and m13Sq_ given m23 and c23.
m23_ = m23; m23Sq_ = m23*m23; c23_ = c23;
Int_t zero(0), one(1), two(2);
m13Sq_ = this->mFromC(m23Sq_, -c23_, m23_, two, one, zero);
m12Sq_ = mParentSq_ - m23Sq_ - m13Sq_ + mSqDTot_;
m13_ = TMath::Sqrt(m13Sq_);
m12_ = TMath::Sqrt(m12Sq_);
//c12_ = cFromM(m12Sq_, m13Sq_, m12_, zero, one, two);
//c13_ = cFromM(m13Sq_, m23Sq_, m13_, two, zero, one);
}
void LauKinematics::updateMassSq_m13(Double_t m13, Double_t c13)
{
// Update the variables m12Sq_ and m23Sq_ given m13 and c13.
m13_ = m13; m13Sq_ = m13*m13; c13_ = c13;
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_);
//c23_ = cFromM(m23Sq_, m12Sq_, m23_, one, two, zero);
//c12_ = cFromM(m12Sq_, m13Sq_, m12_, zero, one, two);
}
void LauKinematics::updateMassSq_m12(Double_t m12, Double_t c12)
{
// Update the variables m23Sq_ and m13Sq_ given m12 and c12.
m12_ = m12; m12Sq_ = m12*m12; c12_ = c12;
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_);
//c23_ = cFromM(m23Sq_, m12Sq_, m23_, one, two, zero);
//c13_ = cFromM(m13Sq_, m23Sq_, m13_, two, zero, one);
}
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;
}
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;
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 6:34 PM (1 d, 14 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805620
Default Alt Text
(77 KB)

Event Timeline