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();
@@ -473,37 +474,37 @@
LauAbsResonance& operator=(const LauAbsResonance& rhs);
//! Information on the resonance
- LauResonanceInfo* resInfo_;
+ LauResonanceInfo* resInfo_{0};
//! Information on the particles
const LauDaughters* daughters_;
//! Parent name
- TString nameParent_;
+ TString nameParent_{""};
//! Daughter 1 name
- TString nameDaug1_;
+ TString nameDaug1_{""};
//! Daughter 2 name
- TString nameDaug2_;
+ TString nameDaug2_{""};
//! Bachelor name
- TString nameBachelor_;
+ TString nameBachelor_{""};
//! Parent charge
- Int_t chargeParent_;
+ Int_t chargeParent_{0};
//! Daughter 1 charge
- Int_t chargeDaug1_;
+ Int_t chargeDaug1_{0};
//! Daughter 2 charge
- Int_t chargeDaug2_;
+ Int_t chargeDaug2_{0};
//! Bachelor charge
- Int_t chargeBachelor_;
+ Int_t chargeBachelor_{0};
//! Parent mass
- Double_t massParent_;
+ Double_t massParent_{0.0};
//! Daughter 1 mass
- Double_t massDaug1_;
+ Double_t massDaug1_{0.0};
//! Daughter 2 mass
- Double_t massDaug2_;
+ Double_t massDaug2_{0.0};
// Bachelor mass
- Double_t massBachelor_;
+ Double_t massBachelor_{0.0};
//! Resonance name
TString resName_;
@@ -512,49 +513,49 @@
TString sanitisedName_;
//! Resonance mass
- LauParameter* resMass_;
+ LauParameter* resMass_{0};
//! Resonance width
- LauParameter* resWidth_;
+ LauParameter* resWidth_{0};
//! All parameters of the resonance
std::vector<LauParameter*> resParameters_;
//! Resonance spin
- Int_t resSpin_;
+ Int_t resSpin_;
//! Resonance charge
- Int_t resCharge_;
+ Int_t resCharge_{0};
//! DP axis identifier
Int_t resPairAmpInt_;
//! Blatt Weisskopf barrier for parent decay
- LauBlattWeisskopfFactor* parBWFactor_;
+ LauBlattWeisskopfFactor* parBWFactor_{0};
//! Blatt Weisskopf barrier for resonance decay
- LauBlattWeisskopfFactor* resBWFactor_;
+ LauBlattWeisskopfFactor* resBWFactor_{0};
//! Spin formalism
- LauSpinType spinType_;
+ LauSpinType spinType_{Zemach_P};
//! Boolean to flip helicity
- Bool_t flipHelicity_;
+ Bool_t flipHelicity_{kFALSE};
//! Boolean to ignore the momentum factors in both the spin factor and the mass-dependent width
- Bool_t ignoreMomenta_;
+ Bool_t ignoreMomenta_{kFALSE};
//! Boolean to set the spinTerm to unity always
- Bool_t ignoreSpin_;
+ Bool_t ignoreSpin_{kFALSE};
//! Boolean to ignore barrier factor scaling in the amplitude numerator, they are still used for the mass-dependent width
- Bool_t ignoreBarrierScaling_;
+ Bool_t ignoreBarrierScaling_{kFALSE};
// Event kinematics information
//! Invariant mass
- Double_t mass_;
+ Double_t mass_{0.0};
//! Helicity angle cosine
- Double_t cosHel_;
+ Double_t cosHel_{0.0};
//! Daughter momentum in resonance rest frame
- Double_t q_;
+ Double_t q_{0.0};
//! Bachelor momentum in resonance rest frame
- Double_t p_;
+ Double_t p_{0.0};
//! Bachelor momentum in parent rest frame
- Double_t pstar_;
+ Double_t pstar_{0.0};
//! Covariant factor
/*!
@@ -566,10 +567,10 @@
\see LauKinematics::getcov13
\see LauKinematics::getcov23
*/
- Double_t erm_;
+ Double_t erm_{1.0};
//! Covariant factor (full spin-dependent expression)
- Double_t covFactor_;
+ Double_t covFactor_{1.0};
ClassDef(LauAbsResonance,0) // Abstract resonance class
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/LauKMatrixProdPole.hh b/inc/LauKMatrixProdPole.hh
--- a/inc/LauKMatrixProdPole.hh
+++ b/inc/LauKMatrixProdPole.hh
@@ -44,7 +44,7 @@
class LauKMatrixProdPole : public LauAbsResonance {
- public:
+ public:
//! Constructor
/*!
\param [in] poleName name of the pole
@@ -54,34 +54,34 @@
\param [in] daughters the daughter particles
\param [in] useProdAdler boolean to turn on/off the production Adler zero factor
*/
- LauKMatrixProdPole(const TString& poleName, Int_t poleIndex, Int_t resPairAmpInt,
- LauKMatrixPropagator* propagator, const LauDaughters* daughters,
- Bool_t useProdAdler = kFALSE);
+ LauKMatrixProdPole( const TString& poleName, Int_t poleIndex, Int_t resPairAmpInt,
+ LauKMatrixPropagator* propagator, const LauDaughters* daughters,
+ Bool_t useProdAdler = kFALSE);
//! Destructor
- virtual ~LauKMatrixProdPole();
+ virtual ~LauKMatrixProdPole();
// 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;}
+ 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] s the invariant-mass squared for the channel
+ \param [in] spinTerm the spin-term for the final channel
+ \return the complex amplitude
+ */
+ virtual LauComplex resAmp(Double_t s, Double_t spinTerm);
- private:
+ private:
//! Copy constructor (not implemented)
LauKMatrixProdPole(const LauKMatrixProdPole& rhs);
@@ -89,14 +89,14 @@
LauKMatrixProdPole& operator=(const LauKMatrixProdPole& rhs);
//! The K-matrix propagator
- LauKMatrixPropagator* thePropagator_;
+ LauKMatrixPropagator* thePropagator_;
//! The number of the pole
Int_t poleIndex_;
- //! 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(LauKMatrixProdPole, 0) // K-matrix production pole
+ ClassDef(LauKMatrixProdPole, 0) // K-matrix production pole
};
diff --git a/inc/LauKMatrixProdSVP.hh b/inc/LauKMatrixProdSVP.hh
--- a/inc/LauKMatrixProdSVP.hh
+++ b/inc/LauKMatrixProdSVP.hh
@@ -64,22 +64,22 @@
//! 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
*/
virtual LauAbsResonance::LauResonanceModel getResonanceModel() const {return LauAbsResonance::KMatrix;}
+
+ 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] s the invariant-mass squared for the channel
+ \param [in] spinTerm the spin-term for the final channel
+ \return the complex amplitude
+ */
+ virtual LauComplex resAmp(Double_t s, Double_t spinTerm);
private:
//! Copy constructor (not implemented)
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,15 +23,15 @@
*/
/*! \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
@@ -49,7 +49,7 @@
class LauKMatrixPropagator {
- public:
+ public:
//! Constructor
/*!
\param [in] name name of the propagator
@@ -59,18 +59,12 @@
\param [in] nPoles the number of poles
\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);
+ LauKMatrixPropagator( const TString& name, const TString& paramFileName,
+ Int_t resPairAmpInt, Int_t nChannels, Int_t nPoles,
+ 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
/*!
@@ -82,25 +76,28 @@
/*!
\param [in] inputFile name of the input file
*/
- void setParameters(const TString& inputFile);
+ void setParameters(const TString& inputFile);
- //! Get the scattering K matrix
- /*!
- \return the real, symmetric scattering K matrix
+ //! 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
*/
- TMatrixD getKMatrix() const {return ScattKMatrix_;}
+ TMatrixD getKMatrix() const {return ScattKMatrix_;}
- //! Get the real part of the propagator full matrix
- /*!
- \return the real part of the propagator full matrix
- */
- TMatrixD getRealPropMatrix() const {return realProp_;}
+ //! Get the real part of the propagator full matrix
+ /*!
+ \return the real part of the propagator full matrix
+ */
+ TMatrixD getRealPropMatrix() const {return realProp_;}
- //! Get the negative imaginary part of the full propagator matrix
- /*!
- \return the negative imaginary part of the full propagator matrix
- */
- TMatrixD getNegImagPropMatrix() const {return negImagProp_;}
+ //! Get the negative imaginary part of the full propagator matrix
+ /*!
+ \return the negative imaginary part of the full propagator matrix
+ */
+ TMatrixD getNegImagPropMatrix() const {return negImagProp_;}
//! Get the real part of the term of the propagator
/*!
@@ -123,6 +120,27 @@
*/
Double_t getPoleDenomTerm(Int_t poleIndex) const;
+ //! 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(Int_t iChannel){return L_[iChannel];}
+
+ //! Get index of final channel
+ /*!
+ \return the index of the channel into which the scattering happens
+ */
+ Int_t getIndex(){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 pole mass
+ */
+ LauParameter& getPoleMassSqParameter(Int_t poleIndex);
+
//! Get coupling constants that were loaded from the input file
/*!
\param [in] poleIndex number of the required pole
@@ -131,6 +149,14 @@
*/
Double_t getCouplingConstant(Int_t poleIndex, Int_t channelIndex) const;
+ //! Get coupling 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
+ */
+ LauParameter& getCouplingParameter(Int_t poleIndex, Int_t channelIndex);
+
//! Get scattering constants that were loaded from the input file
/*!
\param [in] channel1Index number of the first channel index
@@ -139,6 +165,32 @@
*/
Double_t getScatteringConstant(Int_t channel1Index, Int_t channel2Index) const;
+ //! Get scattering parameters, set according to the input file
+ /*!
+ \param [in] channel1Index number of the first channel index
+ \param [in] channel2Index number of the second channel index
+ \return the parameter of the scattering constant
+ */
+ LauParameter& getScatteringParameter(Int_t channel1Index, Int_t channel2Index);
+
+ //! Get s0 production parameter
+ /*!
+ \return the s0Prod parameter
+ */
+ LauParameter& gets0Prod() {return s0Prod_;}
+
+ //! Get sA production parameter
+ /*!
+ \return the sA parameter
+ */
+ LauParameter& getsA() {return sA_;}
+
+ //! Get sA0 production parameter
+ /*!
+ \return the sA0 parameter
+ */
+ LauParameter& getsA0() {return sA0_;}
+
//! Get the "slowly-varying part" term of the amplitude
/*!
\return the svp term
@@ -154,62 +206,59 @@
//! Get the DP axis identifier
/*!
- /return the value to identify the DP axis in question
+ \return the value to identify the DP axis in question
*/
Int_t getResPairAmpInt() const {return resPairAmpInt_;}
//! Get the number of channels
/*!
- /return the number of channels
+ \return the number of channels
*/
Int_t getNChannels() const {return nChannels_;}
//! Get the number of poles
/*!
- /return the number of poles
+ \return the number of poles
*/
Int_t getNPoles() const {return nPoles_;}
//! Get the propagator name
/*!
- /return the name of the propagator
+ \return the name of the propagator
*/
TString getName() const {return name_;}
- //! Get the unitary transition amplitude for the given channel
- /*!
- \param [in] s The invariant mass squared
- \param [in] channel The index number of the channel process
- \return the complex amplitude T
+ //! Get the unitary transition amplitude for the given channel
+ /*!
+ \param [in] s The invariant mass squared
+ \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(Double_t s, 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(Double_t s, 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(Double_t s, Int_t channel);
- protected:
+ protected:
//! Calculate the scattering K-matrix for the given value of s
/*!
\param [in] s the invariant mass squared
@@ -222,12 +271,39 @@
*/
void calcRhoMatrix(Double_t s);
+ //! 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, const Int_t 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);
+ //! Calculate the D0K+ phase space factor
+ /*!
+ \param [in] s the invariant mass squared
+ \return the complex phase space factor
+ */
+ LauComplex calcD0KRho(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(Double_t s) const;
+
//! Calculate the pipi phase space factor
/*!
\param [in] s the invariant mass squared
@@ -318,22 +394,22 @@
Bool_t checkPhaseSpaceType(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(Double_t s);
- //! Get the square root of the phase space matrix
- void getSqrtRhoMatrix();
-
- private:
+ //! Get the square root of the phase space matrix
+ void getSqrtRhoMatrix();
+
+ private:
//! Copy constructor (not implemented)
LauKMatrixPropagator(const LauKMatrixPropagator& rhs);
@@ -341,39 +417,58 @@
LauKMatrixPropagator& operator=(const LauKMatrixPropagator& rhs);
//! Create a map for the K-matrix parameters
- typedef std::map<int, std::vector<LauParameter> > KMatrixParamMap;
+ typedef std::map<int, std::vector<LauParameter> > KMatrixParamMap;
- //! Initialise and set the dimensions for the internal matrices and parameter arrays
- void initialiseMatrices();
+ //! 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
+ //! 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);
+ 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
+ //! 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);
+
+ //! 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 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 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
*/
- void storePole(const std::vector<std::string>& theLine);
+ void storeOrbitalAngularMomenta(const std::vector<std::string>& theLine);
- //! 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
*/
- void storeScattering(const std::vector<std::string>& theLine);
+ void storeBarrierFactorParameter(const std::vector<std::string>& theLine);
- //! 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
+ //! String to store the propagator name
TString name_;
//! Name of the input parameter file
TString paramFileName_;
@@ -397,8 +492,8 @@
// 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};
+ enum KMatrixChannels {Zero, PiPi, KK, FourPi, EtaEta, EtaEtaP,
+ KPi, KEtaP, KThreePi, D0K, Dstar0K, TotChannels};
//! Scattering K-matrix
TMatrixD ScattKMatrix_;
@@ -406,6 +501,8 @@
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
@@ -415,23 +512,33 @@
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_;
+ //! Vector of denominator parameters in centrifugal barrier factor
+ std::vector<Int_t> a_;
+ //! 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_;
+ LauParArray gCouplings_;
//! Array of scattering SVP values
- LauParArray fScattering_;
+ 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_;
@@ -481,6 +588,14 @@
Double_t fourPiFactor1_;
//! Factor used to calculate the pipipipi phase space term
Double_t fourPiFactor2_;
+ //! Defined as (mD0+mK)^2
+ Double_t mD0KSumSq_;
+ //! Defined as (mD0-mK)^2
+ Double_t mD0KDiffSq_;
+ //! Defined as (mD*0+mK)^2
+ Double_t mDstar0KSumSq_;
+ //! Defined as (mD*0-mK)^2
+ Double_t mDstar0KDiffSq_;
//! Multiplicative factor containing various Adler zero constants
Double_t adlerZeroFactor_;
@@ -490,8 +605,8 @@
//! Control the output of the functions
Bool_t verbose_;
- //! Control if scattering constants are channel symmetric: f_ji = f_ij
- Bool_t scattSymmetry_;
+ //! Control if scattering constants are channel symmetric: f_ji = f_ij
+ Bool_t scattSymmetry_;
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
@@ -85,30 +85,13 @@
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),
- spinType_(Zemach_P),
- flipHelicity_(kFALSE),
- ignoreMomenta_(kFALSE),
- ignoreSpin_(kFALSE),
- ignoreBarrierScaling_(kFALSE),
- mass_(0.0),
- cosHel_(0.0),
- q_(0.0),
- p_(0.0),
- pstar_(0.0),
- erm_(1.0),
- covFactor_(1.0)
+ resPairAmpInt_(resPairAmpInt)
{
if ( resInfo == 0 ) {
std::cerr << "ERROR in LauAbsResonance constructor : null LauResonanceInfo object provided" << std::endl;
@@ -142,33 +125,12 @@
}
// Constructor
-LauAbsResonance::LauAbsResonance(const TString& resName, const Int_t resPairAmpInt, const LauDaughters* daughters) :
- resInfo_(0),
+LauAbsResonance::LauAbsResonance(const TString& resName, const Int_t resPairAmpInt, const LauDaughters* daughters, const Int_t resSpin) :
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),
- spinType_(Zemach_P),
- flipHelicity_(kFALSE),
- ignoreMomenta_(kFALSE),
- ignoreSpin_(kFALSE),
- ignoreBarrierScaling_(kFALSE),
- mass_(0.0),
- cosHel_(0.0),
- q_(0.0),
- p_(0.0),
- pstar_(0.0),
- erm_(1.0),
- covFactor_(1.0)
+ resSpin_(resSpin),
+ resPairAmpInt_(resPairAmpInt)
{
if ( daughters_ == 0 ) {
std::cerr << "ERROR in LauAbsResonance constructor : null LauDaughters object provided" << std::endl;
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/LauKMatrixProdPole.cc b/src/LauKMatrixProdPole.cc
--- a/src/LauKMatrixProdPole.cc
+++ b/src/LauKMatrixProdPole.cc
@@ -28,28 +28,31 @@
#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()
@@ -57,13 +60,6 @@
}
LauComplex LauKMatrixProdPole::resAmp(Double_t mass, 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.
@@ -76,7 +72,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);
// Sum the pole denominator terms over all channels j, multiplying by
// the propagator terms. Note that we do not sum over poles, since we
@@ -103,6 +130,41 @@
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;
}
+
+const std::vector<LauParameter*>& LauKMatrixProdPole::getFloatingParameters()
+{
+
+ this->clearFloatingParameters();
+
+ Int_t nChannels = thePropagator_->getNChannels();
+
+ for (int jChannel = 0 ; jChannel < nChannels ; jChannel++)
+ {
+ LauParameter& par_gj_ = thePropagator_->getCouplingParameter(poleIndex_, jChannel);
+ if ( !par_gj_.fixed() )
+ this->addFloatingParameter( &par_gj_ );
+ }
+
+ for (int iPole = 0 ; iPole < thePropagator_->getNPoles() ; iPole++)
+ {
+ LauParameter& par_polemasssq_ = thePropagator_->getPoleMassSqParameter(poleIndex_);
+ 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()
@@ -59,13 +58,6 @@
}
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)
{
// 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,7 +113,59 @@
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;
}
+const std::vector<LauParameter*>& LauKMatrixProdSVP::getFloatingParameters()
+{
+
+ this->clearFloatingParameters();
+
+ Int_t nChannels = thePropagator_->getNChannels();
+
+ for (int jChannel = 0 ; jChannel < nChannels ; jChannel++)
+ {
+ LauParameter& par_f_ = thePropagator_->getScatteringParameter(channelIndex_, jChannel);
+ if ( !par_f_.fixed() )
+ {
+ this->addFloatingParameter( &par_f_ );
+ }
+
+ LauParameter& par_mSq0_ = thePropagator_->getPoleMassSqParameter(jChannel);
+ if ( !par_mSq0_.fixed() )
+ {
+ this->addFloatingParameter( &par_mSq0_ );
+ }
+
+ }
+
+ LauParameter& par_s0Prod_ = thePropagator_->gets0Prod();
+ if ( !par_s0Prod_.fixed() )
+ {
+ this->addFloatingParameter( &par_s0Prod_ );
+ }
+
+ LauParameter& par_sA_ = thePropagator_->getsA();
+ if ( !par_sA_.fixed() )
+ {
+ this->addFloatingParameter( &par_sA_ );
+ }
+
+ LauParameter& par_sA0_ = thePropagator_->getsA0();
+ if ( !par_sA0_.fixed() )
+ {
+ this->addFloatingParameter( &par_sA0_ );
+ }
+
+ return this->getParameters();
+
+}
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,14 @@
*/
/*! \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 "LauResonanceMaker.hh"
+#include "LauResonanceInfo.hh"
#include "TMath.h"
#include "TSystem.h"
@@ -45,8 +47,8 @@
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),
@@ -70,6 +72,10 @@
k3piFactor_(TMath::Power((1.44 - mK3piDiffSq_)/1.44, -2.5)),
fourPiFactor1_(16.0*LauConstants::mPiSq),
fourPiFactor2_(TMath::Sqrt(1.0 - fourPiFactor1_)),
+ mD0KSumSq_((LauConstants::mD0 + LauConstants::mK)*(LauConstants::mD0 + LauConstants::mK)),
+ mD0KDiffSq_((LauConstants::mD0 - LauConstants::mK)*(LauConstants::mD0 - LauConstants::mK)),
+ mDstar0KSumSq_((LauResonanceMaker::get().getResInfo("D*0")->getMass()->value() + LauConstants::mK)*(LauResonanceMaker::get().getResInfo("D*0")->getMass()->value() + LauConstants::mK)),
+ mDstar0KDiffSq_((LauResonanceMaker::get().getResInfo("D*0")->getMass()->value() - LauConstants::mK)*(LauResonanceMaker::get().getResInfo("D*0")->getMass()->value() - LauConstants::mK)),
adlerZeroFactor_(0.0),
parametersSet_(kFALSE),
verbose_(kFALSE),
@@ -77,15 +83,35 @@
{
// 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);
}
this->setParameters(paramFile);
+
+ // 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);
+ }
+ }
}
LauKMatrixPropagator::~LauKMatrixPropagator()
@@ -93,9 +119,10 @@
// 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();
@@ -136,36 +163,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)
{
- // 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 +193,26 @@
// 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)
// 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.
// 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 +222,6 @@
TMatrixD invC(A);
invC += B_invA_B;
-
// Set the real part of the propagator matrix ("C")
realProp_ = TMatrixD(TMatrixD::kInverted, invC);
@@ -226,6 +231,10 @@
negImagProp_ = TMatrixD(invA);
negImagProp_ *= BC;
+ // Pre-multiply by the Gamma matrix:
+ realProp_ = GammaMatrix_ * realProp_;
+ negImagProp_ = GammaMatrix_ * negImagProp_;
+
// Also calculate the production SVP term, since this uses Adler-zero parameters
// defined in the parameter file.
this->updateProdSVPTerm(s);
@@ -261,7 +270,9 @@
// 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
+ // Possible names are mSq0, s0Scatt, s0Prod, sA, sA0
+ // 5) Characteristic radius for each channel. If not set here, defaults to 3.0 GeV^{-1}
+ // "Radii radChannel1 radChannel2 ... radChannelN"
// By default, the scattering constants are symmetrised: f_{ji} = f_{ij}.
// To not assume this use "ScattSymmetry 0" on a line
@@ -272,9 +283,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 +316,21 @@
// Scattering terms
this->storeScattering(theLine);
+ } else if (!keyword.CompareTo("angularmomentum")) {
+
+ // Orbital angular momentum state for each channel
+ this->storeOrbitalAngularMomenta(theLine);
+
+ } else if (!keyword.CompareTo("barrierfactorparameter")) {
+
+ // Value of parameter "a" in denominator of centrifugal barrier factor, gamma
+ this->storeBarrierFactorParameter(theLine);
+
+ } else if (!keyword.CompareTo("radii")) {
+
+ // Values of characteristic radius
+ this->storeRadii(theLine);
+
} else {
// Usually Adler-zero constants
@@ -331,6 +357,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 +396,40 @@
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_.resize(nChannels_);
+ for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) {
+ L_.push_back(0);
+ }
+
+ // Characteristic radius (diagonal) vector (default to 3.0)
+ radii_.clear();
+ radii_.resize(nChannels_);
+ for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) {
+ radii_.push_back(3.0);
+ }
+
+ // Vector to store barrier factor parameter a
+ a_.clear();
+ a_.resize(nChannels_);
+
+ // 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
@@ -413,7 +470,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 +481,11 @@
if (checkChannel == kTRUE) {
cout<<"Adding phase space channel index "<<phaseSpaceInt
- <<" to K-matrix propagator "<<name_<<endl;
+ <<" to K-matrix propagator "<<name_<<endl;
phaseSpaceTypes_[iChannel] = phaseSpaceInt;
} else {
cerr<<"Phase space channel index "<<phaseSpaceInt
- <<" should be between 1 and "<<LauKMatrixPropagator::TotChannels-1<<endl;
+ <<" should be between 1 and "<<LauKMatrixPropagator::TotChannels-1<<endl;
}
}
@@ -452,7 +509,7 @@
Double_t poleMass = std::atof(theLine[2].c_str());
Double_t poleMassSq = poleMass*poleMass;
- LauParameter mPoleParam(poleMassSq);
+ LauParameter mPoleParam(Form("KM_poleMassSq_%i",poleIndex),poleMassSq);
mSqPoles_[poleIndex] = mPoleParam;
cout<<"Added bare pole mass "<<poleMass<<" GeV for pole number "<<poleIndex+1<<endl;
@@ -460,11 +517,11 @@
for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) {
Double_t couplingConst = std::atof(theLine[iChannel+3].c_str());
- LauParameter couplingParam(couplingConst);
+ LauParameter couplingParam(Form("KM_%s_gCoupling_%i_%i",name_.Data(),poleIndex,iChannel),couplingConst);
gCouplings_[poleIndex][iChannel] = couplingParam;
cout<<"Added coupling parameter g^{"<<poleIndex+1<<"}_"
- <<iChannel+1<<" = "<<couplingConst<<endl;
+ <<iChannel+1<<" = "<<couplingConst<<endl;
}
@@ -473,12 +530,11 @@
} else {
cerr<<"Error in LauKMatrixPropagator::storePole. Expecting "<<nExpect
- <<" numbers for pole definition but found "<<nWords
- <<" values instead"<<endl;
+ <<" numbers for pole definition but found "<<nWords
+ <<" values instead"<<endl;
}
-
}
void LauKMatrixPropagator::storeScattering(const std::vector<std::string>& theLine)
@@ -499,11 +555,11 @@
for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) {
Double_t scattConst = std::atof(theLine[iChannel+2].c_str());
- LauParameter scattParam(scattConst);
+ LauParameter scattParam(Form("KM_fScatteringConst_%i_%i",scattIndex,iChannel),scattConst);
fScattering_[scattIndex][iChannel] = scattParam;
cout<<"Added scattering parameter f("<<scattIndex+1<<","
- <<iChannel+1<<") = "<<scattConst<<endl;
+ <<iChannel+1<<") = "<<scattConst<<endl;
}
@@ -512,14 +568,108 @@
} else {
cerr<<"Error in LauKMatrixPropagator::storeScattering. Expecting "<<nExpect
- <<" numbers for scattering constants definition but found "<<nWords
- <<" values instead"<<endl;
+ <<" numbers for radii constants definition but found "<<nWords
+ <<" values instead"<<endl;
}
+}
+
+void LauKMatrixPropagator::storeOrbitalAngularMomenta(const std::vector<std::string>& theLine)
+{
+
+ // 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::atof(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;
+
+ }
}
+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)
+{
+
+ // 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;
+
+ }
+
+ } 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)
{
@@ -533,25 +683,25 @@
Double_t s0ScattValue = std::atof(parString.Data());
cout<<"Adler zero constant s0Scatt = "<<s0ScattValue<<endl;
- s0Scatt_ = LauParameter("s0Scatt", s0ScattValue);
+ s0Scatt_ = LauParameter(Form("KM_%s_s0Scatt",name_.Data()), s0ScattValue);
} else if (!keyword.CompareTo("s0prod")) {
Double_t s0ProdValue = std::atof(parString.Data());
cout<<"Adler zero constant s0Prod = "<<s0ProdValue<<endl;
- s0Prod_ = LauParameter("s0Scatt", s0ProdValue);
+ s0Prod_ = LauParameter(Form("KM_%s_s0Prod",name_.Data()), s0ProdValue);
} else if (!keyword.CompareTo("sa0")) {
Double_t sA0Value = std::atof(parString.Data());
cout<<"Adler zero constant sA0 = "<<sA0Value<<endl;
- sA0_ = LauParameter("sA0", sA0Value);
+ sA0_ = LauParameter(Form("KM_%s_sA0",name_.Data()), sA0Value);
} else if (!keyword.CompareTo("sa")) {
Double_t sAValue = std::atof(parString.Data());
cout<<"Adler zero constant sA = "<<sAValue<<endl;
- sA_ = LauParameter("sA", sAValue);
+ sA_ = LauParameter(Form("KM_%s_sA",name_.Data()), sAValue);
} else if (!keyword.CompareTo("scattsymmetry")) {
@@ -565,7 +715,6 @@
}
-
void LauKMatrixPropagator::calcScattKMatrix(Double_t s)
{
@@ -637,7 +786,6 @@
poleDenomVect_.push_back(invPoleTerm);
}
-
}
Double_t LauKMatrixPropagator::getPoleDenomTerm(Int_t poleIndex) const
@@ -651,6 +799,16 @@
}
+LauParameter& LauKMatrixPropagator::getPoleMassSqParameter(Int_t poleIndex)
+{
+ if ( (parametersSet_ == kFALSE) || (poleIndex < 0 || poleIndex >= nPoles_) ) {
+ std::cerr << "ERROR from LauKMatrixPropagator::getPoleMassSqParameter(). Invalid pole." << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+
+ return mSqPoles_[poleIndex];
+}
+
Double_t LauKMatrixPropagator::getCouplingConstant(Int_t poleIndex, Int_t channelIndex) const
{
@@ -663,6 +821,18 @@
}
+LauParameter& LauKMatrixPropagator::getCouplingParameter(Int_t poleIndex, Int_t channelIndex)
+{
+
+ if ( (parametersSet_ == kFALSE) || (poleIndex < 0 || poleIndex >= nPoles_) || (channelIndex < 0 || channelIndex >= nChannels_) ) {
+ std::cerr << "ERROR from LauKMatrixPropagator::getCouplingParameter(). Invalid coupling." << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+
+ //std::cout << "Minvalue + range for " << poleIndex << ", " << channelIndex << ": " << gCouplings_[poleIndex][channelIndex].minValue() << " => + " << gCouplings_[poleIndex][channelIndex].range() <<
+ // " and init value: " << gCouplings_[poleIndex][channelIndex].initValue() << std::endl;
+ return gCouplings_[poleIndex][channelIndex];
+}
Double_t LauKMatrixPropagator::getScatteringConstant(Int_t channel1Index, Int_t channel2Index) const
{
@@ -676,6 +846,17 @@
}
+LauParameter& LauKMatrixPropagator::getScatteringParameter(Int_t channel1Index, Int_t channel2Index)
+{
+
+ if ( (parametersSet_ == kFALSE) || (channel1Index < 0 || channel1Index >= nChannels_) || (channel2Index < 0 || channel2Index >= nChannels_) ) {
+ std::cerr << "ERROR from LauKMatrixPropagator::getScatteringParameter(). Invalid chanel index." << std::endl;
+ gSystem->Exit(EXIT_FAILURE);
+ }
+
+ return fScattering_[channel1Index][channel2Index];
+}
+
Double_t LauKMatrixPropagator::calcSVPTerm(Double_t s, Double_t s0) const
{
@@ -686,7 +867,7 @@
Double_t deltaS = s - s0;
if (TMath::Abs(deltaS) > 1.0e-6) {
result = (mSq0_.unblindValue() - s0)/deltaS;
- }
+ }
return result;
@@ -717,8 +898,62 @@
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);
+ Int_t phaseSpaceIndex(0);
+
+ for (Int_t iChannel (0); iChannel < nChannels_; ++iChannel) {
+
+ 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;
+
+ }
+
+}
+
+Double_t LauKMatrixPropagator::calcGamma(const Int_t iCh, const Double_t s, const Int_t phaseSpaceIndex) const
+{
+ // Calculate the DK P-wave angular momentum barrier factor
+ Double_t gamma(0.0);
+
+ Double_t sqrtTerm(0);
+ if (phaseSpaceIndex == LauKMatrixPropagator::D0K)
+ sqrtTerm = (s - mD0KSumSq_) * (s - mD0KDiffSq_) / (4*s);
+ else if (phaseSpaceIndex == LauKMatrixPropagator::Dstar0K)
+ sqrtTerm = (s - mDstar0KSumSq_) * (s - mDstar0KDiffSq_) / (4*s);
+ Double_t q(0.0);
+
+ q = TMath::Sqrt( fabs(sqrtTerm) );
+
+ gamma = pow(q,L_[iCh]);
+ if (includeBWBarrierFactor_)
+ {
+ gamma /= pow( q*q + gamAInvRadSq_[iCh] , L_[iCh]/2. );
+ }
+ return gamma;
}
void LauKMatrixPropagator::calcRhoMatrix(Double_t s)
@@ -729,8 +964,8 @@
// 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();
+ // Initialise all entries to zero
+ ReRhoMatrix_.Zero(); ImRhoMatrix_.Zero();
LauComplex rho(0.0, 0.0);
Int_t phaseSpaceIndex(0);
@@ -755,6 +990,10 @@
rho = this->calcKEtaPRho(s);
} else if (phaseSpaceIndex == LauKMatrixPropagator::KThreePi) {
rho = this->calcKThreePiRho(s);
+ } else if (phaseSpaceIndex == LauKMatrixPropagator::D0K) {
+ rho = this->calcD0KRho(s);
+ } else if (phaseSpaceIndex == LauKMatrixPropagator::Dstar0K) {
+ rho = this->calcDstar0KRho(s);
}
if (verbose_) {
@@ -763,12 +1002,48 @@
}
ReRhoMatrix_(iChannel, iChannel) = rho.re();
- ImRhoMatrix_(iChannel, iChannel) = rho.im();
+ ImRhoMatrix_(iChannel, iChannel) = rho.im();
}
}
+LauComplex LauKMatrixPropagator::calcD0KRho(Double_t s) const
+{
+ // Calculate the D0K+ phase space factor
+ LauComplex rho(0.0, 0.0);
+ if (TMath::Abs(s) < 1e-10) {return rho;}
+
+ 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(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
{
// Calculate the pipi phase space factor
@@ -804,30 +1079,30 @@
LauComplex LauKMatrixPropagator::calcFourPiRho(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;}
@@ -858,15 +1133,15 @@
LauComplex LauKMatrixPropagator::calcEtaEtaPRho(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;}
@@ -875,7 +1150,7 @@
if (sqrtTerm < 0.0) {
rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
} else {
- rho.setRealPart( TMath::Sqrt(sqrtTerm) );
+ rho.setRealPart( TMath::Sqrt(sqrtTerm) );
}
return rho;
@@ -957,13 +1232,13 @@
LauComplex LauKMatrixPropagator::getTransitionAmp(Double_t s, 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;}
- this->getTMatrix(s);
+ this->getTMatrix(s);
TAmp.setRealPart(ReTMatrix_[index_][channel]);
TAmp.setImagPart(ImTMatrix_[index_][channel]);
@@ -975,7 +1250,7 @@
LauComplex LauKMatrixPropagator::getPhaseSpaceTerm(Double_t s, 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;
@@ -983,7 +1258,7 @@
// 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]);
@@ -995,11 +1270,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;}
@@ -1023,22 +1298,22 @@
void LauKMatrixPropagator::getTMatrix(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_;
@@ -1075,44 +1350,43 @@
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 THat(0.0, 0.0);
+ LauComplex THat(0.0, 0.0);
channel -= 1;
if (channel < 0 || channel >= nChannels_) {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;