D41.1759117217.diff
No OneTemporary
Actions

Size

57 KB

Referenced Files

None

Subscribers

None

D41.1759117217.diff
View Options

	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,36 @@
	\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;}
	+ 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 LauComplex amplitude(const LauKinematics* kinematics);

	//! 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);
	+ virtual LauComplex resAmp(Double_t mass, Double_t spinTerm);

	- private:
	+ private:
	//! Copy constructor (not implemented)
	LauKMatrixProdPole(const LauKMatrixProdPole& rhs);

	@@ -89,14 +91,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
	@@ -76,6 +76,8 @@
	\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
	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
	@@ -58,13 +58,14 @@
	\param [in] nChannels the number of channels
	\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
	+ \param [in] L this specifies the spin of this K-matrix
	*/
	- 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, Int_t L = 0 );

	//! Destructor
	- virtual ~LauKMatrixPropagator();
	+ virtual ~LauKMatrixPropagator();

	//! Calculate the invariant mass squared s
	/*!
	@@ -82,25 +83,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 +127,12 @@
	*/
	Double_t getPoleDenomTerm(Int_t poleIndex) const;

	+ //! Get spin of K-matrix
	+ /*!
	+ \return the value of the spin, L_
	+ */
	+ Int_t getL(){return L_;}
	+
	//! Get coupling constants that were loaded from the input file
	/*!
	\param [in] poleIndex number of the required pole
	@@ -131,6 +141,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 +157,20 @@
	*/
	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 the "slowly-varying part" term of the amplitude
	/*!
	\return the svp term
	@@ -154,62 +186,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 +251,42 @@
	*/
	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 Xspin matrix
	+ void calcXspinMatrix();
	+
	+ //! 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 +377,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();

	- //! Get the square root of the phase space matrix
	- void getSqrtRhoMatrix();
	-
	- private:
	+ private:
	//! Copy constructor (not implemented)
	LauKMatrixPropagator(const LauKMatrixPropagator& rhs);

	@@ -341,39 +400,52 @@
	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 storePole(const std::vector<std::string>& theLine);
	+ void storeRadii(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 +469,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 +478,10 @@
	TMatrixD ReRhoMatrix_;
	//! Imaginary part of the phase space density diagonal matrix
	TMatrixD ImRhoMatrix_;
	+ //! Gamma angular-momentum barrier matrix
	+ TMatrixD GammaMatrix_;
	+ //! Xspin matrix
	+ TMatrixD XspinMatrix_;
	//! Identity matrix
	TMatrixD IMatrix_;
	//! Null matrix
	@@ -415,23 +491,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_;
	+ //! Spin of K-matrix
	+ Int_t L_;
	+ //! Denominator parameter in centrifugal barrier factor
	+ 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,17 +567,28 @@
	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_;
	//! Tracks if all params have been set
	Bool_t parametersSet_;

	+ //! Helicity angle cosine
	+ Double_t cosHel_;
	+
	//! 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/src/LauKMatrixProdPole.cc b/src/LauKMatrixProdPole.cc
	--- a/src/LauKMatrixProdPole.cc
	+++ b/src/LauKMatrixProdPole.cc
	@@ -38,16 +38,16 @@
	Bool_t useProdAdler) :
	LauAbsResonance(poleName, resPairAmpInt, daughters),
	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;
	}

	}
	@@ -76,6 +76,45 @@
	return amp;
	}

	+ const Int_t resSpin = this->getSpin();
	+ const Double_t q = this->getQ();
	+ const Double_t p = this->getP();
	+ const Double_t pstar = this->getPstar();
	+
	+ // Get barrier scaling factors
	+ Double_t fFactorR(1.0);
	+ Double_t fFactorB(1.0);
	+ if ( thePropagator_->getL() > 0 ) {
	+ const LauBlattWeisskopfFactor* resBWFactor = this->getResBWFactor();
	+ if ( resBWFactor != nullptr ) {
	+ fFactorR = resBWFactor->calcFormFactor(q);
	+ }
	+
	+ const LauBlattWeisskopfFactor* parBWFactor = this->getParBWFactor();
	+ if ( parBWFactor != nullptr ) {
	+ switch ( parBWFactor->getRestFrame() ) {
	+ case LauBlattWeisskopfFactor::ResonanceFrame:
	+ fFactorB = parBWFactor->calcFormFactor(p);
	+ 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(kinematics);

	// Sum the pole denominator terms over all channels j, multiplying by
	@@ -103,6 +142,27 @@

	amp.rescale(poleDenom*adlerZero);

	+ amp.rescale(fFactorR);
	+ amp.rescale(fFactorB);
	+
	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_ );
	+ }
	+
	+ 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
	@@ -42,15 +42,14 @@
	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()
	@@ -78,6 +77,45 @@
	return amp;
	}

	+ const Int_t resSpin = this->getSpin();
	+ const Double_t q = this->getQ();
	+ const Double_t p = this->getP();
	+ const Double_t pstar = this->getPstar();
	+
	+ // Get barrier scaling factors
	+ Double_t fFactorR(1.0);
	+ Double_t fFactorB(1.0);
	+ if ( thePropagator_->getL() > 0 ) {
	+ const LauBlattWeisskopfFactor* resBWFactor = this->getResBWFactor();
	+ if ( resBWFactor != nullptr ) {
	+ fFactorR = resBWFactor->calcFormFactor(q);
	+ }
	+
	+ const LauBlattWeisskopfFactor* parBWFactor = this->getParBWFactor();
	+ if ( parBWFactor != nullptr ) {
	+ switch ( parBWFactor->getRestFrame() ) {
	+ case LauBlattWeisskopfFactor::ResonanceFrame:
	+ fFactorB = parBWFactor->calcFormFactor(p);
	+ 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(kinematics);

	Double_t SVPTerm = thePropagator_->getProdSVPTerm();
	@@ -90,7 +128,33 @@

	amp.rescale(SVPTerm*adlerZero);

	+ amp.rescale(fFactorR);
	+ amp.rescale(fFactorB);
	+
	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_s0Prod_ = thePropagator_->gets0Prod();
	+ if ( !par_s0Prod_.fixed() )
	+ {
	+ this->addFloatingParameter( &par_s0Prod_ );
	+ }
	+
	+ return this->getParameters();
	+
	+}
	\ No newline at end of file
	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, Int_t L) :
	name_(name),
	paramFileName_(paramFile),
	resPairAmpInt_(resPairAmpInt),
	@@ -56,6 +58,7 @@
	prodSVP_(0.0),
	nChannels_(nChannels),
	nPoles_(nPoles),
	+ L_(L),
	sAConst_(0.0),
	m2piSq_(4.0*LauConstants::mPiSq),
	m2KSq_( 4.0*LauConstants::mKSq),
	@@ -70,6 +73,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("D0")->getMass()->value() + LauConstants::mK)(LauResonanceMaker::get().getResInfo("D*0")->getMass()->value() + LauConstants::mK)),
	+ mDstar0KDiffSq_((LauResonanceMaker::get().getResInfo("D0")->getMass()->value() - LauConstants::mK)(LauResonanceMaker::get().getResInfo("D*0")->getMass()->value() - LauConstants::mK)),
	adlerZeroFactor_(0.0),
	parametersSet_(kFALSE),
	verbose_(kFALSE),
	@@ -77,11 +84,28 @@
	{
	// 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;
	+ // Set default value of spin-dependent centrifugal-barrier-factor parameter
	+ switch(L_) {
	+ case 0:
	+ a_ = 0;
	+ break;
	+ case 1:
	+ a_ = 1;
	+ break;
	+ case 2:
	+ a_ = 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);
	+ }
	+
	+ // 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,9 +117,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();
	@@ -139,7 +164,7 @@
	void LauKMatrixPropagator::updatePropagator(const LauKinematics* kinematics)
	{

	- // 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 - iKrho)^-1 P_j, where P is the production K-matrix.
	// i = index for the state (e.g. S-wave index = 0).
	@@ -153,10 +178,13 @@

	if (resPairAmpInt_ == 1) {
	s = kinematics->getm23Sq();
	+ cosHel_ = kinematics->getc23();
	} else if (resPairAmpInt_ == 2) {
	s = kinematics->getm13Sq();
	+ cosHel_ = kinematics->getc13();
	} else if (resPairAmpInt_ == 3) {
	s = kinematics->getm12Sq();
	+ cosHel_ = kinematics->getc12();
	}

	this->updatePropagator(s);
	@@ -165,7 +193,7 @@

	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 - iKrho)^-1 P_j, where P is the production K-matrix.
	// i = index for the state (e.g. S-wave index = 0).
	@@ -193,20 +221,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 the Xspin matrix, which is real and diagonal
	+ this->calcXspinMatrix();
	+
	+ // Calculate Krho(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 + KImag(rho), B = -KReal(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 +253,6 @@

	TMatrixD invC(A);
	invC += B_invA_B;
	-
	// Set the real part of the propagator matrix ("C")
	realProp_ = TMatrixD(TMatrixD::kInverted, invC);

	@@ -226,6 +262,15 @@
	negImagProp_ = TMatrixD(invA);
	negImagProp_ *= BC;

	+ // Pre-multiply by the Gamma matrix:
	+ realProp_ = GammaMatrix_ * realProp_;
	+ negImagProp_ = GammaMatrix_ * negImagProp_;
	+
	+ // Pre-multiply by the Xspin matrix:
	+ realProp_ = XspinMatrix_ * realProp_;
	+ negImagProp_ = XspinMatrix_ * negImagProp_;
	+
	+
	// Also calculate the production SVP term, since this uses Adler-zero parameters
	// defined in the parameter file.
	this->updateProdSVPTerm(s);
	@@ -261,7 +306,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 +319,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 +352,11 @@
	// Scattering terms
	this->storeScattering(theLine);

	+ } else if (!keyword.CompareTo("radii")) {
	+
	+ // Values of characteristic radius
	+ this->storeRadii(theLine);
	+
	} else {

	// Usually Adler-zero constants
	@@ -331,6 +383,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_/(radii_[iChannel]*radii_[iChannel]);
	+ }
	+
	// All required parameters have been set
	parametersSet_ = kTRUE;

	@@ -365,14 +422,34 @@
	ReRhoMatrix_.ResizeTo(nChannels_, nChannels_);
	ImRhoMatrix_.ResizeTo(nChannels_, nChannels_);

	+ // Gamma matrices
	+ GammaMatrix_.Clear();
	+ GammaMatrix_.ResizeTo(nChannels_, nChannels_);
	+
	+ // Xspin matrices
	+ XspinMatrix_.Clear();
	+ XspinMatrix_.ResizeTo(nChannels_, nChannels_);
	+ cosHel_=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 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 +490,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 +501,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;
	}

	}
	@@ -460,11 +537,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_gCoupling_%i_%i",poleIndex,iChannel),couplingConst);
	gCouplings_[poleIndex][iChannel] = couplingParam;

	cout<<"Added coupling parameter g^{"<<poleIndex+1<<"}_"
	- <<iChannel+1<<" = "<<couplingConst<<endl;
	+ <<iChannel+1<<" = "<<couplingConst<<endl;

	}

	@@ -473,12 +550,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 +575,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,11 +588,42 @@
	} else {

	cerr<<"Error in LauKMatrixPropagator::storeScattering. Expecting "<<nExpect
	- <<" numbers for scattering constants definition but found "<<nWords
	- <<" values instead"<<endl;
	+ <<" numbers for scattering constants definition but found "<<nWords
	+ <<" values instead"<<endl;

	}

	+}
	+
	+void LauKMatrixPropagator::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 scattering constants definition but found "<<nWords
	+ <<" values instead"<<endl;
	+
	+ }

	}

	@@ -539,7 +646,7 @@

	Double_t s0ProdValue = std::atof(parString.Data());
	cout<<"Adler zero constant s0Prod = "<<s0ProdValue<<endl;
	- s0Prod_ = LauParameter("s0Scatt", s0ProdValue);
	+ s0Prod_ = LauParameter("s0Prod", s0ProdValue);

	} else if (!keyword.CompareTo("sa0")) {

	@@ -561,11 +668,17 @@
	scattSymmetry_ = kFALSE;
	}

	+ } else if (!keyword.CompareTo("barrierfactorparameter")) {
	+
	+ // Value of parameter "a" in denominator of centrifugal barrier factor, gamma
	+ Double_t aValue = std::atof(parString.Data());
	+ cout<<"K-matrix barrier factor denominator parameter = "<<aValue<<endl;
	+ a_ = aValue;
	+
	}

	}

	-
	void LauKMatrixPropagator::calcScattKMatrix(Double_t s)
	{

	@@ -637,7 +750,6 @@
	poleDenomVect_.push_back(invPoleTerm);

	}
	-
	}

	Double_t LauKMatrixPropagator::getPoleDenomTerm(Int_t poleIndex) const
	@@ -663,6 +775,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 +800,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 +821,7 @@
	Double_t deltaS = s - s0;
	if (TMath::Abs(deltaS) > 1.0e-6) {
	result = (mSq0_.unblindValue() - s0)/deltaS;
	- }
	+ }

	return result;

	@@ -717,7 +852,94 @@
	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_ != 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_);
	+ if (includeBWBarrierFactor_)
	+ {
	+ gamma /= pow( q*q + gamAInvRadSq_[iCh] , L_/2. );
	+ }
	+
	+ return gamma;
	+}
	+
	+void LauKMatrixPropagator::calcXspinMatrix()
	+{
	+ // Calculate the X spin matrix (diagonal: 1 for S wave channel, cos(theta) for P wave channel)
	+
	+ // Initialise all entries to zero
	+ XspinMatrix_.Zero();
	+
	+ Double_t Xspin(0.0);
	+
	+ for (Int_t iChannel (0); iChannel < nChannels_; ++iChannel) {
	+
	+ switch(L_) {
	+ case 0 :
	+ Xspin = 1.0;
	+ break;
	+ case 1 :
	+ Xspin = cosHel_;
	+ break;
	+ case 2 :
	+ Xspin = 3cosHel_cosHel_ - 1.0;
	+ break;
	+ }
	+
	+ if (verbose_) {
	+ cout<<"XspinMatrix("<<iChannel<<", "<<iChannel<<") = "<<Xspin<<endl;
	+ }
	+
	+ XspinMatrix_(iChannel, iChannel) = Xspin;
	+
	+ }

	}

	@@ -729,8 +951,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 +977,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 +989,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 +1066,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 = (2mpi)^2 to (sqrt(s) - 2mpi)^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_0rho1^3(s), where rho1 = sqrt(1.0 - 4mpiSq/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 = (2mpi)^2 to (sqrt(s) - 2mpi)^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_0rho1^3(s), where rho1 = sqrt(1.0 - 4mpiSq/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.07885s + 0.13655)s - 0.29744)*s - 0.20840;
	- rhoTerm = ((rhoTerms + 0.13851)s - 0.01933)*s + 0.00051;
	+ Double_t rhoTerm = ((1.07885s + 0.13655)s - 0.29744)*s - 0.20840;
	+ rhoTerm = ((rhoTerms + 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 +1120,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 +1137,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 +1219,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 +1237,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 +1245,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 +1257,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 +1285,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 +1337,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(realRhorealRho + imagRhoimagRho);
	- 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;}

File Metadata

Mime Type: text/plain
Expires: Mon, Sep 29, 4:40 AM (7 h, 54 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 6550777
Default Alt Text: D41.1759117217.diff (57 KB)

D41.1759117217.diffNo OneTemporaryActions

D41.1759117217.diffView Options

File Metadata

Event Timeline

D41.1759117217.diff
No OneTemporary
Actions

D41.1759117217.diff
View Options