laura is hosted by Hepforge, IPPP Durham
Laura++  v3r0
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 {
63  // Constructor
64 
65  // Check that the index is OK
66  if (index_ < 0 || index_ >= nChannels_) {
67  std::cerr << "ERROR in LauKMatrixPropagator constructor. The rowIndex, which is set to "
68  << rowIndex << ", must be between 1 and the number of channels "
69  << nChannels_ << std::endl;
70  gSystem->Exit(EXIT_FAILURE);
71  }
72 
73  this->setParameters(paramFile);
74 }
75 
77 {
78  // Destructor
79  realProp_.Clear();
80  negImagProp_.Clear();
81  ScattKMatrix_.Clear();
82  ReRhoMatrix_.Clear();
83  ImRhoMatrix_.Clear();
84  ReTMatrix_.Clear();
85  ImTMatrix_.Clear();
86  IMatrix_.Clear();
87  zeroMatrix_.Clear();
88 }
89 
91 {
92  // Get the (i,j) = (index_, channel) term of the propagator
93  // matrix. This allows us not to return the full propagator matrix.
94  Double_t realTerm = this->getRealPropTerm(channel);
95  Double_t imagTerm = this->getImagPropTerm(channel);
96 
97  LauComplex propTerm(realTerm, imagTerm);
98  return propTerm;
99 
100 }
101 
102 Double_t LauKMatrixPropagator::getRealPropTerm(Int_t channel) const
103 {
104  // Get the real part of the (i,j) = (index_, channel) term of the propagator
105  // matrix. This allows us not to return the full propagator matrix.
106  if (parametersSet_ == kFALSE) {return 0.0;}
107 
108  Double_t propTerm = realProp_[index_][channel];
109  return propTerm;
110 
111 }
112 
113 Double_t LauKMatrixPropagator::getImagPropTerm(Int_t channel) const
114 {
115  // Get the imaginary part of the (i,j) = (index_, channel) term of the propagator
116  // matrix. This allows us not to return the full propagator matrix.
117  if (parametersSet_ == kFALSE) {return 0.0;}
118 
119  Double_t imagTerm = -1.0*negImagProp_[index_][channel];
120  return imagTerm;
121 
122 }
123 
125 {
126 
127  // Calculate the K-matrix propagator for the given s value.
128  // The K-matrix amplitude is given by
129  // T_i = sum_{ij} (I - iK*rho)^-1 * P_j, where P is the production K-matrix.
130  // i = index for the state (e.g. S-wave index = 0).
131  // Here, we only find the (I - iK*rho)^-1 matrix part.
132 
133  if (!kinematics) {return;}
134 
135  // Get the invariant mass squared (s) from the kinematics object.
136  // Use the resPairAmpInt to find which mass-squared combination to use.
137  Double_t s(0.0);
138 
139  if (resPairAmpInt_ == 1) {
140  s = kinematics->getm23Sq();
141  } else if (resPairAmpInt_ == 2) {
142  s = kinematics->getm13Sq();
143  } else if (resPairAmpInt_ == 3) {
144  s = kinematics->getm12Sq();
145  }
146 
147  this->updatePropagator(s);
148 
149 }
150 
152 {
153  // Calculate the K-matrix propagator for the given s value.
154  // The K-matrix amplitude is given by
155  // T_i = sum_{ij} (I - iK*rho)^-1 * P_j, where P is the production K-matrix.
156  // i = index for the state (e.g. S-wave index = 0).
157  // Here, we only find the (I - iK*rho)^-1 matrix part.
158 
159  // Check if we have almost the same s value as before. If so, don't re-calculate
160  // the propagator nor any of the pole mass summation terms.
161  if (TMath::Abs(s - previousS_) < 1e-6*s) {
162  //cout<<"Already got propagator for s = "<<s<<endl;
163  return;
164  }
165 
166  if (parametersSet_ == kFALSE) {
167  //cerr<<"ERROR in LauKMatrixPropagator::updatePropagator. Parameters have not been set."<<endl;
168  return;
169  }
170 
171  // Calculate the denominator pole mass terms and Adler zero factor
172  this->calcPoleDenomVect(s);
173  this->updateAdlerZeroFactor(s);
174 
175  // Calculate the scattering K-matrix (real and symmetric)
176  this->calcScattKMatrix(s);
177  // Calculate the phase space density matrix, which is diagonal, but can be complex
178  // if the quantity s is below various threshold values (analytic continuation).
179  this->calcRhoMatrix(s);
180 
181  // Calculate K*rho (real and imaginary parts, since rho can be complex)
182  TMatrixD K_realRho(ScattKMatrix_);
183  K_realRho *= ReRhoMatrix_;
184  TMatrixD K_imagRho(ScattKMatrix_);
185  K_imagRho *= ImRhoMatrix_;
186 
187  // A = I + K*Imag(rho), B = -K*Real(Rho)
188  // Calculate C and D matrices such that (A + iB)*(C + iD) = I,
189  // ie. C + iD = (I - i K*rho)^-1, separated into real and imaginary parts.
190  TMatrixD A(IMatrix_);
191  A += K_imagRho;
192  TMatrixD B(zeroMatrix_);
193  B -= K_realRho;
194 
195  TMatrixD invA(TMatrixD::kInverted, A);
196  TMatrixD invA_B(invA);
197  invA_B *= B;
198  TMatrixD B_invA_B(B);
199  B_invA_B *= invA_B;
200 
201  TMatrixD invC(A);
202  invC += B_invA_B;
203 
204  // Set the real part of the propagator matrix ("C")
205  realProp_ = TMatrixD(TMatrixD::kInverted, invC);
206 
207  // Set the (negative) imaginary part of the propagator matrix ("-D")
208  TMatrixD BC(B);
209  BC *= realProp_;
210  negImagProp_ = TMatrixD(invA);
211  negImagProp_ *= BC;
212 
213  // Also calculate the production SVP term, since this uses Adler-zero parameters
214  // defined in the parameter file.
215  this->updateProdSVPTerm(s);
216 
217  // Finally, keep track of the value of s we just used.
218  previousS_ = s;
219 
220 }
221 
222 void LauKMatrixPropagator::setParameters(const TString& inputFile)
223 {
224 
225  // Read an input file that specifies the values of the coupling constants g_i for
226  // the given number of poles and their (bare) masses. Also provided are the f_{ab}
227  // slow-varying constants. The input file should also provide the Adler zero
228  // constants s_0, s_A and s_A0.
229 
230  // Format of input file:
231  // Indices (nChannels) of phase space channel types (defined in KMatrixChannels enum)
232  // Then the bare pole mass and coupling constants over all channels for each pole:
233  // m_{alpha} g_1^{alpha} g_2^{alpha} ... g_n^{alpha} where n = nChannels
234  // Repeat for N poles. Then define the f_{ab} scattering coefficients
235  // f_{11} f_{12} ... f_{1n}
236  // f_{21} f_{22} ... f_{2n} etc.. where n = nChannels
237  // Note that the K-matrix will be symmetrised (by definition), so it is sufficient
238  // to only have the upper half filled with sensible values.
239  // Then define the Adler-zero and slowly-varying part (SVP) constants
240  // m0^2 s0Scatt s0Prod sA sA0
241  // where m0^2 = mass-squared constant for scattering production f_{ab} numerator
242  parametersSet_ = kFALSE;
243 
244  cout<<"Initialising K-matrix propagator "<<name_<<" parameters from "<<inputFile.Data()<<endl;
245 
246  cout<<"nChannels = "<<nChannels_<<", nPoles = "<<nPoles_<<endl;
247 
248  LauTextFileParser readFile(inputFile);
249  readFile.processFile();
250 
251  // Get the first (non-comment) line
252  std::vector<std::string> channelTypeLine = readFile.getNextLine();
253  Int_t nTypes = static_cast<Int_t>(channelTypeLine.size());
254  if (nTypes != nChannels_) {
255  cerr<<"Error in LauKMatrixPropagator::setParameters. Input file "<<inputFile
256  <<" defines "<<nTypes<<" channels when "<<nChannels_<<" are expected"<<endl;
257  return;
258  }
259 
260  // Get the list of channel indices to specify what phase space factors should be used
261  // e.g. pipi, Kpi, eta eta', 4pi etc..
262  Int_t iChannel(0);
263  phaseSpaceTypes_.clear();
264  for (iChannel = 0; iChannel < nChannels_; iChannel++) {
265 
266  Int_t phaseSpaceInt = atoi(channelTypeLine[iChannel].c_str());
267  Bool_t checkChannel = this->checkPhaseSpaceType(phaseSpaceInt);
268 
269  if (checkChannel == kTRUE) {
270  cout<<"Adding phase space channel index "<<phaseSpaceInt
271  <<" to K-matrix propagator "<<name_<<endl;
272  phaseSpaceTypes_.push_back(phaseSpaceInt);
273  } else {
274  cerr<<"Phase space channel index "<<phaseSpaceInt
275  <<" should be between 1 and "<<LauKMatrixPropagator::TotChannels-1<<endl;
276  cerr<<"Stopping initialisation of the K-matrix propagator "<<name_<<endl;
277  cerr<<"Please correct the channel index on the first line of "<<inputFile<<endl;
278  return;
279  }
280 
281  }
282 
283  // Clear vector of pole masses^2 as well as the map of the coupling constant vectors.
284  mSqPoles_.clear(); gCouplings_.clear();
285 
286  Double_t poleMass(0.0), poleMassSq(0.0), couplingConst(0.0);
287  Int_t iPole(0);
288  Int_t nPoleNumbers(nChannels_+1);
289 
290  for (iPole = 0; iPole < nPoles_; iPole++) {
291 
292  // Read the next line of bare pole mass and its coupling constants
293  // over all channels (nChannels_)
294  std::vector<std::string> poleLine = readFile.getNextLine();
295  Int_t nPoleValues = static_cast<Int_t>(poleLine.size());
296  if (nPoleValues != nPoleNumbers) {
297  cerr<<"Error in LauKMatrixPropagator::setParameters. Expecting "
298  <<nPoleNumbers<<" numbers for the bare pole "<<iPole<<" line (mass and "
299  <<nChannels_<<" coupling constants). Found "<<nPoleValues<<" values instead."
300  <<endl;
301  return;
302  }
303 
304  poleMass = atof(poleLine[0].c_str());
305  poleMassSq = poleMass*poleMass;
306  LauParameter mPoleParam(poleMassSq);
307  mSqPoles_.push_back(mPoleParam);
308 
309  std::vector<LauParameter> couplingVector;
310 
311  cout<<"Added bare pole mass "<<poleMass<<" GeV for pole number "<<iPole+1<<endl;
312 
313  for (iChannel = 0; iChannel < nChannels_; iChannel++) {
314 
315  couplingConst = atof(poleLine[iChannel+1].c_str());
316  LauParameter couplingParam(couplingConst);
317  couplingVector.push_back(couplingParam);
318 
319  cout<<"Added coupling parameter g^"<<iPole+1<<"_"<<iChannel+1<<" = "<<couplingConst<<endl;
320 
321  }
322 
323  gCouplings_[iPole] = couplingVector;
324 
325  }
326 
327  // Scattering (slowly-varying part, or SVP) values
328 
329  fScattering_.clear();
330  Int_t jChannel(0);
331 
332  std::vector<std::string> scatteringLine = readFile.getNextLine();
333  Int_t nScatteringValues = static_cast<Int_t>(scatteringLine.size());
334  if (nScatteringValues != nChannels_) {
335  cerr<<"Error in LauKMatrixPropagator::setParameters. Expecting "
336  <<nChannels_<<" scattering SVP f constants instead of "<<nScatteringValues
337  <<" for the K-matrix row index "<<index_<<endl;
338  return;
339  }
340 
341  std::vector<LauParameter> scatteringVect;
342  for (jChannel = 0; jChannel < nChannels_; jChannel++) {
343 
344  Double_t fScattValue = atof(scatteringLine[jChannel].c_str());
345  LauParameter scatteringParam(fScattValue);
346  scatteringVect.push_back(scatteringParam);
347 
348  cout<<"Added scattering SVP f("<<index_+1<<","<<jChannel+1<<") = "<<fScattValue<<endl;
349 
350  }
351 
352  fScattering_[index_] = scatteringVect;
353 
354  // Now extract the constants for the "Adler-zero" terms
355  std::vector<std::string> constLine = readFile.getNextLine();
356  Int_t nConstValues = static_cast<Int_t>(constLine.size());
357  if (nConstValues != 5) {
358  cerr<<"Error in LauKMatrixPropagator::setParameters. Expecting 5 values for the Adler-zero constants:"
359  <<" m0^2 s0Scatt s0Prod sA sA0"<<endl;
360  return;
361  }
362 
363  Double_t mSq0Value = atof(constLine[0].c_str());
364  Double_t s0ScattValue = atof(constLine[1].c_str());
365  Double_t s0ProdValue = atof(constLine[2].c_str());
366  Double_t sAValue = atof(constLine[3].c_str());
367  Double_t sA0Value = atof(constLine[4].c_str());
368 
369  cout<<"Adler zero constants:"<<endl;
370  cout<<"m0Sq = "<<mSq0Value<<", s0Scattering = "<<s0ScattValue<<", s0Production = "
371  <<s0ProdValue<<", sA = "<<sAValue<<" and sA0 = "<<sA0Value<<endl;
372 
373  mSq0_ = LauParameter("mSq0", mSq0Value);
374  s0Scatt_ = LauParameter("s0Scatt", s0ScattValue);
375  s0Prod_ = LauParameter("s0Prod", s0ProdValue);
376  sA_ = LauParameter("sA", sAValue);
377  sA0_ = LauParameter("sA0", sA0Value);
378  sAConst_ = 0.5*sAValue*LauConstants::mPiSq;
379 
380  // Identity and null matrices
381  IMatrix_.Clear();
382  IMatrix_.ResizeTo(nChannels_, nChannels_);
383  for (iChannel = 0; iChannel < nChannels_; iChannel++) {
384  IMatrix_[iChannel][iChannel] = 1.0;
385  }
386 
387  zeroMatrix_.Clear();
388  zeroMatrix_.ResizeTo(nChannels_, nChannels_);
389 
390  // Real K matrix
391  ScattKMatrix_.Clear();
392  ScattKMatrix_.ResizeTo(nChannels_, nChannels_);
393 
394  // Real and imaginary propagator matrices
395  realProp_.Clear(); negImagProp_.Clear();
396  realProp_.ResizeTo(nChannels_, nChannels_);
397  negImagProp_.ResizeTo(nChannels_, nChannels_);
398 
399  // Phase space matrices
400  ReRhoMatrix_.Clear(); ImRhoMatrix_.Clear();
401  ReRhoMatrix_.ResizeTo(nChannels_, nChannels_);
402  ImRhoMatrix_.ResizeTo(nChannels_, nChannels_);
403 
404  // Square-root phase space matrices
405  ReSqrtRhoMatrix_.Clear(); ImSqrtRhoMatrix_.Clear();
406  ReSqrtRhoMatrix_.ResizeTo(nChannels_, nChannels_);
407  ImSqrtRhoMatrix_.ResizeTo(nChannels_, nChannels_);
408 
409  // T matrices
410  ReTMatrix_.Clear(); ImTMatrix_.Clear();
411  ReTMatrix_.ResizeTo(nChannels_, nChannels_);
412  ImTMatrix_.ResizeTo(nChannels_, nChannels_);
413 
414  // All required parameters have been set
415  parametersSet_ = kTRUE;
416 
417  cout<<"Finished initialising K-matrix propagator "<<name_<<endl;
418 
419 }
420 
422 {
423 
424  // Calculate the scattering K-matrix for the given value of s.
425  // We need to obtain the complete matrix (not just one row/column)
426  // to get the correct inverted (I - i K rho) terms for the propagator.
427 
428  if (verbose_) {cout<<"Within calcScattKMatrix for s = "<<s<<endl;}
429 
430  // Initialise the K matrix to zero
431  ScattKMatrix_.Zero();
432 
433  Int_t iChannel(0), jChannel(0), iPole(0);
434 
435  // The pole denominator 1/(m^2 - s) terms should already be calculated
436  // by the calcPoleDenomVect() function. These same terms are also
437  // used for calculating the production K-matrix elements.
438 
439  // Calculate the "slowly-varying part" (SVP), e.g. (1 GeV - s0)/(s - s0)
440  this->updateScattSVPTerm(s);
441 
442  // Now loop over iChannel, jChannel to calculate Kij = Kji.
443  for (iChannel = 0; iChannel < nChannels_; iChannel++) {
444 
445  std::vector<LauParameter> SVPVect(0);
446 
447  KMatrixParamMap::iterator SVPIter = fScattering_.find(iChannel);
448  if (SVPIter != fScattering_.end()) {SVPVect = SVPIter->second;}
449 
450  Int_t SVPVectSize = static_cast<Int_t>(SVPVect.size());
451 
452  // Scattering matrix is real and symmetric. Start j loop from i.
453  for (jChannel = iChannel; jChannel < nChannels_; jChannel++) {
454 
455  Double_t Kij(0.0);
456 
457  // Calculate pole mass summation term
458  for (iPole = 0; iPole < nPoles_; iPole++) {
459 
460  KMatrixParamMap::iterator couplingIter = gCouplings_.find(iPole);
461 
462  std::vector<LauParameter> couplingVect = couplingIter->second;
463 
464  Double_t g_i = couplingVect[iChannel].value();
465  Double_t g_j = couplingVect[jChannel].value();
466 
467  Kij += poleDenomVect_[iPole]*g_i*g_j;
468  if (verbose_) {cout<<"1: Kij for i = "<<iChannel<<", j = "<<jChannel<<" = "<<Kij<<endl;}
469 
470  }
471 
472  if (SVPVectSize != 0 && jChannel < SVPVectSize) {
473  Double_t fij = SVPVect[jChannel].value();
474  Kij += fij*scattSVP_;
475  }
476 
477  Kij *= adlerZeroFactor_;
478  if (verbose_) {cout<<"2: Kij for i = "<<iChannel<<", j = "<<jChannel<<" = "<<Kij<<endl;}
479 
480  // Assign the TMatrix (i,j) element to the variable Kij and Kji (symmetry)
481  ScattKMatrix_(iChannel, jChannel) = Kij;
482  ScattKMatrix_(jChannel, iChannel) = Kij;
483 
484  } // j loop
485 
486  } // i loop
487 
488 }
489 
491 {
492  // Calculate the 1/(m_pole^2 - s) terms for the scattering
493  // and production K-matrix formulae.
494  poleDenomVect_.clear();
495  Int_t iPole(0);
496  for (iPole = 0; iPole < nPoles_; iPole++) {
497 
498  Double_t poleTerm = mSqPoles_[iPole] - s;
499  Double_t invPoleTerm(0.0);
500  if (TMath::Abs(poleTerm) > 1.0e-6) {invPoleTerm = 1.0/poleTerm;}
501 
502  poleDenomVect_.push_back(invPoleTerm);
503 
504  }
505 
506 }
507 
508 Double_t LauKMatrixPropagator::getPoleDenomTerm(Int_t poleIndex) const
509 {
510 
511  if (parametersSet_ == kFALSE) {return 0.0;}
512 
513  Double_t poleDenom(0.0);
514  poleDenom = poleDenomVect_[poleIndex];
515  return poleDenom;
516 
517 }
518 
519 Double_t LauKMatrixPropagator::getCouplingConstant(Int_t poleIndex, Int_t channelIndex) const
520 {
521 
522  if (parametersSet_ == kFALSE) {return 0.0;}
523 
524  Double_t couplingConst(0.0);
525  KMatrixParamMap::const_iterator couplingIter = gCouplings_.find(poleIndex);
526  std::vector<LauParameter> couplingVect = couplingIter->second;
527  couplingConst = couplingVect[channelIndex].value();
528 
529  return couplingConst;
530 
531 }
532 
533 Double_t LauKMatrixPropagator::calcSVPTerm(Double_t s, Double_t s0) const
534 {
535 
536  if (parametersSet_ == kFALSE) {return 0.0;}
537 
538  // Calculate the "slowly-varying part" (SVP)
539  Double_t result(0.0);
540  Double_t deltaS = s - s0;
541  if (TMath::Abs(deltaS) > 1.0e-6) {
542  result = (mSq0_.value() - s0)/deltaS;
543  }
544 
545  return result;
546 
547 }
548 
550 {
551  // Update the scattering "slowly-varying part" (SVP)
552  Double_t s0Scatt = s0Scatt_.value();
553  scattSVP_ = this->calcSVPTerm(s, s0Scatt);
554 }
555 
557 {
558  // Update the production "slowly-varying part" (SVP)
559  Double_t s0Prod = s0Prod_.value();
560  prodSVP_ = this->calcSVPTerm(s, s0Prod);
561 }
562 
564 {
565 
566  // Calculate the multiplicative factor containing various Adler zero
567  // constants.
568  adlerZeroFactor_ = 0.0;
569 
570  Double_t sA0Val = sA0_.value();
571  Double_t deltaS = s - sA0Val;
572  if (TMath::Abs(deltaS) > 1e-6) {
573  adlerZeroFactor_ = (s - sAConst_)*(1.0 - sA0Val)/deltaS;
574  }
575 
576 }
577 
579 {
580 
581  // Calculate the real and imaginary part of the phase space density
582  // diagonal matrix for the given invariant mass squared quantity, s.
583  // The matrix can be complex if s is below threshold (so that
584  // the amplitude continues analytically).
585 
586  // Initialise all entries to zero
587  ReRhoMatrix_.Zero(); ImRhoMatrix_.Zero();
588 
589  LauComplex rho(0.0, 0.0);
590  Int_t phaseSpaceIndex(0);
591 
592  for (Int_t iChannel (0); iChannel < nChannels_; ++iChannel) {
593 
594  phaseSpaceIndex = phaseSpaceTypes_[iChannel];
595 
596  if (phaseSpaceIndex == LauKMatrixPropagator::PiPi) {
597  rho = this->calcPiPiRho(s);
598  } else if (phaseSpaceIndex == LauKMatrixPropagator::KK) {
599  rho = this->calcKKRho(s);
600  } else if (phaseSpaceIndex == LauKMatrixPropagator::FourPi) {
601  rho = this->calcFourPiRho(s);
602  } else if (phaseSpaceIndex == LauKMatrixPropagator::EtaEta) {
603  rho = this->calcEtaEtaRho(s);
604  } else if (phaseSpaceIndex == LauKMatrixPropagator::EtaEtaP) {
605  rho = this->calcEtaEtaPRho(s);
606  } else if (phaseSpaceIndex == LauKMatrixPropagator::KPi) {
607  rho = this->calcKPiRho(s);
608  } else if (phaseSpaceIndex == LauKMatrixPropagator::KEtaP) {
609  rho = this->calcKEtaPRho(s);
610  } else if (phaseSpaceIndex == LauKMatrixPropagator::KThreePi) {
611  rho = this->calcKThreePiRho(s);
612  }
613 
614  if (verbose_) {
615  cout<<"ReRhoMatrix("<<iChannel<<", "<<iChannel<<") = "<<rho.re()<<endl;
616  cout<<"ImRhoMatrix("<<iChannel<<", "<<iChannel<<") = "<<rho.im()<<endl;
617  }
618 
619  ReRhoMatrix_(iChannel, iChannel) = rho.re();
620  ImRhoMatrix_(iChannel, iChannel) = rho.im();
621 
622  }
623 
624 }
625 
627 {
628  // Calculate the pipi phase space factor
629  LauComplex rho(0.0, 0.0);
630  if (TMath::Abs(s) < 1e-10) {return rho;}
631 
632  Double_t sqrtTerm = (-m2piSq_/s) + 1.0;
633  if (sqrtTerm < 0.0) {
634  rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
635  } else {
636  rho.setRealPart( TMath::Sqrt(sqrtTerm) );
637  }
638 
639  return rho;
640 }
641 
643 {
644  // Calculate the KK phase space factor
645  LauComplex rho(0.0, 0.0);
646  if (TMath::Abs(s) < 1e-10) {return rho;}
647 
648  Double_t sqrtTerm = (-m2KSq_/s) + 1.0;
649  if (sqrtTerm < 0.0) {
650  rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
651  } else {
652  rho.setRealPart( TMath::Sqrt(sqrtTerm) );
653  }
654 
655  return rho;
656 }
657 
659 {
660  // Calculate the 4pi phase space factor. This uses a 6th-order polynomial
661  // parameterisation that approximates the multi-body phase space double integral
662  // defined in Eq 4 of the A&S paper hep-ph/0204328. This form agrees with the
663  // BaBar model (another 6th order polynomial from s^4 down to 1/s^2), but avoids the
664  // exponential increase at small values of s (~< 0.1) arising from 1/s and 1/s^2.
665  // Eq 4 is evaluated for each value of s by assuming incremental steps of 1e-3 GeV^2
666  // for s1 and s2, the invariant energy squared of each of the di-pion states,
667  // with the integration limits of s1 = (2*mpi)^2 to (sqrt(s) - 2*mpi)^2 and
668  // s2 = (2*mpi)^2 to (sqrt(s) - sqrt(s1))^2. The mass M of the rho is taken to be
669  // 0.775 GeV and the energy-dependent width of the 4pi system
670  // Gamma(s) = gamma_0*rho1^3(s), where rho1 = sqrt(1.0 - 4*mpiSq/s) and gamma_0 is
671  // the "width" of the 4pi state at s = 1, which is taken to be 0.3 GeV
672  // (~75% of the total width from PDG estimates of the f0(1370) -> 4pi state).
673  // The normalisation term rho_0 is found by ensuring that the phase space integral
674  // at s = 1 is equal to sqrt(1.0 - 16*mpiSq/s). Note that the exponent for this
675  // factor in hep-ph/0204328 is wrong; it should be 0.5, i.e. sqrt, not n = 1 to 5.
676  // Plotting the value of this double integral as a function of s can then be fitted
677  // to a 6th-order polynomial (for s < 1), which is the result used below
678 
679  LauComplex rho(0.0, 0.0);
680  if (TMath::Abs(s) < 1e-10) {return rho;}
681 
682  if (s <= 1.0) {
683  Double_t rhoTerm = ((1.07885*s + 0.13655)*s - 0.29744)*s - 0.20840;
684  rhoTerm = ((rhoTerm*s + 0.13851)*s - 0.01933)*s + 0.00051;
685  // For some values of s (below 2*mpi), this term is a very small
686  // negative number. Check for this and set the rho term to zero.
687  if (rhoTerm < 0.0) {rhoTerm = 0.0;}
688  rho.setRealPart( rhoTerm );
689  } else {
690  rho.setRealPart( TMath::Sqrt(1.0 - (fourPiFactor1_/s)) );
691  }
692 
693  return rho;
694 }
695 
697 {
698  // Calculate the eta-eta phase space factor
699  LauComplex rho(0.0, 0.0);
700  if (TMath::Abs(s) < 1e-10) {return rho;}
701 
702  Double_t sqrtTerm = (-m2EtaSq_/s) + 1.0;
703  if (sqrtTerm < 0.0) {
704  rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
705  } else {
706  rho.setRealPart( TMath::Sqrt(sqrtTerm) );
707  }
708 
709  return rho;
710 }
711 
713 {
714  // Calculate the eta-eta' phase space factor. Note that the
715  // mass difference term m_eta - m_eta' is not included,
716  // since this corresponds to a "t or u-channel crossing",
717  // which means that we cannot simply analytically continue
718  // this part of the phase space factor below threshold, which
719  // we can do for s-channel contributions. This is actually an
720  // unsolved problem, e.g. see Guo et al 1409.8652, and
721  // Danilkin et al 1409.7708. Anisovich and Sarantsev in
722  // hep-ph/0204328 "solve" this issue by setting the mass
723  // difference term to unity, which is what we do here...
724 
725  LauComplex rho(0.0, 0.0);
726  if (TMath::Abs(s) < 1e-10) {return rho;}
727 
728  Double_t sqrtTerm = (-mEtaEtaPSumSq_/s) + 1.0;
729  if (sqrtTerm < 0.0) {
730  rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
731  } else {
732  rho.setRealPart( TMath::Sqrt(sqrtTerm) );
733  }
734 
735  return rho;
736 }
737 
738 
740 {
741  // Calculate the K-pi phase space factor
742  LauComplex rho(0.0, 0.0);
743  if (TMath::Abs(s) < 1e-10) {return rho;}
744 
745  Double_t sqrtTerm1 = (-mKpiSumSq_/s) + 1.0;
746  Double_t sqrtTerm2 = (-mKpiDiffSq_/s) + 1.0;
747  Double_t sqrtTerm = sqrtTerm1*sqrtTerm2;
748  if (sqrtTerm < 0.0) {
749  rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
750  } else {
751  rho.setRealPart( TMath::Sqrt(sqrtTerm) );
752  }
753 
754  return rho;
755 }
756 
758 {
759  // Calculate the K-eta' phase space factor
760  LauComplex rho(0.0, 0.0);
761  if (TMath::Abs(s) < 1e-10) {return rho;}
762 
763  Double_t sqrtTerm1 = (-mKEtaPSumSq_/s) + 1.0;
764  Double_t sqrtTerm2 = (-mKEtaPDiffSq_/s) + 1.0;
765  Double_t sqrtTerm = sqrtTerm1*sqrtTerm2;
766  if (sqrtTerm < 0.0) {
767  rho.setImagPart( TMath::Sqrt(-sqrtTerm) );
768  } else {
769  rho.setRealPart( TMath::Sqrt(sqrtTerm) );
770  }
771 
772  return rho;
773 }
774 
776 {
777  // Calculate the Kpipipi + multimeson phase space factor.
778  // Use the simplest definition in hep-ph/9705401 (Eq 14), which is the form
779  // used for the rest of that paper (thankfully, the amplitude does not depend
780  // significantly on the form used for the K3pi phase space factor).
781  LauComplex rho(0.0, 0.0);
782  if (TMath::Abs(s) < 1e-10) {return rho;}
783 
784  if (s < 1.44) {
785 
786  Double_t powerTerm = (-mK3piDiffSq_/s) + 1.0;
787  if (powerTerm < 0.0) {
788  rho.setImagPart( k3piFactor_*TMath::Power(-powerTerm, 2.5) );
789  } else {
790  rho.setRealPart( k3piFactor_*TMath::Power(powerTerm, 2.5) );
791  }
792 
793  } else {
794  rho.setRealPart( 1.0 );
795  }
796 
797  return rho;
798 }
799 
800 Bool_t LauKMatrixPropagator::checkPhaseSpaceType(Int_t phaseSpaceInt) const
801 {
802  Bool_t passed(kFALSE);
803 
804  if (phaseSpaceInt >= 1 && phaseSpaceInt < LauKMatrixPropagator::TotChannels) {
805  passed = kTRUE;
806  }
807 
808  return passed;
809 }
810 
812 {
813 
814  // Get the complex (unitary) transition amplitude T for the given channel
815  LauComplex TAmp(0.0, 0.0);
816 
817  channel -= 1;
818  if (channel < 0 || channel >= nChannels_) {return TAmp;}
819 
820  this->getTMatrix(s);
821 
822  TAmp.setRealPart(ReTMatrix_[index_][channel]);
823  TAmp.setImagPart(ImTMatrix_[index_][channel]);
824 
825  return TAmp;
826 
827 }
828 
830 {
831 
832  // Get the complex (unitary) transition amplitude T for the given channel
833  LauComplex rho(0.0, 0.0);
834 
835  channel -= 1;
836  if (channel < 0 || channel >= nChannels_) {return rho;}
837 
838  // If s has changed from the previous value, recalculate rho
839  if (TMath::Abs(s - previousS_) > 1e-6*s) {
840  this->calcRhoMatrix(s);
841  }
842 
843  rho.setRealPart(ReRhoMatrix_[channel][channel]);
844  rho.setImagPart(ImRhoMatrix_[channel][channel]);
845 
846  return rho;
847 
848 }
849 
851 
852  // Find the unitary T matrix, where T = [sqrt(rho)]^{*} T_hat sqrt(rho),
853  // and T_hat = (I - i K rho)^-1 * K is the Lorentz-invariant T matrix,
854  // which has phase-space factors included (rho). This function is not
855  // needed to calculate the K-matrix amplitudes, but allows us
856  // to check the variation of T as a function of s (kinematics)
857 
858  if (!kinematics) {return;}
859 
860  // Get the invariant mass squared (s) from the kinematics object.
861  // Use the resPairAmpInt to find which mass-squared combination to use.
862  Double_t s(0.0);
863 
864  if (resPairAmpInt_ == 1) {
865  s = kinematics->getm23Sq();
866  } else if (resPairAmpInt_ == 2) {
867  s = kinematics->getm13Sq();
868  } else if (resPairAmpInt_ == 3) {
869  s = kinematics->getm12Sq();
870  }
871 
872  this->getTMatrix(s);
873 
874 }
875 
876 
878 {
879 
880  // Find the unitary transition T matrix, where
881  // T = [sqrt(rho)]^{*} T_hat sqrt(rho), and
882  // T_hat = (I - i K rho)^-1 * K is the Lorentz-invariant T matrix,
883  // which has phase-space factors included (rho). Note that the first
884  // sqrt of the rho matrix is complex conjugated.
885 
886  // This function is not needed to calculate the K-matrix amplitudes, but
887  // allows us to check the variation of T as a function of s (kinematics)
888 
889  // Initialse the real and imaginary parts of the T matrix to zero
890  ReTMatrix_.Zero(); ImTMatrix_.Zero();
891 
892  if (parametersSet_ == kFALSE) {return;}
893 
894  // Update K, rho and the propagator (I - i K rho)^-1
895  this->updatePropagator(s);
896 
897  // Find the real and imaginary T_hat matrices
898  TMatrixD THatReal = realProp_*ScattKMatrix_;
899  TMatrixD THatImag(zeroMatrix_);
900  THatImag -= negImagProp_*ScattKMatrix_;
901 
902  // Find the square-root of the phase space matrix
903  this->getSqrtRhoMatrix();
904 
905  // Let sqrt(rho) = A + iB and T_hat = C + iD
906  // => T = A(CA-DB) + B(DA+CB) + i[A(DA+CB) + B(DB-CA)]
907  TMatrixD CA(THatReal);
908  CA *= ReSqrtRhoMatrix_;
909  TMatrixD DA(THatImag);
910  DA *= ReSqrtRhoMatrix_;
911  TMatrixD CB(THatReal);
912  CB *= ImSqrtRhoMatrix_;
913  TMatrixD DB(THatImag);
914  DB *= ImSqrtRhoMatrix_;
915 
916  TMatrixD CAmDB(CA);
917  CAmDB -= DB;
918  TMatrixD DApCB(DA);
919  DApCB += CB;
920  TMatrixD DBmCA(DB);
921  DBmCA -= CA;
922 
923  // Find the real and imaginary parts of the transition matrix T
926 
927 }
928 
930 {
931 
932  // Find the square root of the (current) phase space matrix so that
933  // we can find T = [sqrt(rho)}^{*} T_hat sqrt(rho), where T_hat is the
934  // Lorentz-invariant T matrix = (I - i K rho)^-1 * K; note that the first
935  // sqrt of rho matrix is complex conjugated
936 
937  // If rho = rho_i + i rho_r = a + i b, then sqrt(rho) = c + i d, where
938  // c = sqrt(0.5*(r+a)) and d = sqrt(0.5(r-a)), where r = sqrt(a^2 + b^2).
939  // Since rho is diagonal, then the square root of rho will also be diagonal,
940  // with its real and imaginary matrix elements equal to c and d, respectively
941 
942  // Initialise the real and imaginary parts of the square root of
943  // the rho matrix to zero
944  ReSqrtRhoMatrix_.Zero(); ImSqrtRhoMatrix_.Zero();
945 
946  for (Int_t iChannel (0); iChannel < nChannels_; ++iChannel) {
947 
948  Double_t realRho = ReRhoMatrix_[iChannel][iChannel];
949  Double_t imagRho = ImRhoMatrix_[iChannel][iChannel];
950 
951  Double_t rhoMag = sqrt(realRho*realRho + imagRho*imagRho);
952  Double_t rhoSum = rhoMag + realRho;
953  Double_t rhoDiff = rhoMag - realRho;
954 
955  Double_t reSqrtRho(0.0), imSqrtRho(0.0);
956  if (rhoSum > 0.0) {reSqrtRho = sqrt(0.5*rhoSum);}
957  if (rhoDiff > 0.0) {imSqrtRho = sqrt(0.5*rhoDiff);}
958 
959  ReSqrtRhoMatrix_[iChannel][iChannel] = reSqrtRho;
960  ImSqrtRhoMatrix_[iChannel][iChannel] = imSqrtRho;
961 
962  }
963 
964 
965 }
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.
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.
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
std::vector< std::string > getNextLine()
Retrieve the next line.
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.
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.
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: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
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 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.
Class for defining the fit parameter objects.
Definition: LauParameter.hh:34
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
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 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.
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.
LauComplex getTransitionAmp(Double_t s, Int_t channel)
Get the unitary transition amplitude for the given channel.
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
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;.
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.