diff --git a/doc/ReleaseNotes.md b/doc/ReleaseNotes.md
--- a/doc/ReleaseNotes.md
+++ b/doc/ReleaseNotes.md
@@ -1,5 +1,8 @@
# Laura++ release notes
+4th February 2020 Dan Johnson
+* Extend the K-matrix implementation to handle non-zero spin
+
1st February 2021 Dan Johnson
* Allow floating of parameters in the K-matrix
- see https://phab.hepforge.org/T59
diff --git a/inc/LauAbsResonance.hh b/inc/LauAbsResonance.hh
--- a/inc/LauAbsResonance.hh
+++ b/inc/LauAbsResonance.hh
@@ -64,7 +64,7 @@
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 */
+ 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 */
@@ -117,8 +117,9 @@
\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);
+ LauAbsResonance(const TString& resName, const Int_t resPairAmpInt, const LauDaughters* daughters, const Int_t resSpin);
//! Destructor
virtual ~LauAbsResonance();
@@ -520,7 +521,7 @@
std::vector<LauParameter*> resParameters_;
//! Resonance spin
- Int_t resSpin_{0};
+ Int_t resSpin_;
//! Resonance charge
Int_t resCharge_{0};
//! DP axis identifier
diff --git a/inc/LauBlattWeisskopfFactor.hh b/inc/LauBlattWeisskopfFactor.hh
--- a/inc/LauBlattWeisskopfFactor.hh
+++ b/inc/LauBlattWeisskopfFactor.hh
@@ -90,7 +90,7 @@
/*!
\param newSpin the value of the spin to use for the created instance
*/
- LauBlattWeisskopfFactor* createClone( const UInt_t newSpin );
+ LauBlattWeisskopfFactor* createClone( const UInt_t newSpin , const BarrierType newBarrierType );
//! Retrieve the radius parameter
const LauParameter* getRadiusParameter() const { return radius_; }
@@ -116,7 +116,7 @@
private:
//! Copy constructor
- LauBlattWeisskopfFactor( const LauBlattWeisskopfFactor& other, const UInt_t newSpin );
+ LauBlattWeisskopfFactor( const LauBlattWeisskopfFactor& other, const UInt_t newSpin, const BarrierType newBarrierType );
//! Copy assignment operator (not implemented)
LauBlattWeisskopfFactor& operator=( const LauBlattWeisskopfFactor& other );
diff --git a/inc/LauIsobarDynamics.hh b/inc/LauIsobarDynamics.hh
--- a/inc/LauIsobarDynamics.hh
+++ b/inc/LauIsobarDynamics.hh
@@ -177,8 +177,8 @@
\param [in] nPoles the number of poles
\param [in] rowIndex the index of the row to be used when summing over all amplitude channels: S-wave corresponds to rowIndex = 1.
*/
- void defineKMatrixPropagator(const TString& propName, const TString& paramFileName,
- Int_t resPairAmpInt, Int_t nChannels, Int_t nPoles, Int_t rowIndex = 1);
+ LauKMatrixPropagator* defineKMatrixPropagator( const TString& propName, const TString& paramFileName,
+ Int_t resPairAmpInt, Int_t nChannels, Int_t nPoles, Int_t rowIndex = 1);
//! Add a K-matrix production pole term to the model
/*!
diff --git a/inc/LauKMatrixProdPole.hh b/inc/LauKMatrixProdPole.hh
--- a/inc/LauKMatrixProdPole.hh
+++ b/inc/LauKMatrixProdPole.hh
@@ -64,13 +64,6 @@
// Initialise the model
virtual void initialise() {return;}
- //! The amplitude calculation
- /*!
- \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
@@ -79,13 +72,18 @@
//! Retrieve the resonance parameters, e.g. so that they can be loaded into a fit
/*!
- \return floating parameters of the resonance
+ \return floating parameters of the resonance
*/
virtual const std::vector<LauParameter*>& getFloatingParameters();
protected:
- //! Function not meant to be called, amplitude is called directly in this case
- virtual LauComplex resAmp(Double_t mass, Double_t spinTerm);
+ //! The amplitude calculation
+ /*!
+ \param [in] mass the invariant-mass for the channel
+ \param [in] spinTerm the spin-term for the final channel
+ \return the complex amplitude
+ */
+ virtual LauComplex resAmp(const Double_t mass, const Double_t spinTerm);
private:
//! Copy constructor (not implemented)
diff --git a/inc/LauKMatrixProdSVP.hh b/inc/LauKMatrixProdSVP.hh
--- a/inc/LauKMatrixProdSVP.hh
+++ b/inc/LauKMatrixProdSVP.hh
@@ -44,7 +44,7 @@
class LauKMatrixProdSVP : public LauAbsResonance {
- public:
+ public:
//! Constructor
/*!
\param [in] SVPName name of the slowly varying part (SVP)
@@ -53,39 +53,37 @@
\param [in] propagator a K-matrix propagator
\param [in] daughters the daughter particles
\param [in] useProdAdler boolean to turn on/off the production Adler zero factor
- */
- LauKMatrixProdSVP(const TString& SVPName, Int_t channelIndex, Int_t resPairAmpInt,
- LauKMatrixPropagator* propagator, const LauDaughters* daughters,
- Bool_t useProdAdler = kFALSE);
+ */
+ LauKMatrixProdSVP( const TString& SVPName, Int_t channelIndex, Int_t resPairAmpInt,
+ LauKMatrixPropagator* propagator, const LauDaughters* daughters,
+ Bool_t useProdAdler = kFALSE);
//! Destructor
virtual ~LauKMatrixProdSVP();
//! Initialise the model
- virtual void initialise() {return;}
-
- //! The amplitude calculation
- /*!
- \param [in] kinematics the kinematic variables of the current event
- \return the complex amplitude
- */
- virtual LauComplex amplitude(const LauKinematics* kinematics);
+ virtual void initialise() {return;}
//! Get the resonance model type
- /*!
- \return the resonance model type
- */
+ /*!
+ \return the resonance model type
+ */
virtual LauAbsResonance::LauResonanceModel getResonanceModel() const {return LauAbsResonance::KMatrix;}
//! Retrieve the resonance parameters, e.g. so that they can be loaded into a fit
/*!
- \return floating parameters of the resonance
+ \return floating parameters of the resonance
*/
const std::vector<LauParameter*>& getFloatingParameters();
-
+
protected:
- //! Function not meant to be called
- virtual LauComplex resAmp(Double_t mass, Double_t spinTerm);
+ //! The amplitude calculation
+ /*!
+ \param [in] mass the invariant-mass for the channel
+ \param [in] spinTerm the spin-term for the final channel
+ \return the complex amplitude
+ */
+ virtual LauComplex resAmp(const Double_t mass, const Double_t spinTerm);
private:
//! Copy constructor (not implemented)
@@ -95,15 +93,15 @@
LauKMatrixProdSVP& operator=(const LauKMatrixProdSVP& rhs);
//! The K-matrix propagator
- LauKMatrixPropagator* thePropagator_;
+ LauKMatrixPropagator* thePropagator_;
//! The number of the channel
Int_t channelIndex_;
- //! Boolean to turn on/off the production Adler zero factor
- Bool_t useProdAdler_;
+ //! Boolean to turn on/off the production Adler zero factor
+ Bool_t useProdAdler_;
ClassDef(LauKMatrixProdSVP, 0) // K-matrix production SVP term
};
-#endif
+#endif
\ No newline at end of file
diff --git a/inc/LauKMatrixPropagator.hh b/inc/LauKMatrixPropagator.hh
--- a/inc/LauKMatrixPropagator.hh
+++ b/inc/LauKMatrixPropagator.hh
@@ -6,7 +6,7 @@
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
+ 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,
@@ -23,23 +23,23 @@
*/
/*! \file LauKMatrixPropagator.hh
- \brief File containing declaration of LauKMatrixPropagator class.
+ \brief File containing declaration of LauKMatrixPropagator class.
*/
/*! \class LauKMatrixPropagator
- \brief Class for defining a K-matrix propagator.
+ \brief Class for defining a K-matrix propagator.
- Class used to define a K-matrix propagator.
- See the following papers for info:
- hep-ph/0204328, hep-ex/0312040, [hep-ex]0804.2089 and hep-ph/9705401.
+ Class used to define a K-matrix propagator.
+ See the following papers for info:
+ hep-ph/0204328, hep-ex/0312040, [hep-ex]0804.2089 and hep-ph/9705401.
*/
#ifndef LAU_KMATRIX_PROPAGATOR
#define LAU_KMATRIX_PROPAGATOR
-#include "LauComplex.hh"
-#include "LauKinematics.hh"
-#include "LauParameter.hh"
+#include "LauConstants.hh"
+#include "LauResonanceMaker.hh"
+#include "LauResonanceInfo.hh"
#include "TMatrixD.h"
#include "TString.h"
@@ -47,6 +47,10 @@
#include <map>
#include <vector>
+class LauParameter;
+class LauKinematics;
+class LauComplex;
+
class LauKMatrixPropagator {
public:
@@ -60,23 +64,17 @@
\param [in] rowIndex this specifies which row of the propagator should be used when summing over the amplitude channels
*/
LauKMatrixPropagator( const TString& name, const TString& paramFileName,
- Int_t resPairAmpInt, Int_t nChannels, Int_t nPoles,
- Int_t rowIndex = 1 );
+ const Int_t resPairAmpInt, const Int_t nChannels, const Int_t nPoles,
+ const Int_t rowIndex = 1 );
//! Destructor
- virtual ~LauKMatrixPropagator();
-
- //! Calculate the invariant mass squared s
- /*!
- \param [in] kinematics the kinematics of the current event
- */
- void updatePropagator(const LauKinematics* kinematics);
+ virtual ~LauKMatrixPropagator();
//! Calculate the K-matrix propagator for the given s value
/*!
\param [in] s the invariant mass squared
*/
- void updatePropagator(Double_t s);
+ void updatePropagator(const Double_t s);
//! Read an input file to set parameters
/*!
@@ -84,6 +82,9 @@
*/
void setParameters(const TString& inputFile);
+ //! Set flag to ignore Blatt-Weisskopf-like barrier factor
+ void ignoreBWBarrierFactor() {includeBWBarrierFactor_=kFALSE;}
+
//! Get the scattering K matrix
/*!
\return the real, symmetric scattering K matrix
@@ -107,29 +108,42 @@
\param [in] channelIndex the channel number
\return the real part of the propagator term
*/
- Double_t getRealPropTerm(Int_t channelIndex) const;
+ Double_t getRealPropTerm(const Int_t channelIndex) const;
//! Get the imaginary part of the term of the propagator
/*!
\param [in] channelIndex the channel number
\return the imaginiary part of the propagator term
*/
- Double_t getImagPropTerm(Int_t channelIndex) const;
+ Double_t getImagPropTerm(const Int_t channelIndex) const;
//! Get the 1/(m_pole^2 -s) terms for the scattering and production K-matrix formulae
/*!
\param [in] poleIndex the number of the pole required
\return the value of 1/(m_pole^2 -s)
*/
- Double_t getPoleDenomTerm(Int_t poleIndex) const;
+ Double_t getPoleDenomTerm(const Int_t poleIndex) const;
- //! Get coupling parameters, set according to the input file
+ //! Get spin of K-matrix
+ /*!
+ \param [in] iChannel the index of the channel
+ \return the value of the orbital angular momentum, L_, for this channel
+ */
+ Int_t getL(const Int_t iChannel) const {return L_[iChannel];}
+
+ //! Get index of final channel
+ /*!
+ \return the index of the channel into which the scattering happens
+ */
+ Int_t getIndex() const {return index_;};
+
+ //! Get pole mass parameters, set according to the input file
/*!
\param [in] poleIndex number of the required pole
\param [in] channelIndex number of the required channel
- \return the parameter of the coupling constant
+ \return the parameter of the pole mass
*/
- LauParameter& getPoleMassSqParameter(Int_t poleIndex);
+ LauParameter& getPoleMassSqParameter(const Int_t poleIndex);
//! Get coupling constants that were loaded from the input file
/*!
@@ -137,7 +151,7 @@
\param [in] channelIndex number of the required channel
\return the value of the coupling constant
*/
- Double_t getCouplingConstant(Int_t poleIndex, Int_t channelIndex) const;
+ Double_t getCouplingConstant(const Int_t poleIndex, const Int_t channelIndex) const;
//! Get coupling parameters, set according to the input file
/*!
@@ -145,7 +159,7 @@
\param [in] channelIndex number of the required channel
\return the parameter of the coupling constant
*/
- LauParameter& getCouplingParameter(Int_t poleIndex, Int_t channelIndex);
+ LauParameter& getCouplingParameter(const Int_t poleIndex, const Int_t channelIndex);
//! Get scattering constants that were loaded from the input file
/*!
@@ -153,7 +167,7 @@
\param [in] channel2Index number of the second channel index
\return the value of the scattering constant
*/
- Double_t getScatteringConstant(Int_t channel1Index, Int_t channel2Index) const;
+ Double_t getScatteringConstant(const Int_t channel1Index, const Int_t channel2Index) const;
//! Get scattering parameters, set according to the input file
/*!
@@ -161,7 +175,7 @@
\param [in] channel2Index number of the second channel index
\return the parameter of the scattering constant
*/
- LauParameter& getScatteringParameter(Int_t channel1Index, Int_t channel2Index);
+ LauParameter& getScatteringParameter(const Int_t channel1Index, const Int_t channel2Index);
//! Get mSq0 production parameter
/*!
@@ -185,13 +199,13 @@
/*!
\return the sA parameter
*/
- LauParameter& getsA() {return sA_;}
+ LauParameter& getsA() {return sA_;}
//! Get sA0 production parameter
/*!
\return the sA0 parameter
*/
- LauParameter& getsA0() {return sA0_;}
+ LauParameter& getsA0() {return sA0_;}
//! Get the "slowly-varying part" term of the amplitude
/*!
@@ -204,7 +218,7 @@
\param [in] channelIndex the number of the required channel
\return the complex propagator term
*/
- LauComplex getPropTerm(Int_t channelIndex) const;
+ LauComplex getPropTerm(const Int_t channelIndex) const;
//! Get the DP axis identifier
/*!
@@ -236,105 +250,146 @@
\param [in] channel The index number of the channel process
\return the complex amplitude T
*/
- LauComplex getTransitionAmp(Double_t s, Int_t channel);
+ LauComplex getTransitionAmp(const Double_t s, const Int_t channel);
- //! Get the complex phase space term for the given channel and invariant mass squared
- /*!
- \param [in] s The invariant mass squared
- \param [in] channel The index number of the channel process
- \return the complex phase space term rho(channel, channel)
+ //! Get the complex phase space term for the given channel and invariant mass squared
+ /*!
+ \param [in] s The invariant mass squared
+ \param [in] channel The index number of the channel process
+ \return the complex phase space term rho(channel, channel)
*/
- LauComplex getPhaseSpaceTerm(Double_t s, Int_t channel);
+ LauComplex getPhaseSpaceTerm(const Double_t s, const Int_t channel);
- //! Get the Adler zero factor, which is set when updatePropagator is called
- /*!
- \return the Adler zero factor
+ //! Get the Adler zero factor, which is set when updatePropagator is called
+ /*!
+ \return the Adler zero factor
*/
- Double_t getAdlerZero() const {return adlerZeroFactor_;}
-
+ Double_t getAdlerZero() const {return adlerZeroFactor_;}
- //! Get the THat amplitude for the given s and channel number
- /*!
- \param [in] s The invariant mass squared
+ //! Get the THat amplitude for the given s and channel number
+ /*!
+ \param [in] s The invariant mass squared
\param [in] channel The index number of the channel process
\return the complex THat amplitude
*/
- LauComplex getTHat(Double_t s, Int_t channel);
+ LauComplex getTHat(const Double_t s, const Int_t channel);
protected:
+ // Integers to specify the allowed channels for the phase space calculations.
+ // Please keep Zero at the start and leave TotChannels at the end
+ // whenever more channels are added to this.
+ //! Integers to specify the allowed channels for the phase space calculations
+ enum class KMatrixChannels {Zero, PiPi, KK, FourPi, EtaEta, EtaEtaP,
+ KPi, KEtaP, KThreePi, D0K, Dstar0K, TotChannels};
+
//! Calculate the scattering K-matrix for the given value of s
/*!
\param [in] s the invariant mass squared
*/
- void calcScattKMatrix(Double_t s);
+ void calcScattKMatrix(const Double_t s);
//! Calculate the real and imaginary part of the phase space density diagonal matrix
/*!
\param [in] s the invariant mass squared
*/
- void calcRhoMatrix(Double_t s);
+ void calcRhoMatrix(const Double_t s);
+
+ //! Retrieve the complex phasespace density for a given channel
+ /*!
+ \param [in] s the invariant mass squared
+ \param [in] phaseSpaceIndex the phasespace index of the channel
+ \return the complex phasespace density
+ */
+ LauComplex getRho(const Double_t s, const:LauKMatrixPropagator::KMatrixChannels) const
+
+ //! Calculate the (real) gamma angular-momentum barrier matrix
+ /*!
+ \param [in] s the invariant mass squared
+ */
+ void calcGammaMatrix(const Double_t s);
+
+ //! Calculate the DK P-wave gamma angular-momentum barrier
+ /*!
+ \param [in] s the invariant mass squared
+ \return the centrifugal barrier factor for L=0,1, or 2
+ */
+ Double_t calcGamma(const Int_t iCh, const Double_t s, KMatrixChannels phaseSpaceIndex) const;
//! Calulate the term 1/(m_pole^2 - s) for the scattering and production K-matrix formulae
/*!
\param [in] s the invariant mass squared
*/
- void calcPoleDenomVect(Double_t s);
+ void calcPoleDenomVect(const Double_t s);
+
+ //! Calculate the D0K+ phase space factor
+ /*!
+ \param [in] s the invariant mass squared
+ \return the complex phase space factor
+ */
+ LauComplex calcD0KRho(const Double_t s) const;
+
+ //! Calculate the D*0K+ phase space factor
+ /*!
+ \param [in] s the invariant mass squared
+ \return the complex phase space factor
+ */
+ LauComplex calcDstar0KRho(const Double_t s) const;
//! Calculate the pipi phase space factor
/*!
\param [in] s the invariant mass squared
\return the complex phase space factor
*/
- LauComplex calcPiPiRho(Double_t s) const;
+ LauComplex calcPiPiRho(const Double_t s) const;
//! Calculate the KK phase space factor
/*!
\param [in] s the invariant mass squared
\return the complex phase space factor
*/
- LauComplex calcKKRho(Double_t s) const;
+ LauComplex calcKKRho(const Double_t s) const;
//! Calculate the 4 pi phase space factor
/*!
\param [in] s the invariant mass squared
\return the complex phase space factor
*/
- LauComplex calcFourPiRho(Double_t s) const;
+ LauComplex calcFourPiRho(const Double_t s) const;
//! Calculate the eta-eta phase space factor
/*!
\param [in] s the invariant mass squared
\return the complex phase space factor
*/
- LauComplex calcEtaEtaRho(Double_t s) const;
+ LauComplex calcEtaEtaRho(const Double_t s) const;
//! Calculate the eta-eta' phase space factor
/*!
\param [in] s the invariant mass squared
\return the complex phase space factor
*/
- LauComplex calcEtaEtaPRho(Double_t s) const;
+ LauComplex calcEtaEtaPRho(const Double_t s) const;
//! Calculate the Kpi phase space factor
/*!
\param [in] s the invariant mass squared
\return the complex phase space factor
*/
- LauComplex calcKPiRho(Double_t s) const;
+ LauComplex calcKPiRho(const Double_t s) const;
//! Calculate the K-eta' phase space factor
/*!
\param [in] s the invariant mass squared
\return the complex phase space factor
*/
- LauComplex calcKEtaPRho(Double_t s) const;
+ LauComplex calcKEtaPRho(const Double_t s) const;
//! Calculate the Kpipipi phase space factor
/*!
\param [in] s the invariant mass squared
\return the complex phase space factor
*/
- LauComplex calcKThreePiRho(Double_t s) const;
+ LauComplex calcKThreePiRho(const Double_t s) const;
//! Calculate the "slow-varying part"
/*!
@@ -342,87 +397,104 @@
\param [in] s0 the invariant mass squared at the Adler zero
\return the SVP term
*/
- Double_t calcSVPTerm(Double_t s, Double_t s0) const;
+ Double_t calcSVPTerm(const Double_t s, const Double_t s0) const;
//! Update the scattering "slowly-varying part"
/*!
\param [in] s the invariant mass squared
*/
- void updateScattSVPTerm(Double_t s);
+ void updateScattSVPTerm(const Double_t s);
//! Update the production "slowly-varying part"
/*!
\param [in] s the invariant mass squared
*/
- void updateProdSVPTerm(Double_t s);
+ void updateProdSVPTerm(const Double_t s);
//! Calculate the multiplicative factor containing severa Adler zero constants
/*!
\param [in] s the invariant mass squared
*/
- void updateAdlerZeroFactor(Double_t s);
+ void updateAdlerZeroFactor(const Double_t s);
//! Check the phase space factors that need to be used
/*!
\param [in] phaseSpaceInt phase space types
\return true of false
*/
- Bool_t checkPhaseSpaceType(Int_t phaseSpaceInt) const;
+ Bool_t checkPhaseSpaceType(const Int_t phaseSpaceInt) const;
- //! Get the unitary transition amplitude matrix for the given kinematics
- /*!
- \param [in] kinematics The pointer to the constant kinematics
+ //! Get the unitary transition amplitude matrix for the given kinematics
+ /*!
+ \param [in] kinematics The pointer to the constant kinematics
*/
- void getTMatrix(const LauKinematics* kinematics);
+ void getTMatrix(const LauKinematics* kinematics);
- //! Get the unitary transition amplitude matrix for the given kinematics
- /*!
- \param [in] s The invariant mass squared of the system
+ //! Get the unitary transition amplitude matrix for the given kinematics
+ /*!
+ \param [in] s The invariant mass squared of the system
*/
- void getTMatrix(Double_t s);
+ void getTMatrix(const Double_t s);
+
+ //! Get the square root of the phase space matrix
+ void getSqrtRhoMatrix();
- //! Get the square root of the phase space matrix
- void getSqrtRhoMatrix();
-
- private:
+ private:
//! Copy constructor (not implemented)
- LauKMatrixPropagator(const LauKMatrixPropagator& rhs);
+ LauKMatrixPropagator(const LauKMatrixPropagator& rhs)=delete;
//! Copy assignment operator (not implemented)
- LauKMatrixPropagator& operator=(const LauKMatrixPropagator& rhs);
+ LauKMatrixPropagator& operator=(const LauKMatrixPropagator& rhs)=delete;
- //! Create a map for the K-matrix parameters
- typedef std::map<int, std::vector<LauParameter> > KMatrixParamMap;
+ //! Initialise and set the dimensions for the internal matrices and parameter arrays
+ void initialiseMatrices();
+ //! Store the (phase space) channel indices from a line in the parameter file
+ /*!
+ \param [in] theLine Vector of strings corresponding to the line from the parameter file
+ */
+ void storeChannels(const std::vector<std::string>& theLine);
+
+ //! Store the pole mass and couplings from a line in the parameter file
+ /*!
+ \param [in] theLine Vector of strings corresponding to the line from the parameter file
+ */
+ void storePole(const std::vector<std::string>& theLine);
- //! Initialise and set the dimensions for the internal matrices and parameter arrays
- void initialiseMatrices();
+ //! Store the scattering coefficients from a line in the parameter file
+ /*!
+ \param [in] theLine Vector of strings corresponding to the line from the parameter file
+ */
+ void storeScattering(const std::vector<std::string>& theLine);
- //! Store the (phase space) channel indices from a line in the parameter file
- /*!
- \param [in] theLine Vector of strings corresponding to the line from the parameter file
+ //! Store the channels' characteristic radii from a line in the parameter file
+ /*!
+ \param [in] theLine Vector of strings corresponding to the line from the parameter file
*/
- void storeChannels(const std::vector<std::string>& theLine);
-
- //! Store the pole mass and couplings from a line in the parameter file
- /*!
- \param [in] theLine Vector of strings corresponding to the line from the parameter file
+ void storeRadii(const std::vector<std::string>& theLine);
+
+ //! Store the channels' angular momenta from a line in the parameter file
+ /*!
+ \param [in] theLine Vector of strings corresponding to the line from the parameter file
+ \param [in] a Vector of integers corresponding to parameter in the propagator denominator
*/
- void storePole(const std::vector<std::string>& theLine);
+ void storeOrbitalAngularMomenta(const std::vector<std::string>& theLine, std::vector<Int_t>& a);
- //! Store the scattering coefficients from a line in the parameter file
- /*!
- \param [in] theLine Vector of strings corresponding to the line from the parameter file
+ //! Store the barrier-factor parameter from a line in the parameter file
+ /*!
+ \param [in] theLine Vector of strings corresponding to the line from the parameter file
+ \param [in] a Vector of integers corresponding to parameter in the propagator denominator
*/
- void storeScattering(const std::vector<std::string>& theLine);
+ void storeBarrierFactorParameter(const std::vector<std::string>& theLine, std::vector<Int_t>& a);
- //! Store miscelleanous parameters from a line in the parameter file
- /*!
- \param [in] keyword the name of the parameter to be set
- \param [in] parString the string containing the value of the parameter
+
+ //! Store miscelleanous parameters from a line in the parameter file
+ /*!
+ \param [in] keyword the name of the parameter to be set
+ \param [in] parString the string containing the value of the parameter
*/
- void storeParameter(const TString& keyword, const TString& parString);
+ void storeParameter(const TString& keyword, const TString& parString);
//! String to store the propagator name
@@ -435,29 +507,24 @@
Int_t index_;
//! s value of the previous pole
- Double_t previousS_;
+ Double_t previousS_{0.0};
//! "slowly-varying part" for the scattering K-matrix
- Double_t scattSVP_;
+ Double_t scattSVP_{0.0};
//! "slowly-varying part" for the production K-matrix
- Double_t prodSVP_;
+ Double_t prodSVP_{0.0};
//! Real part of the propagator matrix
TMatrixD realProp_;
//! Imaginary part of the propagator matrix
TMatrixD negImagProp_;
- // Integers to specify the allowed channels for the phase space calculations.
- // Please keep Zero at the start and leave TotChannels at the end
- // whenever more channels are added to this.
- //! Integers to specify the allowed channels for the phase space calculations
- enum KMatrixChannels {Zero, PiPi, KK, FourPi, EtaEta, EtaEtaP,
- KPi, KEtaP, KThreePi, TotChannels};
-
//! Scattering K-matrix
TMatrixD ScattKMatrix_;
//! Real part of the phase space density diagonal matrix
TMatrixD ReRhoMatrix_;
//! Imaginary part of the phase space density diagonal matrix
TMatrixD ImRhoMatrix_;
+ //! Gamma angular-momentum barrier matrix
+ TMatrixD GammaMatrix_;
//! Identity matrix
TMatrixD IMatrix_;
//! Null matrix
@@ -467,26 +534,36 @@
TMatrixD ReSqrtRhoMatrix_;
//! Imaginary part of the square root of the phase space density diagonal matrix
TMatrixD ImSqrtRhoMatrix_;
- //! Real part of the unitary T matrix
- TMatrixD ReTMatrix_;
- //! Imaginary part of the unitary T matrix
- TMatrixD ImTMatrix_;
+ //! Real part of the unitary T matrix
+ TMatrixD ReTMatrix_;
+ //! Imaginary part of the unitary T matrix
+ TMatrixD ImTMatrix_;
//! Number of channels
Int_t nChannels_;
//! Number of poles
Int_t nPoles_;
+ //! Vector of orbital angular momenta
+ std::vector<Int_t> L_;
+ //! Boolean to indicate whether storeBarrierFactorParameter has been called
+ Bool_t haveCalled_storeBarrierFactorParameter{kFALSE};
+ //! Boolean to dictate whether to include Blatt-Weisskopf-like denominator in K-matrix centrifugal barrier factor
+ Bool_t includeBWBarrierFactor_{kTRUE};
//! Vector of squared pole masses
- std::vector<LauParameter> mSqPoles_;
+ std::vector<LauParameter> mSqPoles_;
//! Array of coupling constants
LauParArray gCouplings_;
//! Array of scattering SVP values
LauParArray fScattering_;
+ //! Vector of characteristic radii
+ std::vector<Double_t> radii_;
+ //! Vector of ratio a/R^2
+ std::vector<Double_t> gamAInvRadSq_;
//! Vector of phase space types
- std::vector<Int_t> phaseSpaceTypes_;
+ std::vector<KMatrixChannels> phaseSpaceTypes_;
//! Vector of squared masses
std::vector<Double_t> mSumSq_;
//! Vector of mass differences
@@ -506,44 +583,53 @@
LauParameter sA0_;
//! Defined as 0.5*sA*mPi*mPi
- Double_t sAConst_;
+ Double_t sAConst_{0.0};
+
//! Defined as 4*mPi*mPi
- Double_t m2piSq_;
+ const Double_t m2piSq_{4.0*LauConstants::mPiSq};
//! Defined as 4*mK*mK
- Double_t m2KSq_;
+ const Double_t m2KSq_{4.0*LauConstants::mKSq};
//! Defined as 4*mEta*mEta
- Double_t m2EtaSq_;
+ const Double_t m2EtaSq_{4.0*LauConstants::mEtaSq};
//! Defined as (mEta+mEta')^2
- Double_t mEtaEtaPSumSq_;
+ const Double_t mEtaEtaPSumSq_{(LauConstants::mEta + LauConstants::mEtaPrime)*(LauConstants::mEta + LauConstants::mEtaPrime)};
//! Defined as (mEta-mEta')^2
- Double_t mEtaEtaPDiffSq_;
+ const Double_t mEtaEtaPDiffSq_{(LauConstants::mEta - LauConstants::mEtaPrime)*(LauConstants::mEta - LauConstants::mEtaPrime)};
//! Defined as (mK+mPi)^2
- Double_t mKpiSumSq_;
+ const Double_t mKpiSumSq_{(LauConstants::mK + LauConstants::mPi)*(LauConstants::mK + LauConstants::mPi)};
//! Defined as (mK-mPi)^2
- Double_t mKpiDiffSq_;
+ const Double_t mKpiDiffSq_{(LauConstants::mK - LauConstants::mPi)*(LauConstants::mK - LauConstants::mPi)};
//! Defined as (mK+mEta')^2
- Double_t mKEtaPSumSq_;
+ const Double_t mKEtaPSumSq_{(LauConstants::mK + LauConstants::mEtaPrime)*(LauConstants::mK + LauConstants::mEtaPrime)};
//! Defined as (mK-mEta')^2
- Double_t mKEtaPDiffSq_;
+ const Double_t mKEtaPDiffSq_{(LauConstants::mK - LauConstants::mEtaPrime)*(LauConstants::mK - LauConstants::mEtaPrime)};
//! Defined as (mK-3*mPi)^2
- Double_t mK3piDiffSq_;
+ const Double_t mK3piDiffSq_{(LauConstants::mK - 3.0*LauConstants::mPi)*(LauConstants::mK - 3.0*LauConstants::mPi)};
//! Factor used to calculate the Kpipipi phase space term
- Double_t k3piFactor_;
+ const Double_t k3piFactor_{TMath::Power((1.44 - mK3piDiffSq_)/1.44, -2.5)};
//! Factor used to calculate the pipipipi phase space term
- Double_t fourPiFactor1_;
+ const Double_t fourPiFactor1_{16.0*LauConstants::mPiSq};
//! Factor used to calculate the pipipipi phase space term
- Double_t fourPiFactor2_;
+ const Double_t fourPiFactor2_{TMath::Sqrt(1.0 - fourPiFactor1_)};
+ //! Defined as (mD0+mK)^2
+ const Double_t mD0KSumSq_{(LauConstants::mD0 + LauConstants::mK)*(LauConstants::mD0 + LauConstants::mK)};
+ //! Defined as (mD0-mK)^2
+ const Double_t mD0KDiffSq_{(LauConstants::mD0 - LauConstants::mK)*(LauConstants::mD0 - LauConstants::mK)};
+ //! Defined as (mD*0+mK)^2
+ const Double_t mDstar0KSumSq_{(LauResonanceMaker::get().getResInfo("D*0")->getMass()->value() + LauConstants::mK)*(LauResonanceMaker::get().getResInfo("D*0")->getMass()->value() + LauConstants::mK)};
+ //! Defined as (mD*0-mK)^2
+ const Double_t mDstar0KDiffSq_{(LauResonanceMaker::get().getResInfo("D*0")->getMass()->value() - LauConstants::mK)*(LauResonanceMaker::get().getResInfo("D*0")->getMass()->value() - LauConstants::mK)};
//! Multiplicative factor containing various Adler zero constants
- Double_t adlerZeroFactor_;
+ Double_t adlerZeroFactor_{0.0};
//! Tracks if all params have been set
- Bool_t parametersSet_;
+ Bool_t parametersSet_{kFALSE};
//! Control the output of the functions
- Bool_t verbose_;
+ static constexpr Bool_t verbose_{kFALSE};
//! Control if scattering constants are channel symmetric: f_ji = f_ij
- Bool_t scattSymmetry_;
+ Bool_t scattSymmetry_{kFALSE};
ClassDef(LauKMatrixPropagator,0) // K-matrix amplitude model
};
diff --git a/inc/LauResonanceMaker.hh b/inc/LauResonanceMaker.hh
--- a/inc/LauResonanceMaker.hh
+++ b/inc/LauResonanceMaker.hh
@@ -127,6 +127,9 @@
*/
LauResonanceInfo* getResInfo(const TString& resName) const;
+ //! Retrieve parent Blatt-Weisskopf factor (for use by K-matrix pole/SVP which doesn't have a `resInfo')
+ LauBlattWeisskopfFactor* getParentBWFactor(Int_t newSpin, LauBlattWeisskopfFactor::BarrierType barrierType);
+
protected:
//! Create the list of known resonances
void createResonanceVector();
diff --git a/src/LauAbsResonance.cc b/src/LauAbsResonance.cc
--- a/src/LauAbsResonance.cc
+++ b/src/LauAbsResonance.cc
@@ -125,10 +125,11 @@
}
// Constructor
-LauAbsResonance::LauAbsResonance(const TString& resName, const Int_t resPairAmpInt, const LauDaughters* daughters) :
+LauAbsResonance::LauAbsResonance(const TString& resName, const Int_t resPairAmpInt, const LauDaughters* daughters, const Int_t resSpin) :
daughters_(daughters),
resName_(resName),
sanitisedName_(resName),
+ resSpin_(resSpin),
resPairAmpInt_(resPairAmpInt)
{
if ( daughters_ == 0 ) {
diff --git a/src/LauBlattWeisskopfFactor.cc b/src/LauBlattWeisskopfFactor.cc
--- a/src/LauBlattWeisskopfFactor.cc
+++ b/src/LauBlattWeisskopfFactor.cc
@@ -62,10 +62,10 @@
{
}
-LauBlattWeisskopfFactor::LauBlattWeisskopfFactor( const LauBlattWeisskopfFactor& other, const UInt_t newSpin ) :
+LauBlattWeisskopfFactor::LauBlattWeisskopfFactor( const LauBlattWeisskopfFactor& other, const UInt_t newSpin, const BarrierType newBarrierType ) :
spin_(newSpin),
radius_(other.radius_->createClone()),
- barrierType_(other.barrierType_),
+ barrierType_(newBarrierType),
restFrame_(other.restFrame_)
{
}
@@ -129,9 +129,9 @@
return categoryName;
}
-LauBlattWeisskopfFactor* LauBlattWeisskopfFactor::createClone( const UInt_t newSpin )
+LauBlattWeisskopfFactor* LauBlattWeisskopfFactor::createClone( const UInt_t newSpin, const BarrierType newBarrierType )
{
- LauBlattWeisskopfFactor* clone = new LauBlattWeisskopfFactor( *this, newSpin );
+ LauBlattWeisskopfFactor* clone = new LauBlattWeisskopfFactor( *this, newSpin, newBarrierType );
return clone;
}
diff --git a/src/LauIsobarDynamics.cc b/src/LauIsobarDynamics.cc
--- a/src/LauIsobarDynamics.cc
+++ b/src/LauIsobarDynamics.cc
@@ -789,8 +789,8 @@
return theResonance;
}
-void LauIsobarDynamics::defineKMatrixPropagator(const TString& propName, const TString& paramFileName, Int_t resPairAmpInt,
- Int_t nChannels, Int_t nPoles, Int_t rowIndex)
+LauKMatrixPropagator* LauIsobarDynamics::defineKMatrixPropagator( const TString& propName, const TString& paramFileName, Int_t resPairAmpInt,
+ Int_t nChannels, Int_t nPoles, Int_t rowIndex)
{
// Define the K-matrix propagator. The resPairAmpInt integer specifies which mass combination should be used
// for the invariant mass-squared variable "s". The pole masses and coupling constants are defined in the
@@ -814,6 +814,7 @@
nPoles, rowIndex);
kMatrixPropagators_[propagatorName] = thePropagator;
+ return thePropagator;
}
void LauIsobarDynamics::addKMatrixProdPole(const TString& poleName, const TString& propName, Int_t poleIndex, Bool_t useProdAdler)
@@ -841,6 +842,7 @@
Int_t resPairAmpInt = thePropagator->getResPairAmpInt();
LauAbsResonance* prodPole = new LauKMatrixProdPole(poleName, poleIndex, resPairAmpInt,
thePropagator, daughters_, useProdAdler);
+ prodPole->setSpinType( LauAbsResonance::Legendre );
resTypAmp_.push_back(poleName);
resPairAmp_.push_back(resPairAmpInt);
@@ -889,6 +891,7 @@
Int_t resPairAmpInt = thePropagator->getResPairAmpInt();
LauAbsResonance* prodSVP = new LauKMatrixProdSVP(SVPName, channelIndex, resPairAmpInt,
thePropagator, daughters_, useProdAdler);
+ prodSVP->setSpinType( LauAbsResonance::Legendre );
resTypAmp_.push_back(SVPName);
resPairAmp_.push_back(resPairAmpInt);
@@ -1239,7 +1242,7 @@
width = omegaWidth;
}
- if ( width > narrowWidth_ || width == 0.0 ) {
+ if ( width > narrowWidth_ || width <= 0.0 ) {
continue;
}
diff --git a/src/LauKMatrixProdPole.cc b/src/LauKMatrixProdPole.cc
--- a/src/LauKMatrixProdPole.cc
+++ b/src/LauKMatrixProdPole.cc
@@ -28,44 +28,39 @@
#include "LauKMatrixProdPole.hh"
#include "LauKMatrixPropagator.hh"
+#include "LauResonanceMaker.hh"
#include <iostream>
ClassImp(LauKMatrixProdPole)
-LauKMatrixProdPole::LauKMatrixProdPole(const TString& poleName, Int_t poleIndex, Int_t resPairAmpInt,
- LauKMatrixPropagator* propagator, const LauDaughters* daughters,
- Bool_t useProdAdler) :
- LauAbsResonance(poleName, resPairAmpInt, daughters),
+LauKMatrixProdPole::LauKMatrixProdPole( const TString& poleName, Int_t poleIndex, Int_t resPairAmpInt,
+ LauKMatrixPropagator* propagator, const LauDaughters* daughters,
+ Bool_t useProdAdler) :
+ LauAbsResonance( poleName, resPairAmpInt, daughters, propagator->getL(propagator->getIndex()) ),
thePropagator_(propagator),
- poleIndex_(poleIndex - 1), // poleIndex goes from 1 to nPoles
- useProdAdler_(useProdAdler)
+ poleIndex_(poleIndex - 1), // poleIndex goes from 1 to nPoles
+ useProdAdler_(useProdAdler)
{
- if (useProdAdler_) {
- std::cout<<"Creating K matrix production pole "<<poleName<<" with poleIndex = "
- <<poleIndex<<" with s-dependent production Adler zero term"<<std::endl;
+ if (useProdAdler_) {
+ std::cout <<"Creating K matrix production pole "<<poleName<<" with poleIndex = "
+ <<poleIndex<<" with s-dependent production Adler zero term"<<std::endl;
} else {
- std::cout<<"Creating K matrix production pole "<<poleName<<" with poleIndex = "
- <<poleIndex<<" with production Adler zero factor = 1"<<std::endl;
+ std::cout <<"Creating K matrix production pole "<<poleName<<" with poleIndex = "
+ <<poleIndex<<" with production Adler zero factor = 1"<<std::endl;
}
+ // `Resonance' Blatt-Weisskopf factor is handled internally, but parent must be set here. For other lineshapes, LauResonanceMaker handles this.
+ this->setBarrierRadii( nullptr,LauResonanceMaker::get().getParentBWFactor( propagator->getL(propagator->getIndex()), LauBlattWeisskopfFactor::BWBarrier ) );
}
LauKMatrixProdPole::~LauKMatrixProdPole()
{
}
-LauComplex LauKMatrixProdPole::resAmp(Double_t mass, Double_t spinTerm)
+LauComplex LauKMatrixProdPole::resAmp(const Double_t mass, const Double_t spinTerm)
{
- std::cerr << "ERROR in LauKMatrixProdPole::resAmp : This method shouldn't get called." << std::endl;
- std::cerr << " Returning zero amplitude for mass = " << mass << " and spinTerm = " << spinTerm << "." << std::endl;
- return LauComplex(0.0, 0.0);
-}
-
-LauComplex LauKMatrixProdPole::amplitude(const LauKinematics* kinematics)
-{
-
// Calculate the amplitude for the K-matrix production pole.
LauComplex amp(0.0, 0.0);
@@ -74,7 +69,40 @@
return amp;
}
- thePropagator_->updatePropagator(kinematics);
+ // Get barrier factors ('resonance' factor is already accounted for internally via propagator 'Gamma' matrix)
+ Double_t fFactorB(1.0);
+
+ const Int_t resSpin = this->getSpin();
+ const Double_t pstar = this->getPstar();
+
+ if ( resSpin > 0 ) {
+ const LauBlattWeisskopfFactor* parBWFactor = this->getParBWFactor();
+ if ( parBWFactor != nullptr ) {
+ switch ( parBWFactor->getRestFrame() ) {
+ case LauBlattWeisskopfFactor::ResonanceFrame:
+ fFactorB = parBWFactor->calcFormFactor(this->getP());
+ break;
+ case LauBlattWeisskopfFactor::ParentFrame:
+ fFactorB = parBWFactor->calcFormFactor(pstar);
+ break;
+ case LauBlattWeisskopfFactor::Covariant:
+ {
+ Double_t covFactor = this->getCovFactor();
+ if ( resSpin > 2 ) {
+ covFactor = TMath::Power( covFactor, 1.0/resSpin );
+ } else if ( resSpin == 2 ) {
+ covFactor = TMath::Sqrt( covFactor );
+ }
+ fFactorB = parBWFactor->calcFormFactor(pstar*covFactor);
+ break;
+ }
+ }
+ }
+ }
+
+ // Make sure the K-matrix propagator is up-to-date for
+ // the given centre-of-mass squared value ("s")
+ thePropagator_->updatePropagator(mass*mass);
// Sum the pole denominator terms over all channels j, multiplying by
// the propagator terms. Note that we do not sum over poles, since we
@@ -101,6 +129,14 @@
amp.rescale(poleDenom*adlerZero);
+ // Scale by the spin term
+ Double_t scale = spinTerm;
+
+ // Include Blatt-Weisskopf barrier factor for parent
+ scale *= fFactorB;
+
+ amp.rescale(scale);
+
return amp;
}
@@ -125,8 +161,8 @@
if ( !par_polemasssq_.fixed() )
{
this->addFloatingParameter( &par_polemasssq_ );
- }
+ }
return this->getParameters();
-}
\ No newline at end of file
+}
diff --git a/src/LauKMatrixProdSVP.cc b/src/LauKMatrixProdSVP.cc
--- a/src/LauKMatrixProdSVP.cc
+++ b/src/LauKMatrixProdSVP.cc
@@ -33,24 +33,23 @@
ClassImp(LauKMatrixProdSVP)
-LauKMatrixProdSVP::LauKMatrixProdSVP(const TString& SVPName, Int_t channelIndex, Int_t resPairAmpInt,
- LauKMatrixPropagator* propagator, const LauDaughters* daughters,
- Bool_t useProdAdler) :
- LauAbsResonance(SVPName, resPairAmpInt, daughters),
+LauKMatrixProdSVP::LauKMatrixProdSVP( const TString& SVPName, Int_t channelIndex, Int_t resPairAmpInt,
+ LauKMatrixPropagator* propagator, const LauDaughters* daughters,
+ Bool_t useProdAdler) :
+ LauAbsResonance( SVPName, resPairAmpInt, daughters, propagator->getL(propagator->getIndex()) ),
thePropagator_(propagator),
- channelIndex_(channelIndex - 1), // channelIndex goes from 1 to nChannels.
- useProdAdler_(useProdAdler)
+ channelIndex_(channelIndex - 1), // channelIndex goes from 1 to nChannels.
+ useProdAdler_(useProdAdler)
{
// Constructor
- if (useProdAdler_) {
- std::cout<<"Creating K matrix production SVP "<<SVPName<<" with channelIndex = "
- <<channelIndex<<" with s-dependent production Adler zero term"<<std::endl;
+ if (useProdAdler_) {
+ std::cout <<"Creating K matrix production SVP "<<SVPName<<" with channelIndex = "
+ <<channelIndex<<" with s-dependent production Adler zero term"<<std::endl;
} else {
- std::cout<<"Creating K matrix production SVP "<<SVPName<<" with channelIndex = "
- <<channelIndex<<" with production Adler zero factor = 1"<<std::endl;
+ std::cout <<"Creating K matrix production SVP "<<SVPName<<" with channelIndex = "
+ <<channelIndex<<" with production Adler zero factor = 1"<<std::endl;
}
-
}
LauKMatrixProdSVP::~LauKMatrixProdSVP()
@@ -58,14 +57,7 @@
// Destructor
}
-LauComplex LauKMatrixProdSVP::resAmp(Double_t mass, Double_t spinTerm)
-{
- std::cerr << "ERROR in LauKMatrixProdSVP::resAmp : This method shouldn't get called." << std::endl;
- std::cerr << " Returning zero amplitude for mass = " << mass << " and spinTerm = " << spinTerm << "." << std::endl;
- return LauComplex(0.0, 0.0);
-}
-
-LauComplex LauKMatrixProdSVP::amplitude(const LauKinematics* kinematics)
+LauComplex LauKMatrixProdSVP::resAmp(const Double_t mass, const Double_t spinTerm)
{
// Calculate the amplitude for the K-matrix production pole.
@@ -78,7 +70,38 @@
return amp;
}
- thePropagator_->updatePropagator(kinematics);
+ // Get barrier factors ('resonance' factor is already accounted for internally via propagator 'Gamma' matrix)
+ Double_t fFactorB(1.0);
+
+ const Int_t resSpin = this->getSpin();
+ const Double_t pstar = this->getPstar();
+
+ if ( resSpin > 0 ) {
+ const LauBlattWeisskopfFactor* parBWFactor = this->getParBWFactor();
+ if ( parBWFactor != nullptr ) {
+ switch ( parBWFactor->getRestFrame() ) {
+ case LauBlattWeisskopfFactor::ResonanceFrame:
+ fFactorB = parBWFactor->calcFormFactor(this->getP());
+ break;
+ case LauBlattWeisskopfFactor::ParentFrame:
+ fFactorB = parBWFactor->calcFormFactor(pstar);
+ break;
+ case LauBlattWeisskopfFactor::Covariant:
+ {
+ Double_t covFactor = this->getCovFactor();
+ if ( resSpin > 2 ) {
+ covFactor = TMath::Power( covFactor, 1.0/resSpin );
+ } else if ( resSpin == 2 ) {
+ covFactor = TMath::Sqrt( covFactor );
+ }
+ fFactorB = parBWFactor->calcFormFactor(pstar*covFactor);
+ break;
+ }
+ }
+ }
+ }
+
+ thePropagator_->updatePropagator(mass*mass);
Double_t SVPTerm = thePropagator_->getProdSVPTerm();
@@ -90,6 +113,14 @@
amp.rescale(SVPTerm*adlerZero);
+ // Scale by the spin term
+ Double_t scale = spinTerm;
+
+ // Include Blatt-Weisskopf barrier factor for parent
+ scale *= fFactorB;
+
+ amp.rescale(scale);
+
return amp;
}
@@ -105,9 +136,9 @@
{
LauParameter& par_f_ = thePropagator_->getScatteringParameter(channelIndex_, jChannel);
if ( !par_f_.fixed() )
- {
+ {
this->addFloatingParameter( &par_f_ );
- }
+ }
}
diff --git a/src/LauKMatrixPropagator.cc b/src/LauKMatrixPropagator.cc
--- a/src/LauKMatrixPropagator.cc
+++ b/src/LauKMatrixPropagator.cc
@@ -6,7 +6,7 @@
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
+ 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,
@@ -23,12 +23,13 @@
*/
/*! \file LauKMatrixPropagator.cc
- \brief File containing implementation of LauKMatrixPropagator class.
+ \brief File containing implementation of LauKMatrixPropagator class.
*/
#include "LauKMatrixPropagator.hh"
-#include "LauConstants.hh"
#include "LauTextFileParser.hh"
+#include "LauKinematics.hh"
+#include "LauComplex.hh"
#include "TMath.h"
#include "TSystem.h"
@@ -45,43 +46,22 @@
ClassImp(LauKMatrixPropagator)
LauKMatrixPropagator::LauKMatrixPropagator(const TString& name, const TString& paramFile,
- Int_t resPairAmpInt, Int_t nChannels,
- Int_t nPoles, Int_t rowIndex) :
+ Int_t resPairAmpInt, Int_t nChannels,
+ Int_t nPoles, Int_t rowIndex) :
name_(name),
paramFileName_(paramFile),
resPairAmpInt_(resPairAmpInt),
index_(rowIndex - 1),
- previousS_(0.0),
- scattSVP_(0.0),
- prodSVP_(0.0),
nChannels_(nChannels),
- nPoles_(nPoles),
- sAConst_(0.0),
- m2piSq_(4.0*LauConstants::mPiSq),
- m2KSq_( 4.0*LauConstants::mKSq),
- m2EtaSq_(4.0*LauConstants::mEtaSq),
- mEtaEtaPSumSq_((LauConstants::mEta + LauConstants::mEtaPrime)*(LauConstants::mEta + LauConstants::mEtaPrime)),
- mEtaEtaPDiffSq_((LauConstants::mEta - LauConstants::mEtaPrime)*(LauConstants::mEta - LauConstants::mEtaPrime)),
- mKpiSumSq_((LauConstants::mK + LauConstants::mPi)*(LauConstants::mK + LauConstants::mPi)),
- mKpiDiffSq_((LauConstants::mK - LauConstants::mPi)*(LauConstants::mK - LauConstants::mPi)),
- mKEtaPSumSq_((LauConstants::mK + LauConstants::mEtaPrime)*(LauConstants::mK + LauConstants::mEtaPrime)),
- mKEtaPDiffSq_((LauConstants::mK - LauConstants::mEtaPrime)*(LauConstants::mK - LauConstants::mEtaPrime)),
- mK3piDiffSq_((LauConstants::mK - 3.0*LauConstants::mPi)*(LauConstants::mK - 3.0*LauConstants::mPi)),
- k3piFactor_(TMath::Power((1.44 - mK3piDiffSq_)/1.44, -2.5)),
- fourPiFactor1_(16.0*LauConstants::mPiSq),
- fourPiFactor2_(TMath::Sqrt(1.0 - fourPiFactor1_)),
- adlerZeroFactor_(0.0),
- parametersSet_(kFALSE),
- verbose_(kFALSE),
- scattSymmetry_(kTRUE)
+ nPoles_(nPoles)
{
// Constructor
- // Check that the index is OK
- if (index_ < 0 || index_ >= nChannels_) {
- std::cerr << "ERROR in LauKMatrixPropagator constructor. The rowIndex, which is set to "
- << rowIndex << ", must be between 1 and the number of channels "
- << nChannels_ << std::endl;
+ // Check that the index is OK
+ if (index_ < 0 || index_ >= nChannels_) {
+ std::cerr << "ERROR in LauKMatrixPropagator constructor. The rowIndex, which is set to "
+ << rowIndex << ", must be between 1 and the number of channels "
+ << nChannels_ << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
@@ -93,16 +73,17 @@
// Destructor
realProp_.Clear();
negImagProp_.Clear();
- ScattKMatrix_.Clear();
- ReRhoMatrix_.Clear();
- ImRhoMatrix_.Clear();
+ ScattKMatrix_.Clear();
+ ReRhoMatrix_.Clear();
+ ImRhoMatrix_.Clear();
+ GammaMatrix_.Clear();
ReTMatrix_.Clear();
ImTMatrix_.Clear();
IMatrix_.Clear();
zeroMatrix_.Clear();
}
-LauComplex LauKMatrixPropagator::getPropTerm(Int_t channel) const
+LauComplex LauKMatrixPropagator::getPropTerm(const Int_t channel) const
{
// Get the (i,j) = (index_, channel) term of the propagator
// matrix. This allows us not to return the full propagator matrix.
@@ -114,7 +95,7 @@
}
-Double_t LauKMatrixPropagator::getRealPropTerm(Int_t channel) const
+Double_t LauKMatrixPropagator::getRealPropTerm(const Int_t channel) const
{
// Get the real part of the (i,j) = (index_, channel) term of the propagator
// matrix. This allows us not to return the full propagator matrix.
@@ -125,7 +106,7 @@
}
-Double_t LauKMatrixPropagator::getImagPropTerm(Int_t channel) const
+Double_t LauKMatrixPropagator::getImagPropTerm(const Int_t channel) const
{
// Get the imaginary part of the (i,j) = (index_, channel) term of the propagator
// matrix. This allows us not to return the full propagator matrix.
@@ -136,36 +117,9 @@
}
-void LauKMatrixPropagator::updatePropagator(const LauKinematics* kinematics)
-{
-
- // Calculate the K-matrix propagator for the given s value.
- // The K-matrix amplitude is given by
- // T_i = sum_{ij} (I - iK*rho)^-1 * P_j, where P is the production K-matrix.
- // i = index for the state (e.g. S-wave index = 0).
- // Here, we only find the (I - iK*rho)^-1 matrix part.
-
- if (!kinematics) {return;}
-
- // Get the invariant mass squared (s) from the kinematics object.
- // Use the resPairAmpInt to find which mass-squared combination to use.
- Double_t s(0.0);
-
- if (resPairAmpInt_ == 1) {
- s = kinematics->getm23Sq();
- } else if (resPairAmpInt_ == 2) {
- s = kinematics->getm13Sq();
- } else if (resPairAmpInt_ == 3) {
- s = kinematics->getm12Sq();
- }
-
- this->updatePropagator(s);
-
-}
-
-void LauKMatrixPropagator::updatePropagator(Double_t s)
+void LauKMatrixPropagator::updatePropagator(const Double_t s)
{
- // Calculate the K-matrix propagator for the given s value.
+ // Calculate the K-matrix propagator for the given s value.
// The K-matrix amplitude is given by
// T_i = sum_{ij} (I - iK*rho)^-1 * P_j, where P is the production K-matrix.
// i = index for the state (e.g. S-wave index = 0).
@@ -193,20 +147,29 @@
// if the quantity s is below various threshold values (analytic continuation).
this->calcRhoMatrix(s);
- // Calculate K*rho (real and imaginary parts, since rho can be complex)
- TMatrixD K_realRho(ScattKMatrix_);
- K_realRho *= ReRhoMatrix_;
- TMatrixD K_imagRho(ScattKMatrix_);
- K_imagRho *= ImRhoMatrix_;
+ // Calculate the angular momentum barrier matrix, which is real and diagonal
+ this->calcGammaMatrix(s);
+
+ // Calculate K*rho*(gamma^2) (real and imaginary parts, since rho can be complex)
+ TMatrixD GammaMatrixSq = (GammaMatrix_*GammaMatrix_);
+ TMatrixD K_realRhoGammaSq(ScattKMatrix_);
+
+ K_realRhoGammaSq *= ReRhoMatrix_;
+ K_realRhoGammaSq *= GammaMatrixSq;
+ TMatrixD K_imagRhoGammaSq(ScattKMatrix_);
+
+ K_imagRhoGammaSq *= ImRhoMatrix_;
+ K_imagRhoGammaSq *= GammaMatrixSq;
- // A = I + K*Imag(rho), B = -K*Real(Rho)
+
+ // A = I + K*Imag(rho)Gamma^2, B = -K*Real(Rho)Gamma^2
// Calculate C and D matrices such that (A + iB)*(C + iD) = I,
- // ie. C + iD = (I - i K*rho)^-1, separated into real and imaginary parts.
+ // ie. C + iD = (I - i K*rhoGamma^2)^-1, separated into real and imaginary parts.
// realProp C = (A + B A^-1 B)^-1, imagProp D = -A^-1 B C
TMatrixD A(IMatrix_);
- A += K_imagRho;
+ A += K_imagRhoGammaSq;
TMatrixD B(zeroMatrix_);
- B -= K_realRho;
+ B -= K_realRhoGammaSq;
TMatrixD invA(TMatrixD::kInverted, A);
TMatrixD invA_B(invA);
@@ -216,7 +179,6 @@
TMatrixD invC(A);
invC += B_invA_B;
-
// Set the real part of the propagator matrix ("C")
realProp_ = TMatrixD(TMatrixD::kInverted, invC);
@@ -226,6 +188,29 @@
negImagProp_ = TMatrixD(invA);
negImagProp_ *= BC;
+ // Pre-multiply by the Gamma matrix:
+ realProp_ = GammaMatrix_ * realProp_;
+ negImagProp_ = GammaMatrix_ * negImagProp_;
+
+ if(verbose_)
+ {
+ std::cout << "In LauKMatrixPropagator::updatePropagator(s). D[1-iKrhoD^2]^-1: " << std::endl;
+ TString realOutput("Real part:"), imagOutput("Imag part:");
+ for (int iChannel = 0; iChannel < nChannels_; iChannel++)
+ {
+ for (int jChannel = 0; jChannel < nChannels_; jChannel++)
+ {
+ realOutput += Form("\t%.6f",realProp_[iChannel][jChannel]);
+ imagOutput += Form("\t%.6f",-1*negImagProp_[iChannel][jChannel]);
+ }
+ realOutput += "\n ";
+ imagOutput += "\n ";
+ }
+ std::cout << realOutput << std::endl;
+ std::cout << imagOutput << std::endl;
+ }
+
+
// Also calculate the production SVP term, since this uses Adler-zero parameters
// defined in the parameter file.
this->updateProdSVPTerm(s);
@@ -250,7 +235,8 @@
cout<<"nChannels = "<<nChannels_<<", nPoles = "<<nPoles_<<endl;
// Initialise various matrices
- this->initialiseMatrices();
+ this->initialiseMatrices();
+ std::vector<Int_t> a(nChannels_,0);
// The format of the input file contains lines starting with a keyword followed by the
// appropriate set of parameters. Keywords are case insensitive (treated as lower-case).
@@ -260,8 +246,17 @@
// "Pole poleIndex mass g_1 g_2 ... g_N"
// 3) Definition of scattering f_{ij} constants: scattering index (1 to N), channel values
// "Scatt index f_{i1} f_{i2} ... f_{iN}", where i = index
- // 4) Various Adler zero and scattering constants; each line is "name value".
- // Possible names are mSq0, s0Scatt, s0Prod, sA, sA0
+ // 4) Orbital angular momentu for each channel. If not set here, defaults to 0
+ // "AngularMomentum L[0] L[1] ... L[N]"
+ // 5) Barrier factor parameter, which appears in the denominator and multiplies the term involving
+ // the nominal radius. If not set here, defaults to 0 or values appropriate to non-zero angular
+ // momenta as set in in (4) above.
+ // "BarrierFactorParameter a[0] a[1] ... a[N]"
+ // 6) Characteristic radius for each channel. If not set here, defaults to 3.0 GeV^{-1}
+ // "Radii radChannel1 radChannel2 ... radChannelN"
+ // 7) Various Adler zero and scattering constants; each line is "name value".
+ // Possible names are mSq0, s0Scatt, s0Prod, sA, sA0
+ //
// By default, the scattering constants are symmetrised: f_{ji} = f_{ij}.
// To not assume this use "ScattSymmetry 0" on a line
@@ -272,9 +267,9 @@
UInt_t nTotLines = readFile.getTotalNumLines();
if (nTotLines == 0) {
- std::cerr << "ERROR in LauKMatrixPropagator::setParameters : K-matrix parameter file not present - exiting." << std::endl;
+ std::cerr << "ERROR in LauKMatrixPropagator::setParameters : K-matrix parameter file not present - exiting." << std::endl;
- gSystem->Exit(EXIT_FAILURE);
+ gSystem->Exit(EXIT_FAILURE);
}
UInt_t iLine(0);
@@ -305,6 +300,21 @@
// Scattering terms
this->storeScattering(theLine);
+ } else if (!keyword.CompareTo("angularmomentum")) {
+
+ // Orbital angular momentum state for each channel & set default a values if called before storeBarrierFactorParameter
+ this->storeOrbitalAngularMomenta(theLine, a);
+
+ } else if (!keyword.CompareTo("barrierfactorparameter")) {
+
+ // Value of parameter "a" in denominator of centrifugal barrier factor, gamma
+ this->storeBarrierFactorParameter(theLine, a);
+
+ } else if (!keyword.CompareTo("radii")) {
+
+ // Values of characteristic radius
+ this->storeRadii(theLine);
+
} else {
// Usually Adler-zero constants
@@ -331,6 +341,11 @@
}
}
+ // Now that radii and barrier-factor-denominator parameters have been set, cache the value of "a/(R*R)"
+ for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) {
+ gamAInvRadSq_[iChannel] = a[iChannel]/(radii_[iChannel]*radii_[iChannel]);
+ }
+
// All required parameters have been set
parametersSet_ = kTRUE;
@@ -365,14 +380,30 @@
ReRhoMatrix_.ResizeTo(nChannels_, nChannels_);
ImRhoMatrix_.ResizeTo(nChannels_, nChannels_);
+ // Gamma matrices
+ GammaMatrix_.Clear();
+ GammaMatrix_.ResizeTo(nChannels_, nChannels_);
+
+ // Vector of orbital angular momenta for the channels (default is S-wave everywhere)
+ L_.clear();
+ L_.assign(nChannels_,0);
+
+ // Characteristic radius (diagonal) vector (default to 3.0)
+ radii_.clear();
+ radii_.assign(nChannels_,3.0);
+
+ // Vector to cache ratio a/R^2
+ gamAInvRadSq_.clear();
+ gamAInvRadSq_.resize(nChannels_);
+
// Square-root phase space matrices
ReSqrtRhoMatrix_.Clear(); ImSqrtRhoMatrix_.Clear();
ReSqrtRhoMatrix_.ResizeTo(nChannels_, nChannels_);
ImSqrtRhoMatrix_.ResizeTo(nChannels_, nChannels_);
// T matrices
- ReTMatrix_.Clear(); ImTMatrix_.Clear();
- ReTMatrix_.ResizeTo(nChannels_, nChannels_);
+ ReTMatrix_.Clear(); ImTMatrix_.Clear();
+ ReTMatrix_.ResizeTo(nChannels_, nChannels_);
ImTMatrix_.ResizeTo(nChannels_, nChannels_);
// For the coupling and scattering constants, use LauParArrays instead of TMatrices
@@ -401,9 +432,10 @@
mSqPoles_.clear();
mSqPoles_.resize(nPoles_);
+ haveCalled_storeBarrierFactorParameter = kFALSE;
}
-void LauKMatrixPropagator::storeChannels(const std::vector<std::string>& theLine)
+void LauKMatrixPropagator::storeChannels(const std::vector<std::string>& theLine)
{
// Get the list of channel indices to specify what phase space factors should be used
@@ -413,7 +445,7 @@
Int_t nTypes = static_cast<Int_t>(theLine.size()) - 1;
if (nTypes != nChannels_) {
cerr<<"Error in LauKMatrixPropagator::storeChannels. The input file defines "
- <<nTypes<<" channels when "<<nChannels_<<" are expected"<<endl;
+ <<nTypes<<" channels when "<<nChannels_<<" are expected"<<endl;
return;
}
@@ -424,11 +456,11 @@
if (checkChannel == kTRUE) {
cout<<"Adding phase space channel index "<<phaseSpaceInt
- <<" to K-matrix propagator "<<name_<<endl;
- phaseSpaceTypes_[iChannel] = phaseSpaceInt;
+ <<" to K-matrix propagator "<<name_<<endl;
+ phaseSpaceTypes_[iChannel] = static_cast<LauKMatrixPropagator::KMatrixChannels>(phaseSpaceInt);
} else {
cerr<<"Phase space channel index "<<phaseSpaceInt
- <<" should be between 1 and "<<LauKMatrixPropagator::TotChannels-1<<endl;
+ <<" should be between 1 and "<<static_cast<int>(LauKMatrixPropagator::KMatrixChannels::TotChannels)-1<<endl;
}
}
@@ -464,7 +496,7 @@
gCouplings_[poleIndex][iChannel] = couplingParam;
cout<<"Added coupling parameter g^{"<<poleIndex+1<<"}_"
- <<iChannel+1<<" = "<<couplingConst<<endl;
+ <<iChannel+1<<" = "<<couplingConst<<endl;
}
@@ -478,7 +510,6 @@
}
-
}
void LauKMatrixPropagator::storeScattering(const std::vector<std::string>& theLine)
@@ -499,11 +530,11 @@
for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) {
Double_t scattConst = std::atof(theLine[iChannel+2].c_str());
- LauParameter scattParam(Form("KM_%s_scattConst_%i_%i",name_.Data(),scattIndex,iChannel),scattConst);
+ LauParameter scattParam(Form("KM_%s_fScatteringConst_%i_%i",name_.Data(),scattIndex,iChannel),scattConst);
fScattering_[scattIndex][iChannel] = scattParam;
cout<<"Added scattering parameter f("<<scattIndex+1<<","
- <<iChannel+1<<") = "<<scattConst<<endl;
+ <<iChannel+1<<") = "<<scattConst<<endl;
}
@@ -512,12 +543,132 @@
} else {
cerr<<"Error in LauKMatrixPropagator::storeScattering. Expecting "<<nExpect
- <<" numbers for scattering constants definition but found "<<nWords
- <<" values instead"<<endl;
+ <<" numbers for scattering constants definition but found "<<nWords
+ <<" values instead"<<endl;
+
+ }
+
+}
+
+void LauKMatrixPropagator::storeOrbitalAngularMomenta(const std::vector<std::string>& theLine, std::vector<Int_t>& a)
+{
+
+ // Store the orbital angular momentum for each channel
+ // Each line will contain: angularmomentum OrbitalAngularMomentumPerChannel
+
+ // Check that the line has nChannels_ + 1 strings
+ Int_t nWords = static_cast<Int_t>(theLine.size());
+ Int_t nExpect = nChannels_ + 1;
+
+ if (nWords == nExpect) {
+
+ for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) {
+
+ Int_t angularMomentum = std::atoi(theLine[iChannel+1].c_str());
+ L_[iChannel] = angularMomentum;
+
+ cout<<"Defined K-matrix orbital angular momentum "<<angularMomentum<<" for channel "
+ <<iChannel<<endl;
+
+ }
+
+ } else {
+
+ cerr<<"Error in LauKMatrixPropagator::storeOrbitalAngularMomenta. Expecting "<<nExpect
+ <<" numbers for orbital angular momenta definition but found "<<nWords
+ <<" values instead"<<endl;
+
+ }
+
+ if (!haveCalled_storeBarrierFactorParameter)
+ {
+ // Set default value of spin-dependent centrifugal-barrier-factor parameter
+ for( Int_t iCh = 0; iCh < nChannels_; iCh++ )
+ {
+ switch(L_[iCh]) {
+ case 0:
+ a[iCh] = 0;
+ break;
+ case 1:
+ a[iCh] = 1;
+ break;
+ case 2:
+ a[iCh] = 3;
+ break;
+ default:
+ std::cerr << "ERROR in LauKMatrixPropagator constructor. Centrifugal barrier factor and angular-momentum terms of K-matrix are only defined for S-, P-, or D-wave."
+ << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+ }
+ }
+
+}
+
+void LauKMatrixPropagator::storeRadii(const std::vector<std::string>& theLine)
+{
+
+ // Store the characteristic radii (measured in GeV^{-1})
+ // Each line will contain: Radii RadiusConstantsPerChannel
+
+ // Check that the line has nChannels_ + 1 strings
+ Int_t nWords = static_cast<Int_t>(theLine.size());
+ Int_t nExpect = nChannels_ + 1;
+
+ if (nWords == nExpect) {
+
+ for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) {
+
+ Double_t radiusConst = std::atof(theLine[iChannel+1].c_str());
+ radii_[iChannel] = radiusConst;
+
+ cout<<"Added K-matrix radius "<<radiusConst<<" for channel "
+ <<iChannel<<endl;
+
+ }
+
+ } else {
+
+ cerr<<"Error in LauKMatrixPropagator::storeRadii. Expecting "<<nExpect
+ <<" numbers for radii definition but found "<<nWords
+ <<" values instead"<<endl;
}
+}
+
+void LauKMatrixPropagator::storeBarrierFactorParameter(const std::vector<std::string>& theLine, std::vector<Int_t>& a)
+{
+
+ // Store the parameter of the barrier factor
+ // Each line will contain: barrierfactorparameter ParameterValuePerchannel
+ // Check that the line has nChannels_ + 1 strings
+ Int_t nWords = static_cast<Int_t>(theLine.size());
+ Int_t nExpect = nChannels_ + 1;
+
+ if (nWords == nExpect) {
+
+ for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) {
+
+ Double_t parameterValue = std::atof(theLine[iChannel+1].c_str());
+ a[iChannel] = parameterValue;
+
+ cout<<"Added K-matrix barrier factor parameter value "<<parameterValue<<" for channel "
+ <<iChannel<<endl;
+
+ }
+
+ // Set flag to stop storeOrbitalAngularMomenta overriding these a values
+ haveCalled_storeBarrierFactorParameter = kTRUE;
+
+ } else {
+
+ cerr<<"Error in LauKMatrixPropagator::storeBarrierFactorParameter. Expecting "<<nExpect
+ <<" numbers for barrier factor parameter definition but found "<<nWords
+ <<" values instead"<<endl;
+
+ }
}
void LauKMatrixPropagator::storeParameter(const TString& keyword, const TString& parString)
@@ -565,8 +716,7 @@
}
-
-void LauKMatrixPropagator::calcScattKMatrix(Double_t s)
+void LauKMatrixPropagator::calcScattKMatrix(const Double_t s)
{
// Calculate the scattering K-matrix for the given value of s.
@@ -622,7 +772,7 @@
}
-void LauKMatrixPropagator::calcPoleDenomVect(Double_t s)
+void LauKMatrixPropagator::calcPoleDenomVect(const Double_t s)
{
// Calculate the 1/(m_pole^2 - s) terms for the scattering
// and production K-matrix formulae.
@@ -637,10 +787,9 @@
poleDenomVect_.push_back(invPoleTerm);
}
-
}
-Double_t LauKMatrixPropagator::getPoleDenomTerm(Int_t poleIndex) const
+Double_t LauKMatrixPropagator::getPoleDenomTerm(const Int_t poleIndex) const
{
if (parametersSet_ == kFALSE) {return 0.0;}
@@ -651,7 +800,7 @@
}
-LauParameter& LauKMatrixPropagator::getPoleMassSqParameter(Int_t poleIndex)
+LauParameter& LauKMatrixPropagator::getPoleMassSqParameter(const Int_t poleIndex)
{
if ( (parametersSet_ == kFALSE) || (poleIndex < 0 || poleIndex >= nPoles_) ) {
std::cerr << "ERROR from LauKMatrixPropagator::getPoleMassSqParameter(). Invalid pole." << std::endl;
@@ -661,7 +810,7 @@
return mSqPoles_[poleIndex];
}
-Double_t LauKMatrixPropagator::getCouplingConstant(Int_t poleIndex, Int_t channelIndex) const
+Double_t LauKMatrixPropagator::getCouplingConstant(const Int_t poleIndex, const Int_t channelIndex) const
{
if (parametersSet_ == kFALSE) {return 0.0;}
@@ -673,7 +822,7 @@
}
-LauParameter& LauKMatrixPropagator::getCouplingParameter(Int_t poleIndex, Int_t channelIndex)
+LauParameter& LauKMatrixPropagator::getCouplingParameter(const Int_t poleIndex, const Int_t channelIndex)
{
if ( (parametersSet_ == kFALSE) || (poleIndex < 0 || poleIndex >= nPoles_) || (channelIndex < 0 || channelIndex >= nChannels_) ) {
@@ -686,7 +835,7 @@
return gCouplings_[poleIndex][channelIndex];
}
-Double_t LauKMatrixPropagator::getScatteringConstant(Int_t channel1Index, Int_t channel2Index) const
+Double_t LauKMatrixPropagator::getScatteringConstant(const Int_t channel1Index, const Int_t channel2Index) const
{
if (parametersSet_ == kFALSE) {return 0.0;}
@@ -698,7 +847,7 @@
}
-LauParameter& LauKMatrixPropagator::getScatteringParameter(Int_t channel1Index, Int_t channel2Index)
+LauParameter& LauKMatrixPropagator::getScatteringParameter(const Int_t channel1Index, const Int_t channel2Index)
{
if ( (parametersSet_ == kFALSE) || (channel1Index < 0 || channel1Index >= nChannels_) || (channel2Index < 0 || channel2Index >= nChannels_) ) {
@@ -709,7 +858,7 @@
return fScattering_[channel1Index][channel2Index];
}
-Double_t LauKMatrixPropagator::calcSVPTerm(Double_t s, Double_t s0) const
+Double_t LauKMatrixPropagator::calcSVPTerm(const Double_t s, const Double_t s0) const
{
if (parametersSet_ == kFALSE) {return 0.0;}
@@ -719,27 +868,27 @@
Double_t deltaS = s - s0;
if (TMath::Abs(deltaS) > 1.0e-6) {
result = (mSq0_.unblindValue() - s0)/deltaS;
- }
+ }
return result;
}
-void LauKMatrixPropagator::updateScattSVPTerm(Double_t s)
+void LauKMatrixPropagator::updateScattSVPTerm(const Double_t s)
{
// Update the scattering "slowly-varying part" (SVP)
Double_t s0Scatt = s0Scatt_.unblindValue();
scattSVP_ = this->calcSVPTerm(s, s0Scatt);
}
-void LauKMatrixPropagator::updateProdSVPTerm(Double_t s)
+void LauKMatrixPropagator::updateProdSVPTerm(const Double_t s)
{
// Update the production "slowly-varying part" (SVP)
Double_t s0Prod = s0Prod_.unblindValue();
prodSVP_ = this->calcSVPTerm(s, s0Prod);
}
-void LauKMatrixPropagator::updateAdlerZeroFactor(Double_t s)
+void LauKMatrixPropagator::updateAdlerZeroFactor(const Double_t s)
{
// Calculate the multiplicative factor containing various Adler zero
@@ -750,11 +899,64 @@
Double_t deltaS = s - sA0Val;
if (TMath::Abs(deltaS) > 1e-6) {
adlerZeroFactor_ = (s - sAConst_)*(1.0 - sA0Val)/deltaS;
- }
+ }
+
+}
+
+void LauKMatrixPropagator::calcGammaMatrix(const Double_t s)
+{
+ // Calculate the gamma angular momentum barrier matrix
+ // for the given invariant mass squared quantity, s.
+
+ // Initialise all entries to zero
+ GammaMatrix_.Zero();
+
+ Double_t gamma(0.0);
+
+ for (Int_t iChannel (0); iChannel < nChannels_; ++iChannel) {
+
+ LauKMatrixPropagator::KMatrixChannels phaseSpaceIndex = phaseSpaceTypes_[iChannel];
+
+ if ( L_[iChannel] != 0 ) {
+ gamma = this->calcGamma(iChannel,s,phaseSpaceIndex);
+ } else {
+ gamma = 1.0; // S-wave
+ }
+
+ if (verbose_) {
+ cout<<"GammaMatrix("<<iChannel<<", "<<iChannel<<") = "<<gamma<<endl;
+ }
+
+ GammaMatrix_(iChannel, iChannel) = gamma;
+
+ }
}
-void LauKMatrixPropagator::calcRhoMatrix(Double_t s)
+Double_t LauKMatrixPropagator::calcGamma(const Int_t iCh, const Double_t s, const LauKMatrixPropagator::KMatrixChannels phaseSpaceIndex) const
+{
+ // Calculate the DK P-wave angular momentum barrier factor
+ Double_t gamma(0.0);
+
+ LauComplex rho = getRho(s,phaseSpaceIndex);
+ Double_t q = 0.5 * sqrt(s) * mag(rho);
+
+ gamma = pow(q,L_[iCh]);
+ if (includeBWBarrierFactor_)
+ {
+ gamma /= pow( q*q + gamAInvRadSq_[iCh] , L_[iCh]/2. );
+ }
+
+ if(verbose_)
+ {
+ std::cout << "In LauKMatrixPropagator::calcGamma(iCh=" << iCh << ", s=" << s << ", prop). ";
+ std::cout << "|q(iCh="<<iCh<<")|: " << q << std::endl;
+ }
+
+ return gamma;
+}
+
+void LauKMatrixPropagator::calcRhoMatrix(const Double_t s)
{
// Calculate the real and imaginary part of the phase space density
@@ -762,47 +964,107 @@
// The matrix can be complex if s is below threshold (so that
// the amplitude continues analytically).
- // Initialise all entries to zero
- ReRhoMatrix_.Zero(); ImRhoMatrix_.Zero();
-
- LauComplex rho(0.0, 0.0);
- Int_t phaseSpaceIndex(0);
+ // Initialise all entries to zero
+ ReRhoMatrix_.Zero(); ImRhoMatrix_.Zero();
for (Int_t iChannel (0); iChannel < nChannels_; ++iChannel) {
- phaseSpaceIndex = phaseSpaceTypes_[iChannel];
+ LauKMatrixPropagator::KMatrixChannels phaseSpaceIndex = phaseSpaceTypes_[iChannel];
+
+ LauComplex rho = getRho(s, phaseSpaceIndex)
+
+ if (verbose_) {
+ cout<<"ReRhoMatrix("<<iChannel<<", "<<iChannel<<") = "<<rho.re()<<endl;
+ cout<<"ImRhoMatrix("<<iChannel<<", "<<iChannel<<") = "<<rho.im()<<endl;
+ }
- if (phaseSpaceIndex == LauKMatrixPropagator::PiPi) {
+ ReRhoMatrix_(iChannel, iChannel) = rho.re();
+ ImRhoMatrix_(iChannel, iChannel) = rho.im();
+
+ }
+
+}
+
+LauComplex LauKMatrixPropagator::getRho(const Double_t s, const:LauKMatrixPropagator::KMatrixChannels) const
+{
+ LauComplex rho(0.0, 0.0);
+ switch (phaseSpaceIndex)
+ {
+ case LauKMatrixPropagator::KMatrixChannels::PiPi :
rho = this->calcPiPiRho(s);
- } else if (phaseSpaceIndex == LauKMatrixPropagator::KK) {
+ break;
+ case LauKMatrixPropagator::KMatrixChannels::KK :
rho = this->calcKKRho(s);
- } else if (phaseSpaceIndex == LauKMatrixPropagator::FourPi) {
+ break;
+ case LauKMatrixPropagator::KMatrixChannels::FourPi :
rho = this->calcFourPiRho(s);
- } else if (phaseSpaceIndex == LauKMatrixPropagator::EtaEta) {
+ break;
+ case LauKMatrixPropagator::KMatrixChannels::EtaEta :
rho = this->calcEtaEtaRho(s);
- } else if (phaseSpaceIndex == LauKMatrixPropagator::EtaEtaP) {
+ break;
+ case LauKMatrixPropagator::KMatrixChannels::EtaEtaP :
rho = this->calcEtaEtaPRho(s);
- } else if (phaseSpaceIndex == LauKMatrixPropagator::KPi) {
+ break;
+ case LauKMatrixPropagator::KMatrixChannels::KPi :
rho = this->calcKPiRho(s);
- } else if (phaseSpaceIndex == LauKMatrixPropagator::KEtaP) {
+ break;
+ case LauKMatrixPropagator::KMatrixChannels::KEtaP :
rho = this->calcKEtaPRho(s);
- } else if (phaseSpaceIndex == LauKMatrixPropagator::KThreePi) {
+ break;
+ case LauKMatrixPropagator::KMatrixChannels::KThreePi :
rho = this->calcKThreePiRho(s);
- }
+ break;
+ case LauKMatrixPropagator::KMatrixChannels::D0K :
+ rho = this->calcD0KRho(s);
+ break;
+ case LauKMatrixPropagator::KMatrixChannels::Dstar0K :
+ rho = this->calcDstar0KRho(s);
+ break;
+ default :
+ std::cerr << "ERROR in LauKMatrixPropagator::getRho(...). Phase-space index not recognised for this channel"
+ << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+ return rho;
+}
- if (verbose_) {
- cout<<"ReRhoMatrix("<<iChannel<<", "<<iChannel<<") = "<<rho.re()<<endl;
- cout<<"ImRhoMatrix("<<iChannel<<", "<<iChannel<<") = "<<rho.im()<<endl;
- }
+LauComplex LauKMatrixPropagator::calcD0KRho(const Double_t s) const
+{
+ // Calculate the D0K+ phase space factor
+ LauComplex rho(0.0, 0.0);
+ if (TMath::Abs(s) < 1e-10) {return rho;}
- ReRhoMatrix_(iChannel, iChannel) = rho.re();
- ImRhoMatrix_(iChannel, iChannel) = rho.im();
+ Double_t sqrtTerm1 = (-mD0KSumSq_/s) + 1.0;
+ Double_t sqrtTerm2 = (-mD0KDiffSq_/s) + 1.0;
+ Double_t sqrtTerm = sqrtTerm1*sqrtTerm2;
+ if (sqrtTerm < 0.0) {
+ rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
+ } else {
+ rho.setRealPart( TMath::Sqrt(sqrtTerm) );
+ }
+ return rho;
+}
+
+LauComplex LauKMatrixPropagator::calcDstar0KRho(const Double_t s) const
+{
+ // Calculate the Dstar0K+ phase space factor
+ LauComplex rho(0.0, 0.0);
+ if (TMath::Abs(s) < 1e-10) {return rho;}
+
+ Double_t sqrtTerm1 = (-mDstar0KSumSq_/s) + 1.0;
+ Double_t sqrtTerm2 = (-mDstar0KDiffSq_/s) + 1.0;
+ Double_t sqrtTerm = sqrtTerm1*sqrtTerm2;
+ if (sqrtTerm < 0.0) {
+ rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
+ } else {
+ rho.setRealPart( TMath::Sqrt(sqrtTerm) );
}
+ return rho;
}
-LauComplex LauKMatrixPropagator::calcPiPiRho(Double_t s) const
+LauComplex LauKMatrixPropagator::calcPiPiRho(const Double_t s) const
{
// Calculate the pipi phase space factor
LauComplex rho(0.0, 0.0);
@@ -818,7 +1080,7 @@
return rho;
}
-LauComplex LauKMatrixPropagator::calcKKRho(Double_t s) const
+LauComplex LauKMatrixPropagator::calcKKRho(const Double_t s) const
{
// Calculate the KK phase space factor
LauComplex rho(0.0, 0.0);
@@ -834,33 +1096,33 @@
return rho;
}
-LauComplex LauKMatrixPropagator::calcFourPiRho(Double_t s) const
+LauComplex LauKMatrixPropagator::calcFourPiRho(const Double_t s) const
{
// Calculate the 4pi phase space factor. This uses a 6th-order polynomial
- // parameterisation that approximates the multi-body phase space double integral
- // defined in Eq 4 of the A&S paper hep-ph/0204328. This form agrees with the
- // BaBar model (another 6th order polynomial from s^4 down to 1/s^2), but avoids the
- // exponential increase at small values of s (~< 0.1) arising from 1/s and 1/s^2.
- // Eq 4 is evaluated for each value of s by assuming incremental steps of 1e-3 GeV^2
- // for s1 and s2, the invariant energy squared of each of the di-pion states,
- // with the integration limits of s1 = (2*mpi)^2 to (sqrt(s) - 2*mpi)^2 and
- // s2 = (2*mpi)^2 to (sqrt(s) - sqrt(s1))^2. The mass M of the rho is taken to be
- // 0.775 GeV and the energy-dependent width of the 4pi system
- // Gamma(s) = gamma_0*rho1^3(s), where rho1 = sqrt(1.0 - 4*mpiSq/s) and gamma_0 is
- // the "width" of the 4pi state at s = 1, which is taken to be 0.3 GeV
- // (~75% of the total width from PDG estimates of the f0(1370) -> 4pi state).
- // The normalisation term rho_0 is found by ensuring that the phase space integral
- // at s = 1 is equal to sqrt(1.0 - 16*mpiSq/s). Note that the exponent for this
- // factor in hep-ph/0204328 is wrong; it should be 0.5, i.e. sqrt, not n = 1 to 5.
- // Plotting the value of this double integral as a function of s can then be fitted
- // to a 6th-order polynomial (for s < 1), which is the result used below
+ // parameterisation that approximates the multi-body phase space double integral
+ // defined in Eq 4 of the A&S paper hep-ph/0204328. This form agrees with the
+ // BaBar model (another 6th order polynomial from s^4 down to 1/s^2), but avoids the
+ // exponential increase at small values of s (~< 0.1) arising from 1/s and 1/s^2.
+ // Eq 4 is evaluated for each value of s by assuming incremental steps of 1e-3 GeV^2
+ // for s1 and s2, the invariant energy squared of each of the di-pion states,
+ // with the integration limits of s1 = (2*mpi)^2 to (sqrt(s) - 2*mpi)^2 and
+ // s2 = (2*mpi)^2 to (sqrt(s) - sqrt(s1))^2. The mass M of the rho is taken to be
+ // 0.775 GeV and the energy-dependent width of the 4pi system
+ // Gamma(s) = gamma_0*rho1^3(s), where rho1 = sqrt(1.0 - 4*mpiSq/s) and gamma_0 is
+ // the "width" of the 4pi state at s = 1, which is taken to be 0.3 GeV
+ // (~75% of the total width from PDG estimates of the f0(1370) -> 4pi state).
+ // The normalisation term rho_0 is found by ensuring that the phase space integral
+ // at s = 1 is equal to sqrt(1.0 - 16*mpiSq/s). Note that the exponent for this
+ // factor in hep-ph/0204328 is wrong; it should be 0.5, i.e. sqrt, not n = 1 to 5.
+ // Plotting the value of this double integral as a function of s can then be fitted
+ // to a 6th-order polynomial (for s < 1), which is the result used below
LauComplex rho(0.0, 0.0);
if (TMath::Abs(s) < 1e-10) {return rho;}
if (s <= 1.0) {
- Double_t rhoTerm = ((1.07885*s + 0.13655)*s - 0.29744)*s - 0.20840;
- rhoTerm = ((rhoTerm*s + 0.13851)*s - 0.01933)*s + 0.00051;
+ Double_t rhoTerm = ((1.07885*s + 0.13655)*s - 0.29744)*s - 0.20840;
+ rhoTerm = ((rhoTerm*s + 0.13851)*s - 0.01933)*s + 0.00051;
// For some values of s (below 2*mpi), this term is a very small
// negative number. Check for this and set the rho term to zero.
if (rhoTerm < 0.0) {rhoTerm = 0.0;}
@@ -872,7 +1134,7 @@
return rho;
}
-LauComplex LauKMatrixPropagator::calcEtaEtaRho(Double_t s) const
+LauComplex LauKMatrixPropagator::calcEtaEtaRho(const Double_t s) const
{
// Calculate the eta-eta phase space factor
LauComplex rho(0.0, 0.0);
@@ -888,18 +1150,18 @@
return rho;
}
-LauComplex LauKMatrixPropagator::calcEtaEtaPRho(Double_t s) const
+LauComplex LauKMatrixPropagator::calcEtaEtaPRho(const Double_t s) const
{
// Calculate the eta-eta' phase space factor. Note that the
- // mass difference term m_eta - m_eta' is not included,
- // since this corresponds to a "t or u-channel crossing",
- // which means that we cannot simply analytically continue
- // this part of the phase space factor below threshold, which
- // we can do for s-channel contributions. This is actually an
- // unsolved problem, e.g. see Guo et al 1409.8652, and
- // Danilkin et al 1409.7708. Anisovich and Sarantsev in
- // hep-ph/0204328 "solve" this issue by setting the mass
- // difference term to unity, which is what we do here...
+ // mass difference term m_eta - m_eta' is not included,
+ // since this corresponds to a "t or u-channel crossing",
+ // which means that we cannot simply analytically continue
+ // this part of the phase space factor below threshold, which
+ // we can do for s-channel contributions. This is actually an
+ // unsolved problem, e.g. see Guo et al 1409.8652, and
+ // Danilkin et al 1409.7708. Anisovich and Sarantsev in
+ // hep-ph/0204328 "solve" this issue by setting the mass
+ // difference term to unity, which is what we do here...
LauComplex rho(0.0, 0.0);
if (TMath::Abs(s) < 1e-10) {return rho;}
@@ -908,14 +1170,14 @@
if (sqrtTerm < 0.0) {
rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
} else {
- rho.setRealPart( TMath::Sqrt(sqrtTerm) );
+ rho.setRealPart( TMath::Sqrt(sqrtTerm) );
}
return rho;
}
-LauComplex LauKMatrixPropagator::calcKPiRho(Double_t s) const
+LauComplex LauKMatrixPropagator::calcKPiRho(const Double_t s) const
{
// Calculate the K-pi phase space factor
LauComplex rho(0.0, 0.0);
@@ -933,7 +1195,7 @@
return rho;
}
-LauComplex LauKMatrixPropagator::calcKEtaPRho(Double_t s) const
+LauComplex LauKMatrixPropagator::calcKEtaPRho(const Double_t s) const
{
// Calculate the K-eta' phase space factor
LauComplex rho(0.0, 0.0);
@@ -951,7 +1213,7 @@
return rho;
}
-LauComplex LauKMatrixPropagator::calcKThreePiRho(Double_t s) const
+LauComplex LauKMatrixPropagator::calcKThreePiRho(const Double_t s) const
{
// Calculate the Kpipipi + multimeson phase space factor.
// Use the simplest definition in hep-ph/9705401 (Eq 14), which is the form
@@ -976,51 +1238,49 @@
return rho;
}
-Bool_t LauKMatrixPropagator::checkPhaseSpaceType(Int_t phaseSpaceInt) const
+Bool_t LauKMatrixPropagator::checkPhaseSpaceType(const Int_t phaseSpaceInt) const
{
Bool_t passed(kFALSE);
- if (phaseSpaceInt >= 1 && phaseSpaceInt < LauKMatrixPropagator::TotChannels) {
+ if (phaseSpaceInt >= 1 && phaseSpaceInt < static_cast<Int_t>(LauKMatrixPropagator::KMatrixChannels::TotChannels)) {
passed = kTRUE;
}
return passed;
}
-LauComplex LauKMatrixPropagator::getTransitionAmp(Double_t s, Int_t channel)
+LauComplex LauKMatrixPropagator::getTransitionAmp(const Double_t s, const Int_t channel)
{
- // Get the complex (unitary) transition amplitude T for the given channel
+ // Get the complex (unitary) transition amplitude T for the given channel
LauComplex TAmp(0.0, 0.0);
- channel -= 1;
- if (channel < 0 || channel >= nChannels_) {return TAmp;}
+ if (channel <= 0 || channel > nChannels_) {return TAmp;}
- this->getTMatrix(s);
+ this->getTMatrix(s);
- TAmp.setRealPart(ReTMatrix_[index_][channel]);
- TAmp.setImagPart(ImTMatrix_[index_][channel]);
+ TAmp.setRealPart(ReTMatrix_[index_][channel-1]);
+ TAmp.setImagPart(ImTMatrix_[index_][channel-1]);
return TAmp;
}
-LauComplex LauKMatrixPropagator::getPhaseSpaceTerm(Double_t s, Int_t channel)
+LauComplex LauKMatrixPropagator::getPhaseSpaceTerm(const Double_t s, const Int_t channel)
{
- // Get the complex (unitary) transition amplitude T for the given channel
+ // Get the complex (unitary) transition amplitude T for the given channel
LauComplex rho(0.0, 0.0);
- channel -= 1;
- if (channel < 0 || channel >= nChannels_) {return rho;}
+ if (channel <= 0 || channel > nChannels_) {return rho;}
// If s has changed from the previous value, recalculate rho
if (TMath::Abs(s - previousS_) > 1e-6*s) {
- this->calcRhoMatrix(s);
+ this->calcRhoMatrix(s);
}
- rho.setRealPart(ReRhoMatrix_[channel][channel]);
- rho.setImagPart(ImRhoMatrix_[channel][channel]);
+ rho.setRealPart(ReRhoMatrix_[channel][channel-1]);
+ rho.setImagPart(ImRhoMatrix_[channel][channel-1]);
return rho;
@@ -1028,11 +1288,11 @@
void LauKMatrixPropagator::getTMatrix(const LauKinematics* kinematics) {
- // Find the unitary T matrix, where T = [sqrt(rho)]^{*} T_hat sqrt(rho),
- // and T_hat = (I - i K rho)^-1 * K is the Lorentz-invariant T matrix,
- // which has phase-space factors included (rho). This function is not
- // needed to calculate the K-matrix amplitudes, but allows us
- // to check the variation of T as a function of s (kinematics)
+ // Find the unitary T matrix, where T = [sqrt(rho)]^{*} T_hat sqrt(rho),
+ // and T_hat = (I - i K rho)^-1 * K is the Lorentz-invariant T matrix,
+ // which has phase-space factors included (rho). This function is not
+ // needed to calculate the K-matrix amplitudes, but allows us
+ // to check the variation of T as a function of s (kinematics)
if (!kinematics) {return;}
@@ -1053,25 +1313,25 @@
}
-void LauKMatrixPropagator::getTMatrix(Double_t s)
+void LauKMatrixPropagator::getTMatrix(const Double_t s)
{
- // Find the unitary transition T matrix, where
- // T = [sqrt(rho)]^{*} T_hat sqrt(rho), and
- // T_hat = (I - i K rho)^-1 * K is the Lorentz-invariant T matrix,
- // which has phase-space factors included (rho). Note that the first
- // sqrt of the rho matrix is complex conjugated.
+ // Find the unitary transition T matrix, where
+ // T = [sqrt(rho)]^{*} T_hat sqrt(rho), and
+ // T_hat = (I - i K rho)^-1 * K is the Lorentz-invariant T matrix,
+ // which has phase-space factors included (rho). Note that the first
+ // sqrt of the rho matrix is complex conjugated.
- // This function is not needed to calculate the K-matrix amplitudes, but
- // allows us to check the variation of T as a function of s (kinematics)
+ // This function is not needed to calculate the K-matrix amplitudes, but
+ // allows us to check the variation of T as a function of s (kinematics)
- // Initialse the real and imaginary parts of the T matrix to zero
- ReTMatrix_.Zero(); ImTMatrix_.Zero();
+ // Initialse the real and imaginary parts of the T matrix to zero
+ ReTMatrix_.Zero(); ImTMatrix_.Zero();
- if (parametersSet_ == kFALSE) {return;}
+ if (parametersSet_ == kFALSE) {return;}
- // Update K, rho and the propagator (I - i K rho)^-1
- this->updatePropagator(s);
+ // Update K, rho and the propagator (I - i K rho)^-1
+ this->updatePropagator(s);
// Find the real and imaginary T_hat matrices
TMatrixD THatReal = realProp_*ScattKMatrix_;
@@ -1108,46 +1368,45 @@
void LauKMatrixPropagator::getSqrtRhoMatrix()
{
- // Find the square root of the (current) phase space matrix so that
- // we can find T = [sqrt(rho)}^{*} T_hat sqrt(rho), where T_hat is the
- // Lorentz-invariant T matrix = (I - i K rho)^-1 * K; note that the first
- // sqrt of rho matrix is complex conjugated
+ // Find the square root of the (current) phase space matrix so that
+ // we can find T = [sqrt(rho)}^{*} T_hat sqrt(rho), where T_hat is the
+ // Lorentz-invariant T matrix = (I - i K rho)^-1 * K; note that the first
+ // sqrt of rho matrix is complex conjugated
- // If rho = rho_i + i rho_r = a + i b, then sqrt(rho) = c + i d, where
- // c = sqrt(0.5*(r+a)) and d = sqrt(0.5(r-a)), where r = sqrt(a^2 + b^2).
- // Since rho is diagonal, then the square root of rho will also be diagonal,
- // with its real and imaginary matrix elements equal to c and d, respectively
+ // If rho = rho_i + i rho_r = a + i b, then sqrt(rho) = c + i d, where
+ // c = sqrt(0.5*(r+a)) and d = sqrt(0.5(r-a)), where r = sqrt(a^2 + b^2).
+ // Since rho is diagonal, then the square root of rho will also be diagonal,
+ // with its real and imaginary matrix elements equal to c and d, respectively
- // Initialise the real and imaginary parts of the square root of
- // the rho matrix to zero
- ReSqrtRhoMatrix_.Zero(); ImSqrtRhoMatrix_.Zero();
+ // Initialise the real and imaginary parts of the square root of
+ // the rho matrix to zero
+ ReSqrtRhoMatrix_.Zero(); ImSqrtRhoMatrix_.Zero();
- for (Int_t iChannel (0); iChannel < nChannels_; ++iChannel) {
+ for (Int_t iChannel (0); iChannel < nChannels_; ++iChannel) {
- Double_t realRho = ReRhoMatrix_[iChannel][iChannel];
- Double_t imagRho = ImRhoMatrix_[iChannel][iChannel];
+ Double_t realRho = ReRhoMatrix_[iChannel][iChannel];
+ Double_t imagRho = ImRhoMatrix_[iChannel][iChannel];
Double_t rhoMag = sqrt(realRho*realRho + imagRho*imagRho);
- Double_t rhoSum = rhoMag + realRho;
- Double_t rhoDiff = rhoMag - realRho;
+ Double_t rhoSum = rhoMag + realRho;
+ Double_t rhoDiff = rhoMag - realRho;
- Double_t reSqrtRho(0.0), imSqrtRho(0.0);
- if (rhoSum > 0.0) {reSqrtRho = sqrt(0.5*rhoSum);}
+ Double_t reSqrtRho(0.0), imSqrtRho(0.0);
+ if (rhoSum > 0.0) {reSqrtRho = sqrt(0.5*rhoSum);}
if (rhoDiff > 0.0) {imSqrtRho = sqrt(0.5*rhoDiff);}
- ReSqrtRhoMatrix_[iChannel][iChannel] = reSqrtRho;
- ImSqrtRhoMatrix_[iChannel][iChannel] = imSqrtRho;
-
- }
+ ReSqrtRhoMatrix_[iChannel][iChannel] = reSqrtRho;
+ ImSqrtRhoMatrix_[iChannel][iChannel] = imSqrtRho;
+ }
}
-LauComplex LauKMatrixPropagator::getTHat(Double_t s, Int_t channel) {
+LauComplex LauKMatrixPropagator::getTHat(const Double_t s, const Int_t channel) {
+
+ LauComplex THat(0.0, 0.0);
- LauComplex THat(0.0, 0.0);
- channel -= 1;
- if (channel < 0 || channel >= nChannels_) {return THat;}
+ if (channel <= 0 || channel > nChannels_) {return THat;}
this->updatePropagator(s);
@@ -1157,8 +1416,8 @@
THatImag -= negImagProp_*ScattKMatrix_;
// Return the specific THat component
- THat.setRealPart(THatReal[index_][channel]);
- THat.setImagPart(THatImag[index_][channel]);
+ THat.setRealPart(THatReal[index_][channel-1]);
+ THat.setImagPart(THatImag[index_][channel-1]);
return THat;
diff --git a/src/LauResonanceMaker.cc b/src/LauResonanceMaker.cc
--- a/src/LauResonanceMaker.cc
+++ b/src/LauResonanceMaker.cc
@@ -625,6 +625,29 @@
}
}
+LauBlattWeisskopfFactor* LauResonanceMaker::getParentBWFactor(Int_t resSpin, LauBlattWeisskopfFactor::BarrierType barrierType)
+{
+ LauBlattWeisskopfFactor* bwFactor(0);
+
+ // Look up the category in the category information map
+ BWFactorCategoryMap::iterator factor_iter = bwFactors_.find( LauBlattWeisskopfFactor::Parent );
+
+ if ( factor_iter != bwFactors_.end() ) {
+ // If it exists, we can check if the factor object has been created
+ BlattWeisskopfCategoryInfo& categoryInfo = factor_iter->second;
+
+ if ( categoryInfo.bwFactor_ != 0 ) {
+ // If so, simply clone it
+ bwFactor = categoryInfo.bwFactor_->createClone( resSpin, barrierType );
+ }
+ }
+
+ if ( bwFactor==0 ) {
+ std::cerr<<"ERROR in LauResonanceMaker::getParentBWFactor : No parent Blatt-Weisskopf factor found to be cloned for K-matrix."<<std::endl;
+ }
+ return bwFactor;
+}
+
LauBlattWeisskopfFactor* LauResonanceMaker::getBWFactor( const LauBlattWeisskopfFactor::BlattWeisskopfCategory bwCategory, const LauResonanceInfo* resInfo )
{
LauBlattWeisskopfFactor* bwFactor(0);
@@ -654,7 +677,7 @@
if ( categoryInfo.bwFactor_ != 0 ) {
// If so, simply clone it
const UInt_t resSpin = resInfo->getSpin();
- bwFactor = categoryInfo.bwFactor_->createClone( resSpin );
+ bwFactor = categoryInfo.bwFactor_->createClone( resSpin, categoryInfo.bwFactor_->getBarrierType() );
} else {
// Otherwise we need to create it, using the default value if it has been set
if ( categoryInfo.defaultRadius_ >= 0.0 ) {
@@ -836,7 +859,7 @@
break;
case LauAbsResonance::KMatrix :
- // K-matrix description of S-wave
+ // K-matrix description
std::cerr<<"ERROR in LauResonanceMaker::getResonance : K-matrix type specified, which should be separately handled."<<std::endl;
break;