diff --git a/doc/ReleaseNotes.md b/doc/ReleaseNotes.md --- a/doc/ReleaseNotes.md +++ b/doc/ReleaseNotes.md @@ -1,5 +1,9 @@ # Laura++ release notes +4th February 2020 Dan Johnson +* Extend the K-matrix implementation to handle non-zero spin + - see https://phab.hepforge.org/T135 + 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 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& 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& 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 #include +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 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 > 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& 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& 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& 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& 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& 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& theLine); + void storeOrbitalAngularMomenta(const std::vector& theLine, std::vector& 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& theLine); + void storeBarrierFactorParameter(const std::vector& theLine, std::vector& 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 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 mSqPoles_; + std::vector mSqPoles_; //! Array of coupling constants LauParArray gCouplings_; //! Array of scattering SVP values LauParArray fScattering_; + //! Vector of characteristic radii + std::vector radii_; + //! Vector of ratio a/R^2 + std::vector gamAInvRadSq_; //! Vector of phase space types - std::vector phaseSpaceTypes_; + std::vector phaseSpaceTypes_; //! Vector of squared masses std::vector 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 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 "<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 "<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 = "<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& theLine) +void LauKMatrixPropagator::storeChannels(const std::vector& 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(theLine.size()) - 1; if (nTypes != nChannels_) { cerr<<"Error in LauKMatrixPropagator::storeChannels. The input file defines " - <(phaseSpaceInt); } else { cerr<<"Phase space channel index "<(LauKMatrixPropagator::KMatrixChannels::TotChannels)-1<& 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("<& theLine, std::vector& 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(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 "<Exit(EXIT_FAILURE); + } + } + } + +} + +void LauKMatrixPropagator::storeRadii(const std::vector& 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(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 "<& theLine, std::vector& 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(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 "<= 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("<Exit(EXIT_FAILURE); + } + return rho; +} - if (verbose_) { - cout<<"ReRhoMatrix("< 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(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."<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."<