laura is hosted by Hepforge, IPPP Durham
Laura++  v3r4
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauKMatrixPropagator.cc
Go to the documentation of this file.
1 
2 /*
3 Copyright 2008 University of Warwick
4 
5 Licensed under the Apache License, Version 2.0 (the "License");
6 you may not use this file except in compliance with the License.
7 You may obtain a copy of the License at
8 
9  http://www.apache.org/licenses/LICENSE-2.0
10 
11 Unless required by applicable law or agreed to in writing, software
12 distributed under the License is distributed on an "AS IS" BASIS,
13 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 See the License for the specific language governing permissions and
15 limitations under the License.
16 */
17 
18 /*
19 Laura++ package authors:
20 John Back
21 Paul Harrison
22 Thomas Latham
23 */
24 
29 #include "LauKMatrixPropagator.hh"
30 #include "LauConstants.hh"
31 #include "LauTextFileParser.hh"
32 
33 #include "TMath.h"
34 #include "TSystem.h"
35 
36 #include <iostream>
37 #include <fstream>
38 #include <cmath>
39 #include <cstdlib>
40 
41 using std::cout;
42 using std::endl;
43 using std::cerr;
44 
46 
47 LauKMatrixPropagator::LauKMatrixPropagator(const TString& name, const TString& paramFile,
48  Int_t resPairAmpInt, Int_t nChannels,
49  Int_t nPoles, Int_t rowIndex) :
50  name_(name),
51  paramFileName_(paramFile),
52  resPairAmpInt_(resPairAmpInt),
53  index_(rowIndex - 1),
54  previousS_(0.0),
55  scattSVP_(0.0),
56  prodSVP_(0.0),
57  nChannels_(nChannels),
58  nPoles_(nPoles),
59  sAConst_(0.0),
60  m2piSq_(4.0*LauConstants::mPiSq),
61  m2KSq_( 4.0*LauConstants::mKSq),
62  m2EtaSq_(4.0*LauConstants::mEtaSq),
63  mEtaEtaPSumSq_((LauConstants::mEta + LauConstants::mEtaPrime)*(LauConstants::mEta + LauConstants::mEtaPrime)),
64  mEtaEtaPDiffSq_((LauConstants::mEta - LauConstants::mEtaPrime)*(LauConstants::mEta - LauConstants::mEtaPrime)),
65  mKpiSumSq_((LauConstants::mK + LauConstants::mPi)*(LauConstants::mK + LauConstants::mPi)),
66  mKpiDiffSq_((LauConstants::mK - LauConstants::mPi)*(LauConstants::mK - LauConstants::mPi)),
67  mKEtaPSumSq_((LauConstants::mK + LauConstants::mEtaPrime)*(LauConstants::mK + LauConstants::mEtaPrime)),
68  mKEtaPDiffSq_((LauConstants::mK - LauConstants::mEtaPrime)*(LauConstants::mK - LauConstants::mEtaPrime)),
69  mK3piDiffSq_((LauConstants::mK - 3.0*LauConstants::mPi)*(LauConstants::mK - 3.0*LauConstants::mPi)),
70  k3piFactor_(TMath::Power((1.44 - mK3piDiffSq_)/1.44, -2.5)),
71  fourPiFactor1_(16.0*LauConstants::mPiSq),
72  fourPiFactor2_(TMath::Sqrt(1.0 - fourPiFactor1_)),
73  adlerZeroFactor_(0.0),
74  parametersSet_(kFALSE),
75  verbose_(kFALSE),
76  scattSymmetry_(kTRUE)
77 {
78  // Constructor
79 
80  // Check that the index is OK
81  if (index_ < 0 || index_ >= nChannels_) {
82  std::cerr << "ERROR in LauKMatrixPropagator constructor. The rowIndex, which is set to "
83  << rowIndex << ", must be between 1 and the number of channels "
84  << nChannels_ << std::endl;
85  gSystem->Exit(EXIT_FAILURE);
86  }
87 
88  this->setParameters(paramFile);
89 }
90 
92 {
93  // Destructor
94  realProp_.Clear();
95  negImagProp_.Clear();
96  ScattKMatrix_.Clear();
97  ReRhoMatrix_.Clear();
98  ImRhoMatrix_.Clear();
99  ReTMatrix_.Clear();
100  ImTMatrix_.Clear();
101  IMatrix_.Clear();
102  zeroMatrix_.Clear();
103 }
104 
106 {
107  // Get the (i,j) = (index_, channel) term of the propagator
108  // matrix. This allows us not to return the full propagator matrix.
109  Double_t realTerm = this->getRealPropTerm(channel);
110  Double_t imagTerm = this->getImagPropTerm(channel);
111 
112  LauComplex propTerm(realTerm, imagTerm);
113  return propTerm;
114 
115 }
116 
117 Double_t LauKMatrixPropagator::getRealPropTerm(Int_t channel) const
118 {
119  // Get the real part of the (i,j) = (index_, channel) term of the propagator
120  // matrix. This allows us not to return the full propagator matrix.
121  if (parametersSet_ == kFALSE) {return 0.0;}
122 
123  Double_t propTerm = realProp_[index_][channel];
124  return propTerm;
125 
126 }
127 
128 Double_t LauKMatrixPropagator::getImagPropTerm(Int_t channel) const
129 {
130  // Get the imaginary part of the (i,j) = (index_, channel) term of the propagator
131  // matrix. This allows us not to return the full propagator matrix.
132  if (parametersSet_ == kFALSE) {return 0.0;}
133 
134  Double_t imagTerm = -1.0*negImagProp_[index_][channel];
135  return imagTerm;
136 
137 }
138 
140 {
141 
142  // Calculate the K-matrix propagator for the given s value.
143  // The K-matrix amplitude is given by
144  // T_i = sum_{ij} (I - iK*rho)^-1 * P_j, where P is the production K-matrix.
145  // i = index for the state (e.g. S-wave index = 0).
146  // Here, we only find the (I - iK*rho)^-1 matrix part.
147 
148  if (!kinematics) {return;}
149 
150  // Get the invariant mass squared (s) from the kinematics object.
151  // Use the resPairAmpInt to find which mass-squared combination to use.
152  Double_t s(0.0);
153 
154  if (resPairAmpInt_ == 1) {
155  s = kinematics->getm23Sq();
156  } else if (resPairAmpInt_ == 2) {
157  s = kinematics->getm13Sq();
158  } else if (resPairAmpInt_ == 3) {
159  s = kinematics->getm12Sq();
160  }
161 
162  this->updatePropagator(s);
163 
164 }
165 
167 {
168  // Calculate the K-matrix propagator for the given s value.
169  // The K-matrix amplitude is given by
170  // T_i = sum_{ij} (I - iK*rho)^-1 * P_j, where P is the production K-matrix.
171  // i = index for the state (e.g. S-wave index = 0).
172  // Here, we only find the (I - iK*rho)^-1 matrix part.
173 
174  // Check if we have almost the same s value as before. If so, don't re-calculate
175  // the propagator nor any of the pole mass summation terms.
176  if (TMath::Abs(s - previousS_) < 1e-6*s) {
177  //cout<<"Already got propagator for s = "<<s<<endl;
178  return;
179  }
180 
181  if (parametersSet_ == kFALSE) {
182  //cerr<<"ERROR in LauKMatrixPropagator::updatePropagator. Parameters have not been set."<<endl;
183  return;
184  }
185 
186  // Calculate the denominator pole mass terms and Adler zero factor
187  this->calcPoleDenomVect(s);
188  this->updateAdlerZeroFactor(s);
189 
190  // Calculate the scattering K-matrix (real and symmetric)
191  this->calcScattKMatrix(s);
192  // Calculate the phase space density matrix, which is diagonal, but can be complex
193  // if the quantity s is below various threshold values (analytic continuation).
194  this->calcRhoMatrix(s);
195 
196  // Calculate K*rho (real and imaginary parts, since rho can be complex)
197  TMatrixD K_realRho(ScattKMatrix_);
198  K_realRho *= ReRhoMatrix_;
199  TMatrixD K_imagRho(ScattKMatrix_);
200  K_imagRho *= ImRhoMatrix_;
201 
202  // A = I + K*Imag(rho), B = -K*Real(Rho)
203  // Calculate C and D matrices such that (A + iB)*(C + iD) = I,
204  // ie. C + iD = (I - i K*rho)^-1, separated into real and imaginary parts.
205  // realProp C = (A + B A^-1 B)^-1, imagProp D = -A^-1 B C
206  TMatrixD A(IMatrix_);
207  A += K_imagRho;
208  TMatrixD B(zeroMatrix_);
209  B -= K_realRho;
210 
211  TMatrixD invA(TMatrixD::kInverted, A);
212  TMatrixD invA_B(invA);
213  invA_B *= B;
214  TMatrixD B_invA_B(B);
215  B_invA_B *= invA_B;
216 
217  TMatrixD invC(A);
218  invC += B_invA_B;
219 
220  // Set the real part of the propagator matrix ("C")
221  realProp_ = TMatrixD(TMatrixD::kInverted, invC);
222 
223  // Set the (negative) imaginary part of the propagator matrix ("-D")
224  TMatrixD BC(B);
225  BC *= realProp_;
226  negImagProp_ = TMatrixD(invA);
227  negImagProp_ *= BC;
228 
229  // Also calculate the production SVP term, since this uses Adler-zero parameters
230  // defined in the parameter file.
231  this->updateProdSVPTerm(s);
232 
233  // Finally, keep track of the value of s we just used.
234  previousS_ = s;
235 
236 }
237 
238 void LauKMatrixPropagator::setParameters(const TString& inputFile)
239 {
240 
241  // Read an input file that specifies the values of the coupling constants g_i for
242  // the given number of poles and their (bare) masses. Also provided are the f_{ab}
243  // slow-varying constants. The input file should also provide the Adler zero
244  // constants s_0, s_A and s_A0.
245 
246  parametersSet_ = kFALSE;
247 
248  cout<<"Initialising K-matrix propagator "<<name_<<" parameters from "<<inputFile.Data()<<endl;
249 
250  cout<<"nChannels = "<<nChannels_<<", nPoles = "<<nPoles_<<endl;
251 
252  // Initialise various matrices
253  this->initialiseMatrices();
254 
255  // The format of the input file contains lines starting with a keyword followed by the
256  // appropriate set of parameters. Keywords are case insensitive (treated as lower-case).
257  // 1) Indices (nChannels) of N phase space channel types (defined in KMatrixChannels enum)
258  // "Channels iChannel1 iChannel2 ... iChannelN"
259  // 2) Definition of poles: bare mass (GeV), pole index (1 to NPoles), N channel couplings g_j
260  // "Pole poleIndex mass g_1 g_2 ... g_N"
261  // 3) Definition of scattering f_{ij} constants: scattering index (1 to N), channel values
262  // "Scatt index f_{i1} f_{i2} ... f_{iN}", where i = index
263  // 4) Various Adler zero and scattering constants; each line is "name value".
264  // Possible names are mSq0, s0Scatt, s0Prod, sA, sA0
265  // By default, the scattering constants are symmetrised: f_{ji} = f_{ij}.
266  // To not assume this use "ScattSymmetry 0" on a line
267 
268  LauTextFileParser readFile(inputFile);
269  readFile.processFile();
270 
271  // Loop over the (non-commented) lines
272  UInt_t nTotLines = readFile.getTotalNumLines();
273 
274  if (nTotLines == 0) {
275  std::cerr << "ERROR in LauKMatrixPropagator::setParameters : K-matrix parameter file not present - exiting." << std::endl;
276 
277  gSystem->Exit(EXIT_FAILURE);
278  }
279 
280  UInt_t iLine(0);
281 
282  for (iLine = 1; iLine <= nTotLines; iLine++) {
283 
284  // Get the line of strings
285  std::vector<std::string> theLine = readFile.getLine(iLine);
286 
287  // There should always be at least two strings: a keyword and at least 1 value
288  if (theLine.size() < 2) {continue;}
289 
290  TString keyword(theLine[0].c_str());
291  keyword.ToLower(); // Use lowercase
292 
293  if (!keyword.CompareTo("channels")) {
294 
295  // Channel indices for phase-space factors
296  this->storeChannels(theLine);
297 
298  } else if (!keyword.CompareTo("pole")) {
299 
300  // Pole terms
301  this->storePole(theLine);
302 
303  } else if (!keyword.CompareTo("scatt")) {
304 
305  // Scattering terms
306  this->storeScattering(theLine);
307 
308  } else {
309 
310  // Usually Adler-zero constants
311  TString parString(theLine[1].c_str());
312  this->storeParameter(keyword, parString);
313 
314  }
315 
316  }
317 
319 
320  // Symmetrise scattering parameters if enabled
321  if (scattSymmetry_ == kTRUE) {
322 
323  for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) {
324 
325  for (Int_t jChannel = iChannel; jChannel < nChannels_; jChannel++) {
326 
327  LauParameter fPar = fScattering_[iChannel][jChannel];
328  fScattering_[jChannel][iChannel] = LauParameter(fPar);
329 
330  }
331  }
332  }
333 
334  // All required parameters have been set
335  parametersSet_ = kTRUE;
336 
337  cout<<"Finished initialising K-matrix propagator "<<name_<<endl;
338 
339 }
340 
342 {
343 
344  // Identity and null matrices
345  IMatrix_.Clear();
346  IMatrix_.ResizeTo(nChannels_, nChannels_);
347  for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) {
348  IMatrix_[iChannel][iChannel] = 1.0;
349  }
350 
351  zeroMatrix_.Clear();
352  zeroMatrix_.ResizeTo(nChannels_, nChannels_);
353 
354  // Real K matrix
355  ScattKMatrix_.Clear();
356  ScattKMatrix_.ResizeTo(nChannels_, nChannels_);
357 
358  // Real and imaginary propagator matrices
359  realProp_.Clear(); negImagProp_.Clear();
360  realProp_.ResizeTo(nChannels_, nChannels_);
361  negImagProp_.ResizeTo(nChannels_, nChannels_);
362 
363  // Phase space matrices
364  ReRhoMatrix_.Clear(); ImRhoMatrix_.Clear();
365  ReRhoMatrix_.ResizeTo(nChannels_, nChannels_);
366  ImRhoMatrix_.ResizeTo(nChannels_, nChannels_);
367 
368  // Square-root phase space matrices
369  ReSqrtRhoMatrix_.Clear(); ImSqrtRhoMatrix_.Clear();
370  ReSqrtRhoMatrix_.ResizeTo(nChannels_, nChannels_);
371  ImSqrtRhoMatrix_.ResizeTo(nChannels_, nChannels_);
372 
373  // T matrices
374  ReTMatrix_.Clear(); ImTMatrix_.Clear();
375  ReTMatrix_.ResizeTo(nChannels_, nChannels_);
376  ImTMatrix_.ResizeTo(nChannels_, nChannels_);
377 
378  // For the coupling and scattering constants, use LauParArrays instead of TMatrices
379  // so that the quantities remain LauParameters instead of just doubles.
380  // Each array is an stl vector of another stl vector of LauParameters:
381  // std::vector< std::vector<LauParameter> >.
382  // Set their sizes using the number of poles and channels defined in the constructor
383  gCouplings_.clear();
384  gCouplings_.resize(nPoles_);
385 
386  for (Int_t iPole = 0; iPole < nPoles_; iPole++) {
387  gCouplings_[iPole].resize(nChannels_);
388  }
389 
390  fScattering_.clear();
391  fScattering_.resize(nChannels_);
392 
393  for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) {
394  fScattering_[iChannel].resize(nChannels_);
395  }
396 
397  // Clear other vectors
398  phaseSpaceTypes_.clear();
399  phaseSpaceTypes_.resize(nChannels_);
400 
401  mSqPoles_.clear();
402  mSqPoles_.resize(nPoles_);
403 
404 }
405 
406 void LauKMatrixPropagator::storeChannels(const std::vector<std::string>& theLine)
407 {
408 
409  // Get the list of channel indices to specify what phase space factors should be used
410  // e.g. pipi, Kpi, eta eta', 4pi etc..
411 
412  // Check that the line has nChannels_+1 strings
413  Int_t nTypes = static_cast<Int_t>(theLine.size()) - 1;
414  if (nTypes != nChannels_) {
415  cerr<<"Error in LauKMatrixPropagator::storeChannels. The input file defines "
416  <<nTypes<<" channels when "<<nChannels_<<" are expected"<<endl;
417  return;
418  }
419 
420  for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) {
421 
422  Int_t phaseSpaceInt = std::atoi(theLine[iChannel+1].c_str());
423  Bool_t checkChannel = this->checkPhaseSpaceType(phaseSpaceInt);
424 
425  if (checkChannel == kTRUE) {
426  cout<<"Adding phase space channel index "<<phaseSpaceInt
427  <<" to K-matrix propagator "<<name_<<endl;
428  phaseSpaceTypes_[iChannel] = phaseSpaceInt;
429  } else {
430  cerr<<"Phase space channel index "<<phaseSpaceInt
431  <<" should be between 1 and "<<LauKMatrixPropagator::TotChannels-1<<endl;
432  }
433 
434  }
435 
436 }
437 
438 void LauKMatrixPropagator::storePole(const std::vector<std::string>& theLine)
439 {
440 
441  // Store the pole mass and its coupling constants for each channel.
442  // Each line will contain: Pole poleNumber poleMass poleCouplingsPerChannel
443 
444  // Check that the line has nChannels_ + 3 strings
445  Int_t nWords = static_cast<Int_t>(theLine.size());
446  Int_t nExpect = nChannels_ + 3;
447 
448  if (nWords == nExpect) {
449 
450  Int_t poleIndex = std::atoi(theLine[1].c_str()) - 1;
451  if (poleIndex >= 0 && poleIndex < nPoles_) {
452 
453  Double_t poleMass = std::atof(theLine[2].c_str());
454  Double_t poleMassSq = poleMass*poleMass;
455  LauParameter mPoleParam(poleMassSq);
456  mSqPoles_[poleIndex] = mPoleParam;
457 
458  cout<<"Added bare pole mass "<<poleMass<<" GeV for pole number "<<poleIndex+1<<endl;
459 
460  for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) {
461 
462  Double_t couplingConst = std::atof(theLine[iChannel+3].c_str());
463  LauParameter couplingParam(couplingConst);
464  gCouplings_[poleIndex][iChannel] = couplingParam;
465 
466  cout<<"Added coupling parameter g^{"<<poleIndex+1<<"}_"
467  <<iChannel+1<<" = "<<couplingConst<<endl;
468 
469  }
470 
471  }
472 
473  } else {
474 
475  cerr<<"Error in LauKMatrixPropagator::storePole. Expecting "<<nExpect
476  <<" numbers for pole definition but found "<<nWords
477  <<" values instead"<<endl;
478 
479  }
480 
481 
482 }
483 
484 void LauKMatrixPropagator::storeScattering(const std::vector<std::string>& theLine)
485 {
486 
487  // Store the scattering constants (along one of the channel rows).
488  // Each line will contain: Scatt ScattIndex ScattConstantsPerChannel
489 
490  // Check that the line has nChannels_ + 2 strings
491  Int_t nWords = static_cast<Int_t>(theLine.size());
492  Int_t nExpect = nChannels_ + 2;
493 
494  if (nWords == nExpect) {
495 
496  Int_t scattIndex = std::atoi(theLine[1].c_str()) - 1;
497  if (scattIndex >= 0 && scattIndex < nChannels_) {
498 
499  for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) {
500 
501  Double_t scattConst = std::atof(theLine[iChannel+2].c_str());
502  LauParameter scattParam(scattConst);
503  fScattering_[scattIndex][iChannel] = scattParam;
504 
505  cout<<"Added scattering parameter f("<<scattIndex+1<<","
506  <<iChannel+1<<") = "<<scattConst<<endl;
507 
508  }
509 
510  }
511 
512  } else {
513 
514  cerr<<"Error in LauKMatrixPropagator::storeScattering. Expecting "<<nExpect
515  <<" numbers for scattering constants definition but found "<<nWords
516  <<" values instead"<<endl;
517 
518  }
519 
520 
521 }
522 
523 void LauKMatrixPropagator::storeParameter(const TString& keyword, const TString& parString)
524 {
525 
526  if (!keyword.CompareTo("msq0")) {
527 
528  Double_t mSq0Value = std::atof(parString.Data());
529  cout<<"Adler zero constant m0Sq = "<<mSq0Value<<endl;
530  mSq0_ = LauParameter("mSq0", mSq0Value);
531 
532  } else if (!keyword.CompareTo("s0scatt")) {
533 
534  Double_t s0ScattValue = std::atof(parString.Data());
535  cout<<"Adler zero constant s0Scatt = "<<s0ScattValue<<endl;
536  s0Scatt_ = LauParameter("s0Scatt", s0ScattValue);
537 
538  } else if (!keyword.CompareTo("s0prod")) {
539 
540  Double_t s0ProdValue = std::atof(parString.Data());
541  cout<<"Adler zero constant s0Prod = "<<s0ProdValue<<endl;
542  s0Prod_ = LauParameter("s0Scatt", s0ProdValue);
543 
544  } else if (!keyword.CompareTo("sa0")) {
545 
546  Double_t sA0Value = std::atof(parString.Data());
547  cout<<"Adler zero constant sA0 = "<<sA0Value<<endl;
548  sA0_ = LauParameter("sA0", sA0Value);
549 
550  } else if (!keyword.CompareTo("sa")) {
551 
552  Double_t sAValue = std::atof(parString.Data());
553  cout<<"Adler zero constant sA = "<<sAValue<<endl;
554  sA_ = LauParameter("sA", sAValue);
555 
556  } else if (!keyword.CompareTo("scattsymmetry")) {
557 
558  Int_t flag = std::atoi(parString.Data());
559  if (flag == 0) {
560  cout<<"Turning off scattering parameter symmetry: f_ji = f_ij will not be assumed"<<endl;
561  scattSymmetry_ = kFALSE;
562  }
563 
564  }
565 
566 }
567 
568 
570 {
571 
572  // Calculate the scattering K-matrix for the given value of s.
573  // We need to obtain the complete matrix (not just one row/column)
574  // to get the correct inverted (I - i K rho) terms for the propagator.
575 
576  if (verbose_) {cout<<"Within calcScattKMatrix for s = "<<s<<endl;}
577 
578  // Initialise the K matrix to zero
579  ScattKMatrix_.Zero();
580 
581  Int_t iChannel(0), jChannel(0), iPole(0);
582 
583  // The pole denominator 1/(m^2 - s) terms should already be calculated
584  // by the calcPoleDenomVect() function. These same terms are also
585  // used for calculating the production K-matrix elements.
586 
587  // Calculate the "slowly-varying part" (SVP), e.g. (1 GeV - s0)/(s - s0)
588  this->updateScattSVPTerm(s);
589 
590  // Now loop over iChannel, jChannel to calculate Kij = Kji.
591  for (iChannel = 0; iChannel < nChannels_; iChannel++) {
592 
593  // Scattering matrix is real and symmetric. Start j loop from i.
594  for (jChannel = iChannel; jChannel < nChannels_; jChannel++) {
595 
596  Double_t Kij(0.0);
597 
598  // Calculate pole mass summation term
599  for (iPole = 0; iPole < nPoles_; iPole++) {
600 
601  Double_t g_i = this->getCouplingConstant(iPole, iChannel);
602  Double_t g_j = this->getCouplingConstant(iPole, jChannel);
603 
604  Kij += poleDenomVect_[iPole]*g_i*g_j;
605  if (verbose_) {cout<<"1: Kij for i = "<<iChannel<<", j = "<<jChannel<<" = "<<Kij<<endl;}
606 
607  }
608 
609  Double_t fij = this->getScatteringConstant(iChannel, jChannel);
610  Kij += fij*scattSVP_;
611 
612  Kij *= adlerZeroFactor_;
613  if (verbose_) {cout<<"2: Kij for i = "<<iChannel<<", j = "<<jChannel<<" = "<<Kij<<endl;}
614 
615  // Assign the TMatrix (i,j) element to the variable Kij and Kji (symmetry)
616  ScattKMatrix_(iChannel, jChannel) = Kij;
617  ScattKMatrix_(jChannel, iChannel) = Kij;
618 
619  } // j loop
620 
621  } // i loop
622 
623 }
624 
626 {
627  // Calculate the 1/(m_pole^2 - s) terms for the scattering
628  // and production K-matrix formulae.
629  poleDenomVect_.clear();
630  Int_t iPole(0);
631  for (iPole = 0; iPole < nPoles_; iPole++) {
632 
633  Double_t poleTerm = mSqPoles_[iPole].unblindValue() - s;
634  Double_t invPoleTerm(0.0);
635  if (TMath::Abs(poleTerm) > 1.0e-6) {invPoleTerm = 1.0/poleTerm;}
636 
637  poleDenomVect_.push_back(invPoleTerm);
638 
639  }
640 
641 }
642 
643 Double_t LauKMatrixPropagator::getPoleDenomTerm(Int_t poleIndex) const
644 {
645 
646  if (parametersSet_ == kFALSE) {return 0.0;}
647 
648  Double_t poleDenom(0.0);
649  poleDenom = poleDenomVect_[poleIndex];
650  return poleDenom;
651 
652 }
653 
654 Double_t LauKMatrixPropagator::getCouplingConstant(Int_t poleIndex, Int_t channelIndex) const
655 {
656 
657  if (parametersSet_ == kFALSE) {return 0.0;}
658  if (poleIndex < 0 || poleIndex >= nPoles_) {return 0.0;}
659  if (channelIndex < 0 || channelIndex >= nChannels_) {return 0.0;}
660 
661  Double_t couplingConst = gCouplings_[poleIndex][channelIndex].unblindValue();
662  return couplingConst;
663 
664 }
665 
666 
667 Double_t LauKMatrixPropagator::getScatteringConstant(Int_t channel1Index, Int_t channel2Index) const
668 {
669 
670  if (parametersSet_ == kFALSE) {return 0.0;}
671  if (channel1Index < 0 || channel1Index >= nChannels_) {return 0.0;}
672  if (channel2Index < 0 || channel2Index >= nChannels_) {return 0.0;}
673 
674  Double_t scatteringConst = fScattering_[channel1Index][channel2Index].unblindValue();
675  return scatteringConst;
676 
677 }
678 
679 Double_t LauKMatrixPropagator::calcSVPTerm(Double_t s, Double_t s0) const
680 {
681 
682  if (parametersSet_ == kFALSE) {return 0.0;}
683 
684  // Calculate the "slowly-varying part" (SVP)
685  Double_t result(0.0);
686  Double_t deltaS = s - s0;
687  if (TMath::Abs(deltaS) > 1.0e-6) {
688  result = (mSq0_.unblindValue() - s0)/deltaS;
689  }
690 
691  return result;
692 
693 }
694 
696 {
697  // Update the scattering "slowly-varying part" (SVP)
698  Double_t s0Scatt = s0Scatt_.unblindValue();
699  scattSVP_ = this->calcSVPTerm(s, s0Scatt);
700 }
701 
703 {
704  // Update the production "slowly-varying part" (SVP)
705  Double_t s0Prod = s0Prod_.unblindValue();
706  prodSVP_ = this->calcSVPTerm(s, s0Prod);
707 }
708 
710 {
711 
712  // Calculate the multiplicative factor containing various Adler zero
713  // constants.
714  adlerZeroFactor_ = 0.0;
715 
716  Double_t sA0Val = sA0_.unblindValue();
717  Double_t deltaS = s - sA0Val;
718  if (TMath::Abs(deltaS) > 1e-6) {
719  adlerZeroFactor_ = (s - sAConst_)*(1.0 - sA0Val)/deltaS;
720  }
721 
722 }
723 
725 {
726 
727  // Calculate the real and imaginary part of the phase space density
728  // diagonal matrix for the given invariant mass squared quantity, s.
729  // The matrix can be complex if s is below threshold (so that
730  // the amplitude continues analytically).
731 
732  // Initialise all entries to zero
733  ReRhoMatrix_.Zero(); ImRhoMatrix_.Zero();
734 
735  LauComplex rho(0.0, 0.0);
736  Int_t phaseSpaceIndex(0);
737 
738  for (Int_t iChannel (0); iChannel < nChannels_; ++iChannel) {
739 
740  phaseSpaceIndex = phaseSpaceTypes_[iChannel];
741 
742  if (phaseSpaceIndex == LauKMatrixPropagator::PiPi) {
743  rho = this->calcPiPiRho(s);
744  } else if (phaseSpaceIndex == LauKMatrixPropagator::KK) {
745  rho = this->calcKKRho(s);
746  } else if (phaseSpaceIndex == LauKMatrixPropagator::FourPi) {
747  rho = this->calcFourPiRho(s);
748  } else if (phaseSpaceIndex == LauKMatrixPropagator::EtaEta) {
749  rho = this->calcEtaEtaRho(s);
750  } else if (phaseSpaceIndex == LauKMatrixPropagator::EtaEtaP) {
751  rho = this->calcEtaEtaPRho(s);
752  } else if (phaseSpaceIndex == LauKMatrixPropagator::KPi) {
753  rho = this->calcKPiRho(s);
754  } else if (phaseSpaceIndex == LauKMatrixPropagator::KEtaP) {
755  rho = this->calcKEtaPRho(s);
756  } else if (phaseSpaceIndex == LauKMatrixPropagator::KThreePi) {
757  rho = this->calcKThreePiRho(s);
758  }
759 
760  if (verbose_) {
761  cout<<"ReRhoMatrix("<<iChannel<<", "<<iChannel<<") = "<<rho.re()<<endl;
762  cout<<"ImRhoMatrix("<<iChannel<<", "<<iChannel<<") = "<<rho.im()<<endl;
763  }
764 
765  ReRhoMatrix_(iChannel, iChannel) = rho.re();
766  ImRhoMatrix_(iChannel, iChannel) = rho.im();
767 
768  }
769 
770 }
771 
773 {
774  // Calculate the pipi phase space factor
775  LauComplex rho(0.0, 0.0);
776  if (TMath::Abs(s) < 1e-10) {return rho;}
777 
778  Double_t sqrtTerm = (-m2piSq_/s) + 1.0;
779  if (sqrtTerm < 0.0) {
780  rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
781  } else {
782  rho.setRealPart( TMath::Sqrt(sqrtTerm) );
783  }
784 
785  return rho;
786 }
787 
789 {
790  // Calculate the KK phase space factor
791  LauComplex rho(0.0, 0.0);
792  if (TMath::Abs(s) < 1e-10) {return rho;}
793 
794  Double_t sqrtTerm = (-m2KSq_/s) + 1.0;
795  if (sqrtTerm < 0.0) {
796  rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
797  } else {
798  rho.setRealPart( TMath::Sqrt(sqrtTerm) );
799  }
800 
801  return rho;
802 }
803 
805 {
806  // Calculate the 4pi phase space factor. This uses a 6th-order polynomial
807  // parameterisation that approximates the multi-body phase space double integral
808  // defined in Eq 4 of the A&S paper hep-ph/0204328. This form agrees with the
809  // BaBar model (another 6th order polynomial from s^4 down to 1/s^2), but avoids the
810  // exponential increase at small values of s (~< 0.1) arising from 1/s and 1/s^2.
811  // Eq 4 is evaluated for each value of s by assuming incremental steps of 1e-3 GeV^2
812  // for s1 and s2, the invariant energy squared of each of the di-pion states,
813  // with the integration limits of s1 = (2*mpi)^2 to (sqrt(s) - 2*mpi)^2 and
814  // s2 = (2*mpi)^2 to (sqrt(s) - sqrt(s1))^2. The mass M of the rho is taken to be
815  // 0.775 GeV and the energy-dependent width of the 4pi system
816  // Gamma(s) = gamma_0*rho1^3(s), where rho1 = sqrt(1.0 - 4*mpiSq/s) and gamma_0 is
817  // the "width" of the 4pi state at s = 1, which is taken to be 0.3 GeV
818  // (~75% of the total width from PDG estimates of the f0(1370) -> 4pi state).
819  // The normalisation term rho_0 is found by ensuring that the phase space integral
820  // at s = 1 is equal to sqrt(1.0 - 16*mpiSq/s). Note that the exponent for this
821  // factor in hep-ph/0204328 is wrong; it should be 0.5, i.e. sqrt, not n = 1 to 5.
822  // Plotting the value of this double integral as a function of s can then be fitted
823  // to a 6th-order polynomial (for s < 1), which is the result used below
824 
825  LauComplex rho(0.0, 0.0);
826  if (TMath::Abs(s) < 1e-10) {return rho;}
827 
828  if (s <= 1.0) {
829  Double_t rhoTerm = ((1.07885*s + 0.13655)*s - 0.29744)*s - 0.20840;
830  rhoTerm = ((rhoTerm*s + 0.13851)*s - 0.01933)*s + 0.00051;
831  // For some values of s (below 2*mpi), this term is a very small
832  // negative number. Check for this and set the rho term to zero.
833  if (rhoTerm < 0.0) {rhoTerm = 0.0;}
834  rho.setRealPart( rhoTerm );
835  } else {
836  rho.setRealPart( TMath::Sqrt(1.0 - (fourPiFactor1_/s)) );
837  }
838 
839  return rho;
840 }
841 
843 {
844  // Calculate the eta-eta phase space factor
845  LauComplex rho(0.0, 0.0);
846  if (TMath::Abs(s) < 1e-10) {return rho;}
847 
848  Double_t sqrtTerm = (-m2EtaSq_/s) + 1.0;
849  if (sqrtTerm < 0.0) {
850  rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
851  } else {
852  rho.setRealPart( TMath::Sqrt(sqrtTerm) );
853  }
854 
855  return rho;
856 }
857 
859 {
860  // Calculate the eta-eta' phase space factor. Note that the
861  // mass difference term m_eta - m_eta' is not included,
862  // since this corresponds to a "t or u-channel crossing",
863  // which means that we cannot simply analytically continue
864  // this part of the phase space factor below threshold, which
865  // we can do for s-channel contributions. This is actually an
866  // unsolved problem, e.g. see Guo et al 1409.8652, and
867  // Danilkin et al 1409.7708. Anisovich and Sarantsev in
868  // hep-ph/0204328 "solve" this issue by setting the mass
869  // difference term to unity, which is what we do here...
870 
871  LauComplex rho(0.0, 0.0);
872  if (TMath::Abs(s) < 1e-10) {return rho;}
873 
874  Double_t sqrtTerm = (-mEtaEtaPSumSq_/s) + 1.0;
875  if (sqrtTerm < 0.0) {
876  rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
877  } else {
878  rho.setRealPart( TMath::Sqrt(sqrtTerm) );
879  }
880 
881  return rho;
882 }
883 
884 
886 {
887  // Calculate the K-pi phase space factor
888  LauComplex rho(0.0, 0.0);
889  if (TMath::Abs(s) < 1e-10) {return rho;}
890 
891  Double_t sqrtTerm1 = (-mKpiSumSq_/s) + 1.0;
892  Double_t sqrtTerm2 = (-mKpiDiffSq_/s) + 1.0;
893  Double_t sqrtTerm = sqrtTerm1*sqrtTerm2;
894  if (sqrtTerm < 0.0) {
895  rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
896  } else {
897  rho.setRealPart( TMath::Sqrt(sqrtTerm) );
898  }
899 
900  return rho;
901 }
902 
904 {
905  // Calculate the K-eta' phase space factor
906  LauComplex rho(0.0, 0.0);
907  if (TMath::Abs(s) < 1e-10) {return rho;}
908 
909  Double_t sqrtTerm1 = (-mKEtaPSumSq_/s) + 1.0;
910  Double_t sqrtTerm2 = (-mKEtaPDiffSq_/s) + 1.0;
911  Double_t sqrtTerm = sqrtTerm1*sqrtTerm2;
912  if (sqrtTerm < 0.0) {
913  rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
914  } else {
915  rho.setRealPart( TMath::Sqrt(sqrtTerm) );
916  }
917 
918  return rho;
919 }
920 
922 {
923  // Calculate the Kpipipi + multimeson phase space factor.
924  // Use the simplest definition in hep-ph/9705401 (Eq 14), which is the form
925  // used for the rest of that paper (thankfully, the amplitude does not depend
926  // significantly on the form used for the K3pi phase space factor).
927  LauComplex rho(0.0, 0.0);
928  if (TMath::Abs(s) < 1e-10) {return rho;}
929 
930  if (s < 1.44) {
931 
932  Double_t powerTerm = (-mK3piDiffSq_/s) + 1.0;
933  if (powerTerm < 0.0) {
934  rho.setImagPart( k3piFactor_*TMath::Power(-powerTerm, 2.5) );
935  } else {
936  rho.setRealPart( k3piFactor_*TMath::Power(powerTerm, 2.5) );
937  }
938 
939  } else {
940  rho.setRealPart( 1.0 );
941  }
942 
943  return rho;
944 }
945 
946 Bool_t LauKMatrixPropagator::checkPhaseSpaceType(Int_t phaseSpaceInt) const
947 {
948  Bool_t passed(kFALSE);
949 
950  if (phaseSpaceInt >= 1 && phaseSpaceInt < LauKMatrixPropagator::TotChannels) {
951  passed = kTRUE;
952  }
953 
954  return passed;
955 }
956 
958 {
959 
960  // Get the complex (unitary) transition amplitude T for the given channel
961  LauComplex TAmp(0.0, 0.0);
962 
963  channel -= 1;
964  if (channel < 0 || channel >= nChannels_) {return TAmp;}
965 
966  this->getTMatrix(s);
967 
968  TAmp.setRealPart(ReTMatrix_[index_][channel]);
969  TAmp.setImagPart(ImTMatrix_[index_][channel]);
970 
971  return TAmp;
972 
973 }
974 
976 {
977 
978  // Get the complex (unitary) transition amplitude T for the given channel
979  LauComplex rho(0.0, 0.0);
980 
981  channel -= 1;
982  if (channel < 0 || channel >= nChannels_) {return rho;}
983 
984  // If s has changed from the previous value, recalculate rho
985  if (TMath::Abs(s - previousS_) > 1e-6*s) {
986  this->calcRhoMatrix(s);
987  }
988 
989  rho.setRealPart(ReRhoMatrix_[channel][channel]);
990  rho.setImagPart(ImRhoMatrix_[channel][channel]);
991 
992  return rho;
993 
994 }
995 
997 
998  // Find the unitary T matrix, where T = [sqrt(rho)]^{*} T_hat sqrt(rho),
999  // and T_hat = (I - i K rho)^-1 * K is the Lorentz-invariant T matrix,
1000  // which has phase-space factors included (rho). This function is not
1001  // needed to calculate the K-matrix amplitudes, but allows us
1002  // to check the variation of T as a function of s (kinematics)
1003 
1004  if (!kinematics) {return;}
1005 
1006  // Get the invariant mass squared (s) from the kinematics object.
1007  // Use the resPairAmpInt to find which mass-squared combination to use.
1008  Double_t s(0.0);
1009 
1010  if (resPairAmpInt_ == 1) {
1011  s = kinematics->getm23Sq();
1012  } else if (resPairAmpInt_ == 2) {
1013  s = kinematics->getm13Sq();
1014  } else if (resPairAmpInt_ == 3) {
1015  s = kinematics->getm12Sq();
1016  }
1017 
1018  this->getTMatrix(s);
1019 
1020 }
1021 
1022 
1024 {
1025 
1026  // Find the unitary transition T matrix, where
1027  // T = [sqrt(rho)]^{*} T_hat sqrt(rho), and
1028  // T_hat = (I - i K rho)^-1 * K is the Lorentz-invariant T matrix,
1029  // which has phase-space factors included (rho). Note that the first
1030  // sqrt of the rho matrix is complex conjugated.
1031 
1032  // This function is not needed to calculate the K-matrix amplitudes, but
1033  // allows us to check the variation of T as a function of s (kinematics)
1034 
1035  // Initialse the real and imaginary parts of the T matrix to zero
1036  ReTMatrix_.Zero(); ImTMatrix_.Zero();
1037 
1038  if (parametersSet_ == kFALSE) {return;}
1039 
1040  // Update K, rho and the propagator (I - i K rho)^-1
1041  this->updatePropagator(s);
1042 
1043  // Find the real and imaginary T_hat matrices
1044  TMatrixD THatReal = realProp_*ScattKMatrix_;
1045  TMatrixD THatImag(zeroMatrix_);
1046  THatImag -= negImagProp_*ScattKMatrix_;
1047 
1048  // Find the square-root of the phase space matrix
1049  this->getSqrtRhoMatrix();
1050 
1051  // Let sqrt(rho) = A + iB and T_hat = C + iD
1052  // => T = A(CA-DB) + B(DA+CB) + i[A(DA+CB) + B(DB-CA)]
1053  TMatrixD CA(THatReal);
1054  CA *= ReSqrtRhoMatrix_;
1055  TMatrixD DA(THatImag);
1056  DA *= ReSqrtRhoMatrix_;
1057  TMatrixD CB(THatReal);
1058  CB *= ImSqrtRhoMatrix_;
1059  TMatrixD DB(THatImag);
1060  DB *= ImSqrtRhoMatrix_;
1061 
1062  TMatrixD CAmDB(CA);
1063  CAmDB -= DB;
1064  TMatrixD DApCB(DA);
1065  DApCB += CB;
1066  TMatrixD DBmCA(DB);
1067  DBmCA -= CA;
1068 
1069  // Find the real and imaginary parts of the transition matrix T
1072 
1073 }
1074 
1076 {
1077 
1078  // Find the square root of the (current) phase space matrix so that
1079  // we can find T = [sqrt(rho)}^{*} T_hat sqrt(rho), where T_hat is the
1080  // Lorentz-invariant T matrix = (I - i K rho)^-1 * K; note that the first
1081  // sqrt of rho matrix is complex conjugated
1082 
1083  // If rho = rho_i + i rho_r = a + i b, then sqrt(rho) = c + i d, where
1084  // c = sqrt(0.5*(r+a)) and d = sqrt(0.5(r-a)), where r = sqrt(a^2 + b^2).
1085  // Since rho is diagonal, then the square root of rho will also be diagonal,
1086  // with its real and imaginary matrix elements equal to c and d, respectively
1087 
1088  // Initialise the real and imaginary parts of the square root of
1089  // the rho matrix to zero
1090  ReSqrtRhoMatrix_.Zero(); ImSqrtRhoMatrix_.Zero();
1091 
1092  for (Int_t iChannel (0); iChannel < nChannels_; ++iChannel) {
1093 
1094  Double_t realRho = ReRhoMatrix_[iChannel][iChannel];
1095  Double_t imagRho = ImRhoMatrix_[iChannel][iChannel];
1096 
1097  Double_t rhoMag = sqrt(realRho*realRho + imagRho*imagRho);
1098  Double_t rhoSum = rhoMag + realRho;
1099  Double_t rhoDiff = rhoMag - realRho;
1100 
1101  Double_t reSqrtRho(0.0), imSqrtRho(0.0);
1102  if (rhoSum > 0.0) {reSqrtRho = sqrt(0.5*rhoSum);}
1103  if (rhoDiff > 0.0) {imSqrtRho = sqrt(0.5*rhoDiff);}
1104 
1105  ReSqrtRhoMatrix_[iChannel][iChannel] = reSqrtRho;
1106  ImSqrtRhoMatrix_[iChannel][iChannel] = imSqrtRho;
1107 
1108  }
1109 
1110 
1111 }
1112 
1113 LauComplex LauKMatrixPropagator::getTHat(Double_t s, Int_t channel) {
1114 
1115  LauComplex THat(0.0, 0.0);
1116  channel -= 1;
1117  if (channel < 0 || channel >= nChannels_) {return THat;}
1118 
1119  this->updatePropagator(s);
1120 
1121  // Find the real and imaginary T_hat matrices
1122  TMatrixD THatReal = realProp_*ScattKMatrix_;
1123  TMatrixD THatImag(zeroMatrix_);
1124  THatImag -= negImagProp_*ScattKMatrix_;
1125 
1126  // Return the specific THat component
1127  THat.setRealPart(THatReal[index_][channel]);
1128  THat.setImagPart(THatImag[index_][channel]);
1129 
1130  return THat;
1131 
1132 }
1133 
void calcPoleDenomVect(Double_t s)
Calulate the term 1/(m_pole^2 - s) for the scattering and production K-matrix formulae.
TMatrixD ImTMatrix_
Imaginary part of the unitary T matrix.
const Double_t mEta
Mass of eta (GeV/c^2)
Definition: LauConstants.hh:62
File containing declaration of LauTextFileParser class.
Double_t previousS_
s value of the previous pole
LauComplex calcKThreePiRho(Double_t s) const
Calculate the Kpipipi phase space factor.
Bool_t scattSymmetry_
Control if scattering constants are channel symmetric: f_ji = f_ij.
Class for parsing text files.
std::vector< Int_t > phaseSpaceTypes_
Vector of phase space types.
Double_t mKEtaPSumSq_
Defined as (mK+mEta&#39;)^2.
TMatrixD IMatrix_
Identity matrix.
virtual ~LauKMatrixPropagator()
Destructor.
void storeParameter(const TString &keyword, const TString &parString)
Store miscelleanous parameters from a line in the parameter file.
Double_t calcSVPTerm(Double_t s, Double_t s0) const
Calculate the &quot;slow-varying part&quot;.
void calcRhoMatrix(Double_t s)
Calculate the real and imaginary part of the phase space density diagonal matrix. ...
LauComplex calcFourPiRho(Double_t s) const
Calculate the 4 pi phase space factor.
LauComplex calcEtaEtaRho(Double_t s) const
Calculate the eta-eta phase space factor.
LauComplex calcKEtaPRho(Double_t s) const
Calculate the K-eta&#39; phase space factor.
void storeScattering(const std::vector< std::string > &theLine)
Store the scattering coefficients from a line in the parameter file.
TString name_
String to store the propagator name.
ClassImp(LauAbsCoeffSet)
TMatrixD ReTMatrix_
Real part of the unitary T matrix.
Double_t mKEtaPDiffSq_
Defined as (mK-mEta&#39;)^2.
LauParameter()
Default constructor.
Definition: LauParameter.cc:44
LauParameter sA0_
Constant from input file.
const TString & name() const
The parameter name.
LauComplex calcEtaEtaPRho(Double_t s) const
Calculate the eta-eta&#39; phase space factor.
Double_t re() const
Get the real part.
Definition: LauComplex.hh:219
Double_t im() const
Get the imaginary part.
Definition: LauComplex.hh:228
Double_t getRealPropTerm(Int_t channelIndex) const
Get the real part of the term of the propagator.
void getSqrtRhoMatrix()
Get the square root of the phase space matrix.
TMatrixD ReSqrtRhoMatrix_
Real part of the square root of the phase space density diagonal matrix.
LauParameter s0Prod_
Constant from input file.
TMatrixD ScattKMatrix_
Scattering K-matrix.
Double_t m2piSq_
Defined as 4*mPi*mPi.
std::vector< LauParameter > mSqPoles_
Vector of squared pole masses.
Double_t getScatteringConstant(Int_t channel1Index, Int_t channel2Index) const
Get scattering constants that were loaded from the input file.
std::vector< Double_t > poleDenomVect_
Vector of 1/(m_pole^2 - s) terms for scattering and production K-matrix formulae. ...
Bool_t checkPhaseSpaceType(Int_t phaseSpaceInt) const
Check the phase space factors that need to be used.
Double_t mEtaEtaPSumSq_
Defined as (mEta+mEta&#39;)^2.
Double_t getm23Sq() const
Get the m23 invariant mass square.
LauParArray fScattering_
Array of scattering SVP values.
File containing declaration of LauKMatrixPropagator class.
const Double_t mPi
Mass of pi+- (GeV/c^2)
Definition: LauConstants.hh:54
void storePole(const std::vector< std::string > &theLine)
Store the pole mass and couplings from a line in the parameter file.
TMatrixD zeroMatrix_
Null matrix.
void storeChannels(const std::vector< std::string > &theLine)
Store the (phase space) channel indices from a line in the parameter file.
void updateProdSVPTerm(Double_t s)
Update the production &quot;slowly-varying part&quot;.
Int_t nChannels_
Number of channels.
LauParameter mSq0_
Constant from input file.
LauComplex calcPiPiRho(Double_t s) const
Calculate the pipi phase space factor.
void setImagPart(Double_t imagpart)
Set the imaginary part.
Definition: LauComplex.hh:318
Double_t mKpiSumSq_
Defined as (mK+mPi)^2.
LauComplex calcKPiRho(Double_t s) const
Calculate the Kpi phase space factor.
Double_t getPoleDenomTerm(Int_t poleIndex) const
Get the 1/(m_pole^2 -s) terms for the scattering and production K-matrix formulae.
LauParameter s0Scatt_
Constant from input file.
const Double_t mEtaSq
Square of eta mass.
Definition: LauConstants.hh:87
void updateAdlerZeroFactor(Double_t s)
Calculate the multiplicative factor containing severa Adler zero constants.
TMatrixD ImRhoMatrix_
Imaginary part of the phase space density diagonal matrix.
Double_t mK3piDiffSq_
Defined as (mK-3*mPi)^2.
Double_t sAConst_
Defined as 0.5*sA*mPi*mPi.
const Double_t mPiSq
Square of pi+- mass.
Definition: LauConstants.hh:79
TMatrixD ImSqrtRhoMatrix_
Imaginary part of the square root of the phase space density diagonal matrix.
Double_t prodSVP_
&quot;slowly-varying part&quot; for the production K-matrix
Bool_t parametersSet_
Tracks if all params have been set.
std::vector< std::string > getLine(UInt_t lineNo)
Retrieve the specified line.
Class for defining the fit parameter objects.
Definition: LauParameter.hh:49
const Double_t mEtaPrime
Mass of eta&#39; (GeV/c^2)
Definition: LauConstants.hh:64
void setRealPart(Double_t realpart)
Set the real part.
Definition: LauComplex.hh:309
Double_t getCouplingConstant(Int_t poleIndex, Int_t channelIndex) const
Get coupling constants that were loaded from the input file.
LauParameter sA_
Constant from input file.
void updatePropagator(const LauKinematics *kinematics)
Calculate the invariant mass squared s.
LauComplex getTHat(Double_t s, Int_t channel)
Get the THat amplitude for the given s and channel number.
Int_t index_
Row index - 1.
Double_t mKpiDiffSq_
Defined as (mK-mPi)^2.
Double_t getImagPropTerm(Int_t channelIndex) const
Get the imaginary part of the term of the propagator.
Double_t getm12Sq() const
Get the m12 invariant mass square.
Bool_t verbose_
Control the output of the functions.
void getTMatrix(const LauKinematics *kinematics)
Get the unitary transition amplitude matrix for the given kinematics.
void setParameters(const TString &inputFile)
Read an input file to set parameters.
Double_t m2KSq_
Defined as 4*mK*mK.
TString name_
The parameter name.
Double_t getm13Sq() const
Get the m13 invariant mass square.
Double_t k3piFactor_
Factor used to calculate the Kpipipi phase space term.
Int_t nPoles_
Number of poles.
File containing LauConstants namespace.
Double_t unblindValue() const
The unblinded value of the parameter.
Class for defining a complex number.
Definition: LauComplex.hh:61
void initialiseMatrices()
Initialise and set the dimensions for the internal matrices and parameter arrays. ...
Double_t adlerZeroFactor_
Multiplicative factor containing various Adler zero constants.
LauParArray gCouplings_
Array of coupling constants.
Class for calculating 3-body kinematic quantities.
LauComplex getTransitionAmp(Double_t s, Int_t channel)
Get the unitary transition amplitude for the given channel.
TMatrixD realProp_
Real part of the propagator matrix.
Int_t resPairAmpInt_
Number to identify the DP axis in question.
Double_t scattSVP_
&quot;slowly-varying part&quot; for the scattering K-matrix
LauComplex getPropTerm(Int_t channelIndex) const
Get the full complex propagator term for a given channel.
LauComplex calcKKRho(Double_t s) const
Calculate the KK phase space factor.
void calcScattKMatrix(Double_t s)
Calculate the scattering K-matrix for the given value of s.
void processFile()
Parse the file.
Double_t m2EtaSq_
Defined as 4*mEta*mEta.
Double_t fourPiFactor1_
Factor used to calculate the pipipipi phase space term.
const Double_t mK
Mass of K+- (GeV/c^2)
Definition: LauConstants.hh:58
LauComplex getPhaseSpaceTerm(Double_t s, Int_t channel)
Get the complex phase space term for the given channel and invariant mass squared.
const Double_t mKSq
Square of K+- mass.
Definition: LauConstants.hh:83
void updateScattSVPTerm(Double_t s)
Update the scattering &quot;slowly-varying part&quot;.
UInt_t getTotalNumLines() const
Get the total number of lines that are not comments.
TMatrixD negImagProp_
Imaginary part of the propagator matrix.
Class for defining a K-matrix propagator.
TMatrixD ReRhoMatrix_
Real part of the phase space density diagonal matrix.