ROOTPWA
massDep.cc
Go to the documentation of this file.
1 #include "massDep.h"
2 
3 
4 #define MAXPRECISION(val) setprecision(numeric_limits<double>::digits10 + 1) << val
5 
6 
7 using namespace std;
8 
9 
10 extern particleDataTable PDGtable;
11 
12 
13 complex<double> flat::val(const particle& p)
14 { return complex<double>(1,0); }
15 
16 
17 complex<double> flatRange::val(const particle& p)
18 {
19 
20  const double m = ~(p.get4P ());
21  //cerr << m << " " << mlow << " " << mhigh << endl;
22  if(m>mlow && m<mhigh) return complex<double>(1,0);
23  else return complex<double>(0,0);
24 }
25 
26 
27 complex<double> breitWigner::val(const particle& p)
28 {
29  const double m0 = p.Mass ();
30  const double Gamma0 = p.Width ();
31  const double q = p.q ();
32  const double q0 = p.q0 ();
33  const double m = ~(p.get4P ());
34  const int l = p.Decay ()->L ();
35 
36  const double GammaV = Gamma0 * (m0 / m) * (q / q0) * (pow(F(l, q), 2) / pow(F(l, q0), 2));
37  complex<double> ret = (m0 * Gamma0) / (m0 * m0 - m * m - complex<double>(0, 1) * m0 * GammaV);
38  return ret;
39 }
40 
41 
42 // rho(1450)/rho(1700)
43 // see DOI: 10.1007/BF01552547
44 complex<double> rhoPrime::val(const particle& p)
45 {
46 
47  const double m1 = 1.465;
48  const double Gamma1 = 0.235;
49  const double m2 = 1.720;
50  const double Gamma2 = 0.220;
51 
52  const double q = p.q ();
53  const double q0 = p.q0 ();
54  const double m = ~(p.get4P ());
55  const int l = p.Decay ()->L ();
56 
57  const double GammaV1 = Gamma1 * (m1 / m) * (q / q0) * (pow(F(l, q), 2) / pow(F(l, q0), 2));
58  const double GammaV2 = Gamma2 * (m2 / m) * (q / q0) * (pow(F(l, q), 2) / pow(F(l, q0), 2));
59 
60  //cout << m << " " << GammaV1 << " " << GammaV2 << endl;
61 
62  complex<double> amp1 = (m1 * Gamma1) / (m1 * m1 - m * m - complex<double>(0, 1) * m1 * GammaV1);
63  complex<double> amp2 = (m2 * Gamma2) / (m2 * m2 - m * m - complex<double>(0, 1) * m2 * GammaV2);
64 
65 
66  complex<double> ret= 1./7 * ( 4. * amp1 - 3. * amp2 );
67  return ret;
68 
69 }
70 
71 
72 AMP_M::AMP_M() {
73  ves_sheet = 0;
74  _Pmax = 2;
75  _Nmax = 5;
76 
77  _rho = matrix<complex<double> >(2, 2);
78  _M = matrix<complex<double> >(2, 2);
79  _T = matrix<complex<double> >(2, 2);
80 
81  _f = matrix<complex<double> >(1, 2);
82  _f.el(0, 0) = 0.1968;
83  _f.el(0, 1) = -0.0154;
84 
85  _a = vector<matrix<complex<double> > >(2, matrix<complex<double> >(2, 2));
86  _a[0].el(0, 0) = 0.1131;
87  _a[0].el(0, 1) = 0.0150;
88  _a[0].el(1, 0) = 0.0150;
89  _a[0].el(1, 1) = -0.3216;
90  _a[1].el(0, 0) = _f.el(0, 0) * _f.el(0, 0);
91  _a[1].el(0, 1) = _f.el(0, 0) * _f.el(0, 1);
92  _a[1].el(1, 0) = _f.el(0, 1) * _f.el(0, 0);
93  _a[1].el(1, 1) = _f.el(0, 1) * _f.el(0, 1);
94 
95 
96  _c = vector<matrix<complex<double> > >(5, matrix<complex<double> >(2, 2));
97  _c[0].el(0, 0) = 0.0337;
98  _c[1].el(0, 0) = -0.3185;
99  _c[2].el(0, 0) = -0.0942;
100  _c[3].el(0, 0) = -0.5927;
101  _c[4].el(0, 0) = 0.1957;
102  _c[0].el(0, 1) = _c[0].el(1, 0) = -0.2826;
103  _c[1].el(0, 1) = _c[1].el(1, 0) = 0.0918;
104  _c[2].el(0, 1) = _c[2].el(1, 0) = 0.1669;
105  _c[3].el(0, 1) = _c[3].el(1, 0) = -0.2082;
106  _c[4].el(0, 1) = _c[4].el(1, 0) = -0.1386;
107  _c[0].el(1, 1) = 0.3010;
108  _c[1].el(1, 1) = -0.5140;
109  _c[2].el(1, 1) = 0.1176;
110  _c[3].el(1, 1) = 0.5204;
111  _c[4].el(1, 1) = -0.3977;
112 
113  _sP = matrix<double>(1, 2);
114  _sP.el(0, 0) = -0.0074;
115  _sP.el(0, 1) = 0.9828;
116 }
117 
118 
119 complex<double> AMP_M::val(const particle& p)
120 {
121  extern particleDataTable PDGtable;
122 
123  _M.el(0, 0) = 0;
124  _M.el(0, 1) = 0;
125  _M.el(1, 0) = 0;
126  _M.el(0, 1) = 0;
127 
128  complex <double> ret;
129  complex <double> i(0, 1);
130 
131  double pi_mass = PDGtable.get("pi").Mass();
132  double pi0_mass = PDGtable.get("pi0").Mass();
133  double K_mass = PDGtable.get("K").Mass();
134  double K0_mass = PDGtable.get("K0").Mass();
135  double Kmean_mass = 0.5 * (K_mass + K0_mass);
136  double My_mass = ~(p.get4P());
137  double s = My_mass * My_mass;
138  if (fabs(s - _sP.el(0, 1)) < 1e-6) {
139  My_mass += 1e-6;
140  s = My_mass * My_mass;
141  }
142 
143  complex<double> q_pipi = q(My_mass, pi_mass, pi_mass);
144  complex<double> q_pi0pi0 = q(My_mass, pi0_mass, pi0_mass);
145  complex<double> q_KK = q(My_mass, K_mass, K_mass);
146  complex<double> q_K0K0 = q(My_mass, K0_mass, K0_mass);
147  complex<double> q_KmKm = q(My_mass, Kmean_mass, Kmean_mass);
148  _rho.el(0, 0) = 0.5 * ((2.0 * q_pipi) / My_mass + (2.0 * q_pi0pi0) / My_mass);
149  _rho.el(1, 1) = 0.5 * ((2.0 * q_KK) / My_mass + (2.0 * q_K0K0) / My_mass);
150 
151  if (ves_sheet) {
152  if (q_KmKm.imag() > 0.0)
153  q_KmKm *= -1;
154  _rho.el(0, 0) = (2.0 * q_pipi) / My_mass;
155  _rho.el(1, 1) = (2.0 * q_KmKm) / My_mass;
156  }
157  // _rho.print();
158 
159  double scale = (s / (4.0 * Kmean_mass * Kmean_mass)) - 1;
160  //cout << "scale=" << scale << endl;
161 
162  for (int p = 0; p < _Pmax; ++p) {
163  complex<double> fa = 1. / (s - _sP.el(0, p));
164  // cout << "fa[" << p <<"] = "<< fa << "; a[" << p << "] = " << _a[p] << endl;
165  _M += fa * _a[p];
166  }
167  for (int n = 0; n < _Nmax; ++n) {
168  complex<double> sc = pow(scale, n);
169  // cout << "sc[" << n <<"] = "<< sc << "; c[" << n << "] = " << _c[n] << endl;
170  _M += sc *_c[n];
171  }
172 
173  // Modification: Here we set off-diagonal terms to 0
174  _M.el(0, 1) = 0;
175  _M.el(1, 0) = 0;
176  // _M.print();
177 
178  _T = (_M - i * _rho).inv();
179  // _T.print();
180 
181  ret = _T.el(0, 0);
182  return ret;
183 }
184 
185 
186 complex<double> AMP_ves::val(const particle& p) {
187  complex<double> bw, amp_m;
188  static complex<double> coupling(-0.3743, 0.3197);
189  static double massf0 = 0.9837, widthf0 = 0.0376;
190 
191  double pi_mass = PDGtable.get("pi").Mass();
192  double My_mass = ~(p.get4P());
193 
194  ves_sheet = 1;
195  amp_m = AMP_M::val(p);
196 
197  if (My_mass > 2.0 * pi_mass) {
198  double p, p0, gam, denom, A, B, C;
199  p = q(My_mass, pi_mass, pi_mass).real();
200  p0 = q(massf0, pi_mass, pi_mass).real();
201  gam = widthf0 * (p / p0);
202  A = massf0 * massf0 - My_mass * My_mass;
203  B = massf0 * gam;
204  C = B * (My_mass / p);
205  denom = C / (A * A + B * B);
206  bw = denom * complex<double>(A, B);
207  }
208 
209  return amp_m - coupling * bw;
210 }
211 
212 
213 AMP_kach::AMP_kach(): AMP_M()
214 {
215  // Change parameters according to Kachaev's prescription
216  _c[4].el(0, 0) = 0; // was 0.1957;
217  _c[4].el(1, 1) = 0; // was -0.3977;
218 
219  _a[0].el(0, 1) = 0; // setting off-diagonal values to 0
220  _a[0].el(1, 0) = 0; //
221 
222  // _a[1] are the f's from the AMP paper! Setting to 0!
223  _a[1].el(0, 0) = 0;
224  _a[1].el(0, 1) = 0;
225  _a[1].el(1, 0) = 0;
226  _a[1].el(1, 1) = 0;
227 }
228 
229 
230 AMP_kach::AMP_kach(const AMP_kach&): AMP_M()
231 {
232  // Change parameters according to Kachaev's prescription
233  _c[4].el(0, 0) = 0; // was 0.1957;
234  _c[4].el(1, 1) = 0; // was -0.3977;
235 
236  _a[0].el(0, 1) = 0; // setting off-diagonal values to 0
237  _a[0].el(1, 0) = 0; //
238 
239  // _a[1] are the f's from the AMP paper! Setting to 0!
240  _a[1].el(0, 0) = 0;
241  _a[1].el(0, 1) = 0;
242  _a[1].el(1, 0) = 0;
243  _a[1].el(1, 1) = 0;
244 }
245 
246 
247 /*
248  SUBROUTINE KPIAMP(BW, W , W0 , G0 , P , P0) ! K+ pi- solution
249  c-------
250  c K-PI S-WAVE PARAMETRISATION
251  c SOURCE: NP B296 493
252  c Here used formula exp(-i*a)*sin(a) + BW*exp(-2.*i*a)
253  c with BW parameters M=1.437 and W=0.355 .
254  c That gives amplitude and phase which coincide rouphly
255  c with pictures from article
256  c-------
257  REAL W0, G0 ! MASS AND WIDTH OF RES.
258  REAL P, P0 ! CURRENT MOM, RESONANCE MOM
259  REAL W, S ! DIMESON MASS, MASS SQUARED (GEV**2)
260  COMPLEX BW ! AMPLITUDE
261  COMPLEX BG ! BACKGROUND
262  PARAMETER ( WPI = 0.1395675) ! PI+- MASS
263  PARAMETER ( WKC = 0.493646 ) ! K+- MASS
264  PARAMETER ( WETAP= 0.95775 ) ! ETA_PRIME MASS
265  PARAMETER ( WLOW = WPI+WKC ) ! LOW LIMIT FOR PARAMETRISATION
266  PARAMETER ( WUP = 1.7 ) ! UPPER LIMIT FOR PARAMETRISATION
267  C---- Model parameters
268  PARAMETER ( EPS = 1.00 ) ! ratio
269  PARAMETER ( A = 1.93 ) !
270  PARAMETER ( B = 2.66 ) !
271  PARAMETER ( PHI = -16.9 ) !
272  C----
273  BW = (0.,0.)
274  S = W*W
275  IF ( S.LE.WLOW**2 ) RETURN
276  IF ( S.GE. WUP**2 ) RETURN
277  C---- BW-PARAMETRISATION
278  GK = G0*W0/W*P/P0 ! WIDTH FOR K-PI CHANNEL
279  BWR = W0**2 - W**2 ! REAL PART
280  BWI = W0*GK ! IMAGE PART
281  C-- C/(A-i*B) = (C*A/(A**2+B**2), C*B/(A**2+B**2))
282  DENOM = EPS*W0*GK / (BWR**2+BWI**2)
283  BW = CMPLX(BWR*DENOM, BWI*DENOM)
284  C---- Background parametrisation
285  C---- sin(a)*exp(i*a)=(ctg(a)+i)/(ctg(a)**2+1.)
286  CTG = 1.0/(A*P) + B*P/2.0
287  C-- COS(ARCCTG(X) = X/SQRT(1+X**2)
288  C-- SIN(ARCCTG(X) = 1/SQRT(1+X**2)
289  XX = SQRT(1.+CTG**2)
290  CS = CTG/XX
291  SN = 1. /XX
292  BG = CMPLX( SN*CS , SN**2 )
293  C---- AMP = BG + BW*EXP(I*2.*A)
294  BW = BG + BW*CMPLX(CS**2-SN**2,2.*SN*CS)
295  RETURN
296  END
297 */
298 complex<double> AMP_LASS::val(const particle& p)
299 {
300  /*
301  const double m0 = p.Mass ();
302  const double Gamma0 = p.Width ();
303  const double q = p.q ();
304  const double q0 = p.q0 ();
305  const double m = ~(p.get4P ());
306  const int l = p.Decay ()->L ();
307 
308  const double GammaV = Gamma0 * (m0 / m) * (q / q0) * (pow(F(l, q), 2) / pow(F(l, q0), 2));
309  complex<double> ret = (m0 * Gamma0) / (m0 * m0 - m * m - complex<double>(0, 1) * m0 * GammaV);
310  return ret;*/
311 
312  //BW, W , W0 , G0 , P , P0) ! K+ pi- solution
313  double W0 = 1.437;//p.Mass();
314  double G0 = 0.355;//p.Width(); // MASS AND WIDTH OF RES.
315  double P = p.q();
316  double P0 = p.q0(); // CURRENT MOM, RESONANCE MOM
317  double W = ~(p.get4P ());
318  double S ; // DIMESON MASS, MASS SQUARED (GEV**2)
319  complex<double> BW; // AMPLITUDE
320  complex<double> BG; // BACKGROUND
321  double WPI = 0.1395675; // PI+- MASS
322  double WKC = 0.493646 ; // K+- MASS
323  //double WETAP= 0.95775 ; // ETA_PRIME MASS
324  double WLOW = WPI+WKC ; // LOW LIMIT FOR PARAMETRISATION
325  double WUP = 1.7 ; // UPPER LIMIT FOR PARAMETRISATION
326  //---- Model parameters
327  double EPS = 1.00 ; // ratio
328  double A = 1.93 ;
329  double B = 2.66 ;
330  //double PHI = -16.9 ;
331  //----
332  BW = complex<double>(0.,0.);
333  S = W*W;
334  if ( S <= WLOW*WLOW ) return BW;
335  if ( S >= WUP*WUP ) return BW;
336  //---- BW-PARAMETRISATION
337  double GK = G0*W0/W*P/P0; // WIDTH FOR K-PI CHANNEL
338  BW = complex<double>(W0*W0 - W*W, W0*GK);
339  // C/(A-i*B) = (C*A/(A**2+B**2), C*B/(A**2+B**2))
340  double DENOM = EPS*W0*GK / (real(BW)*real(BW)+imag(BW)*imag(BW));
341  BW = complex<double>(real(BW)*DENOM, imag(BW)*DENOM);
342  //---- Background parametrisation
343  //---- sin(a)*exp(i*a)=(ctg(a)+i)/(ctg(a)**2+1.)
344  double CTG = 1.0/(A*P) + B*P/2.0;
345  //-- COS(ARCCTG(X) = X/SQRT(1+X**2)
346  //-- SIN(ARCCTG(X) = 1/SQRT(1+X**2)
347  double XX = sqrt(1.+CTG*CTG);
348  double CS = CTG/XX;
349  double SN = 1. /XX;
350  BG = complex<double>( SN*CS , SN*SN );
351  //---- AMP = BG + BW*EXP(I*2.*A)
352  BW = BG + BW*complex<double>(CS*CS-SN*SN,2.*SN*CS);
353  return BW;
354 }
Generated by   doxygen 1.8.1.2