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& 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& 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 /*! @@ -131,6 +135,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 +151,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 +180,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 +245,42 @@ */ void calcRhoMatrix(Double_t s); + //! Calculate the (real) 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 + \return the centrifugal barrier factor for L=0,1, or 2 + */ + Double_t calcGamma(Int_t iCh, Double_t s, 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 +371,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 +394,52 @@ LauKMatrixPropagator& operator=(const LauKMatrixPropagator& rhs); //! Create a map for the K-matrix parameters - typedef std::map > KMatrixParamMap; + typedef std::map > 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& theLine); + 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 + //! 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); + void storePole(const std::vector& 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& theLine); + void storeScattering(const std::vector& 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 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 storeParameter(const TString& keyword, const TString& parString); + void storeRadii(const std::vector& theLine); + + //! 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 storeBarrierFactorParameter(const std::vector& theLine); - //! String to store the propagator name + //! 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); + + + //! String to store the propagator name TString name_; //! Name of the input parameter file TString paramFileName_; @@ -397,8 +463,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 +472,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 +485,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 mSqPoles_; + std::vector mSqPoles_; //! Array of coupling constants - LauParArray gCouplings_; + LauParArray gCouplings_; //! Array of scattering SVP values - LauParArray fScattering_; + 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_; @@ -481,17 +561,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 @@ -106,3 +106,21 @@ return amp; } + +const std::vector& 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& 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("D*0")->getMass()->value() + LauConstants::mK)*(LauResonanceMaker::get().getResInfo("D*0")->getMass()->value() + LauConstants::mK)), // TODO: Can we avoid hard-coding the D*0 mass? Can't see how to access LauResonanceMaker info without creating a resonance + mDstar0KDiffSq_((LauResonanceMaker::get().getResInfo("D*0")->getMass()->value() - LauConstants::mK)*(LauResonanceMaker::get().getResInfo("D*0")->getMass()->value() - LauConstants::mK)), adlerZeroFactor_(0.0), parametersSet_(kFALSE), verbose_(kFALSE), @@ -77,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 - iK*rho)^-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 - 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 +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 K*rho*(gamma^2) (real and imaginary parts, since rho can be complex) + TMatrixD GammaMatrixSq = (GammaMatrix_*GammaMatrix_); + TMatrixD K_realRhoGammaSq(ScattKMatrix_); + K_realRhoGammaSq *= ReRhoMatrix_; + K_realRhoGammaSq *= GammaMatrixSq; + TMatrixD K_imagRhoGammaSq(ScattKMatrix_); + K_imagRhoGammaSq *= ImRhoMatrix_; + K_imagRhoGammaSq *= GammaMatrixSq; // A = I + K*Imag(rho), B = -K*Real(Rho) // Calculate C and D matrices such that (A + iB)*(C + iD) = I, // ie. C + iD = (I - i K*rho)^-1, separated into real and imaginary parts. // realProp C = (A + B A^-1 B)^-1, imagProp D = -A^-1 B C TMatrixD A(IMatrix_); - A += K_imagRho; + A += K_imagRhoGammaSq; TMatrixD B(zeroMatrix_); - B -= K_realRho; + B -= K_realRhoGammaSq; TMatrixD invA(TMatrixD::kInverted, A); TMatrixD invA_B(invA); @@ -216,7 +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(theLine.size()) - 1; if (nTypes != nChannels_) { cerr<<"Error in LauKMatrixPropagator::storeChannels. The input file defines " - <& 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("<& 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 "<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,93 @@ Double_t deltaS = s - sA0Val; if (TMath::Abs(deltaS) > 1e-6) { adlerZeroFactor_ = (s - sAConst_)*(1.0 - sA0Val)/deltaS; - } + } + +} + +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 ( L_ != 0 ) { + gamma = this->calcGamma(iChannel,s,phaseSpaceIndex); + } else { + gamma = 1.0; // S-wave + } + + if (verbose_) { + cout<<"GammaMatrix("< 4pi state). - // The normalisation term rho_0 is found by ensuring that the phase space integral - // at s = 1 is equal to sqrt(1.0 - 16*mpiSq/s). Note that the exponent for this - // factor in hep-ph/0204328 is wrong; it should be 0.5, i.e. sqrt, not n = 1 to 5. - // Plotting the value of this double integral as a function of s can then be fitted - // to a 6th-order polynomial (for s < 1), which is the result used below + // parameterisation that approximates the multi-body phase space double integral + // defined in Eq 4 of the A&S paper hep-ph/0204328. This form agrees with the + // BaBar model (another 6th order polynomial from s^4 down to 1/s^2), but avoids the + // exponential increase at small values of s (~< 0.1) arising from 1/s and 1/s^2. + // Eq 4 is evaluated for each value of s by assuming incremental steps of 1e-3 GeV^2 + // for s1 and s2, the invariant energy squared of each of the di-pion states, + // with the integration limits of s1 = (2*mpi)^2 to (sqrt(s) - 2*mpi)^2 and + // s2 = (2*mpi)^2 to (sqrt(s) - sqrt(s1))^2. The mass M of the rho is taken to be + // 0.775 GeV and the energy-dependent width of the 4pi system + // Gamma(s) = gamma_0*rho1^3(s), where rho1 = sqrt(1.0 - 4*mpiSq/s) and gamma_0 is + // the "width" of the 4pi state at s = 1, which is taken to be 0.3 GeV + // (~75% of the total width from PDG estimates of the f0(1370) -> 4pi state). + // The normalisation term rho_0 is found by ensuring that the phase space integral + // at s = 1 is equal to sqrt(1.0 - 16*mpiSq/s). Note that the exponent for this + // factor in hep-ph/0204328 is wrong; it should be 0.5, i.e. sqrt, not n = 1 to 5. + // Plotting the value of this double integral as a function of s can then be fitted + // to a 6th-order polynomial (for s < 1), which is the result used below LauComplex rho(0.0, 0.0); if (TMath::Abs(s) < 1e-10) {return rho;} if (s <= 1.0) { - Double_t rhoTerm = ((1.07885*s + 0.13655)*s - 0.29744)*s - 0.20840; - rhoTerm = ((rhoTerm*s + 0.13851)*s - 0.01933)*s + 0.00051; + Double_t rhoTerm = ((1.07885*s + 0.13655)*s - 0.29744)*s - 0.20840; + rhoTerm = ((rhoTerm*s + 0.13851)*s - 0.01933)*s + 0.00051; // For some values of s (below 2*mpi), this term is a very small // negative number. Check for this and set the rho term to zero. if (rhoTerm < 0.0) {rhoTerm = 0.0;} @@ -858,15 +1119,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 +1136,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 +1218,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 +1236,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 +1244,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 +1256,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 +1284,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 +1336,43 @@ void LauKMatrixPropagator::getSqrtRhoMatrix() { - // Find the square root of the (current) phase space matrix so that - // we can find T = [sqrt(rho)}^{*} T_hat sqrt(rho), where T_hat is the - // Lorentz-invariant T matrix = (I - i K rho)^-1 * K; note that the first - // sqrt of rho matrix is complex conjugated + // Find the square root of the (current) phase space matrix so that + // we can find T = [sqrt(rho)}^{*} T_hat sqrt(rho), where T_hat is the + // Lorentz-invariant T matrix = (I - i K rho)^-1 * K; note that the first + // sqrt of rho matrix is complex conjugated - // If rho = rho_i + i rho_r = a + i b, then sqrt(rho) = c + i d, where - // c = sqrt(0.5*(r+a)) and d = sqrt(0.5(r-a)), where r = sqrt(a^2 + b^2). - // Since rho is diagonal, then the square root of rho will also be diagonal, - // with its real and imaginary matrix elements equal to c and d, respectively + // If rho = rho_i + i rho_r = a + i b, then sqrt(rho) = c + i d, where + // c = sqrt(0.5*(r+a)) and d = sqrt(0.5(r-a)), where r = sqrt(a^2 + b^2). + // Since rho is diagonal, then the square root of rho will also be diagonal, + // with its real and imaginary matrix elements equal to c and d, respectively - // Initialise the real and imaginary parts of the square root of - // the rho matrix to zero - ReSqrtRhoMatrix_.Zero(); ImSqrtRhoMatrix_.Zero(); + // Initialise the real and imaginary parts of the square root of + // the rho matrix to zero + ReSqrtRhoMatrix_.Zero(); ImSqrtRhoMatrix_.Zero(); - for (Int_t iChannel (0); iChannel < nChannels_; ++iChannel) { + for (Int_t iChannel (0); iChannel < nChannels_; ++iChannel) { - Double_t realRho = ReRhoMatrix_[iChannel][iChannel]; - Double_t imagRho = ImRhoMatrix_[iChannel][iChannel]; + Double_t realRho = ReRhoMatrix_[iChannel][iChannel]; + Double_t imagRho = ImRhoMatrix_[iChannel][iChannel]; Double_t rhoMag = sqrt(realRho*realRho + imagRho*imagRho); - Double_t rhoSum = rhoMag + realRho; - Double_t rhoDiff = rhoMag - realRho; + Double_t rhoSum = rhoMag + realRho; + Double_t rhoDiff = rhoMag - realRho; - Double_t reSqrtRho(0.0), imSqrtRho(0.0); - if (rhoSum > 0.0) {reSqrtRho = sqrt(0.5*rhoSum);} + Double_t reSqrtRho(0.0), imSqrtRho(0.0); + if (rhoSum > 0.0) {reSqrtRho = sqrt(0.5*rhoSum);} if (rhoDiff > 0.0) {imSqrtRho = sqrt(0.5*rhoDiff);} - ReSqrtRhoMatrix_[iChannel][iChannel] = reSqrtRho; - ImSqrtRhoMatrix_[iChannel][iChannel] = imSqrtRho; - - } + ReSqrtRhoMatrix_[iChannel][iChannel] = reSqrtRho; + ImSqrtRhoMatrix_[iChannel][iChannel] = imSqrtRho; + } } LauComplex LauKMatrixPropagator::getTHat(Double_t s, Int_t channel) { - LauComplex THat(0.0, 0.0); + LauComplex THat(0.0, 0.0); channel -= 1; if (channel < 0 || channel >= nChannels_) {return THat;}