ROOTPWA
massDependence.cc
Go to the documentation of this file.
1 
2 //
3 // Copyright 2010
4 //
5 // This file is part of rootpwa
6 //
7 // rootpwa is free software: you can redistribute it and/or modify
8 // it under the terms of the GNU General Public License as published by
9 // the Free Software Foundation, either version 3 of the License, or
10 // (at your option) any later version.
11 //
12 // rootpwa is distributed in the hope that it will be useful,
13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 // GNU General Public License for more details.
16 //
17 // You should have received a copy of the GNU General Public License
18 // along with rootpwa. If not, see <http://www.gnu.org/licenses/>.
19 //
21 //-------------------------------------------------------------------------
22 // File and Version Information:
23 // $Rev:: $: revision of last commit
24 // $Author:: $: author of last commit
25 // $Date:: $: date of last commit
26 //
27 // Description:
28 // functor class hierarchy for mass-dependent part of the amplitude
29 //
30 //
31 // Author List:
32 // Boris Grube TUM (original author)
33 //
34 //
35 //-------------------------------------------------------------------------
36 
37 
38 #include <boost/numeric/ublas/io.hpp>
39 
40 #include "physUtils.hpp"
41 #include "isobarDecayVertex.h"
42 #include "particleDataTable.h"
43 #include "massDependence.h"
44 
45 
46 using namespace std;
47 using namespace boost::numeric::ublas;
48 using namespace rpwa;
49 
50 
52 bool massDependence::_debug = false;
53 
54 
55 ostream&
56 massDependence::print(ostream& out) const
57 {
58  out << name();
59  return out;
60 }
61 
62 
64 complex<double>
65 flatMassDependence::amp(const isobarDecayVertex&)
66 {
67  if (_debug)
68  printDebug << name() << " = 1" << endl;
69  return 1;
70 }
71 
72 
74 complex<double>
75 relativisticBreitWigner::amp(const isobarDecayVertex& v)
76 {
77  const particlePtr& parent = v.parent();
78 
79  // get Breit-Wigner parameters
80  const double M = parent->lzVec().M(); // parent mass
81  const double m1 = v.daughter1()->lzVec().M(); // daughter 1 mass
82  const double m2 = v.daughter2()->lzVec().M(); // daughter 2 mass
83  const double q = breakupMomentum(M, m1, m2);
84  const double M0 = parent->mass(); // resonance peak position
85  const double q02 = breakupMomentumSquared(M0, m1, m2, true);
86  // !NOTE! the following is incorrect but this is how it was done in PWA2000
87  const double q0 = sqrt(fabs(q02));
88  const double Gamma0 = parent->width(); // resonance peak width
89  const unsigned int L = v.L();
90 
91  const complex<double> bw = breitWigner(M, M0, Gamma0, L, q, q0);
92  if (_debug)
93  printDebug << name() << "(m = " << maxPrecision(M) << " GeV/c^2, m_0 = " << maxPrecision(M0)
94  << " GeV/c^2, Gamma_0 = " << maxPrecision(Gamma0) << " GeV/c^2, L = " << spinQn(L)
95  << ", q = " << maxPrecision(q) << " GeV/c, q0 = "
96  << maxPrecision(q0) << " GeV/c) = " << maxPrecisionDouble(bw) << endl;
97  return bw;
98 }
99 
100 
102 complex<double>
103 constWidthBreitWigner::amp(const isobarDecayVertex& v)
104 {
105  const particlePtr& parent = v.parent();
106 
107  // get Breit-Wigner parameters
108  const double M = parent->lzVec().M(); // parent mass
109  const double M0 = parent->mass(); // resonance peak position
110  const double Gamma0 = parent->width(); // resonance peak width
111 
112  // A / (B - iA) = (A / (B^2 + A^2)) * (B + iA)
113  const double A = M0 * Gamma0;
114  const double B = M0 * M0 - M * M;
115  const complex<double> bw = (A / (B * B + A * A)) * complex<double>(B, A);
116  // const complex<double> bw = (M0 * Gamma0) / (M0 * M0 - M * M - imag * M0 * Gamma0);
117  if (_debug)
118  printDebug << name() << "(m = " << maxPrecision(M) << " GeV/c^2, m_0 = " << maxPrecision(M0)
119  << " GeV/c^2, Gamma_0 = " << maxPrecision(Gamma0) << " GeV/c^2) = "
120  << maxPrecisionDouble(bw) << endl;
121  return bw;
122 }
123 
124 
126 complex<double>
127 rhoBreitWigner::amp(const isobarDecayVertex& v)
128 {
129  const particlePtr& parent = v.parent();
130 
131  // get Breit-Wigner parameters
132  const double M = parent->lzVec().M(); // parent mass
133  const double m1 = v.daughter1()->lzVec().M(); // daughter 1 mass
134  const double m2 = v.daughter2()->lzVec().M(); // daughter 2 mass
135  const double q2 = breakupMomentumSquared(M, m1, m2);
136  const double q = sqrt(q2);
137  const double M0 = parent->mass(); // resonance peak position
138  const double q02 = breakupMomentumSquared(M0, m1, m2);
139  const double q0 = sqrt(q02);
140  const double Gamma0 = parent->width(); // resonance peak width
141 
142  const double F = 2 * q2 / (q02 + q2);
143  const double Gamma = Gamma0 * (M0 / M) * (q / q0) * F;
144  // in the original publication the width reads
145  // Gamma = Gamma0 * (q / q0) * F
146 
147  // A / (B - iC) = (A / (B^2 + C^2)) * (B + iC)
148  const double A = M0 * Gamma0 * sqrt(F);
149  // in the original publication A reads
150  // A = sqrt(M0 * Gamma0 * (m / q0) * F)
151  const double B = M0 * M0 - M * M;
152  const double C = M0 * Gamma;
153  const complex<double> bw = (A / (B * B + C * C)) * std::complex<double>(B, C);
154  // return (M0 * Gamma0 * sqrt(F)) / (M0 * M0 - M * M - imag * M0 * Gamma);
155  if (_debug)
156  printDebug << name() << "(m = " << maxPrecision(M) << " GeV/c^2, m_0 = " << maxPrecision(M0)
157  << " GeV/c^2, Gamma_0 = " << maxPrecision(Gamma0) << " GeV/c^2, "
158  << "q = " << maxPrecision(q) << " GeV/c, q0 = " << maxPrecision(q0) << " GeV/c) "
159  << "= " << maxPrecisionDouble(bw) << endl;
160  return bw;
161 }
162 
163 
165 complex<double>
166 f0980BreitWigner::amp(const isobarDecayVertex& v)
167 {
168  const particlePtr& parent = v.parent();
169 
170  // get Breit-Wigner parameters
171  const double M = parent->lzVec().M(); // parent mass
172  const double m1 = v.daughter1()->lzVec().M(); // daughter 1 mass
173  const double m2 = v.daughter2()->lzVec().M(); // daughter 2 mass
174  const double q = breakupMomentum(M, m1, m2);
175  const double M0 = parent->mass(); // resonance peak position
176  const double q0 = breakupMomentum(M0, m1, m2);
177  const double Gamma0 = parent->width(); // resonance peak width
178 
179  const double Gamma = Gamma0 * (q / q0);
180 
181  // A / (B - iC) = (A / (B^2 + C^2)) * (B + iC)
182  const double C = M0 * Gamma;
183  const double A = C * M / q;
184  const double B = M0 * M0 - M * M;
185  const complex<double> bw = (A / (B * B + C * C)) * std::complex<double>(B, C);
186  // return ((M0 * Gamma0 * M / q) / (M0 * M0 - M * M - imag * M0 * Gamma);
187  if (_debug)
188  printDebug << name() << "(m = " << maxPrecision(M) << " GeV/c^2, m_0 = " << maxPrecision(M0)
189  << " GeV/c^2, Gamma_0 = " << maxPrecision(Gamma0) << " GeV/c^2, "
190  << "q = " << maxPrecision(q) << " GeV/c, q0 = " << maxPrecision(q0) << " GeV/c) "
191  << "= " << maxPrecisionDouble(bw) << endl;
192  return bw;
193 }
194 
195 
197 piPiSWaveAuMorganPenningtonM::piPiSWaveAuMorganPenningtonM()
198  : massDependence(),
199  _T (2, 2),
200  _a (2, matrix<complex<double> >(2, 2)),
201  _c (5, matrix<complex<double> >(2, 2)),
202  _sP (1, 2),
203  _vesSheet(0)
204 {
205  const double f[2] = {0.1968, -0.0154}; // AMP Table 1, M solution: f_1^1 and f_2^1
206 
207  _a[0](0, 0) = 0.1131; // AMP Table 1, M solution: f_2^2
208  _a[0](0, 1) = 0.0150; // AMP Table 1, M solution: f_1^3
209  _a[0](1, 0) = 0.0150; // AMP Table 1, M solution: f_1^3
210  _a[0](1, 1) = -0.3216; // AMP Table 1, M solution: f_2^3
211  _a[1](0, 0) = f[0] * f[0];
212  _a[1](0, 1) = f[0] * f[1];
213  _a[1](1, 0) = f[1] * f[0];
214  _a[1](1, 1) = f[1] * f[1];
215 
216  _c[0](0, 0) = 0.0337; // AMP Table 1, M solution: c_11^0
217  _c[1](0, 0) = -0.3185; // AMP Table 1, M solution: c_11^1
218  _c[2](0, 0) = -0.0942; // AMP Table 1, M solution: c_11^2
219  _c[3](0, 0) = -0.5927; // AMP Table 1, M solution: c_11^3
220  _c[4](0, 0) = 0.1957; // AMP Table 1, M solution: c_11^4
221  _c[0](0, 1) = _c[0](1, 0) = -0.2826; // AMP Table 1, M solution: c_12^0
222  _c[1](0, 1) = _c[1](1, 0) = 0.0918; // AMP Table 1, M solution: c_12^1
223  _c[2](0, 1) = _c[2](1, 0) = 0.1669; // AMP Table 1, M solution: c_12^2
224  _c[3](0, 1) = _c[3](1, 0) = -0.2082; // AMP Table 1, M solution: c_12^3
225  _c[4](0, 1) = _c[4](1, 0) = -0.1386; // AMP Table 1, M solution: c_12^4
226  _c[0](1, 1) = 0.3010; // AMP Table 1, M solution: c_22^0
227  _c[1](1, 1) = -0.5140; // AMP Table 1, M solution: c_12^1
228  _c[2](1, 1) = 0.1176; // AMP Table 1, M solution: c_12^2
229  _c[3](1, 1) = 0.5204; // AMP Table 1, M solution: c_12^3
230  _c[4](1, 1) = -0.3977; // AMP Table 1, M solution: c_12^4
231 
232  _sP(0, 0) = -0.0074; // AMP Table 1, M solution: s_0
233  _sP(0, 1) = 0.9828; // AMP Table 1, M solution: s_1
234 
235  particleDataTable& pdt = particleDataTable::instance();
236  const string partList[] = {"pi+", "pi0", "K+", "K0"};
237  for (unsigned int i = 0; i < sizeof(partList) / sizeof(partList[0]); ++i)
238  if (not pdt.isInTable(partList[i])) {
239  printErr << "cannot find particle " << partList[i] << " in particle data table. "
240  << "aborting." << endl;
241  throw;
242  }
243  _piChargedMass = pdt.entry("pi+")->mass();
244  _piNeutralMass = pdt.entry("pi0")->mass();
245  _kaonChargedMass = pdt.entry("K+" )->mass();
246  _kaonNeutralMass = pdt.entry("K0" )->mass();
247  _kaonMeanMass = (_kaonChargedMass + _kaonNeutralMass) / 2;
248 }
249 
250 
251 complex<double>
252 piPiSWaveAuMorganPenningtonM::amp(const isobarDecayVertex& v)
253 {
254  const complex<double> imag(0, 1);
255 
256  double mass = v.parent()->lzVec().M();
257  double s = mass * mass;
258  if (fabs(s - _sP(0, 1)) < 1e-6) {
259  mass += 1e-6;
260  s = mass * mass;
261  }
262 
263  const complex<double> qPiPi = breakupMomentumComplex(mass, _piChargedMass, _piChargedMass );
264  const complex<double> qPi0Pi0 = breakupMomentumComplex(mass, _piNeutralMass, _piNeutralMass );
265  const complex<double> qKK = breakupMomentumComplex(mass, _kaonChargedMass, _kaonChargedMass);
266  const complex<double> qK0K0 = breakupMomentumComplex(mass, _kaonNeutralMass, _kaonNeutralMass);
267  complex<double> qKmKm = breakupMomentumComplex(mass, _kaonMeanMass, _kaonMeanMass );
268 
269  matrix<complex<double> > rho(2, 2);
270  if (_vesSheet) {
271  if (qKmKm.imag() > 0)
272  qKmKm *= -1;
273  rho(0, 0) = (2. * qPiPi) / mass;
274  rho(1, 1) = (2. * qKmKm) / mass;
275  } else {
276  rho(0, 0) = ((2. * qPiPi) / mass + (2. * qPi0Pi0) / mass) / 2.;
277  rho(1, 1) = ((2. * qKK) / mass + (2. * qK0K0) / mass) / 2.;
278  }
279  rho(0, 1) = rho(1, 0) = 0;
280 
281  const double scale = (s / (4 * _kaonMeanMass * _kaonMeanMass)) - 1;
282 
283  matrix<complex<double> > M(zero_matrix<complex<double> >(2, 2));
284  for (unsigned int i = 0; i < _sP.size2(); ++i) {
285  const complex<double> fa = 1. / (s - _sP(0, i));
286  M += fa * _a[i];
287  }
288  for (unsigned int i = 0; i < _c.size(); ++i) {
289  const complex<double> sc = pow(scale, (int)i);
290  M += sc *_c[i];
291  }
292 
293  // modification: off-diagonal terms set to 0
294  M(0, 1) = 0;
295  M(1, 0) = 0;
296 
297  invertMatrix<complex<double> >(M - imag * rho, _T);
298  const complex<double> amp = _T(0, 0);
299  if (_debug)
300  printDebug << name() << "(m = " << maxPrecision(mass) << " GeV) = "
301  << maxPrecisionDouble(amp) << endl;
302 
303  return amp;
304 }
305 
306 
308 piPiSWaveAuMorganPenningtonVes::piPiSWaveAuMorganPenningtonVes()
309  : piPiSWaveAuMorganPenningtonM()
310 {
311  _vesSheet = 1;
312 }
313 
314 
315 complex<double>
316 piPiSWaveAuMorganPenningtonVes::amp(const isobarDecayVertex& v)
317 {
318  const double M = v.parent()->lzVec().M();
319 
320  const double f0Mass = 0.9837; // [GeV]
321  const double f0Width = 0.0376; // [GeV]
322  const complex<double> coupling(-0.3743, 0.3197);
323 
324  const complex<double> ampM = piPiSWaveAuMorganPenningtonM::amp(v);
325 
326  complex<double> bw = 0;
327  if (M > 2 * _piChargedMass) {
328  const double q = breakupMomentum(M, _piChargedMass, _piChargedMass);
329  const double q0 = breakupMomentum(f0Mass, _piChargedMass, _piChargedMass);
330  const double Gamma = f0Width * (q / q0);
331  // A / (B - iC) = (A / (B^2 + C^2)) * (B + iC)
332  const double C = f0Mass * Gamma;
333  const double A = C * (M / q);
334  const double B = f0Mass * f0Mass - M * M;
335  bw = (A / (B * B + C * C)) * complex<double>(B, C);
336  }
337 
338  const complex<double> amp = ampM - coupling * bw;
339  if (_debug)
340  printDebug << name() << "(m = " << maxPrecision(M) << " GeV) = "
341  << maxPrecisionDouble(amp) << endl;
342 
343  return amp;
344 }
345 
346 
348 piPiSWaveAuMorganPenningtonKachaev::piPiSWaveAuMorganPenningtonKachaev()
349  : piPiSWaveAuMorganPenningtonM()
350 {
351  // change parameters according to Kachaev's prescription
352  _c[4](0, 0) = 0; // was 0.1957;
353  _c[4](1, 1) = 0; // was -0.3977;
354 
355  _a[0](0, 1) = 0; // was 0.0150
356  _a[0](1, 0) = 0; // was 0.0150
357 
358  // _a[1] are the f's from the AMP paper
359  _a[1](0, 0) = 0;
360  _a[1](0, 1) = 0;
361  _a[1](1, 0) = 0;
362  _a[1](1, 1) = 0;
363 }
364 
365 
367 complex<double>
368 rhoPrimeMassDep::amp(const isobarDecayVertex& v)
369 {
370  const particlePtr& parent = v.parent();
371 
372  // get Breit-Wigner parameters
373  const double M = parent->lzVec().M(); // parent mass
374  const double m1 = v.daughter1()->lzVec().M(); // daughter 1 measured/assumed mass
375  const double m2 = v.daughter2()->lzVec().M(); // daughter 2 measured/assumed mass
376  const double q = breakupMomentum(M, m1, m2);
377  // const double M0 = parent->mass(); // resonance peak position
378  // const double q02 = breakupMomentumSquared(M0, m1, m2, true);
379  // const double q0 = sqrt(fabs(q02)); // !NOTE! this is incorrect but this is how it was done in PWA2000
380  const unsigned int L = v.L();
381 
382  // rho' parameters
383  const double M01 = 1.465; // rho(1450) mass [GeV/c^]
384  const double Gamma01 = 0.235; // rho(1450) width [GeV/c^]
385  const double M02 = 1.700; // rho(1700) mass [GeV/c^]
386  const double Gamma02 = 0.220; // rho(1700) width [GeV/c^]
387  // const double M02 = 1.720; // rho(1700) mass; PDG12 [GeV/c^]
388  // const double Gamma02 = 0.250; // rho(1700) width; PDG12 [GeV/c^]
389  const double q102 = breakupMomentumSquared(M01, m1, m2, true);
390  const double q10 = sqrt(fabs(q102));
391  const double q202 = breakupMomentumSquared(M02, m1, m2, true);
392  const double q20 = sqrt(fabs(q202));
393 
394  // const complex<double> bw1 = breitWigner(M, M01, Gamma01, L, q, q0);
395  // const complex<double> bw2 = breitWigner(M, M02, Gamma02, L, q, q0);
396  const complex<double> bw1 = breitWigner(M, M01, Gamma01, L, q, q10);
397  const complex<double> bw2 = breitWigner(M, M02, Gamma02, L, q, q20);
398  const complex<double> amp = (4 * bw1 - 3 * bw2) / 7;
399 
400  if (_debug)
401  // printDebug << name() << "(m = " << maxPrecision(M) << " GeV/c^2, m_0 = " << maxPrecision(M0)
402  // << " GeV/c^2, L = " << spinQn(L) << ", q = " << maxPrecision(q) << " GeV/c, "
403  // << "q_0 = " << maxPrecision(q0) << " GeV/c) = " << maxPrecisionDouble(amp) << endl;
404  printDebug << name() << "(m = " << maxPrecision(M) << " GeV/c^2, "
405  << "L = " << spinQn(L) << ", q = " << maxPrecision(q) << " GeV/c, "
406  << "q1_0 = " << maxPrecision(q10) << " GeV/c, "
407  << "q2_0 = " << maxPrecision(q20) << " GeV/c) = " << maxPrecisionDouble(amp) << endl;
408  return amp;
409 }
Generated by   doxygen 1.8.1.2