Page MenuHomeHEPForge
Differential D41 Diff 198 src/LauKMatrixPropagator.cc

Changeset View

Standalone View

View Options

src/LauKMatrixPropagator.cc

/* /*
Copyright 2008 University of Warwick Copyright 2008 University of Warwick
Licensed under the Apache License, Version 2.0 (the "License"); Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License. you may not use this file except in compliance with the License.
You may obtain a copy of the License at 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 Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS, distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and See the License for the specific language governing permissions and
limitations under the License. limitations under the License.
*/ */
/* /*
Laura++ package authors: Laura++ package authors:
John Back John Back
Paul Harrison Paul Harrison
Thomas Latham Thomas Latham
*/ */
/*! \file LauKMatrixPropagator.cc /*! \file LauKMatrixPropagator.cc
\brief File containing implementation of LauKMatrixPropagator class. \brief File containing implementation of LauKMatrixPropagator class.
*/ */
#include "LauKMatrixPropagator.hh" #include "LauKMatrixPropagator.hh"
#include "LauConstants.hh"
#include "LauTextFileParser.hh" #include "LauTextFileParser.hh"
#include "LauKinematics.hh"
#include "LauComplex.hh"
#include "TMath.h" #include "TMath.h"
#include "TSystem.h" #include "TSystem.h"
#include <iostream> #include <iostream>
#include <fstream> #include <fstream>
#include <cmath> #include <cmath>
#include <cstdlib> #include <cstdlib>
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) :
johndanAuthorUnsubmitted
Done
Inline Actions

@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.

johndan: @jback I do not include the J^P for the K-matrix, I simply removed L. I think it's reasonable…
jbackUnsubmitted
Done
Inline Actions

OK, sounds reasonable.

jback: OK, sounds reasonable.
name_(name), name_(name),
paramFileName_(paramFile), paramFileName_(paramFile),
resPairAmpInt_(resPairAmpInt), resPairAmpInt_(resPairAmpInt),
index_(rowIndex - 1), index_(rowIndex - 1),
previousS_(0.0),
scattSVP_(0.0),
prodSVP_(0.0),
nChannels_(nChannels), nChannels_(nChannels),
nPoles_(nPoles), nPoles_(nPoles)
sAConst_(0.0),
m2piSq_(4.0*LauConstants::mPiSq),
m2KSq_( 4.0*LauConstants::mKSq),
m2EtaSq_(4.0*LauConstants::mEtaSq),
mEtaEtaPSumSq_((LauConstants::mEta + LauConstants::mEtaPrime)*(LauConstants::mEta + LauConstants::mEtaPrime)),
mEtaEtaPDiffSq_((LauConstants::mEta - LauConstants::mEtaPrime)*(LauConstants::mEta - LauConstants::mEtaPrime)),
mKpiSumSq_((LauConstants::mK + LauConstants::mPi)*(LauConstants::mK + LauConstants::mPi)),
mKpiDiffSq_((LauConstants::mK - LauConstants::mPi)*(LauConstants::mK - LauConstants::mPi)),
mKEtaPSumSq_((LauConstants::mK + LauConstants::mEtaPrime)*(LauConstants::mK + LauConstants::mEtaPrime)),
mKEtaPDiffSq_((LauConstants::mK - LauConstants::mEtaPrime)*(LauConstants::mK - LauConstants::mEtaPrime)),
mK3piDiffSq_((LauConstants::mK - 3.0*LauConstants::mPi)*(LauConstants::mK - 3.0*LauConstants::mPi)),
k3piFactor_(TMath::Power((1.44 - mK3piDiffSq_)/1.44, -2.5)),
fourPiFactor1_(16.0*LauConstants::mPiSq),
fourPiFactor2_(TMath::Sqrt(1.0 - fourPiFactor1_)),
adlerZeroFactor_(0.0),
parametersSet_(kFALSE),
verbose_(kFALSE),
scattSymmetry_(kTRUE)
{ {
// Constructor // Constructor
tlathamUnsubmitted
Done
Inline Actions

You 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()…
// 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);
// Set default value of spin-dependent centrifugal-barrier-factor parameter
for( Int_t iCh = 0; iCh < nChannels_; iCh++ )
johndanAuthorUnsubmitted
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…
jbackUnsubmitted
Done
Inline Actions

No, 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.
johndanAuthorUnsubmitted
Done
Inline Actions

Hmmm not sure. If a = 3 then calcGamma gives:
pow(q,L) / pow( q*q + gamAInvRadSq_[iCh] , L_[iCh]/2. );
where gamAInvRadSq_ = a/(R*R)

i.e. the return is:

L = 2 : sqrt[ (qR)^2L / ( (qR)^2 + 3 )^L ]
L = 2 : sqrt[ (qR)^4 / ( (qR)^2 + 3 )^2 ]
L = 2 : sqrt[ (qR)^4 / ( (qR)^4 + 6 (qR)^2 + 9 ) ]

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_…
jbackUnsubmitted
Done
Inline Actions

You 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…
johndanAuthorUnsubmitted
Done
Inline Actions

OK! Agreed, the dependence is what matters.

johndan: OK! Agreed, the dependence is what matters.
{
switch(L_[iCh]) {
case 0:
a_[iCh] = 0;
break;
case 1:
a_[iCh] = 1;
break;
case 2:
a_[iCh] = 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);
}
}
} }
tlathamUnsubmitted
Done
Inline Actions

