D41.1759196425.diff
No OneTemporary
Actions

Size

29 KB

Referenced Files

None

Subscribers

None

D41.1759196425.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
	@@ -49,7 +49,7 @@

	class LauKMatrixPropagator {

	- public:
	+ public:
	//! Constructor
	/*!
	\param [in] name name of the propagator
	@@ -59,12 +59,12 @@
	\param [in] nPoles the number of poles
	\param [in] rowIndex this specifies which row of the propagator should be used when summing over the amplitude channels
	*/
	- LauKMatrixPropagator(const TString& name, const TString& paramFileName,
	- Int_t resPairAmpInt, Int_t nChannels, Int_t nPoles,
	- Int_t rowIndex = 1);
	+ LauKMatrixPropagator( const TString& name, const TString& paramFileName,
	+ Int_t resPairAmpInt, Int_t nChannels, Int_t nPoles,
	+ Int_t rowIndex = 1 );

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

	//! Calculate the invariant mass squared s
	/*!
	@@ -82,25 +82,25 @@
	/*!
	\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
	+ //! 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
	/*!
	@@ -131,6 +131,12 @@
	*/
	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
	+ */s =
	//! Get scattering constants that were loaded from the input file
	/*!
	\param [in] channel1Index number of the first channel index
	@@ -143,6 +149,21 @@
	/*!
	\return the svp term
	*/
	+
	+ //! 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_;}
	+
	Double_t getProdSVPTerm() const {return prodSVP_;}

	//! Get the full complex propagator term for a given channel
	@@ -154,62 +175,60 @@

	//! 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 +241,41 @@
	*/
	void calcRhoMatrix(Double_t s);

	+ //! Calculate the real and imaginary part of the gamma angular-momentum barrier matrix
	+ /*!
	+ \param [in] s the invariant mass squared
	+ */
	+ void calcGammaMatrix(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
	+ */
	+ Double_t calcD0KPwaveGamma(Double_t s) 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 +366,22 @@
	Bool_t checkPhaseSpaceType(Int_t phaseSpaceInt) const;


	- //! Get the unitary transition amplitude matrix for the given kinematics
	- /*!
	- \param [in] kinematics The pointer to the constant kinematics
	+ //! Get the unitary transition amplitude matrix for the given kinematics
	+ /*!
	+ \param [in] kinematics The pointer to the constant kinematics
	*/
	- void getTMatrix(const LauKinematics* kinematics);
	+ void getTMatrix(const LauKinematics* kinematics);

	- //! Get the unitary transition amplitude matrix for the given kinematics
	- /*!
	- \param [in] s The invariant mass squared of the system
	+ //! Get the unitary transition amplitude matrix for the given kinematics
	+ /*!
	+ \param [in] s The invariant mass squared of the system
	*/
	- void getTMatrix(Double_t s);
	+ void getTMatrix(Double_t s);

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

	@@ -341,39 +389,39 @@
	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);
	+ 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
	+ //! 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);
	+ void storeScattering(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_;
	@@ -398,7 +446,7 @@
	// 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};
	+ KPi, KEtaP, KThreePi, D0K_Pwave, Dstar0K_Swave, TotChannels};

	//! Scattering K-matrix
	TMatrixD ScattKMatrix_;
	@@ -406,6 +454,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,10 +467,10 @@
	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_;
	@@ -426,12 +478,12 @@
	Int_t nPoles_;

	//! 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 phase space types
	std::vector<Int_t> phaseSpaceTypes_;
	@@ -481,17 +533,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/LauBelleNR.cc b/src/LauBelleNR.cc
	--- a/src/LauBelleNR.cc
	+++ b/src/LauBelleNR.cc
	@@ -49,7 +49,7 @@
	parName += "_alpha";
	alpha_ = resInfo->getExtraParameter( parName );
	if ( alpha_ == 0 ) {
	- alpha_ = new LauParameter( parName, 0.0, 0.0, 10.0, kTRUE );
	+ alpha_ = new LauParameter( parName, 0.0, -2.0, 10.0, kTRUE );
	alpha_->secondStage(kTRUE);
	resInfo->addExtraParameter( alpha_ );
	}
	diff --git a/src/LauKMatrixProdPole.cc b/src/LauKMatrixProdPole.cc
	--- a/src/LauKMatrixProdPole.cc
	+++ b/src/LauKMatrixProdPole.cc
	@@ -106,3 +106,21 @@
	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
	@@ -94,3 +94,26 @@

	}

	+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
	@@ -45,7 +45,7 @@
	ClassImp(LauKMatrixPropagator)

	LauKMatrixPropagator::LauKMatrixPropagator(const TString& name, const TString& paramFile,
	- Int_t resPairAmpInt, Int_t nChannels,
	+ Int_t resPairAmpInt, Int_t nChannels,
	Int_t nPoles, Int_t rowIndex) :
	name_(name),
	paramFileName_(paramFile),
	@@ -70,6 +70,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_((2.00696 + LauConstants::mK)*(2.00696 + LauConstants::mK)),
	+ mDstar0KDiffSq_((2.00696 - LauConstants::mK)*(2.00696 - LauConstants::mK)),
	adlerZeroFactor_(0.0),
	parametersSet_(kFALSE),
	verbose_(kFALSE),
	@@ -93,9 +97,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();
	@@ -153,10 +158,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);
	@@ -193,20 +201,28 @@
	// 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 diagonal
	+ this->calcGammaMatrix(s);
	+
	+ // Calculate the Xspin matrix, which is diagonal
	+ this->calcXspinMatrix();
	+
	+ // Calculate Krho(gamma^2) (real and imaginary parts, since rho can be complex)
	+ TMatrixD K_realRhoGammasq(ScattKMatrix_);
	+ K_realRhoGammasq *= ReRhoMatrix_;
	+ K_realRhoGammasq = (GammaMatrix_GammaMatrix_);
	+ TMatrixD K_imagRhoGammasq(ScattKMatrix_);
	+ K_imagRhoGammasq *= ImRhoMatrix_;
	+ K_imagRhoGammasq = (GammaMatrix_GammaMatrix_);

	// 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 +232,6 @@

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

	@@ -226,6 +241,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);
	@@ -365,14 +389,23 @@
	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.;
	+
	// 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
	@@ -460,7 +493,7 @@
	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<<"}_"
	@@ -663,6 +696,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 +721,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
	{

	@@ -755,7 +811,11 @@
	rho = this->calcKEtaPRho(s);
	} else if (phaseSpaceIndex == LauKMatrixPropagator::KThreePi) {
	rho = this->calcKThreePiRho(s);
	- }
	+ } else if (phaseSpaceIndex == LauKMatrixPropagator::D0K_Pwave) {
	+ rho = this->calcD0KRho(s);
	+ } else if (phaseSpaceIndex == LauKMatrixPropagator::Dstar0K_Swave) {
	+ rho = this->calcDstar0KRho(s);
	+ }

	if (verbose_) {
	cout<<"ReRhoMatrix("<<iChannel<<", "<<iChannel<<") = "<<rho.re()<<endl;
	@@ -769,6 +829,123 @@

	}

	+void LauKMatrixPropagator::calcGammaMatrix(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 (phaseSpaceIndex == LauKMatrixPropagator::D0K_Pwave) {
	+ gamma = this->calcD0KPwaveGamma(s);
	+ } else {
	+ gamma = 1.0; // S-wave
	+ }
	+
	+ if (verbose_) {
	+ cout<<"GammaMatrix("<<iChannel<<", "<<iChannel<<") = "<<gamma<<endl;
	+ }
	+
	+ GammaMatrix_(iChannel, iChannel) = 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);
	+ Int_t phaseSpaceIndex(0);
	+
	+ for (Int_t iChannel (0); iChannel < nChannels_; ++iChannel) {
	+
	+ phaseSpaceIndex = phaseSpaceTypes_[iChannel];
	+
	+ if (phaseSpaceIndex == LauKMatrixPropagator::D0K_Pwave) {
	+
	+ Xspin = cosHel_;
	+ } else {
	+ Xspin = 1.0; // S-wave
	+ }
	+
	+ if (verbose_) {
	+ cout<<"XspinMatrix("<<iChannel<<", "<<iChannel<<") = "<<Xspin<<endl;
	+ }
	+
	+ XspinMatrix_(iChannel, iChannel) = v;
	+
	+ }
	+
	+}
	+
	+
	+Double_t LauKMatrixPropagator::calcD0KPwaveGamma(Double_t s) const
	+{
	+ // Calculate the DK P-wave angular momentum barrier factor
	+ Double_t gamma(0.0);
	+
	+ Double_t sqrtTerm = (s - mD0KSumSq_) * (s - mD0KDiffSq_) / (4*s);
	+
	+ if (sqrtTerm < 0.0) {
	+ std::cerr << "ERROR in LauKMatrixPropagator::calcD0KPwaveGamma(). The centre of mass energy for the D0K system is below the m(D0K) threshold. "
	+ << "This should not happen" << std::endl;
	+ gSystem->Exit(EXIT_FAILURE);
	+ } else {
	+ gamma = TMath::Sqrt(sqrtTerm);
	+ }
	+
	+ return gamma;
	+}
	+
	+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

File Metadata

Mime Type: text/plain
Expires: Tue, Sep 30, 2:40 AM (22 h, 23 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 6567260
Default Alt Text: D41.1759196425.diff (29 KB)

D41.1759196425.diffNo OneTemporaryActions

D41.1759196425.diffView Options

File Metadata

Event Timeline

D41.1759196425.diff
No OneTemporary
Actions

D41.1759196425.diff
View Options