Changeset View
Standalone View
src/LauKMatrixPropagator.cc
Show All 39 Lines | |||||
using std::cout; | using std::cout; | ||||
using std::endl; | using std::endl; | ||||
using std::cerr; | using std::cerr; | ||||
ClassImp(LauKMatrixPropagator) | ClassImp(LauKMatrixPropagator) | ||||
LauKMatrixPropagator::LauKMatrixPropagator(const TString& name, const TString& paramFile, | 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) : | Int_t nPoles, Int_t rowIndex) : | ||||
name_(name), | name_(name), | ||||
johndan: @jback I do not include the J^P for the K-matrix, I simply removed L. I think it's reasonable… | |||||
Done Inline ActionsOK, sounds reasonable. jback: OK, sounds reasonable. | |||||
paramFileName_(paramFile), | paramFileName_(paramFile), | ||||
resPairAmpInt_(resPairAmpInt), | resPairAmpInt_(resPairAmpInt), | ||||
index_(rowIndex - 1), | index_(rowIndex - 1), | ||||
previousS_(0.0), | previousS_(0.0), | ||||
scattSVP_(0.0), | scattSVP_(0.0), | ||||
prodSVP_(0.0), | prodSVP_(0.0), | ||||
nChannels_(nChannels), | nChannels_(nChannels), | ||||
nPoles_(nPoles), | nPoles_(nPoles), | ||||
sAConst_(0.0), | sAConst_(0.0), | ||||
m2piSq_(4.0*LauConstants::mPiSq), | m2piSq_(4.0*LauConstants::mPiSq), | ||||
m2KSq_( 4.0*LauConstants::mKSq), | m2KSq_( 4.0*LauConstants::mKSq), | ||||
m2EtaSq_(4.0*LauConstants::mEtaSq), | m2EtaSq_(4.0*LauConstants::mEtaSq), | ||||
mEtaEtaPSumSq_((LauConstants::mEta + LauConstants::mEtaPrime)*(LauConstants::mEta + LauConstants::mEtaPrime)), | mEtaEtaPSumSq_((LauConstants::mEta + LauConstants::mEtaPrime)*(LauConstants::mEta + LauConstants::mEtaPrime)), | ||||
mEtaEtaPDiffSq_((LauConstants::mEta - LauConstants::mEtaPrime)*(LauConstants::mEta - LauConstants::mEtaPrime)), | mEtaEtaPDiffSq_((LauConstants::mEta - LauConstants::mEtaPrime)*(LauConstants::mEta - LauConstants::mEtaPrime)), | ||||
mKpiSumSq_((LauConstants::mK + LauConstants::mPi)*(LauConstants::mK + LauConstants::mPi)), | mKpiSumSq_((LauConstants::mK + LauConstants::mPi)*(LauConstants::mK + LauConstants::mPi)), | ||||
mKpiDiffSq_((LauConstants::mK - LauConstants::mPi)*(LauConstants::mK - LauConstants::mPi)), | mKpiDiffSq_((LauConstants::mK - LauConstants::mPi)*(LauConstants::mK - LauConstants::mPi)), | ||||
mKEtaPSumSq_((LauConstants::mK + LauConstants::mEtaPrime)*(LauConstants::mK + LauConstants::mEtaPrime)), | mKEtaPSumSq_((LauConstants::mK + LauConstants::mEtaPrime)*(LauConstants::mK + LauConstants::mEtaPrime)), | ||||
mKEtaPDiffSq_((LauConstants::mK - LauConstants::mEtaPrime)*(LauConstants::mK - LauConstants::mEtaPrime)), | mKEtaPDiffSq_((LauConstants::mK - LauConstants::mEtaPrime)*(LauConstants::mK - LauConstants::mEtaPrime)), | ||||
mK3piDiffSq_((LauConstants::mK - 3.0*LauConstants::mPi)*(LauConstants::mK - 3.0*LauConstants::mPi)), | mK3piDiffSq_((LauConstants::mK - 3.0*LauConstants::mPi)*(LauConstants::mK - 3.0*LauConstants::mPi)), | ||||
k3piFactor_(TMath::Power((1.44 - mK3piDiffSq_)/1.44, -2.5)), | k3piFactor_(TMath::Power((1.44 - mK3piDiffSq_)/1.44, -2.5)), | ||||
fourPiFactor1_(16.0*LauConstants::mPiSq), | fourPiFactor1_(16.0*LauConstants::mPiSq), | ||||
fourPiFactor2_(TMath::Sqrt(1.0 - fourPiFactor1_)), | 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)), | |||||
Done Inline ActionsYou can use something like: LauResonanceMaker::get().getResInfo("D*0")->getMass()->value() Not very pretty but it should work I think. tlatham: You can use something like:
```
LauResonanceMaker::get().getResInfo("D*0")->getMass()->value()… | |||||
adlerZeroFactor_(0.0), | adlerZeroFactor_(0.0), | ||||
parametersSet_(kFALSE), | parametersSet_(kFALSE), | ||||
verbose_(kFALSE), | verbose_(kFALSE), | ||||
scattSymmetry_(kTRUE) | scattSymmetry_(kTRUE) | ||||
{ | { | ||||
// Constructor | // Constructor | ||||
// Check that the index is OK | // Check that the index is OK | ||||
if (index_ < 0 || index_ >= nChannels_) { | if (index_ < 0 || index_ >= nChannels_) { | ||||
std::cerr << "ERROR in LauKMatrixPropagator constructor. The rowIndex, which is set to " | std::cerr << "ERROR in LauKMatrixPropagator constructor. The rowIndex, which is set to " | ||||
<< rowIndex << ", must be between 1 and the number of channels " | << rowIndex << ", must be between 1 and the number of channels " | ||||
<< nChannels_ << std::endl; | << nChannels_ << std::endl; | ||||
gSystem->Exit(EXIT_FAILURE); | gSystem->Exit(EXIT_FAILURE); | ||||
} | } | ||||
this->setParameters(paramFile); | this->setParameters(paramFile); | ||||
} | } | ||||
Done Inline ActionsI don't think that this loop does anything because the values of a_ should already have been set and then used to set gamAInvRadSq_ in setParameters above. If they were used anywhere else, this would be bad because they're being overwritten. But actually they don't seem to be used anywhere else, I'm not sure that these even need to be a member variable, they could be local to setParameters. tlatham: I don't think that this loop does anything because the values of a_ should already have been… | |||||
LauKMatrixPropagator::~LauKMatrixPropagator() | LauKMatrixPropagator::~LauKMatrixPropagator() | ||||
{ | { | ||||
// Destructor | // Destructor | ||||
realProp_.Clear(); | realProp_.Clear(); | ||||
negImagProp_.Clear(); | negImagProp_.Clear(); | ||||
Done Inline Actions@jback I was thinking about this in passing: this will give a L=2 barrier-factor-squared denominator of the form q^4 + 6q^2/R^2 + 9/R^4 which is, I think, not what you should have here. I think we are wanting the second term to be 3q^2R^2 aren't we? johndan: @jback I was thinking about this in passing: this will give a L=2 barrier-factor-squared… | |||||
Done Inline ActionsNo, the denominator factor uses L/2 for the power not L, so we get q^2 + 3/R^2 as expected. jback: No, the denominator factor uses L/2 for the power not L, so we get q^2 + 3/R^2 as expected. | |||||
Done Inline ActionsHmmm not sure. If a = 3 then calcGamma gives: i.e. the return is: L = 2 : sqrt[ (qR)^2L / ( (qR)^2 + 3 )^L ] But looking at equation 98 in the Laura++ paper I would expect the coefficient of (qR)^2 in the denominator to be 3 not 6. What am I missing? johndan: Hmmm not sure. If a = 3 then `calcGamma` gives:
pow(q,L) / pow( q*q + gamAInvRadSq_[iCh] , L_… | |||||
Done Inline ActionsYou are not missing anything; it's just that we can't exactly reproduce the z^4 + 3z^2 + 9 factor using a perfect square product between two D terms without introducing complex terms. Remember, this barrier term is an ad-hoc "best guess" at what the momentum dependence should be. The important thing is that we have a quadratic dependence on z^2, and we can change "a" and "R" at initialisation for systematic variations. jback: You are not missing anything; it's just that we can't exactly reproduce the z^4 + 3z^2 + 9… | |||||
Done Inline ActionsOK! Agreed, the dependence is what matters. johndan: OK! Agreed, the dependence is what matters. | |||||
ScattKMatrix_.Clear(); | ScattKMatrix_.Clear(); | ||||
ReRhoMatrix_.Clear(); | ReRhoMatrix_.Clear(); | ||||
ImRhoMatrix_.Clear(); | ImRhoMatrix_.Clear(); | ||||
GammaMatrix_.Clear(); | |||||
ReTMatrix_.Clear(); | ReTMatrix_.Clear(); | ||||
ImTMatrix_.Clear(); | ImTMatrix_.Clear(); | ||||
IMatrix_.Clear(); | IMatrix_.Clear(); | ||||
zeroMatrix_.Clear(); | zeroMatrix_.Clear(); | ||||
} | } | ||||
LauComplex LauKMatrixPropagator::getPropTerm(Int_t channel) const | LauComplex LauKMatrixPropagator::getPropTerm(Int_t channel) const | ||||
{ | { | ||||
▲ Show 20 Lines • Show All 41 Lines • ▼ Show 20 Lines | void LauKMatrixPropagator::updatePropagator(const LauKinematics* kinematics) | ||||
if (!kinematics) {return;} | if (!kinematics) {return;} | ||||
// Get the invariant mass squared (s) from the kinematics object. | // Get the invariant mass squared (s) from the kinematics object. | ||||
// Use the resPairAmpInt to find which mass-squared combination to use. | // Use the resPairAmpInt to find which mass-squared combination to use. | ||||
Double_t s(0.0); | Double_t s(0.0); | ||||
if (resPairAmpInt_ == 1) { | if (resPairAmpInt_ == 1) { | ||||
s = kinematics->getm23Sq(); | s = kinematics->getm23Sq(); | ||||
cosHel_ = kinematics->getc23(); | |||||
} else if (resPairAmpInt_ == 2) { | } else if (resPairAmpInt_ == 2) { | ||||
s = kinematics->getm13Sq(); | s = kinematics->getm13Sq(); | ||||
cosHel_ = kinematics->getc13(); | |||||
} else if (resPairAmpInt_ == 3) { | } else if (resPairAmpInt_ == 3) { | ||||
s = kinematics->getm12Sq(); | s = kinematics->getm12Sq(); | ||||
cosHel_ = kinematics->getc12(); | |||||
} | } | ||||
this->updatePropagator(s); | this->updatePropagator(s); | ||||
} | } | ||||
void LauKMatrixPropagator::updatePropagator(Double_t s) | void LauKMatrixPropagator::updatePropagator(Double_t s) | ||||
{ | { | ||||
Show All 20 Lines | void LauKMatrixPropagator::updatePropagator(Double_t s) | ||||
this->updateAdlerZeroFactor(s); | this->updateAdlerZeroFactor(s); | ||||
// Calculate the scattering K-matrix (real and symmetric) | // Calculate the scattering K-matrix (real and symmetric) | ||||
this->calcScattKMatrix(s); | this->calcScattKMatrix(s); | ||||
// Calculate the phase space density matrix, which is diagonal, but can be complex | // Calculate the phase space density matrix, which is diagonal, but can be complex | ||||
// if the quantity s is below various threshold values (analytic continuation). | // if the quantity s is below various threshold values (analytic continuation). | ||||
this->calcRhoMatrix(s); | this->calcRhoMatrix(s); | ||||
// Calculate K*rho (real and imaginary parts, since rho can be complex) | // Calculate the angular momentum barrier matrix, which is diagonal | ||||
TMatrixD K_realRho(ScattKMatrix_); | this->calcGammaMatrix(s); | ||||
K_realRho *= ReRhoMatrix_; | |||||
TMatrixD K_imagRho(ScattKMatrix_); | // Calculate the Xspin matrix, which is diagonal | ||||
K_imagRho *= ImRhoMatrix_; | this->calcXspinMatrix(); | ||||
// Calculate K*rho*(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_); | |||||
Done Inline ActionsWe should reuse the gamma matrix product a few lines earlier (set it as another matrix), since we need to reduce matrix multiplications as much as possible! jback: We should reuse the gamma matrix product a few lines earlier (set it as another matrix), since… | |||||
// A = I + K*Imag(rho), B = -K*Real(Rho) | // A = I + K*Imag(rho), B = -K*Real(Rho) | ||||
// Calculate C and D matrices such that (A + iB)*(C + iD) = I, | // Calculate C and D matrices such that (A + iB)*(C + iD) = I, | ||||
// ie. C + iD = (I - i K*rho)^-1, separated into real and imaginary parts. | // ie. C + iD = (I - i K*rho)^-1, separated into real and imaginary parts. | ||||
// realProp C = (A + B A^-1 B)^-1, imagProp D = -A^-1 B C | // realProp C = (A + B A^-1 B)^-1, imagProp D = -A^-1 B C | ||||
TMatrixD A(IMatrix_); | TMatrixD A(IMatrix_); | ||||
A += K_imagRho; | A += K_imagRhoGammasq; | ||||
TMatrixD B(zeroMatrix_); | TMatrixD B(zeroMatrix_); | ||||
B -= K_realRho; | B -= K_realRhoGammasq; | ||||
TMatrixD invA(TMatrixD::kInverted, A); | TMatrixD invA(TMatrixD::kInverted, A); | ||||
TMatrixD invA_B(invA); | TMatrixD invA_B(invA); | ||||
invA_B *= B; | invA_B *= B; | ||||
TMatrixD B_invA_B(B); | TMatrixD B_invA_B(B); | ||||
B_invA_B *= invA_B; | B_invA_B *= invA_B; | ||||
TMatrixD invC(A); | TMatrixD invC(A); | ||||
invC += B_invA_B; | invC += B_invA_B; | ||||
// Set the real part of the propagator matrix ("C") | // Set the real part of the propagator matrix ("C") | ||||
realProp_ = TMatrixD(TMatrixD::kInverted, invC); | realProp_ = TMatrixD(TMatrixD::kInverted, invC); | ||||
// Set the (negative) imaginary part of the propagator matrix ("-D") | // Set the (negative) imaginary part of the propagator matrix ("-D") | ||||
TMatrixD BC(B); | TMatrixD BC(B); | ||||
BC *= realProp_; | BC *= realProp_; | ||||
negImagProp_ = TMatrixD(invA); | negImagProp_ = TMatrixD(invA); | ||||
negImagProp_ *= BC; | 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 | // Also calculate the production SVP term, since this uses Adler-zero parameters | ||||
// defined in the parameter file. | // defined in the parameter file. | ||||
this->updateProdSVPTerm(s); | this->updateProdSVPTerm(s); | ||||
// Finally, keep track of the value of s we just used. | // Finally, keep track of the value of s we just used. | ||||
previousS_ = s; | previousS_ = s; | ||||
} | } | ||||
Show All 19 Lines | void LauKMatrixPropagator::setParameters(const TString& inputFile) | ||||
// appropriate set of parameters. Keywords are case insensitive (treated as lower-case). | // appropriate set of parameters. Keywords are case insensitive (treated as lower-case). | ||||
// 1) Indices (nChannels) of N phase space channel types (defined in KMatrixChannels enum) | // 1) Indices (nChannels) of N phase space channel types (defined in KMatrixChannels enum) | ||||
// "Channels iChannel1 iChannel2 ... iChannelN" | // "Channels iChannel1 iChannel2 ... iChannelN" | ||||
// 2) Definition of poles: bare mass (GeV), pole index (1 to NPoles), N channel couplings g_j | // 2) Definition of poles: bare mass (GeV), pole index (1 to NPoles), N channel couplings g_j | ||||
// "Pole poleIndex mass g_1 g_2 ... g_N" | // "Pole poleIndex mass g_1 g_2 ... g_N" | ||||
// 3) Definition of scattering f_{ij} constants: scattering index (1 to N), channel values | // 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 | // "Scatt index f_{i1} f_{i2} ... f_{iN}", where i = index | ||||
// 4) Various Adler zero and scattering constants; each line is "name value". | // 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 | ||||
// By default, the scattering constants are symmetrised: f_{ji} = f_{ij}. | // By default, the scattering constants are symmetrised: f_{ji} = f_{ij}. | ||||
Done Inline ActionsNeed to mention the angular momentum and the barrier factor parameters here too. tlatham: Need to mention the angular momentum and the barrier factor parameters here too. | |||||
// To not assume this use "ScattSymmetry 0" on a line | // To not assume this use "ScattSymmetry 0" on a line | ||||
LauTextFileParser readFile(inputFile); | LauTextFileParser readFile(inputFile); | ||||
readFile.processFile(); | readFile.processFile(); | ||||
// Loop over the (non-commented) lines | // Loop over the (non-commented) lines | ||||
UInt_t nTotLines = readFile.getTotalNumLines(); | UInt_t nTotLines = readFile.getTotalNumLines(); | ||||
▲ Show 20 Lines • Show All 86 Lines • ▼ Show 20 Lines | void LauKMatrixPropagator::initialiseMatrices() | ||||
realProp_.ResizeTo(nChannels_, nChannels_); | realProp_.ResizeTo(nChannels_, nChannels_); | ||||
negImagProp_.ResizeTo(nChannels_, nChannels_); | negImagProp_.ResizeTo(nChannels_, nChannels_); | ||||
// Phase space matrices | // Phase space matrices | ||||
ReRhoMatrix_.Clear(); ImRhoMatrix_.Clear(); | ReRhoMatrix_.Clear(); ImRhoMatrix_.Clear(); | ||||
ReRhoMatrix_.ResizeTo(nChannels_, nChannels_); | ReRhoMatrix_.ResizeTo(nChannels_, nChannels_); | ||||
ImRhoMatrix_.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 | // Square-root phase space matrices | ||||
ReSqrtRhoMatrix_.Clear(); ImSqrtRhoMatrix_.Clear(); | ReSqrtRhoMatrix_.Clear(); ImSqrtRhoMatrix_.Clear(); | ||||
ReSqrtRhoMatrix_.ResizeTo(nChannels_, nChannels_); | ReSqrtRhoMatrix_.ResizeTo(nChannels_, nChannels_); | ||||
ImSqrtRhoMatrix_.ResizeTo(nChannels_, nChannels_); | ImSqrtRhoMatrix_.ResizeTo(nChannels_, nChannels_); | ||||
// T matrices | // T matrices | ||||
ReTMatrix_.Clear(); ImTMatrix_.Clear(); | ReTMatrix_.Clear(); ImTMatrix_.Clear(); | ||||
ReTMatrix_.ResizeTo(nChannels_, nChannels_); | ReTMatrix_.ResizeTo(nChannels_, nChannels_); | ||||
Done Inline ActionsThese two vectors will have size 2*nChannels_ after this. L_.clear(); L_.assign( nChannels_, 0 ); and similarly for radii_. tlatham: These two vectors will have size `2*nChannels_` after this.
I would suggest this formulation… | |||||
ImTMatrix_.ResizeTo(nChannels_, nChannels_); | ImTMatrix_.ResizeTo(nChannels_, nChannels_); | ||||
// For the coupling and scattering constants, use LauParArrays instead of TMatrices | // For the coupling and scattering constants, use LauParArrays instead of TMatrices | ||||
// so that the quantities remain LauParameters instead of just doubles. | // so that the quantities remain LauParameters instead of just doubles. | ||||
Done Inline ActionsAs discussed above, probably this shouldn't be a member variable. However, if it is decided to keep it as one then can also use assign here. tlatham: As discussed above, probably this shouldn't be a member variable. However, if it is decided to… | |||||
// Each array is an stl vector of another stl vector of LauParameters: | // Each array is an stl vector of another stl vector of LauParameters: | ||||
// std::vector< std::vector<LauParameter> >. | // std::vector< std::vector<LauParameter> >. | ||||
// Set their sizes using the number of poles and channels defined in the constructor | // Set their sizes using the number of poles and channels defined in the constructor | ||||
gCouplings_.clear(); | gCouplings_.clear(); | ||||
gCouplings_.resize(nPoles_); | gCouplings_.resize(nPoles_); | ||||
for (Int_t iPole = 0; iPole < nPoles_; iPole++) { | for (Int_t iPole = 0; iPole < nPoles_; iPole++) { | ||||
gCouplings_[iPole].resize(nChannels_); | gCouplings_[iPole].resize(nChannels_); | ||||
▲ Show 20 Lines • Show All 59 Lines • ▼ Show 20 Lines | void LauKMatrixPropagator::storePole(const std::vector<std::string>& theLine) | ||||
if (nWords == nExpect) { | if (nWords == nExpect) { | ||||
Int_t poleIndex = std::atoi(theLine[1].c_str()) - 1; | Int_t poleIndex = std::atoi(theLine[1].c_str()) - 1; | ||||
if (poleIndex >= 0 && poleIndex < nPoles_) { | if (poleIndex >= 0 && poleIndex < nPoles_) { | ||||
Double_t poleMass = std::atof(theLine[2].c_str()); | Double_t poleMass = std::atof(theLine[2].c_str()); | ||||
Double_t poleMassSq = poleMass*poleMass; | Double_t poleMassSq = poleMass*poleMass; | ||||
LauParameter mPoleParam(poleMassSq); | LauParameter mPoleParam(poleMassSq); | ||||
Done Inline ActionsLooks like this change from D39 got reverted by mistake here. tlatham: Looks like this change from D39 got reverted by mistake here. | |||||
mSqPoles_[poleIndex] = mPoleParam; | mSqPoles_[poleIndex] = mPoleParam; | ||||
cout<<"Added bare pole mass "<<poleMass<<" GeV for pole number "<<poleIndex+1<<endl; | cout<<"Added bare pole mass "<<poleMass<<" GeV for pole number "<<poleIndex+1<<endl; | ||||
for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) { | for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) { | ||||
Double_t couplingConst = std::atof(theLine[iChannel+3].c_str()); | 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; | gCouplings_[poleIndex][iChannel] = couplingParam; | ||||
cout<<"Added coupling parameter g^{"<<poleIndex+1<<"}_" | cout<<"Added coupling parameter g^{"<<poleIndex+1<<"}_" | ||||
<<iChannel+1<<" = "<<couplingConst<<endl; | <<iChannel+1<<" = "<<couplingConst<<endl; | ||||
} | } | ||||
} | } | ||||
Show All 22 Lines | void LauKMatrixPropagator::storeScattering(const std::vector<std::string>& theLine) | ||||
if (nWords == nExpect) { | if (nWords == nExpect) { | ||||
Int_t scattIndex = std::atoi(theLine[1].c_str()) - 1; | Int_t scattIndex = std::atoi(theLine[1].c_str()) - 1; | ||||
if (scattIndex >= 0 && scattIndex < nChannels_) { | if (scattIndex >= 0 && scattIndex < nChannels_) { | ||||
for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) { | for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) { | ||||
Double_t scattConst = std::atof(theLine[iChannel+2].c_str()); | Double_t scattConst = std::atof(theLine[iChannel+2].c_str()); | ||||
LauParameter scattParam(scattConst); | LauParameter scattParam(scattConst); | ||||
Done Inline ActionsAnother accidental reversion tlatham: Another accidental reversion | |||||
fScattering_[scattIndex][iChannel] = scattParam; | fScattering_[scattIndex][iChannel] = scattParam; | ||||
cout<<"Added scattering parameter f("<<scattIndex+1<<"," | cout<<"Added scattering parameter f("<<scattIndex+1<<"," | ||||
<<iChannel+1<<") = "<<scattConst<<endl; | <<iChannel+1<<") = "<<scattConst<<endl; | ||||
} | } | ||||
} | } | ||||
} else { | } else { | ||||
cerr<<"Error in LauKMatrixPropagator::storeScattering. Expecting "<<nExpect | cerr<<"Error in LauKMatrixPropagator::storeScattering. Expecting "<<nExpect | ||||
<<" numbers for scattering constants definition but found "<<nWords | <<" numbers for scattering constants definition but found "<<nWords | ||||
<<" values instead"<<endl; | <<" values instead"<<endl; | ||||
Done Inline ActionsNot sure why this wording has changed tlatham: Not sure why this wording has changed | |||||
Done Inline Actionsthat's an error, thanks. johndan: that's an error, thanks. | |||||
} | } | ||||
} | } | ||||
void LauKMatrixPropagator::storeParameter(const TString& keyword, const TString& parString) | void LauKMatrixPropagator::storeParameter(const TString& keyword, const TString& parString) | ||||
{ | { | ||||
if (!keyword.CompareTo("msq0")) { | if (!keyword.CompareTo("msq0")) { | ||||
Double_t mSq0Value = std::atof(parString.Data()); | Double_t mSq0Value = std::atof(parString.Data()); | ||||
cout<<"Adler zero constant m0Sq = "<<mSq0Value<<endl; | cout<<"Adler zero constant m0Sq = "<<mSq0Value<<endl; | ||||
mSq0_ = LauParameter("mSq0", mSq0Value); | mSq0_ = LauParameter("mSq0", mSq0Value); | ||||
} else if (!keyword.CompareTo("s0scatt")) { | } else if (!keyword.CompareTo("s0scatt")) { | ||||
Double_t s0ScattValue = std::atof(parString.Data()); | Double_t s0ScattValue = std::atof(parString.Data()); | ||||
cout<<"Adler zero constant s0Scatt = "<<s0ScattValue<<endl; | cout<<"Adler zero constant s0Scatt = "<<s0ScattValue<<endl; | ||||
s0Scatt_ = LauParameter("s0Scatt", s0ScattValue); | s0Scatt_ = LauParameter("s0Scatt", s0ScattValue); | ||||
Done Inline ActionsShould be atoi here rather than atof tlatham: Should be `atoi` here rather than `atof` | |||||
} else if (!keyword.CompareTo("s0prod")) { | } else if (!keyword.CompareTo("s0prod")) { | ||||
Double_t s0ProdValue = std::atof(parString.Data()); | Double_t s0ProdValue = std::atof(parString.Data()); | ||||
cout<<"Adler zero constant s0Prod = "<<s0ProdValue<<endl; | cout<<"Adler zero constant s0Prod = "<<s0ProdValue<<endl; | ||||
s0Prod_ = LauParameter("s0Scatt", s0ProdValue); | s0Prod_ = LauParameter("s0Scatt", s0ProdValue); | ||||
} else if (!keyword.CompareTo("sa0")) { | } else if (!keyword.CompareTo("sa0")) { | ||||
▲ Show 20 Lines • Show All 112 Lines • ▼ Show 20 Lines | Double_t LauKMatrixPropagator::getCouplingConstant(Int_t poleIndex, Int_t channelIndex) const | ||||
if (poleIndex < 0 || poleIndex >= nPoles_) {return 0.0;} | if (poleIndex < 0 || poleIndex >= nPoles_) {return 0.0;} | ||||
if (channelIndex < 0 || channelIndex >= nChannels_) {return 0.0;} | if (channelIndex < 0 || channelIndex >= nChannels_) {return 0.0;} | ||||
Double_t couplingConst = gCouplings_[poleIndex][channelIndex].unblindValue(); | Double_t couplingConst = gCouplings_[poleIndex][channelIndex].unblindValue(); | ||||
return couplingConst; | return couplingConst; | ||||
} | } | ||||
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 | Double_t LauKMatrixPropagator::getScatteringConstant(Int_t channel1Index, Int_t channel2Index) const | ||||
{ | { | ||||
if (parametersSet_ == kFALSE) {return 0.0;} | if (parametersSet_ == kFALSE) {return 0.0;} | ||||
if (channel1Index < 0 || channel1Index >= nChannels_) {return 0.0;} | if (channel1Index < 0 || channel1Index >= nChannels_) {return 0.0;} | ||||
if (channel2Index < 0 || channel2Index >= nChannels_) {return 0.0;} | if (channel2Index < 0 || channel2Index >= nChannels_) {return 0.0;} | ||||
Double_t scatteringConst = fScattering_[channel1Index][channel2Index].unblindValue(); | Double_t scatteringConst = fScattering_[channel1Index][channel2Index].unblindValue(); | ||||
return scatteringConst; | return scatteringConst; | ||||
} | } | ||||
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 | Double_t LauKMatrixPropagator::calcSVPTerm(Double_t s, Double_t s0) const | ||||
{ | { | ||||
if (parametersSet_ == kFALSE) {return 0.0;} | if (parametersSet_ == kFALSE) {return 0.0;} | ||||
// Calculate the "slowly-varying part" (SVP) | // Calculate the "slowly-varying part" (SVP) | ||||
Double_t result(0.0); | Double_t result(0.0); | ||||
Double_t deltaS = s - s0; | Double_t deltaS = s - s0; | ||||
▲ Show 20 Lines • Show All 44 Lines • ▼ Show 20 Lines | // the amplitude continues analytically). | ||||
// Initialise all entries to zero | // Initialise all entries to zero | ||||
ReRhoMatrix_.Zero(); ImRhoMatrix_.Zero(); | ReRhoMatrix_.Zero(); ImRhoMatrix_.Zero(); | ||||
LauComplex rho(0.0, 0.0); | LauComplex rho(0.0, 0.0); | ||||
Int_t phaseSpaceIndex(0); | Int_t phaseSpaceIndex(0); | ||||
for (Int_t iChannel (0); iChannel < nChannels_; ++iChannel) { | for (Int_t iChannel (0); iChannel < nChannels_; ++iChannel) { | ||||
Done Inline ActionsI think we can remove the if statement and just use: gamma = this->calcGamma(iChannel, s, L_); where L_ is the spin integer from the constructor. This code will then be general for any channel. jback: I think we can remove the if statement and just use:
gamma = this->calcGamma(iChannel, s, L_)… | |||||
Done Inline ActionsI agree that calcGamma should be general for any channel, so no need to check the phasespace index here. But shouldn't we retain the if statement in order to avoid calling "calcGamma" (and it's internal "if" statements) if the spin is zero? johndan: I agree that calcGamma should be general for any channel, so no need to check the phasespace… | |||||
Done Inline ActionsYes, OK, this will avoid unnecessary calculations for spin zero. jback: Yes, OK, this will avoid unnecessary calculations for spin zero. | |||||
phaseSpaceIndex = phaseSpaceTypes_[iChannel]; | phaseSpaceIndex = phaseSpaceTypes_[iChannel]; | ||||
if (phaseSpaceIndex == LauKMatrixPropagator::PiPi) { | if (phaseSpaceIndex == LauKMatrixPropagator::PiPi) { | ||||
rho = this->calcPiPiRho(s); | rho = this->calcPiPiRho(s); | ||||
} else if (phaseSpaceIndex == LauKMatrixPropagator::KK) { | } else if (phaseSpaceIndex == LauKMatrixPropagator::KK) { | ||||
rho = this->calcKKRho(s); | rho = this->calcKKRho(s); | ||||
} else if (phaseSpaceIndex == LauKMatrixPropagator::FourPi) { | } else if (phaseSpaceIndex == LauKMatrixPropagator::FourPi) { | ||||
rho = this->calcFourPiRho(s); | rho = this->calcFourPiRho(s); | ||||
} else if (phaseSpaceIndex == LauKMatrixPropagator::EtaEta) { | } else if (phaseSpaceIndex == LauKMatrixPropagator::EtaEta) { | ||||
rho = this->calcEtaEtaRho(s); | rho = this->calcEtaEtaRho(s); | ||||
} else if (phaseSpaceIndex == LauKMatrixPropagator::EtaEtaP) { | } else if (phaseSpaceIndex == LauKMatrixPropagator::EtaEtaP) { | ||||
rho = this->calcEtaEtaPRho(s); | rho = this->calcEtaEtaPRho(s); | ||||
} else if (phaseSpaceIndex == LauKMatrixPropagator::KPi) { | } else if (phaseSpaceIndex == LauKMatrixPropagator::KPi) { | ||||
rho = this->calcKPiRho(s); | rho = this->calcKPiRho(s); | ||||
} else if (phaseSpaceIndex == LauKMatrixPropagator::KEtaP) { | } else if (phaseSpaceIndex == LauKMatrixPropagator::KEtaP) { | ||||
rho = this->calcKEtaPRho(s); | rho = this->calcKEtaPRho(s); | ||||
} else if (phaseSpaceIndex == LauKMatrixPropagator::KThreePi) { | } else if (phaseSpaceIndex == LauKMatrixPropagator::KThreePi) { | ||||
rho = this->calcKThreePiRho(s); | rho = this->calcKThreePiRho(s); | ||||
} else if (phaseSpaceIndex == LauKMatrixPropagator::D0K_Pwave) { | |||||
Done Inline Actions"D0K" and "Dstar0K" channel enumerations jback: "D0K" and "Dstar0K" channel enumerations | |||||
rho = this->calcD0KRho(s); | |||||
} else if (phaseSpaceIndex == LauKMatrixPropagator::Dstar0K_Swave) { | |||||
rho = this->calcDstar0KRho(s); | |||||
Done Inline ActionsWe need to change this to find gamma for all channels, not just DK. I think we can do: q = sqrt(s)*mag(rho)/2 where rho is the (complex) phase space factor depending on the channel, similar to what is done in calcRhoMatrix(). We should also probably move the if statements within calcRhoMatrix() to another function like getRho(s, channel), which can also be called here to avoid code repetition. jback: We need to change this to find gamma for all channels, not just DK. I think we can do:
```
q… | |||||
Done Inline ActionsI think all bases were covered in that calcGamma(...) was not called (by calcGammaMatrix(...)) for any case other than DK because only the two DK channels had L!=0 and the default 'gamma' was trivially returned in other cases. But I see the general usefulness of what you propose, and that it nicely avoids redetermining the break-up momentum in calcGamma so I've implemented that as you proposed. Thanks! johndan: I think all bases were covered in that `calcGamma(...)` was not called (by `calcGammaMatrix(... | |||||
} | } | ||||
Done Inline ActionsMight be nicer to convert this to a switch - that then has the compile-time check that there is a case for each state of the enum. tlatham: Might be nicer to convert this to a switch - that then has the compile-time check that there is… | |||||
if (verbose_) { | if (verbose_) { | ||||
cout<<"ReRhoMatrix("<<iChannel<<", "<<iChannel<<") = "<<rho.re()<<endl; | cout<<"ReRhoMatrix("<<iChannel<<", "<<iChannel<<") = "<<rho.re()<<endl; | ||||
cout<<"ImRhoMatrix("<<iChannel<<", "<<iChannel<<") = "<<rho.im()<<endl; | cout<<"ImRhoMatrix("<<iChannel<<", "<<iChannel<<") = "<<rho.im()<<endl; | ||||
} | } | ||||
ReRhoMatrix_(iChannel, iChannel) = rho.re(); | ReRhoMatrix_(iChannel, iChannel) = rho.re(); | ||||
ImRhoMatrix_(iChannel, iChannel) = rho.im(); | ImRhoMatrix_(iChannel, iChannel) = rho.im(); | ||||
Done Inline ActionsWe could store "gamAInvRadSq_[iCh] = a_/(radii_[iCh]*radii_[iCh])" at initialisation to avoid repeating these calculations every time calcGamma is called. jback: We could store "gamAInvRadSq_[iCh] = a_/(radii_[iCh]*radii_[iCh])" at initialisation to avoid… | |||||
} | } | ||||
} | } | ||||
Done Inline ActionsWe could use TMath::Sqrt( fabs(sqrtTerm) ) to allow cases below threshold. jback: We could use TMath::Sqrt( fabs(sqrtTerm) ) to allow cases below threshold.
| |||||
void LauKMatrixPropagator::calcGammaMatrix(Double_t s) | |||||
{ | |||||
// Calculate the gamma angular momentum barrier matrix | |||||
Done Inline ActionsWe need (a/R^2) instead of (1/R) here. We should set default values of a = 0, 1 and 3 (which could be an optional parameter that can be varied) for L = 0, 1, and 2, respectively. Some analyses don't include the denominator term, so we could make this optional (perhaps via a boolean passed to the function) to allow systematic comparisons, but will be enabled by default. Don't have to do this now. jback: We need (a/R^2) instead of (1/R) here. We should set default values of a = 0, 1 and 3 (which… | |||||
Done Inline ActionsAbsolutely. By 'parameter' and 'varied', do you mean a constant (that can be set by the user, but is fixed during the fit) or a LauParameter? Concerning the denominator, I've included a function to ignore this if needed. johndan: Absolutely. By 'parameter' and 'varied', do you mean a constant (that can be set by the user… | |||||
Done Inline ActionsFixed during the fit, but can be initialised to a different constant value for systematic checks. jback: Fixed during the fit, but can be initialised to a different constant value for systematic… | |||||
// 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_) { | |||||
Done Inline ActionsCheck spin enumeration integer again. jback: Check spin enumeration integer again. | |||||
Done Inline ActionsYes. For L==2 should the function be simply "3 cos^2 - 1"? Or are other normalisation factors required? johndan: Yes. For L==2 should the function be simply "3 cos^2 - 1"? Or are other normalisation factors… | |||||
Done Inline ActionsLet us use "3cos^2 - 1". jback: Let us use "3cos^2 - 1". | |||||
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 | LauComplex LauKMatrixPropagator::calcPiPiRho(Double_t s) const | ||||
{ | { | ||||
// Calculate the pipi phase space factor | // Calculate the pipi phase space factor | ||||
LauComplex rho(0.0, 0.0); | LauComplex rho(0.0, 0.0); | ||||
if (TMath::Abs(s) < 1e-10) {return rho;} | if (TMath::Abs(s) < 1e-10) {return rho;} | ||||
Double_t sqrtTerm = (-m2piSq_/s) + 1.0; | Double_t sqrtTerm = (-m2piSq_/s) + 1.0; | ||||
if (sqrtTerm < 0.0) { | if (sqrtTerm < 0.0) { | ||||
▲ Show 20 Lines • Show All 354 Lines • Show Last 20 Lines |
@jback I do not include the J^P for the K-matrix, I simply removed L. I think it's reasonable to think that the user has to be trusted to get that right as it's equally important that they check the validity of the coupled channels in this regard and we don't have the information about the J^P of the coupled channel particles that we would need in order to check for them.