I 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.
I would say that they should be all initialised to 0 at the start of that function, then there should be a loop like this one in storeOrbitalAngularMomenta but that should only be used if storeBarrierFactorParameter has not been called already - so need some sort of flag for that.

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();
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(const Int_t channel) const noexcept
{ {
// Get the (i,j) = (index_, channel) term of the propagator // Get the (i,j) = (index_, channel) term of the propagator
// matrix. This allows us not to return the full propagator matrix. // matrix. This allows us not to return the full propagator matrix.
Double_t realTerm = this->getRealPropTerm(channel); Double_t realTerm = this->getRealPropTerm(channel);
Double_t imagTerm = this->getImagPropTerm(channel); Double_t imagTerm = this->getImagPropTerm(channel);
LauComplex propTerm(realTerm, imagTerm); LauComplex propTerm(realTerm, imagTerm);
return propTerm; return propTerm;
} }
Double_t LauKMatrixPropagator::getRealPropTerm(Int_t channel) const Double_t LauKMatrixPropagator::getRealPropTerm(const Int_t channel) const noexcept
{ {
// Get the real part of the (i,j) = (index_, channel) term of the propagator // Get the real part of the (i,j) = (index_, channel) term of the propagator
// matrix. This allows us not to return the full propagator matrix. // matrix. This allows us not to return the full propagator matrix.
if (parametersSet_ == kFALSE) {return 0.0;} if (parametersSet_ == kFALSE) {return 0.0;}
Double_t propTerm = realProp_[index_][channel]; Double_t propTerm = realProp_[index_][channel];
return propTerm; return propTerm;
} }
Double_t LauKMatrixPropagator::getImagPropTerm(Int_t channel) const Double_t LauKMatrixPropagator::getImagPropTerm(const Int_t channel) const noexcept
{ {
// Get the imaginary part of the (i,j) = (index_, channel) term of the propagator // Get the imaginary part of the (i,j) = (index_, channel) term of the propagator
// matrix. This allows us not to return the full propagator matrix. // matrix. This allows us not to return the full propagator matrix.
if (parametersSet_ == kFALSE) {return 0.0;} if (parametersSet_ == kFALSE) {return 0.0;}
Double_t imagTerm = -1.0*negImagProp_[index_][channel]; Double_t imagTerm = -1.0*negImagProp_[index_][channel];
return imagTerm; return imagTerm;
} }
void LauKMatrixPropagator::updatePropagator(const LauKinematics* kinematics) void LauKMatrixPropagator::updatePropagator(const Double_t s) noexcept
{
// Calculate the K-matrix propagator for the given s value.
// The K-matrix amplitude is given by
// T_i = sum_{ij} (I - iK*rho)^-1 * P_j, where P is the production K-matrix.
// i = index for the state (e.g. S-wave index = 0).
// Here, we only find the (I - iK*rho)^-1 matrix part.
if (!kinematics) {return;}
// Get the invariant mass squared (s) from the kinematics object.
// Use the resPairAmpInt to find which mass-squared combination to use.
Double_t s(0.0);
if (resPairAmpInt_ == 1) {
s = kinematics->getm23Sq();
} else if (resPairAmpInt_ == 2) {
s = kinematics->getm13Sq();
} else if (resPairAmpInt_ == 3) {
s = kinematics->getm12Sq();
}
this->updatePropagator(s);
}
void LauKMatrixPropagator::updatePropagator(Double_t s)
{ {
// 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 // The K-matrix amplitude is given by
// T_i = sum_{ij} (I - iK*rho)^-1 * P_j, where P is the production K-matrix. // 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). // i = index for the state (e.g. S-wave index = 0).
// Here, we only find the (I - iK*rho)^-1 matrix part. // Here, we only find the (I - iK*rho)^-1 matrix part.
// Check if we have almost the same s value as before. If so, don't re-calculate // Check if we have almost the same s value as before. If so, don't re-calculate
// the propagator nor any of the pole mass summation terms. // the propagator nor any of the pole mass summation terms.
if (TMath::Abs(s - previousS_) < 1e-6*s) { if (TMath::Abs(s - previousS_) < 1e-6*s) {
Show All 11 Linesvoid LauKMatrixPropagator::updatePropagator(const Double_t s) noexcept
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 real and diagonal
TMatrixD K_realRho(ScattKMatrix_); this->calcGammaMatrix(s);
K_realRho *= ReRhoMatrix_;
TMatrixD K_imagRho(ScattKMatrix_); // Calculate K*rho*(gamma^2) (real and imaginary parts, since rho can be complex)
K_imagRho *= ImRhoMatrix_; TMatrixD GammaMatrixSq = (GammaMatrix_*GammaMatrix_);
TMatrixD K_realRhoGammaSq(ScattKMatrix_);
K_realRhoGammaSq *= ReRhoMatrix_;
K_realRhoGammaSq *= GammaMatrixSq;
TMatrixD K_imagRhoGammaSq(ScattKMatrix_);
K_imagRhoGammaSq *= ImRhoMatrix_;
K_imagRhoGammaSq *= GammaMatrixSq;
jbackUnsubmitted
Done
Inline Actions

We 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)Gamma^2, B = -K*Real(Rho)Gamma^2
// 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*rhoGamma^2)^-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_;
if(verbose_)
{
std::cout << "In LauKMatrixPropagator::updatePropagator(s). D[1-iKrhoD^2]^-1: " << std::endl;
TString realOutput("Real part:"), imagOutput("Imag part:");
for (int iChannel = 0; iChannel < nChannels_; iChannel++)
{
for (int jChannel = 0; jChannel < nChannels_; jChannel++)
{
realOutput += Form("\t%.2f",realProp_[iChannel][jChannel]);
imagOutput += Form("\t%.2f",-1*negImagProp_[iChannel][jChannel]);
}
realOutput += "\n ";
imagOutput += "\n ";
}
std::cout << realOutput << std::endl;
std::cout << imagOutput << std::endl;
}
// Also calculate the production SVP term, since this uses Adler-zero parameters // 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;
} }
void LauKMatrixPropagator::setParameters(const TString& inputFile) void LauKMatrixPropagator::setParameters(const TString& inputFile) noexcept
{ {
// Read an input file that specifies the values of the coupling constants g_i for // Read an input file that specifies the values of the coupling constants g_i for
// the given number of poles and their (bare) masses. Also provided are the f_{ab} // the given number of poles and their (bare) masses. Also provided are the f_{ab}
// slow-varying constants. The input file should also provide the Adler zero // slow-varying constants. The input file should also provide the Adler zero
// constants s_0, s_A and s_A0. // constants s_0, s_A and s_A0.
parametersSet_ = kFALSE; parametersSet_ = kFALSE;
Show All 9 Linesvoid LauKMatrixPropagator::setParameters(const TString& inputFile) noexcept
// 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
// 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}. // By default, the scattering constants are symmetrised: f_{ji} = f_{ij}.
tlathamUnsubmitted
Done
Inline Actions

Need 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();
if (nTotLines == 0) { 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); UInt_t iLine(0);
for (iLine = 1; iLine <= nTotLines; iLine++) { for (iLine = 1; iLine <= nTotLines; iLine++) {
// Get the line of strings // Get the line of strings
std::vector<std::string> theLine = readFile.getLine(iLine); std::vector<std::string> theLine = readFile.getLine(iLine);
Show All 14 Lines if (!keyword.CompareTo("channels")) {
// Pole terms // Pole terms
this->storePole(theLine); this->storePole(theLine);
} else if (!keyword.CompareTo("scatt")) { } else if (!keyword.CompareTo("scatt")) {
// Scattering terms // Scattering terms
this->storeScattering(theLine); this->storeScattering(theLine);
} else if (!keyword.CompareTo("angularmomentum")) {
// Orbital angular momentum state for each channel
this->storeOrbitalAngularMomenta(theLine);
} else if (!keyword.CompareTo("barrierfactorparameter")) {
// Value of parameter "a" in denominator of centrifugal barrier factor, gamma
this->storeBarrierFactorParameter(theLine);
} else if (!keyword.CompareTo("radii")) {
// Values of characteristic radius
this->storeRadii(theLine);
} else { } else {
// Usually Adler-zero constants // Usually Adler-zero constants
TString parString(theLine[1].c_str()); TString parString(theLine[1].c_str());
this->storeParameter(keyword, parString); this->storeParameter(keyword, parString);
} }
Show All 10 Lines for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) {
LauParameter fPar = fScattering_[iChannel][jChannel]; LauParameter fPar = fScattering_[iChannel][jChannel];
fScattering_[jChannel][iChannel] = LauParameter(fPar); fScattering_[jChannel][iChannel] = LauParameter(fPar);
} }
} }
} }
// Now that radii and barrier-factor-denominator parameters have been set, cache the value of "a/(R*R)"
for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) {
gamAInvRadSq_[iChannel] = a_[iChannel]/(radii_[iChannel]*radii_[iChannel]);
}
// All required parameters have been set // All required parameters have been set
parametersSet_ = kTRUE; parametersSet_ = kTRUE;
cout<<"Finished initialising K-matrix propagator "<<name_<<endl; cout<<"Finished initialising K-matrix propagator "<<name_<<endl;
} }
void LauKMatrixPropagator::initialiseMatrices() void LauKMatrixPropagator::initialiseMatrices() noexcept
{ {
// Identity and null matrices // Identity and null matrices
IMatrix_.Clear(); IMatrix_.Clear();
IMatrix_.ResizeTo(nChannels_, nChannels_); IMatrix_.ResizeTo(nChannels_, nChannels_);
for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) { for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) {
IMatrix_[iChannel][iChannel] = 1.0; IMatrix_[iChannel][iChannel] = 1.0;
} }
Show All 10 Linesvoid LauKMatrixPropagator::initialiseMatrices() noexcept
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_);
// Vector of orbital angular momenta for the channels (default is S-wave everywhere)
L_.clear();
L_.resize(nChannels_);
for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) {
L_.push_back(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);
}
tlathamUnsubmitted
Done
Inline Actions

