Page MenuHomeHEPForge
Differential D41 Diff 185 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 "LauConstants.hh"
#include "LauTextFileParser.hh" #include "LauTextFileParser.hh"
#include "LauResonanceMaker.hh"
#include "LauResonanceInfo.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), 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_((LauResonanceMaker::get().getResInfo("D*0")->getMass()->value() + LauConstants::mK)*(LauResonanceMaker::get().getResInfo("D*0")->getMass()->value() + LauConstants::mK)),
mDstar0KDiffSq_((LauResonanceMaker::get().getResInfo("D*0")->getMass()->value() - LauConstants::mK)*(LauResonanceMaker::get().getResInfo("D*0")->getMass()->value() - LauConstants::mK)),
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()…
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);
// 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(Int_t channel) const
{ {
Show All 24 LinesDouble_t LauKMatrixPropagator::getImagPropTerm(Int_t channel) const
// 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)
{
// 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) 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(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 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;
// A = I + K*Imag(rho), B = -K*Real(Rho) // A = I + K*Imag(rho), B = -K*Real(Rho)
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…
// 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_;
// 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 Linesvoid 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
// 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()
Show All 18 Linesvoid 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_);
// 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 24 Linesvoid LauKMatrixPropagator::storeChannels(const std::vector<std::string>& theLine)
// 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)
{ {
// 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(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_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;
} }
} }
} 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)
{ {
// 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)
{
// 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)
{
// 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)
{
// 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;
}
} 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) 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);
} 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("s0Prod", 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("sA0", sA0Value);
} else if (!keyword.CompareTo("sa")) { } else if (!keyword.CompareTo("sa")) {
Show All 9 Lines 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(Double_t s) 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 20 Lines • Show All 55 Lines • ▼ Show 20 Lines 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(Int_t poleIndex) const
{ {
if (parametersSet_ == kFALSE) {return 0.0;} if (parametersSet_ == kFALSE) {return 0.0;}
Double_t poleDenom(0.0); Double_t poleDenom(0.0);
Show All 9 LinesDouble_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;
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(Double_t s)
{ {
// Update the scattering "slowly-varying part" (SVP) // Update the scattering "slowly-varying part" (SVP)
Show All 14 Linesvoid LauKMatrixPropagator::updateAdlerZeroFactor(Double_t s)
// 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::calcGammaMatrix(const Double_t s)
{
// Calculate the gamma angular momentum barrier matrix
// for the given invariant mass squared quantity, s.
// Initialise all entries to zero
GammaMatrix_.Zero();
Double_t gamma(0.0);
Int_t phaseSpaceIndex(0);
for (Int_t iChannel (0); iChannel < nChannels_; ++iChannel) {
phaseSpaceIndex = phaseSpaceTypes_[iChannel];
if ( L_[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 Int_t phaseSpaceIndex) const
{
// Calculate the DK P-wave angular momentum barrier factor
Double_t gamma(0.0);
Double_t sqrtTerm(0);
if (phaseSpaceIndex == LauKMatrixPropagator::D0K)
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.
return gamma;
}
void LauKMatrixPropagator::calcRhoMatrix(Double_t s) void LauKMatrixPropagator::calcRhoMatrix(Double_t s)
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…
{ {
// 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);
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);
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".
} 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::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 All 19 LinesLauComplex LauKMatrixPropagator::calcKKRho(Double_t s) const
} }
return rho; return rho;
} }
LauComplex LauKMatrixPropagator::calcFourPiRho(Double_t s) const LauComplex LauKMatrixPropagator::calcFourPiRho(Double_t s) const
{ {
// 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)) );
} }
Show All 14 LinesLauComplex LauKMatrixPropagator::calcEtaEtaRho(Double_t s) const
} }
return rho; return rho;
} }
LauComplex LauKMatrixPropagator::calcEtaEtaPRho(Double_t s) const LauComplex LauKMatrixPropagator::calcEtaEtaPRho(Double_t s) const
{ {
// 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(Double_t s) const
{ {
▲ Show 20 Lines • Show All 65 Lines • ▼ Show 20 LinesBool_t LauKMatrixPropagator::checkPhaseSpaceType(Int_t phaseSpaceInt) const
} }
return passed; return passed;
} }
LauComplex LauKMatrixPropagator::getTransitionAmp(Double_t s, Int_t channel) LauComplex LauKMatrixPropagator::getTransitionAmp(Double_t s, Int_t channel)
{ {
// Get the complex (unitary) transition amplitude T for the given channel // Get the complex (unitary) transition amplitude T for the given channel
LauComplex TAmp(0.0, 0.0); LauComplex TAmp(0.0, 0.0);
channel -= 1; 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]);
TAmp.setImagPart(ImTMatrix_[index_][channel]); TAmp.setImagPart(ImTMatrix_[index_][channel]);
return TAmp; return TAmp;
} }
LauComplex LauKMatrixPropagator::getPhaseSpaceTerm(Double_t s, Int_t channel) LauComplex LauKMatrixPropagator::getPhaseSpaceTerm(Double_t s, Int_t channel)
{ {
// Get the complex (unitary) transition amplitude T for the given channel // Get the complex (unitary) transition amplitude T for the given channel
LauComplex rho(0.0, 0.0); LauComplex rho(0.0, 0.0);
channel -= 1; 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]);
rho.setImagPart(ImRhoMatrix_[channel][channel]); rho.setImagPart(ImRhoMatrix_[channel][channel]);
return rho; return rho;
} }
void LauKMatrixPropagator::getTMatrix(const LauKinematics* kinematics) { void LauKMatrixPropagator::getTMatrix(const LauKinematics* kinematics) {
// 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(Double_t s)
{ {
// 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 20 Linesvoid LauKMatrixPropagator::getTMatrix(Double_t s)
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()
{ {
// 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(Double_t s, Int_t channel) {
LauComplex THat(0.0, 0.0); LauComplex THat(0.0, 0.0);
channel -= 1; 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_);
Show All 10 Lines