laura is hosted by Hepforge, IPPP Durham
Laura++  v1r0
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 
21 #include <iostream>
22 #include <fstream>
23 #include <cmath>
24 #include <cstdlib>
25 
26 using std::cout;
27 using std::endl;
28 using std::cerr;
29 
30 ClassImp(LauKMatrixPropagator)
31 
32 LauKMatrixPropagator::LauKMatrixPropagator(const TString& name, const TString& paramFile,
33  Int_t resPairAmpInt, Int_t nChannels,
34  Int_t nPoles, Int_t rowIndex) :
35  name_(name),
36  paramFileName_(paramFile),
37  resPairAmpInt_(resPairAmpInt),
38  index_(rowIndex - 1),
39  previousS_(0.0),
40  scattSVP_(0.0),
41  prodSVP_(0.0),
42  nChannels_(nChannels),
43  nPoles_(nPoles),
44  sAConst_(0.0),
45  m2piSq_(4.0*LauConstants::mPiSq),
46  m2KSq_( 4.0*LauConstants::mKSq),
47  m2EtaSq_(4.0*LauConstants::mEtaSq),
48  mEtaEtaPSumSq_((LauConstants::mEta + LauConstants::mEtaPrime)*(LauConstants::mEta + LauConstants::mEtaPrime)),
49  mEtaEtaPDiffSq_((LauConstants::mEta - LauConstants::mEtaPrime)*(LauConstants::mEta - LauConstants::mEtaPrime)),
50  mKpiSumSq_((LauConstants::mK + LauConstants::mPi)*(LauConstants::mK + LauConstants::mPi)),
51  mKpiDiffSq_((LauConstants::mK - LauConstants::mPi)*(LauConstants::mK - LauConstants::mPi)),
52  mKEtaPSumSq_((LauConstants::mK + LauConstants::mEtaPrime)*(LauConstants::mK + LauConstants::mEtaPrime)),
53  mKEtaPDiffSq_((LauConstants::mK - LauConstants::mEtaPrime)*(LauConstants::mK - LauConstants::mEtaPrime)),
54  mK3piDiffSq_((LauConstants::mK - 3.0*LauConstants::mPi)*(LauConstants::mK - 3.0*LauConstants::mPi)),
55  k3piFactor_(TMath::Power((1.44 - mK3piDiffSq_)/1.44, -2.5)),
56  fourPiFactor1_(16.0*LauConstants::mPiSq),
57  fourPiFactor2_(TMath::Sqrt(1.0 - fourPiFactor1_)),
58  adlerZeroFactor_(0.0),
59  parametersSet_(kFALSE),
60  verbose_(kFALSE)
61 {
62  // Constructor
63  realProp_.Clear(); negImagProp_.Clear();
64  this->setParameters(paramFile);
65 }
66 
68 {
69  // Destructor
70  realProp_.Clear();
71  negImagProp_.Clear();
72 }
73 
75 {
76  // Get the (i,j) = (index_, channel) term of the propagator
77  // matrix. This allows us not to return the full propagator matrix.
78  Double_t realTerm = this->getRealPropTerm(channel);
79  Double_t imagTerm = this->getImagPropTerm(channel);
80 
81  LauComplex propTerm(realTerm, imagTerm);
82  return propTerm;
83 
84 }
85 
86 Double_t LauKMatrixPropagator::getRealPropTerm(Int_t channel) const
87 {
88  // Get the real part of the (i,j) = (index_, channel) term of the propagator
89  // matrix. This allows us not to return the full propagator matrix.
90  if (parametersSet_ == kFALSE) {return 0.0;}
91 
92  Double_t propTerm = realProp_[index_][channel];
93  return propTerm;
94 
95 }
96 
97 Double_t LauKMatrixPropagator::getImagPropTerm(Int_t channel) const
98 {
99  // Get the imaginary part of the (i,j) = (index_, channel) term of the propagator
100  // matrix. This allows us not to return the full propagator matrix.
101  if (parametersSet_ == kFALSE) {return 0.0;}
102 
103  Double_t imagTerm = -1.0*negImagProp_[index_][channel];
104  return imagTerm;
105 
106 }
107 
109 {
110 
111  // Calculate the K-matrix propagator for the given s value.
112  // The K-matrix amplitude is given by
113  // T_i = sum_{ij} (I - iK*rho)^-1 * P_j, where P is the production K-matrix.
114  // i = index for the state (e.g. S-wave index = 0).
115  // Here, we only find the (I - iK*rho)^-1 matrix part.
116 
117  // Get the invariant mass squared (s) from the kinematics object.
118  // Use the resPairAmpInt to find which mass-squared combination to use.
119  Double_t s(0.0);
120  if (resPairAmpInt_ == 1) {
121  s = kinematics->getm23Sq();
122  } else if (resPairAmpInt_ == 2) {
123  s = kinematics->getm13Sq();
124  } else if (resPairAmpInt_ == 3) {
125  s = kinematics->getm12Sq();
126  }
127 
128  this->updatePropagator(s);
129 
130 }
131 
133 {
134  // Calculate the K-matrix propagator for the given s value.
135  // The K-matrix amplitude is given by
136  // T_i = sum_{ij} (I - iK*rho)^-1 * P_j, where P is the production K-matrix.
137  // i = index for the state (e.g. S-wave index = 0).
138  // Here, we only find the (I - iK*rho)^-1 matrix part.
139 
140  // Check if we have almost the same s value as before. If so, don't re-calculate
141  // the propagator nor any of the pole mass summation terms.
142  if (TMath::Abs(s - previousS_) < 1e-6) {
143  //cout<<"Already got propagator for s = "<<s<<endl;
144  return;
145  }
146 
147  if (parametersSet_ == kFALSE) {
148  //cerr<<"ERROR in LauKMatrixPropagator::updatePropagator. Parameters have not been set."<<endl;
149  return;
150  }
151 
152  // Calculate the denominator pole mass terms and Adler zero factor
153  this->calcPoleDenomVect(s);
154  this->updateAdlerZeroFactor(s);
155 
156  // Calculate the scattering K-matrix (real and symmetric)
157  this->calcScattKMatrix(s);
158  // Calculate the phase space density matrix, which is diagonal, but can be complex
159  // if the quantity s is below various threshold values (analytic continuation).
160  this->calcRhoMatrix(s);
161 
162  // Calculate K*rho (real and imaginary parts, since rho can be complex)
163  TMatrixD K_realRho(ScattKMatrix_);
164  K_realRho *= ReRhoMatrix_;
165  TMatrixD K_imagRho(ScattKMatrix_);
166  K_imagRho *= ImRhoMatrix_;
167 
168  // A = I + K*Imag(rho), B = -K*Real(Rho)
169  // Calculate C and D matrices such that (A + iB)*(C + iD) = I,
170  // ie. C + iD = (I - i K*rho)^-1, separated into real and imaginary parts.
171  TMatrixD A(IMatrix_);
172  A += K_imagRho;
173  TMatrixD B(zeroMatrix_);
174  B -= K_realRho;
175 
176  TMatrixD invA(TMatrixD::kInverted, A);
177  TMatrixD invA_B(invA);
178  invA_B *= B;
179  TMatrixD B_invA_B(B);
180  B_invA_B *= invA_B;
181 
182  TMatrixD invC(A);
183  invC += B_invA_B;
184 
185  // Set the real part of the propagator matrix ("C")
186  realProp_ = TMatrixD(TMatrixD::kInverted, invC);
187 
188  // Set the (negative) imaginary part of the propagator matrix ("-D")
189  TMatrixD BC(B);
190  BC *= realProp_;
191  negImagProp_ = TMatrixD(invA);
192  negImagProp_ *= BC;
193 
194  // Also calculate the production SVP term, since this uses Adler-zero parameters
195  // defined in the parameter file.
196  this->updateProdSVPTerm(s);
197 
198  // Finally, keep track of the value of s we just used.
199  previousS_ = s;
200 
201 }
202 
203 void LauKMatrixPropagator::setParameters(const TString& inputFile)
204 {
205 
206  // Read an input file that specifies the values of the coupling constants g_i for
207  // the given number of poles and their (bare) masses. Also provided are the f_{ab}
208  // slow-varying constants. The input file should also provide the Adler zero
209  // constants s_0, s_A and s_A0.
210 
211  // Format of input file:
212  // Indices (nChannels) of phase space channel types (defined in KMatrixChannels enum)
213  // Then the bare pole mass and coupling constants over all channels for each pole:
214  // m_{alpha} g_1^{alpha} g_2^{alpha} ... g_n^{alpha} where n = nChannels
215  // Repeat for N poles. Then define the f_{ab} scattering coefficients
216  // f_{11} f_{12} ... f_{1n}
217  // f_{21} f_{22} ... f_{2n} etc.. where n = nChannels
218  // Note that the K-matrix will be symmetrised (by definition), so it is sufficient
219  // to only have the upper half filled with sensible values.
220  // Then define the Adler-zero and slowly-varying part (SVP) constants
221  // m0^2 s0Scatt s0Prod sA sA0
222  // where m0^2 = mass-squared constant for scattering production f_{ab} numerator
223  parametersSet_ = kFALSE;
224 
225  cout<<"Initialising K-matrix propagator "<<name_<<" parameters from "<<inputFile.Data()<<endl;
226 
227  cout<<"nChannels = "<<nChannels_<<", nPoles = "<<nPoles_<<endl;
228 
229  LauTextFileParser readFile(inputFile);
230  readFile.processFile();
231 
232  // Get the first (non-comment) line
233  std::vector<std::string> channelTypeLine = readFile.getNextLine();
234  Int_t nTypes = static_cast<Int_t>(channelTypeLine.size());
235  if (nTypes != nChannels_) {
236  cerr<<"Error in LauKMatrixPropagator::setParameters. Input file "<<inputFile
237  <<" defines "<<nTypes<<" channels when "<<nChannels_<<" are expected"<<endl;
238  return;
239  }
240 
241  // Get the list of channel indices to specify what phase space factors should be used
242  // e.g. pipi, Kpi, eta eta', 4pi etc..
243  Int_t iChannel(0);
244  phaseSpaceTypes_.clear();
245  for (iChannel = 0; iChannel < nChannels_; iChannel++) {
246 
247  Int_t phaseSpaceInt = atoi(channelTypeLine[iChannel].c_str());
248  Bool_t checkChannel = this->checkPhaseSpaceType(phaseSpaceInt);
249 
250  if (checkChannel == kTRUE) {
251  cout<<"Adding phase space channel index "<<phaseSpaceInt
252  <<" to K-matrix propagator "<<name_<<endl;
253  phaseSpaceTypes_.push_back(phaseSpaceInt);
254  } else {
255  cerr<<"Phase space channel index "<<phaseSpaceInt
256  <<" should be between 1 and "<<LauKMatrixPropagator::TotChannels-1<<endl;
257  cerr<<"Stopping initialisation of the K-matrix propagator "<<name_<<endl;
258  cerr<<"Please correct the channel index on the first line of "<<inputFile<<endl;
259  return;
260  }
261 
262  }
263 
264  // Clear vector of pole masses^2 as well as the map of the coupling constant vectors.
265  mSqPoles_.clear(); gCouplings_.clear();
266 
267  Double_t poleMass(0.0), poleMassSq(0.0), couplingConst(0.0);
268  Int_t iPole(0);
269  Int_t nPoleNumbers(nChannels_+1);
270 
271  for (iPole = 0; iPole < nPoles_; iPole++) {
272 
273  // Read the next line of bare pole mass and its coupling constants
274  // over all channels (nChannels_)
275  std::vector<std::string> poleLine = readFile.getNextLine();
276  Int_t nPoleValues = static_cast<Int_t>(poleLine.size());
277  if (nPoleValues != nPoleNumbers) {
278  cerr<<"Error in LauKMatrixPropagator::setParameters. Expecting "
279  <<nPoleNumbers<<" numbers for the bare pole "<<iPole<<" line (mass and "
280  <<nChannels_<<" coupling constants). Found "<<nPoleValues<<" values instead."
281  <<endl;
282  return;
283  }
284 
285  poleMass = atof(poleLine[0].c_str());
286  poleMassSq = poleMass*poleMass;
287  LauParameter mPoleParam(poleMassSq);
288  mSqPoles_.push_back(mPoleParam);
289 
290  std::vector<LauParameter> couplingVector;
291 
292  cout<<"Added bare pole mass "<<poleMass<<" GeV for pole number "<<iPole+1<<endl;
293 
294  for (iChannel = 0; iChannel < nChannels_; iChannel++) {
295 
296  couplingConst = atof(poleLine[iChannel+1].c_str());
297  LauParameter couplingParam(couplingConst);
298  couplingVector.push_back(couplingParam);
299 
300  cout<<"Added coupling parameter g^"<<iPole+1<<"_"<<iChannel+1<<" = "<<couplingConst<<endl;
301 
302  }
303 
304  gCouplings_[iPole] = couplingVector;
305 
306  }
307 
308  // Scattering (slowly-varying part, or SVP) values
309 
310  fScattering_.clear();
311  Int_t jChannel(0);
312 
313  std::vector<std::string> scatteringLine = readFile.getNextLine();
314  Int_t nScatteringValues = static_cast<Int_t>(scatteringLine.size());
315  if (nScatteringValues != nChannels_) {
316  cerr<<"Error in LauKMatrixPropagator::setParameters. Expecting "
317  <<nChannels_<<" scattering SVP f constants instead of "<<nScatteringValues
318  <<" for the K-matrix row index "<<index_<<endl;
319  return;
320  }
321 
322  std::vector<LauParameter> scatteringVect;
323  for (jChannel = 0; jChannel < nChannels_; jChannel++) {
324 
325  Double_t fScattValue = atof(scatteringLine[jChannel].c_str());
326  LauParameter scatteringParam(fScattValue);
327  scatteringVect.push_back(scatteringParam);
328 
329  cout<<"Added scattering SVP f("<<index_+1<<","<<jChannel+1<<") = "<<fScattValue<<endl;
330 
331  }
332 
333  fScattering_[index_] = scatteringVect;
334 
335  // Now extract the constants for the "Adler-zero" terms
336  std::vector<std::string> constLine = readFile.getNextLine();
337  Int_t nConstValues = static_cast<Int_t>(constLine.size());
338  if (nConstValues != 5) {
339  cerr<<"Error in LauKMatrixPropagator::setParameters. Expecting 5 values for the Adler-zero constants:"
340  <<" m0^2 s0Scatt s0Prod sA sA0"<<endl;
341  return;
342  }
343 
344  Double_t mSq0Value = atof(constLine[0].c_str());
345  Double_t s0ScattValue = atof(constLine[1].c_str());
346  Double_t s0ProdValue = atof(constLine[2].c_str());
347  Double_t sAValue = atof(constLine[3].c_str());
348  Double_t sA0Value = atof(constLine[4].c_str());
349 
350  cout<<"Adler zero constants:"<<endl;
351  cout<<"m0Sq = "<<mSq0Value<<", s0Scattering = "<<s0ScattValue<<", s0Production = "
352  <<s0ProdValue<<", sA = "<<sAValue<<" and sA0 = "<<sA0Value<<endl;
353 
354  mSq0_ = LauParameter("mSq0", mSq0Value);
355  s0Scatt_ = LauParameter("s0Scatt", s0ScattValue);
356  s0Prod_ = LauParameter("s0Prod", s0ProdValue);
357  sA_ = LauParameter("sA", sAValue);
358  sA0_ = LauParameter("sA0", sA0Value);
359  sAConst_ = 0.5*sAValue*LauConstants::mPiSq;
360 
361  // Identity and null matrices
362  IMatrix_.Clear();
363  IMatrix_.ResizeTo(nChannels_, nChannels_);
364  for (iChannel = 0; iChannel < nChannels_; iChannel++) {
365  IMatrix_[iChannel][iChannel] = 1.0;
366  }
367 
368  zeroMatrix_.Clear();
369  zeroMatrix_.ResizeTo(nChannels_, nChannels_);
370 
371  realProp_.Clear(); negImagProp_.Clear();
372  realProp_.ResizeTo(nChannels_, nChannels_);
373  negImagProp_.ResizeTo(nChannels_, nChannels_);
374 
375  // All required parameters have been set
376  parametersSet_ = kTRUE;
377 
378  cout<<"Finished initialising K-matrix propagator "<<name_<<endl;
379 
380 }
381 
383 {
384 
385  // Calculate the scattering K-matrix for the given value of s.
386  // We need to obtain the complete matrix (not just one row/column)
387  // to get the correct inverted (I - i K rho) terms for the propagator.
388 
389  if (verbose_) {cout<<"Within calcScattKMatrix for s = "<<s<<endl;}
390 
391  // Clear the internal matrix and resize.
392  ScattKMatrix_.Clear();
394 
395  Int_t iChannel(0), jChannel(0), iPole(0);
396 
397  // The pole denominator 1/(m^2 - s) terms should already be calculated
398  // by the calcPoleDenomVect() function. These same terms are also
399  // used for calculating the production K-matrix elements.
400 
401  // Calculate the "slowly-varying part" (SVP), e.g. (1 GeV - s0)/(s - s0)
402  this->updateScattSVPTerm(s);
403 
404  // Now loop over iChannel, jChannel to calculate Kij = Kji.
405  for (iChannel = 0; iChannel < nChannels_; iChannel++) {
406 
407  std::vector<LauParameter> SVPVect(0);
408 
409  KMatrixParamMap::iterator SVPIter = fScattering_.find(iChannel);
410  if (SVPIter != fScattering_.end()) {SVPVect = SVPIter->second;}
411 
412  Int_t SVPVectSize = static_cast<Int_t>(SVPVect.size());
413 
414  // Scattering matrix is real and symmetric. Start j loop from i.
415  for (jChannel = iChannel; jChannel < nChannels_; jChannel++) {
416 
417  Double_t Kij(0.0);
418 
419  // Calculate pole mass summation term
420  for (iPole = 0; iPole < nPoles_; iPole++) {
421 
422  KMatrixParamMap::iterator couplingIter = gCouplings_.find(iPole);
423 
424  std::vector<LauParameter> couplingVect = couplingIter->second;
425 
426  Double_t g_i = couplingVect[iChannel].value();
427  Double_t g_j = couplingVect[jChannel].value();
428 
429  Kij += poleDenomVect_[iPole]*g_i*g_j;
430  if (verbose_) {cout<<"1: Kij for i = "<<iChannel<<", j = "<<jChannel<<" = "<<Kij<<endl;}
431 
432  }
433 
434  if (SVPVectSize != 0 && jChannel < SVPVectSize) {
435  Double_t fij = SVPVect[jChannel].value();
436  Kij += fij*scattSVP_;
437  }
438 
439  Kij *= adlerZeroFactor_;
440  if (verbose_) {cout<<"2: Kij for i = "<<iChannel<<", j = "<<jChannel<<" = "<<Kij<<endl;}
441 
442  // Assign the TMatrix (i,j) element to the variable Kij and Kji (symmetry)
443  ScattKMatrix_(iChannel, jChannel) = Kij;
444  ScattKMatrix_(jChannel, iChannel) = Kij;
445 
446  } // j loop
447 
448  } // i loop
449 
450 }
451 
453 {
454  // Calculate the 1/(m_pole^2 - s) terms for the scattering
455  // and production K-matrix formulae.
456  poleDenomVect_.clear();
457  Int_t iPole(0);
458  for (iPole = 0; iPole < nPoles_; iPole++) {
459 
460  Double_t poleTerm = mSqPoles_[iPole] - s;
461  Double_t invPoleTerm(0.0);
462  if (TMath::Abs(poleTerm) > 1.0e-6) {invPoleTerm = 1.0/poleTerm;}
463 
464  poleDenomVect_.push_back(invPoleTerm);
465 
466  }
467 
468 }
469 
470 Double_t LauKMatrixPropagator::getPoleDenomTerm(Int_t poleIndex) const
471 {
472 
473  if (parametersSet_ == kFALSE) {return 0.0;}
474 
475  Double_t poleDenom(0.0);
476  poleDenom = poleDenomVect_[poleIndex];
477  return poleDenom;
478 
479 }
480 
481 Double_t LauKMatrixPropagator::getCouplingConstant(Int_t poleIndex, Int_t channelIndex) const
482 {
483 
484  if (parametersSet_ == kFALSE) {return 0.0;}
485 
486  Double_t couplingConst(0.0);
487  KMatrixParamMap::const_iterator couplingIter = gCouplings_.find(poleIndex);
488  std::vector<LauParameter> couplingVect = couplingIter->second;
489  couplingConst = couplingVect[channelIndex].value();
490 
491  return couplingConst;
492 
493 }
494 
495 Double_t LauKMatrixPropagator::calcSVPTerm(Double_t s, Double_t s0) const
496 {
497 
498  if (parametersSet_ == kFALSE) {return 0.0;}
499 
500  // Calculate the "slowly-varying part" (SVP)
501  Double_t result(0.0);
502  Double_t deltaS = s - s0;
503  if (TMath::Abs(deltaS) > 1.0e-6) {
504  result = (mSq0_.value() - s0)/deltaS;
505  }
506 
507  return result;
508 
509 }
510 
512 {
513  // Update the scattering "slowly-varying part" (SVP)
514  Double_t s0Scatt = s0Scatt_.value();
515  scattSVP_ = this->calcSVPTerm(s, s0Scatt);
516 }
517 
519 {
520  // Update the production "slowly-varying part" (SVP)
521  Double_t s0Prod = s0Prod_.value();
522  prodSVP_ = this->calcSVPTerm(s, s0Prod);
523 }
524 
526 {
527 
528  // Calculate the multiplicative factor containing various Adler zero
529  // constants.
530  adlerZeroFactor_ = 0.0;
531 
532  Double_t sA0Val = sA0_.value();
533  Double_t deltaS = s - sA0Val;
534  if (TMath::Abs(deltaS) > 1e-6) {
535  adlerZeroFactor_ = (s - sAConst_)*(1.0 - sA0Val)/deltaS;
536  }
537 
538 }
539 
541 {
542 
543  // Calculate the real and imaginary part of the phase space density
544  // diagonal matrix for the given invariant mass squared quantity, s.
545  // The matrix can be complex if s is below threshold (so that
546  // the amplitude continues analytically).
547  ReRhoMatrix_.Clear(); ImRhoMatrix_.Clear();
550 
551  LauComplex rho(0.0, 0.0);
552  Int_t phaseSpaceIndex(0);
553 
554  for (Int_t iChannel (0); iChannel < nChannels_; ++iChannel) {
555 
556  phaseSpaceIndex = phaseSpaceTypes_[iChannel];
557 
558  if (phaseSpaceIndex == LauKMatrixPropagator::PiPi) {
559  rho = this->calcPiPiRho(s);
560  } else if (phaseSpaceIndex == LauKMatrixPropagator::KK) {
561  rho = this->calcKKRho(s);
562  } else if (phaseSpaceIndex == LauKMatrixPropagator::FourPi) {
563  rho = this->calcFourPiRho(s);
564  } else if (phaseSpaceIndex == LauKMatrixPropagator::EtaEta) {
565  rho = this->calcEtaEtaRho(s);
566  } else if (phaseSpaceIndex == LauKMatrixPropagator::EtaEtaP) {
567  rho = this->calcEtaEtaPRho(s);
568  } else if (phaseSpaceIndex == LauKMatrixPropagator::KPi) {
569  rho = this->calcKPiRho(s);
570  } else if (phaseSpaceIndex == LauKMatrixPropagator::KEtaP) {
571  rho = this->calcKEtaPRho(s);
572  } else if (phaseSpaceIndex == LauKMatrixPropagator::KThreePi) {
573  rho = this->calcKThreePiRho(s);
574  }
575 
576  if (verbose_) {
577  cout<<"ReRhoMatrix("<<iChannel<<", "<<iChannel<<") = "<<rho.re()<<endl;
578  cout<<"ImRhoMatrix("<<iChannel<<", "<<iChannel<<") = "<<rho.im()<<endl;
579  }
580 
581  ReRhoMatrix_(iChannel, iChannel) = rho.re();
582  ImRhoMatrix_(iChannel, iChannel) = rho.im();
583 
584  }
585 
586 }
587 
589 {
590  // Calculate the pipi phase space factor
591  LauComplex rho(0.0, 0.0);
592  if (TMath::Abs(s) < 1e-10) {return rho;}
593 
594  Double_t sqrtTerm = (-m2piSq_/s) + 1.0;
595  if (sqrtTerm < 0.0) {
596  rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
597  } else {
598  rho.setRealPart( TMath::Sqrt(sqrtTerm) );
599  }
600 
601  return rho;
602 }
603 
605 {
606  // Calculate the KK phase space factor
607  LauComplex rho(0.0, 0.0);
608  if (TMath::Abs(s) < 1e-10) {return rho;}
609 
610  Double_t sqrtTerm = (-m2KSq_/s) + 1.0;
611  if (sqrtTerm < 0.0) {
612  rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
613  } else {
614  rho.setRealPart( TMath::Sqrt(sqrtTerm) );
615  }
616 
617  return rho;
618 }
619 
621 {
622  // Calculate the 4pi phase space factor. This uses a parameterisation that approximates
623  // the multi-body phase space double integral.
624  LauComplex rho(0.0, 0.0);
625  if (TMath::Abs(s) < 1e-10) {return rho;}
626 
627  if (s <= 1.0) {
628  Double_t term1 = (((11.3153*s - 21.8845)*s + 16.8358)*s - 6.39017)*s + 1.2274;
629  term1 += 0.00370909/(s*s);
630  term1 -= 0.111203/s;
631  rho.setRealPart( term1*fourPiFactor2_ );
632  } else {
633  rho.setRealPart( TMath::Sqrt(1.0 - (fourPiFactor1_/s)) );
634  }
635 
636  return rho;
637 }
638 
640 {
641  // Calculate the eta-eta phase space factor
642  LauComplex rho(0.0, 0.0);
643  if (TMath::Abs(s) < 1e-10) {return rho;}
644 
645  Double_t sqrtTerm = (-m2EtaSq_/s) + 1.0;
646  if (sqrtTerm < 0.0) {
647  rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
648  } else {
649  rho.setRealPart( TMath::Sqrt(sqrtTerm) );
650  }
651 
652  return rho;
653 }
654 
656 {
657  // Calculate the eta-eta' phase space factor
658  LauComplex rho(0.0, 0.0);
659  if (TMath::Abs(s) < 1e-10) {return rho;}
660 
661  Double_t sqrtTerm1 = (-mEtaEtaPSumSq_/s) + 1.0;
662  Double_t sqrtTerm2 = (-mEtaEtaPDiffSq_/s) + 1.0;
663  Double_t sqrtTerm = sqrtTerm1*sqrtTerm2;
664  if (sqrtTerm < 0.0) {
665  rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
666  } else {
667  rho.setRealPart( TMath::Sqrt(sqrtTerm) );
668  }
669 
670  return rho;
671 }
672 
673 
675 {
676  // Calculate the K-pi phase space factor
677  LauComplex rho(0.0, 0.0);
678  if (TMath::Abs(s) < 1e-10) {return rho;}
679 
680  Double_t sqrtTerm1 = (-mKpiSumSq_/s) + 1.0;
681  Double_t sqrtTerm2 = (-mKpiDiffSq_/s) + 1.0;
682  Double_t sqrtTerm = sqrtTerm1*sqrtTerm2;
683  if (sqrtTerm < 0.0) {
684  rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
685  } else {
686  rho.setRealPart( TMath::Sqrt(sqrtTerm) );
687  }
688 
689  return rho;
690 }
691 
693 {
694  // Calculate the K-eta' phase space factor
695  LauComplex rho(0.0, 0.0);
696  if (TMath::Abs(s) < 1e-10) {return rho;}
697 
698  Double_t sqrtTerm1 = (-mKEtaPSumSq_/s) + 1.0;
699  Double_t sqrtTerm2 = (-mKEtaPDiffSq_/s) + 1.0;
700  Double_t sqrtTerm = sqrtTerm1*sqrtTerm2;
701  if (sqrtTerm < 0.0) {
702  rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
703  } else {
704  rho.setRealPart( TMath::Sqrt(sqrtTerm) );
705  }
706 
707  return rho;
708 }
709 
711 {
712  // Calculate the Kpipipi + multimeson phase space factor.
713  // Use the simplest definition in hep-ph/9705401 (Eq 14), which is the form
714  // used for the rest of that paper (thankfully, the amplitude does not depend
715  // significantly on the form used for the K3pi phase space factor).
716  LauComplex rho(0.0, 0.0);
717  if (TMath::Abs(s) < 1e-10) {return rho;}
718 
719  if (s < 1.44) {
720 
721  Double_t powerTerm = (-mK3piDiffSq_/s) + 1.0;
722  if (powerTerm < 0.0) {
723  rho.setImagPart( k3piFactor_*TMath::Power(-powerTerm, 2.5) );
724  } else {
725  rho.setRealPart( k3piFactor_*TMath::Power(powerTerm, 2.5) );
726  }
727 
728  } else {
729  rho.setRealPart( 1.0 );
730  }
731 
732  return rho;
733 }
734 
735 Bool_t LauKMatrixPropagator::checkPhaseSpaceType(Int_t phaseSpaceInt) const
736 {
737  Bool_t passed(kFALSE);
738 
739  if (phaseSpaceInt >= 1 && phaseSpaceInt < LauKMatrixPropagator::TotChannels) {
740  passed = kTRUE;
741  }
742 
743  return passed;
744 }
void calcPoleDenomVect(Double_t s)
Calulate the term 1/(m_pole^2 - s) for the scattering and production K-matrix formulae.
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.
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.
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.
TString name_
String to store the propagator name.
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:202
std::vector< std::string > getNextLine()
Retrieve the next line.
Double_t im() const
Get the imaginary part.
Definition: LauComplex.hh:211
Double_t getRealPropTerm(Int_t channelIndex) const
Get the real part of the term of the propagator.
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.
std::vector< Double_t > poleDenomVect_
Vector of 1/(m_pole^2 - s) terms for scattering and production K-matrix formulae. ...
Double_t mEtaEtaPDiffSq_
Defined as (mEta-mEta&#39;)^2.
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.
File containing declaration of LauKMatrixPropagator class.
const Double_t mPi
Mass of pi+- (GeV/c^2)
Definition: LauConstants.hh:40
TMatrixD zeroMatrix_
Null matrix.
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:301
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
KMatrixParamMap gCouplings_
Map of coupling constants.
void updateAdlerZeroFactor(Double_t s)
Calculate the multiplicative factor containing severa Adler zero constants.
TMatrixD ImRhoMatrix_
Imaginary part of the pahse 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
Double_t fourPiFactor2_
Factor used to calculate the pipipipi phase space term.
Double_t prodSVP_
&quot;slowly-varying part&quot; for the production K-matrix
Bool_t parametersSet_
Tracks if all params have been set.
Class for defining the fit parameter objects.
Definition: LauParameter.hh:31
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:292
KMatrixParamMap fScattering_
Map of scattering SVP values.
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.
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 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.
Class for defining a complex number.
Definition: LauComplex.hh:47
Double_t adlerZeroFactor_
Multiplicative factor containing various Adler zero constants.
Class for calculating 3-body kinematic quantities.
Double_t value() const
The value of the parameter.
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
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;.
TMatrixD negImagProp_
Imaginary part of the propagator matrix.
Class for defining a K-matrix propagator.
TMatrixD ReRhoMatrix_
Real part of the pahse space density diagonal matrix.