These two vectors will have size 2*nChannels_ after this.
I would suggest this formulation instead:

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…
// Vector to store barrier factor parameter a
a_.clear();
a_.resize(nChannels_);
tlathamUnsubmitted
Done
Inline Actions

As 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…
// Vector to cache ratio a/R^2
gamAInvRadSq_.clear();
gamAInvRadSq_.resize(nChannels_);
// 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_);
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.
// 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();
Show All 14 Linesvoid LauKMatrixPropagator::initialiseMatrices() noexcept
phaseSpaceTypes_.clear(); phaseSpaceTypes_.clear();
phaseSpaceTypes_.resize(nChannels_); phaseSpaceTypes_.resize(nChannels_);
mSqPoles_.clear(); mSqPoles_.clear();
mSqPoles_.resize(nPoles_); mSqPoles_.resize(nPoles_);
} }
void LauKMatrixPropagator::storeChannels(const std::vector<std::string>& theLine) void LauKMatrixPropagator::storeChannels(const std::vector<std::string>& theLine) noexcept
{ {
// Get the list of channel indices to specify what phase space factors should be used // Get the list of channel indices to specify what phase space factors should be used
// e.g. pipi, Kpi, eta eta', 4pi etc.. // e.g. pipi, Kpi, eta eta', 4pi etc..
// Check that the line has nChannels_+1 strings // Check that the line has nChannels_+1 strings
Int_t nTypes = static_cast<Int_t>(theLine.size()) - 1; Int_t nTypes = static_cast<Int_t>(theLine.size()) - 1;
if (nTypes != nChannels_) { if (nTypes != nChannels_) {
cerr<<"Error in LauKMatrixPropagator::storeChannels. The input file defines " cerr<<"Error in LauKMatrixPropagator::storeChannels. The input file defines "
<<nTypes<<" channels when "<<nChannels_<<" are expected"<<endl; <<nTypes<<" channels when "<<nChannels_<<" are expected"<<endl;
return; return;
} }
for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) { for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) {
Int_t phaseSpaceInt = std::atoi(theLine[iChannel+1].c_str()); Int_t phaseSpaceInt = std::atoi(theLine[iChannel+1].c_str());
Bool_t checkChannel = this->checkPhaseSpaceType(phaseSpaceInt); Bool_t checkChannel = this->checkPhaseSpaceType(phaseSpaceInt);
if (checkChannel == kTRUE) { if (checkChannel == kTRUE) {
cout<<"Adding phase space channel index "<<phaseSpaceInt cout<<"Adding phase space channel index "<<phaseSpaceInt
<<" to K-matrix propagator "<<name_<<endl; <<" to K-matrix propagator "<<name_<<endl;
phaseSpaceTypes_[iChannel] = phaseSpaceInt; phaseSpaceTypes_[iChannel] = phaseSpaceInt;
} else { } else {
cerr<<"Phase space channel index "<<phaseSpaceInt cerr<<"Phase space channel index "<<phaseSpaceInt
<<" should be between 1 and "<<LauKMatrixPropagator::TotChannels-1<<endl; <<" should be between 1 and "<<LauKMatrixPropagator::TotChannels-1<<endl;
} }
} }
} }
void LauKMatrixPropagator::storePole(const std::vector<std::string>& theLine) void LauKMatrixPropagator::storePole(const std::vector<std::string>& theLine) noexcept
{ {
// Store the pole mass and its coupling constants for each channel. // Store the pole mass and its coupling constants for each channel.
// Each line will contain: Pole poleNumber poleMass poleCouplingsPerChannel // Each line will contain: Pole poleNumber poleMass poleCouplingsPerChannel
// Check that the line has nChannels_ + 3 strings // Check that the line has nChannels_ + 3 strings
Int_t nWords = static_cast<Int_t>(theLine.size()); Int_t nWords = static_cast<Int_t>(theLine.size());
Int_t nExpect = nChannels_ + 3; Int_t nExpect = nChannels_ + 3;
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(Form("KM_poleMassSq_%i",poleIndex),poleMassSq);
tlathamUnsubmitted
Done
Inline Actions

Looks 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_%s_gCoupling_%i_%i",name_.Data(),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;
} }
} }
} else { } else {
cerr<<"Error in LauKMatrixPropagator::storePole. Expecting "<<nExpect cerr<<"Error in LauKMatrixPropagator::storePole. Expecting "<<nExpect
<<" numbers for pole definition but found "<<nWords <<" numbers for pole definition but found "<<nWords
<<" values instead"<<endl; <<" values instead"<<endl;
} }
} }
void LauKMatrixPropagator::storeScattering(const std::vector<std::string>& theLine) void LauKMatrixPropagator::storeScattering(const std::vector<std::string>& theLine) noexcept
{ {
// Store the scattering constants (along one of the channel rows). // Store the scattering constants (along one of the channel rows).
// Each line will contain: Scatt ScattIndex ScattConstantsPerChannel // Each line will contain: Scatt ScattIndex ScattConstantsPerChannel
// Check that the line has nChannels_ + 2 strings // Check that the line has nChannels_ + 2 strings
Int_t nWords = static_cast<Int_t>(theLine.size()); Int_t nWords = static_cast<Int_t>(theLine.size());
Int_t nExpect = nChannels_ + 2; Int_t nExpect = nChannels_ + 2;
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(Form("KM_fScatteringConst_%i_%i",scattIndex,iChannel),scattConst);
tlathamUnsubmitted
Done
Inline Actions

Another 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 radii constants definition but found "<<nWords
<<" values instead"<<endl;
}
}
void LauKMatrixPropagator::storeOrbitalAngularMomenta(const std::vector<std::string>& theLine) noexcept
{
// Store the orbital angular momentum for each channel
// Each line will contain: angularmomentum OrbitalAngularMomentumPerChannel
// Check that the line has nChannels_ + 1 strings
Int_t nWords = static_cast<Int_t>(theLine.size());
Int_t nExpect = nChannels_ + 1;
if (nWords == nExpect) {
for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) {
Int_t angularMomentum = std::atof(theLine[iChannel+1].c_str());
tlathamUnsubmitted
Done
Inline Actions

Should be atoi here rather than atof

tlatham: Should be `atoi` here rather than `atof`
L_[iChannel] = angularMomentum;
cout<<"Defined K-matrix orbital angular momentum "<<angularMomentum<<" for channel "
<<iChannel<<endl;
}
} else {
cerr<<"Error in LauKMatrixPropagator::storeOrbitalAngularMomenta. Expecting "<<nExpect
<<" numbers for orbital angular momenta definition but found "<<nWords
<<" values instead"<<endl;
tlathamUnsubmitted
Done
Inline Actions

Not sure why this wording has changed

tlatham: Not sure why this wording has changed
johndanAuthorUnsubmitted
Done
Inline Actions

that's an error, thanks.

johndan: that's an error, thanks.
}
}
void LauKMatrixPropagator::storeRadii(const std::vector<std::string>& theLine) noexcept
{
// Store the characteristic radii (measured in GeV^{-1})
// Each line will contain: Radii RadiusConstantsPerChannel
// Check that the line has nChannels_ + 1 strings
Int_t nWords = static_cast<Int_t>(theLine.size());
Int_t nExpect = nChannels_ + 1;
if (nWords == nExpect) {
for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) {
Double_t radiusConst = std::atof(theLine[iChannel+1].c_str());
radii_[iChannel] = radiusConst;
cout<<"Added K-matrix radius "<<radiusConst<<" for channel "
<<iChannel<<endl;
}
} else {
cerr<<"Error in LauKMatrixPropagator::storeRadii. Expecting "<<nExpect
<<" numbers for radii definition but found "<<nWords
<<" values instead"<<endl; <<" values instead"<<endl;
} }
}
void LauKMatrixPropagator::storeBarrierFactorParameter(const std::vector<std::string>& theLine) noexcept
{
// Store the parameter of the barrier factor
// Each line will contain: barrierfactorparameter ParameterValuePerchannel
// Check that the line has nChannels_ + 1 strings
Int_t nWords = static_cast<Int_t>(theLine.size());
Int_t nExpect = nChannels_ + 1;
if (nWords == nExpect) {
for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) {
Double_t parameterValue = std::atof(theLine[iChannel+1].c_str());
a_[iChannel] = parameterValue;
cout<<"Added K-matrix barrier factor parameter value "<<parameterValue<<" for channel "
<<iChannel<<endl;
} }
void LauKMatrixPropagator::storeParameter(const TString& keyword, const TString& parString) } else {
cerr<<"Error in LauKMatrixPropagator::storeBarrierFactorParameter. Expecting "<<nExpect
<<" numbers for barrier factor parameter definition but found "<<nWords
<<" values instead"<<endl;
}
}
void LauKMatrixPropagator::storeParameter(const TString& keyword, const TString& parString) noexcept
{ {
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(Form("KM_%s_s0Scatt",name_.Data()), s0ScattValue);
} 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(Form("KM_%s_s0Prod",name_.Data()), s0ProdValue);
} else if (!keyword.CompareTo("sa0")) { } else if (!keyword.CompareTo("sa0")) {
Double_t sA0Value = std::atof(parString.Data()); Double_t sA0Value = std::atof(parString.Data());
cout<<"Adler zero constant sA0 = "<<sA0Value<<endl; cout<<"Adler zero constant sA0 = "<<sA0Value<<endl;
sA0_ = LauParameter("sA0", sA0Value); sA0_ = LauParameter(Form("KM_%s_sA0",name_.Data()), sA0Value);
} else if (!keyword.CompareTo("sa")) { } else if (!keyword.CompareTo("sa")) {
Double_t sAValue = std::atof(parString.Data()); Double_t sAValue = std::atof(parString.Data());
cout<<"Adler zero constant sA = "<<sAValue<<endl; cout<<"Adler zero constant sA = "<<sAValue<<endl;
sA_ = LauParameter("sA", sAValue); sA_ = LauParameter(Form("KM_%s_sA",name_.Data()), sAValue);
} else if (!keyword.CompareTo("scattsymmetry")) { } else if (!keyword.CompareTo("scattsymmetry")) {
Int_t flag = std::atoi(parString.Data()); Int_t flag = std::atoi(parString.Data());
if (flag == 0) { if (flag == 0) {
cout<<"Turning off scattering parameter symmetry: f_ji = f_ij will not be assumed"<<endl; cout<<"Turning off scattering parameter symmetry: f_ji = f_ij will not be assumed"<<endl;
scattSymmetry_ = kFALSE; scattSymmetry_ = kFALSE;
} }
} }
} }
void LauKMatrixPropagator::calcScattKMatrix(const Double_t s) noexcept
void LauKMatrixPropagator::calcScattKMatrix(Double_t s)
{ {
// Calculate the scattering K-matrix for the given value of s. // Calculate the scattering K-matrix for the given value of s.
// We need to obtain the complete matrix (not just one row/column) // We need to obtain the complete matrix (not just one row/column)
// to get the correct inverted (I - i K rho) terms for the propagator. // to get the correct inverted (I - i K rho) terms for the propagator.
if (verbose_) {cout<<"Within calcScattKMatrix for s = "<<s<<endl;} if (verbose_) {cout<<"Within calcScattKMatrix for s = "<<s<<endl;}
Show All 17 Lines for (jChannel = iChannel; jChannel < nChannels_; jChannel++) {
Double_t Kij(0.0); Double_t Kij(0.0);
// Calculate pole mass summation term // Calculate pole mass summation term
for (iPole = 0; iPole < nPoles_; iPole++) { for (iPole = 0; iPole < nPoles_; iPole++) {
Double_t g_i = this->getCouplingConstant(iPole, iChannel); Double_t g_i = this->getCouplingConstant(iPole, iChannel);
Double_t g_j = this->getCouplingConstant(iPole, jChannel); Double_t g_j = this->getCouplingConstant(iPole, jChannel);
std::cout << "For pole " << iPole << " and channel " << iChannel << " got g_i: " << g_i << std::endl;
std::cout << "For pole " << iPole << " and channel " << jChannel << " got g_j: " << g_j << std::endl;
std::cout << "... and poleDenomVect[iPole]: " << poleDenomVect_[iPole] << std::endl;
Kij += poleDenomVect_[iPole]*g_i*g_j; Kij += poleDenomVect_[iPole]*g_i*g_j;
if (verbose_) {cout<<"1: Kij for i = "<<iChannel<<", j = "<<jChannel<<" = "<<Kij<<endl;} if (verbose_) {cout<<"1: Kij for i = "<<iChannel<<", j = "<<jChannel<<" = "<<Kij<<endl;}
} }
Double_t fij = this->getScatteringConstant(iChannel, jChannel); Double_t fij = this->getScatteringConstant(iChannel, jChannel);
Kij += fij*scattSVP_; Kij += fij*scattSVP_;
Kij *= adlerZeroFactor_; Kij *= adlerZeroFactor_;
if (verbose_) {cout<<"2: Kij for i = "<<iChannel<<", j = "<<jChannel<<" = "<<Kij<<endl;} if (verbose_) {cout<<"2: Kij for i = "<<iChannel<<", j = "<<jChannel<<" = "<<Kij<<endl;}
// Assign the TMatrix (i,j) element to the variable Kij and Kji (symmetry) // Assign the TMatrix (i,j) element to the variable Kij and Kji (symmetry)
ScattKMatrix_(iChannel, jChannel) = Kij; ScattKMatrix_(iChannel, jChannel) = Kij;
ScattKMatrix_(jChannel, iChannel) = Kij; ScattKMatrix_(jChannel, iChannel) = Kij;
} // j loop } // j loop
} // i loop } // i loop
} }
void LauKMatrixPropagator::calcPoleDenomVect(Double_t s) void LauKMatrixPropagator::calcPoleDenomVect(const Double_t s) noexcept
{ {
// Calculate the 1/(m_pole^2 - s) terms for the scattering // Calculate the 1/(m_pole^2 - s) terms for the scattering
// and production K-matrix formulae. // and production K-matrix formulae.
poleDenomVect_.clear(); poleDenomVect_.clear();
Int_t iPole(0); Int_t iPole(0);
for (iPole = 0; iPole < nPoles_; iPole++) { for (iPole = 0; iPole < nPoles_; iPole++) {
Double_t poleTerm = mSqPoles_[iPole].unblindValue() - s; Double_t poleTerm = mSqPoles_[iPole].unblindValue() - s;
Double_t invPoleTerm(0.0); Double_t invPoleTerm(0.0);
if (TMath::Abs(poleTerm) > 1.0e-6) {invPoleTerm = 1.0/poleTerm;} if (TMath::Abs(poleTerm) > 1.0e-6) {invPoleTerm = 1.0/poleTerm;}
poleDenomVect_.push_back(invPoleTerm); poleDenomVect_.push_back(invPoleTerm);
} }
} }
Double_t LauKMatrixPropagator::getPoleDenomTerm(Int_t poleIndex) const Double_t LauKMatrixPropagator::getPoleDenomTerm(const Int_t poleIndex) const noexcept
{ {
if (parametersSet_ == kFALSE) {return 0.0;} if (parametersSet_ == kFALSE) {return 0.0;}
Double_t poleDenom(0.0); Double_t poleDenom(0.0);
poleDenom = poleDenomVect_[poleIndex]; poleDenom = poleDenomVect_[poleIndex];
return poleDenom; return poleDenom;
} }
Double_t LauKMatrixPropagator::getCouplingConstant(Int_t poleIndex, Int_t channelIndex) const LauParameter& LauKMatrixPropagator::getPoleMassSqParameter(const Int_t poleIndex) noexcept
{
if ( (parametersSet_ == kFALSE) || (poleIndex < 0 || poleIndex >= nPoles_) ) {
std::cerr << "ERROR from LauKMatrixPropagator::getPoleMassSqParameter(). Invalid pole." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
return mSqPoles_[poleIndex];
}
Double_t LauKMatrixPropagator::getCouplingConstant(const Int_t poleIndex, const Int_t channelIndex) const noexcept
{ {
if (parametersSet_ == kFALSE) {return 0.0;} if (parametersSet_ == kFALSE) {return 0.0;}
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(const Int_t poleIndex, const Int_t channelIndex) noexcept
{
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);
}
Double_t LauKMatrixPropagator::getScatteringConstant(Int_t channel1Index, Int_t channel2Index) const //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(const Int_t channel1Index, const Int_t channel2Index) const noexcept
{ {
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;
} }
Double_t LauKMatrixPropagator::calcSVPTerm(Double_t s, Double_t s0) const LauParameter& LauKMatrixPropagator::getScatteringParameter(const Int_t channel1Index, const Int_t channel2Index) noexcept
{
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(const Double_t s, const Double_t s0) const noexcept
{ {
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;
if (TMath::Abs(deltaS) > 1.0e-6) { if (TMath::Abs(deltaS) > 1.0e-6) {
result = (mSq0_.unblindValue() - s0)/deltaS; result = (mSq0_.unblindValue() - s0)/deltaS;
} }
return result; return result;
} }
void LauKMatrixPropagator::updateScattSVPTerm(Double_t s) void LauKMatrixPropagator::updateScattSVPTerm(const Double_t s) noexcept
{ {
// Update the scattering "slowly-varying part" (SVP) // Update the scattering "slowly-varying part" (SVP)
Double_t s0Scatt = s0Scatt_.unblindValue(); Double_t s0Scatt = s0Scatt_.unblindValue();
scattSVP_ = this->calcSVPTerm(s, s0Scatt); scattSVP_ = this->calcSVPTerm(s, s0Scatt);
} }
void LauKMatrixPropagator::updateProdSVPTerm(Double_t s) void LauKMatrixPropagator::updateProdSVPTerm(const Double_t s) noexcept
{ {
// Update the production "slowly-varying part" (SVP) // Update the production "slowly-varying part" (SVP)
Double_t s0Prod = s0Prod_.unblindValue(); Double_t s0Prod = s0Prod_.unblindValue();
prodSVP_ = this->calcSVPTerm(s, s0Prod); prodSVP_ = this->calcSVPTerm(s, s0Prod);
} }
void LauKMatrixPropagator::updateAdlerZeroFactor(Double_t s) void LauKMatrixPropagator::updateAdlerZeroFactor(const Double_t s) noexcept
{ {
// Calculate the multiplicative factor containing various Adler zero // Calculate the multiplicative factor containing various Adler zero
// constants. // constants.
adlerZeroFactor_ = 0.0; adlerZeroFactor_ = 0.0;
Double_t sA0Val = sA0_.unblindValue(); Double_t sA0Val = sA0_.unblindValue();
Double_t deltaS = s - sA0Val; Double_t deltaS = s - sA0Val;
if (TMath::Abs(deltaS) > 1e-6) { if (TMath::Abs(deltaS) > 1e-6) {
adlerZeroFactor_ = (s - sAConst_)*(1.0 - sA0Val)/deltaS; adlerZeroFactor_ = (s - sAConst_)*(1.0 - sA0Val)/deltaS;
} }
} }
void LauKMatrixPropagator::calcRhoMatrix(Double_t s) void LauKMatrixPropagator::calcGammaMatrix(const Double_t s) noexcept
{
// 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_[iChannel] != 0 ) {
jbackUnsubmitted
Done
Inline Actions

I 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_)…
johndanAuthorUnsubmitted
Done
Inline Actions

I 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…
jbackUnsubmitted
Done
Inline Actions

Yes, OK, this will avoid unnecessary calculations for spin zero.

jback: Yes, OK, this will avoid unnecessary calculations for spin zero.
gamma = this->calcGamma(iChannel,s,phaseSpaceIndex);
} else {
gamma = 1.0; // S-wave
}
if (verbose_) {
cout<<"GammaMatrix("<<iChannel<<", "<<iChannel<<") = "<<gamma<<endl;
}
GammaMatrix_(iChannel, iChannel) = gamma;
}
}
Double_t LauKMatrixPropagator::calcGamma(const Int_t iCh, const Double_t s, const LauKMatrixPropagator::KMatrixChannels phaseSpaceIndex) const noexcept
{
// Calculate the DK P-wave angular momentum barrier factor
Double_t gamma(0.0);
Double_t sqrtTerm(0);
if (phaseSpaceIndex == LauKMatrixPropagator::D0K)
jbackUnsubmitted
Done
Inline Actions

We 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…
johndanAuthorUnsubmitted
Done
Inline Actions

I 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(...
sqrtTerm = (s - mD0KSumSq_) * (s - mD0KDiffSq_) / (4*s);
else if (phaseSpaceIndex == LauKMatrixPropagator::Dstar0K)
sqrtTerm = (s - mDstar0KSumSq_) * (s - mDstar0KDiffSq_) / (4*s);
Double_t q(0.0);
q = TMath::Sqrt( fabs(sqrtTerm) );
gamma = pow(q,L_[iCh]);
if (includeBWBarrierFactor_)
{
jbackUnsubmitted
Done
Inline Actions

We 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…
gamma /= pow( q*q + gamAInvRadSq_[iCh] , L_[iCh]/2. );
}
jbackUnsubmitted
Done
Inline Actions

We could use TMath::Sqrt( fabs(sqrtTerm) ) to allow cases below threshold.

jback: We could use TMath::Sqrt( fabs(sqrtTerm) ) to allow cases below threshold.
if(verbose_)
{
std::cout << "In LauKMatrixPropagator::calcGamma(iCh=" << iCh << ", s=" << s << ", prop). ";
std::cout << "|q(iCh="<<iCh<<")|: " << q << std::endl;
jbackUnsubmitted
Done
Inline Actions

We 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…
johndanAuthorUnsubmitted
Done
Inline Actions

Absolutely. 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…
jbackUnsubmitted
Done
Inline Actions

Fixed 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…
}
return gamma;
}
void LauKMatrixPropagator::calcRhoMatrix(const Double_t s) noexcept
{ {
// Calculate the real and imaginary part of the phase space density // Calculate the real and imaginary part of the phase space density
// diagonal matrix for the given invariant mass squared quantity, s. // diagonal matrix for the given invariant mass squared quantity, s.
// The matrix can be complex if s is below threshold (so that // The matrix can be complex if s is below threshold (so that
// the amplitude continues analytically). // 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);
jbackUnsubmitted
Done
Inline Actions

Check spin enumeration integer again.

jback: Check spin enumeration integer again.
johndanAuthorUnsubmitted
Done
Inline Actions

Yes. 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…
jbackUnsubmitted
Done
Inline Actions

Let us use "3cos^2 - 1".

jback: Let us use "3cos^2 - 1".
for (Int_t iChannel (0); iChannel < nChannels_; ++iChannel) { for (Int_t iChannel (0); iChannel < nChannels_; ++iChannel) {
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) {
jbackUnsubmitted
Done
Inline Actions

"D0K" and "Dstar0K" channel enumerations

jback: "D0K" and "Dstar0K" channel enumerations
rho = this->calcD0KRho(s);
} else if (phaseSpaceIndex == LauKMatrixPropagator::Dstar0K) {
rho = this->calcDstar0KRho(s);
} }
tlathamUnsubmitted
Done
Inline Actions

Might 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();
} }
} }
LauComplex LauKMatrixPropagator::calcPiPiRho(Double_t s) const LauComplex LauKMatrixPropagator::calcD0KRho(const Double_t s) const noexcept
{
// 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(const Double_t s) const noexcept
{
// 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(const Double_t s) const noexcept
{ {
// 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) {
rho.setImagPart( TMath::Sqrt(-sqrtTerm) ); rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
} else { } else {
rho.setRealPart( TMath::Sqrt(sqrtTerm) ); rho.setRealPart( TMath::Sqrt(sqrtTerm) );
} }
return rho; return rho;
} }
LauComplex LauKMatrixPropagator::calcKKRho(Double_t s) const LauComplex LauKMatrixPropagator::calcKKRho(const Double_t s) const noexcept
{ {
// Calculate the KK phase space factor // Calculate the KK 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 = (-m2KSq_/s) + 1.0; Double_t sqrtTerm = (-m2KSq_/s) + 1.0;
if (sqrtTerm < 0.0) { if (sqrtTerm < 0.0) {
rho.setImagPart( TMath::Sqrt(-sqrtTerm) ); rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
} else { } else {
rho.setRealPart( TMath::Sqrt(sqrtTerm) ); rho.setRealPart( TMath::Sqrt(sqrtTerm) );
} }
return rho; return rho;
} }
LauComplex LauKMatrixPropagator::calcFourPiRho(Double_t s) const LauComplex LauKMatrixPropagator::calcFourPiRho(const Double_t s) const noexcept
{ {
// Calculate the 4pi phase space factor. This uses a 6th-order polynomial // Calculate the 4pi phase space factor. This uses a 6th-order polynomial
// parameterisation that approximates the multi-body phase space double integral // 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 // 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 // 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. // 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 // 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, // 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 // 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 // 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 // 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 // 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 // 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). // (~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 // 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 // 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. // 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 // 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 // to a 6th-order polynomial (for s < 1), which is the result used below
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;}
if (s <= 1.0) { if (s <= 1.0) {
Double_t rhoTerm = ((1.07885*s + 0.13655)*s - 0.29744)*s - 0.20840; 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; 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 // 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. // negative number. Check for this and set the rho term to zero.
if (rhoTerm < 0.0) {rhoTerm = 0.0;} if (rhoTerm < 0.0) {rhoTerm = 0.0;}
rho.setRealPart( rhoTerm ); rho.setRealPart( rhoTerm );
} else { } else {
rho.setRealPart( TMath::Sqrt(1.0 - (fourPiFactor1_/s)) ); rho.setRealPart( TMath::Sqrt(1.0 - (fourPiFactor1_/s)) );
} }
return rho; return rho;
} }
LauComplex LauKMatrixPropagator::calcEtaEtaRho(Double_t s) const LauComplex LauKMatrixPropagator::calcEtaEtaRho(const Double_t s) const noexcept
{ {
// Calculate the eta-eta phase space factor // Calculate the eta-eta 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 = (-m2EtaSq_/s) + 1.0; Double_t sqrtTerm = (-m2EtaSq_/s) + 1.0;
if (sqrtTerm < 0.0) { if (sqrtTerm < 0.0) {
rho.setImagPart( TMath::Sqrt(-sqrtTerm) ); rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
} else { } else {
rho.setRealPart( TMath::Sqrt(sqrtTerm) ); rho.setRealPart( TMath::Sqrt(sqrtTerm) );
} }
return rho; return rho;
} }
LauComplex LauKMatrixPropagator::calcEtaEtaPRho(Double_t s) const LauComplex LauKMatrixPropagator::calcEtaEtaPRho(const Double_t s) const noexcept
{ {
// Calculate the eta-eta' phase space factor. Note that the // Calculate the eta-eta' phase space factor. Note that the
// mass difference term m_eta - m_eta' is not included, // mass difference term m_eta - m_eta' is not included,
// since this corresponds to a "t or u-channel crossing", // since this corresponds to a "t or u-channel crossing",
// which means that we cannot simply analytically continue // which means that we cannot simply analytically continue
// this part of the phase space factor below threshold, which // this part of the phase space factor below threshold, which
// we can do for s-channel contributions. This is actually an // we can do for s-channel contributions. This is actually an
// unsolved problem, e.g. see Guo et al 1409.8652, and // unsolved problem, e.g. see Guo et al 1409.8652, and
// Danilkin et al 1409.7708. Anisovich and Sarantsev in // Danilkin et al 1409.7708. Anisovich and Sarantsev in
// hep-ph/0204328 "solve" this issue by setting the mass // hep-ph/0204328 "solve" this issue by setting the mass
// difference term to unity, which is what we do here... // difference term to unity, which is what we do here...
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 = (-mEtaEtaPSumSq_/s) + 1.0; Double_t sqrtTerm = (-mEtaEtaPSumSq_/s) + 1.0;
if (sqrtTerm < 0.0) { if (sqrtTerm < 0.0) {
rho.setImagPart( TMath::Sqrt(-sqrtTerm) ); rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
} else { } else {
rho.setRealPart( TMath::Sqrt(sqrtTerm) ); rho.setRealPart( TMath::Sqrt(sqrtTerm) );
} }
return rho; return rho;
} }
LauComplex LauKMatrixPropagator::calcKPiRho(Double_t s) const LauComplex LauKMatrixPropagator::calcKPiRho(const Double_t s) const noexcept
{ {
// Calculate the K-pi phase space factor // Calculate the K-pi 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 sqrtTerm1 = (-mKpiSumSq_/s) + 1.0; Double_t sqrtTerm1 = (-mKpiSumSq_/s) + 1.0;
Double_t sqrtTerm2 = (-mKpiDiffSq_/s) + 1.0; Double_t sqrtTerm2 = (-mKpiDiffSq_/s) + 1.0;
Double_t sqrtTerm = sqrtTerm1*sqrtTerm2; Double_t sqrtTerm = sqrtTerm1*sqrtTerm2;
if (sqrtTerm < 0.0) { if (sqrtTerm < 0.0) {
rho.setImagPart( TMath::Sqrt(-sqrtTerm) ); rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
} else { } else {
rho.setRealPart( TMath::Sqrt(sqrtTerm) ); rho.setRealPart( TMath::Sqrt(sqrtTerm) );
} }
return rho; return rho;
} }
LauComplex LauKMatrixPropagator::calcKEtaPRho(Double_t s) const LauComplex LauKMatrixPropagator::calcKEtaPRho(const Double_t s) const noexcept
{ {
// Calculate the K-eta' phase space factor // Calculate the K-eta' 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 sqrtTerm1 = (-mKEtaPSumSq_/s) + 1.0; Double_t sqrtTerm1 = (-mKEtaPSumSq_/s) + 1.0;
Double_t sqrtTerm2 = (-mKEtaPDiffSq_/s) + 1.0; Double_t sqrtTerm2 = (-mKEtaPDiffSq_/s) + 1.0;
Double_t sqrtTerm = sqrtTerm1*sqrtTerm2; Double_t sqrtTerm = sqrtTerm1*sqrtTerm2;
if (sqrtTerm < 0.0) { if (sqrtTerm < 0.0) {
rho.setImagPart( TMath::Sqrt(-sqrtTerm) ); rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
} else { } else {
rho.setRealPart( TMath::Sqrt(sqrtTerm) ); rho.setRealPart( TMath::Sqrt(sqrtTerm) );
} }
return rho; return rho;
} }
LauComplex LauKMatrixPropagator::calcKThreePiRho(Double_t s) const LauComplex LauKMatrixPropagator::calcKThreePiRho(const Double_t s) const noexcept
{ {
// Calculate the Kpipipi + multimeson phase space factor. // Calculate the Kpipipi + multimeson phase space factor.
// Use the simplest definition in hep-ph/9705401 (Eq 14), which is the form // Use the simplest definition in hep-ph/9705401 (Eq 14), which is the form
// used for the rest of that paper (thankfully, the amplitude does not depend // used for the rest of that paper (thankfully, the amplitude does not depend
// significantly on the form used for the K3pi phase space factor). // significantly on the form used for the K3pi 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;}
if (s < 1.44) { if (s < 1.44) {
Double_t powerTerm = (-mK3piDiffSq_/s) + 1.0; Double_t powerTerm = (-mK3piDiffSq_/s) + 1.0;
if (powerTerm < 0.0) { if (powerTerm < 0.0) {
rho.setImagPart( k3piFactor_*TMath::Power(-powerTerm, 2.5) ); rho.setImagPart( k3piFactor_*TMath::Power(-powerTerm, 2.5) );
} else { } else {
rho.setRealPart( k3piFactor_*TMath::Power(powerTerm, 2.5) ); rho.setRealPart( k3piFactor_*TMath::Power(powerTerm, 2.5) );
} }
} else { } else {
rho.setRealPart( 1.0 ); rho.setRealPart( 1.0 );
} }
return rho; return rho;
} }
Bool_t LauKMatrixPropagator::checkPhaseSpaceType(Int_t phaseSpaceInt) const Bool_t LauKMatrixPropagator::checkPhaseSpaceType(const Int_t phaseSpaceInt) const noexcept
{ {
Bool_t passed(kFALSE); Bool_t passed(kFALSE);
if (phaseSpaceInt >= 1 && phaseSpaceInt < LauKMatrixPropagator::TotChannels) { if (phaseSpaceInt >= 1 && phaseSpaceInt < LauKMatrixPropagator::TotChannels) {
passed = kTRUE; passed = kTRUE;
} }
return passed; return passed;
} }
LauComplex LauKMatrixPropagator::getTransitionAmp(Double_t s, Int_t channel) LauComplex LauKMatrixPropagator::getTransitionAmp(const Double_t s, const Int_t channel) noexcept
{ {
// 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); LauComplex TAmp(0.0, 0.0);
channel -= 1; if (channel <= 0 || channel > nChannels_) {return TAmp;}
if (channel < 0 || channel >= nChannels_) {return TAmp;}
this->getTMatrix(s); this->getTMatrix(s);
TAmp.setRealPart(ReTMatrix_[index_][channel]); TAmp.setRealPart(ReTMatrix_[index_][channel-1]);
TAmp.setImagPart(ImTMatrix_[index_][channel]); TAmp.setImagPart(ImTMatrix_[index_][channel-1]);
return TAmp; return TAmp;
} }
LauComplex LauKMatrixPropagator::getPhaseSpaceTerm(Double_t s, Int_t channel) LauComplex LauKMatrixPropagator::getPhaseSpaceTerm(const Double_t s, const Int_t channel) noexcept
{ {
// 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); LauComplex rho(0.0, 0.0);
channel -= 1; if (channel <= 0 || channel > nChannels_) {return rho;}
if (channel < 0 || channel >= nChannels_) {return rho;}
// If s has changed from the previous value, recalculate rho // If s has changed from the previous value, recalculate rho
if (TMath::Abs(s - previousS_) > 1e-6*s) { if (TMath::Abs(s - previousS_) > 1e-6*s) {
this->calcRhoMatrix(s); this->calcRhoMatrix(s);
} }
rho.setRealPart(ReRhoMatrix_[channel][channel]); rho.setRealPart(ReRhoMatrix_[channel][channel-1]);
rho.setImagPart(ImRhoMatrix_[channel][channel]); rho.setImagPart(ImRhoMatrix_[channel][channel-1]);
return rho; return rho;
} }
void LauKMatrixPropagator::getTMatrix(const LauKinematics* kinematics) { void LauKMatrixPropagator::getTMatrix(const LauKinematics* kinematics) noexcept {
// Find the unitary T matrix, where T = [sqrt(rho)]^{*} T_hat sqrt(rho), // 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, // 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 // which has phase-space factors included (rho). This function is not
// needed to calculate the K-matrix amplitudes, but allows us // needed to calculate the K-matrix amplitudes, but allows us
// to check the variation of T as a function of s (kinematics) // to check the variation of T as a function of s (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();
} else if (resPairAmpInt_ == 2) { } else if (resPairAmpInt_ == 2) {
s = kinematics->getm13Sq(); s = kinematics->getm13Sq();
} else if (resPairAmpInt_ == 3) { } else if (resPairAmpInt_ == 3) {
s = kinematics->getm12Sq(); s = kinematics->getm12Sq();
} }
this->getTMatrix(s); this->getTMatrix(s);
} }
void LauKMatrixPropagator::getTMatrix(Double_t s) void LauKMatrixPropagator::getTMatrix(const Double_t s) noexcept
{ {
// Find the unitary transition T matrix, where // Find the unitary transition T matrix, where
// T = [sqrt(rho)]^{*} T_hat sqrt(rho), and // T = [sqrt(rho)]^{*} T_hat sqrt(rho), and
// T_hat = (I - i K rho)^-1 * K is the Lorentz-invariant T matrix, // 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 // which has phase-space factors included (rho). Note that the first
// sqrt of the rho matrix is complex conjugated. // sqrt of the rho matrix is complex conjugated.
// This function is not needed to calculate the K-matrix amplitudes, but // 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) // 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 // Initialse the real and imaginary parts of the T matrix to zero
ReTMatrix_.Zero(); ImTMatrix_.Zero(); ReTMatrix_.Zero(); ImTMatrix_.Zero();
if (parametersSet_ == kFALSE) {return;} if (parametersSet_ == kFALSE) {return;}
// Update K, rho and the propagator (I - i K rho)^-1 // Update K, rho and the propagator (I - i K rho)^-1
this->updatePropagator(s); this->updatePropagator(s);
// Find the real and imaginary T_hat matrices // Find the real and imaginary T_hat matrices
TMatrixD THatReal = realProp_*ScattKMatrix_; TMatrixD THatReal = realProp_*ScattKMatrix_;
TMatrixD THatImag(zeroMatrix_); TMatrixD THatImag(zeroMatrix_);
THatImag -= negImagProp_*ScattKMatrix_; THatImag -= negImagProp_*ScattKMatrix_;
// Find the square-root of the phase space matrix // Find the square-root of the phase space matrix
this->getSqrtRhoMatrix(); this->getSqrtRhoMatrix();
Show All 17 Linesvoid LauKMatrixPropagator::getTMatrix(const Double_t s) noexcept
DBmCA -= CA; DBmCA -= CA;
// Find the real and imaginary parts of the transition matrix T // Find the real and imaginary parts of the transition matrix T
ReTMatrix_ = ReSqrtRhoMatrix_*CAmDB + ImSqrtRhoMatrix_*DApCB; ReTMatrix_ = ReSqrtRhoMatrix_*CAmDB + ImSqrtRhoMatrix_*DApCB;
ImTMatrix_ = ReSqrtRhoMatrix_*DApCB + ImSqrtRhoMatrix_*DBmCA; ImTMatrix_ = ReSqrtRhoMatrix_*DApCB + ImSqrtRhoMatrix_*DBmCA;
} }
void LauKMatrixPropagator::getSqrtRhoMatrix() void LauKMatrixPropagator::getSqrtRhoMatrix() noexcept
{ {
// Find the square root of the (current) phase space matrix so that // 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 // 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 // Lorentz-invariant T matrix = (I - i K rho)^-1 * K; note that the first
// sqrt of rho matrix is complex conjugated // sqrt of rho matrix is complex conjugated
// If rho = rho_i + i rho_r = a + i b, then sqrt(rho) = c + i d, where // 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). // 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, // 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 // with its real and imaginary matrix elements equal to c and d, respectively
// Initialise the real and imaginary parts of the square root of // Initialise the real and imaginary parts of the square root of
// the rho matrix to zero // the rho matrix to zero
ReSqrtRhoMatrix_.Zero(); ImSqrtRhoMatrix_.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 realRho = ReRhoMatrix_[iChannel][iChannel];
Double_t imagRho = ImRhoMatrix_[iChannel][iChannel]; Double_t imagRho = ImRhoMatrix_[iChannel][iChannel];
Double_t rhoMag = sqrt(realRho*realRho + imagRho*imagRho); Double_t rhoMag = sqrt(realRho*realRho + imagRho*imagRho);
Double_t rhoSum = rhoMag + realRho; Double_t rhoSum = rhoMag + realRho;
Double_t rhoDiff = rhoMag - realRho; Double_t rhoDiff = rhoMag - realRho;
Double_t reSqrtRho(0.0), imSqrtRho(0.0); Double_t reSqrtRho(0.0), imSqrtRho(0.0);
if (rhoSum > 0.0) {reSqrtRho = sqrt(0.5*rhoSum);} if (rhoSum > 0.0) {reSqrtRho = sqrt(0.5*rhoSum);}
if (rhoDiff > 0.0) {imSqrtRho = sqrt(0.5*rhoDiff);} if (rhoDiff > 0.0) {imSqrtRho = sqrt(0.5*rhoDiff);}
ReSqrtRhoMatrix_[iChannel][iChannel] = reSqrtRho; ReSqrtRhoMatrix_[iChannel][iChannel] = reSqrtRho;
ImSqrtRhoMatrix_[iChannel][iChannel] = imSqrtRho; ImSqrtRhoMatrix_[iChannel][iChannel] = imSqrtRho;
} }
} }
LauComplex LauKMatrixPropagator::getTHat(Double_t s, Int_t channel) { LauComplex LauKMatrixPropagator::getTHat(const Double_t s, const Int_t channel) noexcept {
LauComplex THat(0.0, 0.0); LauComplex THat(0.0, 0.0);
channel -= 1;
if (channel < 0 || channel >= nChannels_) {return THat;} if (channel <= 0 || channel > nChannels_) {return THat;}
this->updatePropagator(s); this->updatePropagator(s);
// Find the real and imaginary T_hat matrices // Find the real and imaginary T_hat matrices
TMatrixD THatReal = realProp_*ScattKMatrix_; TMatrixD THatReal = realProp_*ScattKMatrix_;
TMatrixD THatImag(zeroMatrix_); TMatrixD THatImag(zeroMatrix_);
THatImag -= negImagProp_*ScattKMatrix_; THatImag -= negImagProp_*ScattKMatrix_;
// Return the specific THat component // Return the specific THat component
THat.setRealPart(THatReal[index_][channel]); THat.setRealPart(THatReal[index_][channel-1]);
THat.setImagPart(THatImag[index_][channel]); THat.setImagPart(THatImag[index_][channel-1]);
return THat; return THat;
} }