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