BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
Dalitz.cxx
Go to the documentation of this file.
1// -*- C++ -*-
2
3//
4// Author: Dan Ambrose
5// Created: Mon Mar 5 2012
6// based on Cleo package by Eric White and Warner Sun
7
8//
9// Revision history
10//
11//
12
13#include "Dalitz.h"
14// #include "EvtGenBase/EvtGenKine.hh"
15#include <complex>
16
17using std::string;
18
19double Dalitz::PI = 3.1415926; // pi
20
21Dalitz::Dalitz() {
22 N = 8; // Default bins if not specified
23}
24
25// Constructor
26Dalitz::Dalitz( int binNum ) {
27 N = binNum; // Set bin number
28} // end constructor
29
30TComplex Dalitz::Amplitude( double x, double y, double z ) {
31
32 // PRD 73, 112009(2006) Belle
33 // for D0 particle 1: K particle 2: pi- particle 3: pi+
34 // the right order is: (1,2), (1,3), (3,2)
35 double m_mass[4] = { 1.86450, 0.497648, 0.139570, 0.139570 }; // mass
36
37 TComplex D0( 0.0, 0.0 );
38
39 // x, y, z already squared, need to get back the mass!!!!
40 x = sqrt( x );
41 y = sqrt( y );
42 z = sqrt( z );
43
44 // Array for resonances
45 TComplex DK2piRes[19];
46
47 // x->3 y->2 z->1
48 DK2piRes[0] = sakurai( x, y, z, 1.00, 0.0, 0.1503, 0.7758 ); // RHO(770)
49 DK2piRes[1] = resAmp( x, y, z, m_mass[3], m_mass[2], m_mass[1], 0.0314, 110.8, 0.00849,
50 0.78259, 1 ); // OMEGA
51 DK2piRes[2] = f_980( z, 0.980, 0.365, 201.9 ); // F_0(980)
52 DK2piRes[3] = resAmp( x, y, z, m_mass[3], m_mass[2], m_mass[1], 1.32, 348, 0.1851, 1.2754,
53 2 ); // F_2(1270)
54 DK2piRes[4] = resAmp( x, y, z, m_mass[3], m_mass[2], m_mass[1], 1.44, 82, 0.173, 1.434,
55 0 ); // F_0(1370) hep-ph 0009168
56 DK2piRes[5] = sakurai( x, y, z, 0.66, 9, 0.400, 1.465 ); // RHO(1450)
57 DK2piRes[6] = resAmp( x, y, z, m_mass[3], m_mass[2], m_mass[1], 1.43, 212, 0.454, 0.519,
58 0 ); // Sigma(600)
59 DK2piRes[7] =
60 resAmp( x, y, z, m_mass[3], m_mass[2], m_mass[1], 0.23, 237, 0.101, 1.050, 0 ); // Sigma
61 DK2piRes[8] = resAmp( x, z, y, m_mass[3], m_mass[1], m_mass[2], 1.644, 132.1, 0.0508,
62 0.89166, 1 ); // K*(892)-
63 DK2piRes[9] = resAmp( x, z, y, m_mass[3], m_mass[1], m_mass[2], 0.61, 113, 0.232, 1.414,
64 1 ); // K*(1410)-
65 DK2piRes[10] = resAmp( x, z, y, m_mass[3], m_mass[1], m_mass[2], 2.15, 353.6, 0.294, 1.412,
66 0 ); // K*_0(1430)-
67 DK2piRes[11] = resAmp( x, z, y, m_mass[3], m_mass[1], m_mass[2], 0.88, 318.7, 0.0985, 1.4256,
68 2 ); // K*_2(1430)-
69 DK2piRes[12] = resAmp( x, z, y, m_mass[3], m_mass[1], m_mass[2], 1.39, 103, 0.322, 1.717,
70 1 ); // K*(1680)-
71 DK2piRes[13] = resAmp( y, z, x, m_mass[2], m_mass[1], m_mass[3], 0.144, 320.3, 0.0508,
72 0.89166, 1 ); // K*(892)+
73 DK2piRes[14] = resAmp( y, z, x, m_mass[2], m_mass[1], m_mass[3], 0.45, 254, 0.232, 1.414,
74 1 ); // K*(1410)+
75 DK2piRes[15] = resAmp( y, z, x, m_mass[2], m_mass[1], m_mass[3], 0.47, 88, 0.294, 1.412,
76 0 ); // K*_0(1430)+
77 DK2piRes[16] = resAmp( y, z, x, m_mass[2], m_mass[1], m_mass[3], 0.25, 265, 0.0985, 1.4256,
78 2 ); // K*_2(1430)+
79 DK2piRes[17] = resAmp( y, z, x, m_mass[2], m_mass[1], m_mass[3], 1.2, 118, 0.322, 1.717,
80 1 ); // K*(1680)+
81
82 double pi180inv = 3.1415926 / 180.;
83 DK2piRes[18] = TComplex( 3.0 * cos( 164 * pi180inv ), 3.0 * sin( 164 * pi180inv ) );
84
85 for ( int i = 0; i < 19; i++ ) { D0 += DK2piRes[i]; }
86
87 return D0;
88} // Amplitude
89
90TComplex Dalitz::CLEO_Amplitude( double x, double y, double z ) {
91
92 double m_mass[4] = { 1.86450, 0.497648, 0.139570, 0.139570 }; // mass
93
94 TComplex D0( 0.0, 0.0 );
95
96 // x, y, z already squared, need to get back the mass!!!!
97 x = sqrt( x );
98 y = sqrt( y );
99 z = sqrt( z );
100
101 // Array for resonances
102 TComplex DK2piRes[3];
103
104 // x->3 y->2 z->1
105 DK2piRes[0] = CLEO_resAmp( x, z, y, m_mass[3], m_mass[1], m_mass[2], 2.31, 109.0, 0.0498,
106 0.89610, 1 ); // K*(892)
107 DK2piRes[1] = CLEO_resAmp( x, y, z, m_mass[3], m_mass[2], m_mass[1], 1.59, -123.0, 0.1491,
108 0.7683, 1 ); // RHO(770)
109 DK2piRes[2] = TComplex( 1.0, 0.0 );
110
111 for ( int i = 0; i < 3; i++ ) { D0 += DK2piRes[i]; }
112
113 return D0;
114}
115
116double Dalitz::Phase( double x, double y, double z, int Babar ) {
117
118 TComplex D0( 0, 0 );
119 TComplex D0bar( 0, 0 );
120
121 if ( Babar == 1 )
122 {
123 D0 = Babar_Amplitude( x, y, z );
124 D0bar = Babar_Amplitude( y, x, z );
125 }
126 else
127 {
128 D0 = Amplitude( x, y, z );
129 D0bar = Amplitude( y, x, z );
130 }
131
132 // Compute phase difference
133 double deltaD = arg( D0 ) - arg( D0bar );
134 if ( x < y ) deltaD = -deltaD; // make it symmetric to the lower half
135
136 if ( deltaD < -PI / N ) deltaD += 2 * PI; // shift deltaD to [-PI/8,15*PI/8]
137 if ( deltaD > ( 2 * N - 1 ) * PI / N ) deltaD -= 2 * PI; // shift deltaD to [-PI/8,15*PI/8]
138
139 return deltaD;
140
141} // end Phase
142
143bool Dalitz::Point_on_DP( double x, double y ) {
144
145 double m_mass[4] = { 1.86450, 0.497648, 0.139570, 0.139570 }; // mass
146 double m_mass2[4] = { 1.86450 * 1.86450, 0.497648 * 0.497648, 0.139570 * 0.139570,
147 0.139570 * 0.139570 }; // mass square
148
149 double m_XmaxDP = m_mass[0] - m_mass[3];
150 m_XmaxDP *= m_XmaxDP;
151 double m_XminDP = m_mass[1] + m_mass[2];
152 m_XminDP *= m_XminDP;
153
154 if ( ( x > m_XmaxDP ) || ( x < m_XminDP ) ) return false;
155
156 double Low = 0;
157 double Up = 0;
158 double HInv_m12 = 0.5 / sqrt( x );
159 double E1 = HInv_m12 * ( x + m_mass2[1] - m_mass2[2] );
160 double E3 = HInv_m12 * ( m_mass2[0] - m_mass2[3] - x );
161 double E1_2 = E1 * E1;
162 double E3_2 = E3 * E3;
163
164 if ( E1 < m_mass[1] )
165 {
166 E1 = m_mass[1];
167 E1_2 = m_mass2[1];
168 }
169 if ( E3 < m_mass[3] )
170 {
171 E3 = m_mass[3];
172 E3_2 = m_mass2[3];
173 }
174
175 double temp = E1_2 - m_mass2[1];
176 if ( temp < 0 ) temp = 0;
177 double P1 = sqrt( temp );
178 temp = E3_2 - m_mass2[3];
179 if ( temp < 0 ) temp = 0;
180 double P3 = sqrt( temp );
181 double E13_2 = ( E1 + E3 ) * ( E1 + E3 );
182
183 // Compute hi and lo y-coord
184 Low = E13_2 - ( P1 + P3 ) * ( P1 + P3 );
185 Up = E13_2 - ( P1 - P3 ) * ( P1 - P3 );
186
187 if ( ( y > Up ) || ( y < Low ) ) return false;
188
189 return true;
190} // end Point_on_DP
191
192// Alternate, for points outside DP
193bool Dalitz::Point_on_DP2( double x, double y ) {
194
195 double m_mass[4] = { 1.86450, 0.497648, 0.139570, 0.139570 }; // mass
196 double m_mass2[4] = { 1.86450 * 1.86450, 0.497648 * 0.497648, 0.139570 * 0.139570,
197 0.139570 * 0.139570 }; // mass square
198
199 double m_XmaxDP = m_mass[0] - m_mass[3];
200 m_XmaxDP *= m_XmaxDP;
201 double m_XminDP = m_mass[1] + m_mass[2];
202 m_XminDP *= m_XminDP;
203
204 if ( ( x > m_XmaxDP ) || ( x < m_XminDP ) ) return false;
205
206 double Low = 0;
207 double Up = 0;
208 double HInv_m12 = 0.5 / sqrt( x );
209 double E1 = HInv_m12 * ( x + m_mass2[1] - m_mass2[2] );
210 double E3 = HInv_m12 * ( m_mass2[0] - m_mass2[3] - x );
211 double E1_2 = E1 * E1;
212 double E3_2 = E3 * E3;
213
214 if ( E1 < m_mass[1] )
215 {
216 E1 = m_mass[1];
217 E1_2 = m_mass2[1];
218 }
219 if ( E3 < m_mass[3] )
220 {
221 E3 = m_mass[3];
222 E3_2 = m_mass2[3];
223 }
224
225 double temp = E1_2 - m_mass2[1];
226 if ( temp < 0 ) temp = 0;
227 double P1 = sqrt( temp );
228 temp = E3_2 - m_mass2[3];
229 if ( temp < 0 ) temp = 0;
230 double P3 = sqrt( temp );
231 double E13_2 = ( E1 + E3 ) * ( E1 + E3 );
232
233 Low = E13_2 - ( P1 + P3 ) * ( P1 + P3 );
234 Up = E13_2 - ( P1 - P3 ) * ( P1 - P3 );
235
236 if ( ( y > ( Up + 0.05 ) ) || ( y < ( Low - 0.05 ) ) ) return false; // make it larger
237
238 return true;
239} // end Point_on_DP2
240
241TComplex Dalitz::CLEO_resAmp( double mAC, double mBC, double mAB, double mB, double mA,
242 double mC, double _ampl, double _theta, double _gamma,
243 double _bwm, int _spin ) {
244
245 double pi180inv = 3.1415926 / 180.;
246
247 TComplex ampl;
248
249 // EvtVector4R _p4_d3 = _p4_p-_p4_d1-_p4_d2;
250
251 // get cos of the angle between the daughters from their 4-momenta
252 // and the 4-momentum of the parent
253
254 // in general, EvtDecayAngle(parent, part1+part2, part1) gives the angle
255 // the missing particle (not listed in the arguments) makes
256 // with part2 in the rest frame of both
257 // listed particles (12)
258
259 // angle 3 makes with 2 in rest frame of 12 (CS3)
260 // double cos_phi_0 = EvtDecayAngle(_p4_p, _p4_d1+_p4_d2, _p4_d1);
261 // double EvtDecayAngle(const EvtVector4R& p,const EvtVector4R& q,
262 // const EvtVector4R& d) {
263
264 // double pd=p*d;
265 // double pq=p*q;
266 // double qd=q*d;
267 // double mp2=p.mass2();
268 // double mq2=q.mass2();
269 // double md2=d.mass2();
270
271 // double cost=(pd*mq2-pq*qd)/sqrt((pq*pq-mq2*mp2)*(qd*qd-mq2*md2));
272
273 // return cost;
274
275 // }
276
277 // double mD = 1.86450 ;
278 double mD = 1.86484;
279 double eA = ( mD * mD - mBC * mBC + mA * mA ) / ( 2. * mD );
280 double eAB = ( mD * mD - mC * mC + mAB * mAB ) / ( 2. * mD );
281
282 // Take D to be at rest
283 double pd = mD * eA;
284 double pq = mD * eAB;
285 double qd = mA * mA + 0.5 * ( mAB * mAB - mA * mA - mB * mB );
286 double mp2 = mD * mD;
287 double mq2 = mAB * mAB;
288 double md2 = mA * mA;
289 double cos_phi_0 =
290 ( pd * mq2 - pq * qd ) / sqrt( ( pq * pq - mq2 * mp2 ) * ( qd * qd - mq2 * md2 ) );
291
292 // angle 3 makes with 1 in 12 is, of course, -cos_phi_0
293
294 // double mAB=(_p4_d1+_p4_d2).mass();
295 // double mBC=(_p4_d2+p4_d3).mass();
296 // double mAC=(_p4_d1+p4_d3).mass();
297 // double mA=_p4_d1.mass();
298 // double mB=_p4_d2.mass();
299 // double mD=_p4_p.mass();
300 // double mC=p4_d3.mass();
301
302 switch ( _spin )
303 {
304
305 case 0:
306 ampl = ( _ampl * TComplex( cos( _theta * pi180inv ), sin( _theta * pi180inv ) ) *
307 sqrt( _gamma / ( 2. * PI ) ) *
308 ( 1.0 / ( mAB - _bwm - TComplex( 0.0, 0.5 * _gamma ) ) ) );
309 break;
310
311 case 1:
312 ampl = ( _ampl * TComplex( cos( _theta * pi180inv ), sin( _theta * pi180inv ) ) *
313 sqrt( _gamma / ( 2. * PI ) ) *
314 ( cos_phi_0 / ( mAB - _bwm - TComplex( 0.0, 0.5 * _gamma ) ) ) );
315 break;
316
317 // case 2:
318 // ampl=(_ampl*TComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
319 // sqrt(_gamma/(2.*PI))*
320 // ((1.5*cos_phi_0*cos_phi_0-0.5)/((_p4_d1+_p4_d2).mass()-_bwm-EvtComplex(0.0,
321 // 0.5*_gamma)))); break;
322
323 // case 3:
324 // ampl=(_ampl*EvtComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
325 // sqrt(_gamma/EvtConst::twoPi)*
326 // ((2.5*cos_phi_0*cos_phi_0*cos_phi_0-1.5*cos_phi_0)/((_p4_d1+_p4_d2).mass()-_bwm-EvtComplex(0.0,
327 // 0.5*_gamma)))); break;
328
329 default:
330 // cout << "EvtGen: wrong spin in CLEO_resAmp()" << endl;
331 ampl = TComplex( 0.0 );
332 break;
333 }
334
335 return ampl;
336}
337
338TComplex Dalitz::resAmp( double mAC, double mBC, double mAB, double mB, double mA, double mC,
339 double _ampl, double _theta, double _gamma, double _bwm, int _spin ) {
340
341 double pi180inv = 3.1415926 / 180.;
342
343 TComplex ampl;
344
345 double mD = 1.86450;
346
347 double mR = _bwm;
348 double gammaR = _gamma;
349
350 double temp = ( mAB * mAB - mA * mA - mB * mB ) * ( mAB * mAB - mA * mA - mB * mB ) / 4.0 -
351 mA * mA * mB * mB;
352 if ( temp < 0 ) temp = 0;
353 double pAB = sqrt( temp / ( mAB * mAB ) );
354
355 temp = ( mR * mR - mA * mA - mB * mB ) * ( mR * mR - mA * mA - mB * mB ) / 4.0 -
356 mA * mA * mB * mB;
357 if ( temp < 0 ) temp = 0;
358 double pR = sqrt( temp / ( mR * mR ) );
359
360 temp = ( mD * mD - mR * mR - mC * mC ) * ( mD * mD - mR * mR - mC * mC ) / 4.0 -
361 mR * mR * mC * mC;
362 if ( temp < 0 ) temp = 0;
363 double pD = sqrt( temp / ( mD * mD ) );
364
365 temp = ( mD * mD - mAB * mAB - mC * mC ) * ( mD * mD - mAB * mAB - mC * mC ) / 4.0 -
366 mAB * mAB * mC * mC;
367 if ( temp < 0 ) temp = 0;
368 double pDAB = sqrt( temp / ( mD * mD ) );
369
370 double fR = 1;
371 double fD = 1;
372 int power = 0;
373 switch ( _spin )
374 {
375 case 0:
376 fR = 1.0;
377 fD = 1.0;
378 power = 1;
379 break;
380 case 1:
381 fR = sqrt( 1.0 + 1.5 * 1.5 * pR * pR ) / sqrt( 1.0 + 1.5 * 1.5 * pAB * pAB );
382 fD = sqrt( 1.0 + 5.0 * 5.0 * pD * pD ) / sqrt( 1.0 + 5.0 * 5.0 * pDAB * pDAB );
383 power = 3;
384 break;
385 case 2:
386 fR = sqrt( ( 9 + 3 * pow( ( 1.5 * pR ), 2 ) + pow( ( 1.5 * pR ), 4 ) ) /
387 ( 9 + 3 * pow( ( 1.5 * pAB ), 2 ) + pow( ( 1.5 * pAB ), 4 ) ) );
388 fD = sqrt( ( 9 + 3 * pow( ( 5.0 * pD ), 2 ) + pow( ( 5.0 * pD ), 4 ) ) /
389 ( 9 + 3 * pow( ( 5.0 * pDAB ), 2 ) + pow( ( 5.0 * pDAB ), 4 ) ) );
390 power = 5;
391 break;
392 }
393
394 double gammaAB = gammaR * pow( pAB / pR, power ) * ( mR / mAB ) * fR * fR;
395
396 switch ( _spin )
397 {
398 case 0:
399 ampl = _ampl * TComplex( cos( _theta * pi180inv ), sin( _theta * pi180inv ) ) * fR * fD /
400 ( mR * mR - mAB * mAB - TComplex( 0.0, mR * gammaAB ) );
401 break;
402 case 1:
403 ampl = _ampl * TComplex( cos( _theta * pi180inv ), sin( _theta * pi180inv ) ) * fR * fD *
404 ( mAC * mAC - mBC * mBC +
405 ( mD * mD - mC * mC ) * ( mB * mB - mA * mA ) / ( mR * mR ) ) /
406 ( mR * mR - mAB * mAB - TComplex( 0.0, mR * gammaAB ) );
407 break;
408 case 2:
409 ampl = _ampl * TComplex( cos( _theta * pi180inv ), sin( _theta * pi180inv ) ) * fR * fD /
410 ( mR * mR - mAB * mAB - TComplex( 0.0, mR * gammaAB ) ) *
411 ( pow( ( mBC * mBC - mAC * mAC +
412 ( mD * mD - mC * mC ) * ( mA * mA - mB * mB ) / ( mR * mR ) ),
413 2 ) -
414 ( 1.0 / 3.0 ) *
415 ( mAB * mAB - 2 * mD * mD - 2 * mC * mC +
416 pow( ( mD * mD - mC * mC ) / mR, 2 ) ) *
417 ( mAB * mAB - 2 * mA * mA - 2 * mB * mB +
418 pow( ( mA * mA - mB * mB ) / mR, 2 ) ) );
419 break;
420 }
421
422 return ampl;
423} // end resAmp
424
425// PRL 86, 765 (2001)
426TComplex Dalitz::f_980( double mPP, double mR, double _ampl, double _theta ) {
427
428 double pi180inv = 3.1415926 / 180.;
429 double mK = 0.493677;
430 double mK0 = 0.497648;
431 double mP = 0.13957;
432
433 double m2_PP = mPP * mPP;
434 double gamma = 0.09 * sqrt( m2_PP / 4. - mP * mP );
435 if ( m2_PP / 4. > mK * mK ) gamma = gamma + 0.02 / 2. * sqrt( m2_PP / 4. - mK * mK );
436 if ( m2_PP / 4. > mK0 * mK0 ) gamma = gamma + 0.02 / 2. * sqrt( m2_PP / 4. - mK0 * mK0 );
437
438 // form factor both equal 1
439 TComplex ampl;
440 ampl = _ampl * TComplex( cos( _theta * pi180inv ), sin( _theta * pi180inv ) ) * 1. /
441 TComplex( mR * mR - m2_PP, -mR * gamma );
442
443 return ampl;
444
445} // end f_980
446
447// PRL 21, 244 (1968)
448TComplex Dalitz::sakurai( double mkp, double mkm, double mpp, double _ampl, double _theta,
449 double gamma_r, double m_r ) {
450
451 double pi180inv = 3.1415926 / 180.;
452 double m_pi = 0.139570;
453 double m_k = 0.497648;
454 double mD = 1.86450;
455 double num, m_a, m_b, m_c, m2_ab, m2_ac, m2_bc;
456 double m_ab, m_ac, m_bc, m2_a, m2_b, m2_c, m2_d;
457 m_a = m_pi;
458 m_b = m_pi;
459 m_c = m_k;
460 m_ab = mpp;
461 m_ac = mkp;
462 m_bc = mkm;
463
464 m2_ab = m_ab * m_ab;
465 m2_ac = m_ac * m_ac;
466 m2_bc = m_bc * m_bc;
467 m2_a = m_a * m_a;
468 m2_b = m_b * m_b;
469 m2_c = m_c * m_c;
470 m2_d = mD * mD;
471
472 // for spin 1 angular term
473 num = m2_ac - m2_bc + ( m2_d - m2_c ) * ( m2_b - m2_a ) / ( m_r * m_r );
474
475 double pi, m2, m_pi2, ss, ppi2, p02, ppi, p0;
476 double d, hs, hm, dhdq, f, gamma_s, dr, di;
477
478 pi = 3.14159265358979;
479
480 m2 = m_r * m_r;
481 m_pi2 = m_pi * m_pi;
482 ss = sqrt( m2_ab );
483
484 ppi2 = ( m2_ab - 4. * m_pi2 ) / 4.;
485 p02 = ( m2 - 4. * m_pi2 ) / 4.;
486 p0 = sqrt( p02 );
487 ppi = sqrt( ppi2 );
488
489 d = 3. * m_pi2 / pi / p02 * log( ( m_r + 2. * p0 ) / 2. / m_pi ) + m_r / 2. / pi / p0 -
490 m_pi2 * m_r / pi / ( p0 * p0 * p0 );
491
492 hs = 2. * ppi / pi / ss * log( ( ss + 2. * ppi ) / 2. / m_pi );
493 hm = 2. * p0 / pi / m_r * log( ( m_r + 2. * p0 ) / 2. / m_pi );
494
495 dhdq = hm * ( 1. / 8. / p02 - 1. / 2. / m2 ) + 1. / 2. / pi / m2;
496
497 f = gamma_r * m_r * m_r / ( p0 * p0 * p0 ) *
498 ( ppi2 * ( hs - hm ) - p02 * ( m2_ab - m2 ) * dhdq );
499
500 gamma_s = gamma_r * m2 * ppi * ppi * ppi / ss / ( p0 * p0 * p0 );
501
502 dr = m2 - m2_ab + f;
503 di = gamma_s;
504
505 TComplex ampl = num / TComplex( dr, -di );
506 ampl = _ampl * TComplex( cos( _theta * pi180inv ), sin( _theta * pi180inv ) ) * ampl;
507
508 return ampl;
509}
510
511TComplex Dalitz::Babar_sakurai( double mkp, double mkm, double mpp, double _ampl,
512 double _theta, double gamma_r, double m_r ) {
513
514 double pi180inv = 3.1415926 / 180.;
515 double m_pi = 0.139570;
516 double m_k = 0.497648;
517 double mD = 1.86450;
518 double num, m_a, m_b, m_c, m2_ab, m2_ac, m2_bc;
519 double m_ab, m_ac, m_bc, m2_a, m2_b, m2_c, m2_d;
520 m_a = m_pi;
521 m_b = m_pi;
522 m_c = m_k;
523 m_ab = mpp;
524 m_ac = mkp;
525 m_bc = mkm;
526
527 m2_ab = m_ab * m_ab;
528 m2_ac = m_ac * m_ac;
529 m2_bc = m_bc * m_bc;
530 m2_a = m_a * m_a;
531 m2_b = m_b * m_b;
532 m2_c = m_c * m_c;
533 m2_d = mD * mD;
534
535 // for spin 1 angular term
536 num = m2_ac - m2_bc + ( m2_d - m2_c ) * ( m2_b - m2_a ) / ( m_r * m_r );
537
538 // form factor ---------------------------------------------------
539 double mAB = m_ab;
540 double mA = m_a;
541 double mB = m_b;
542 double mC = m_c;
543 double mR = m_r;
544 double temp = ( mAB * mAB - mA * mA - mB * mB ) * ( mAB * mAB - mA * mA - mB * mB ) / 4.0 -
545 mA * mA * mB * mB;
546 if ( temp < 0 ) temp = 0;
547 double pAB = sqrt( temp / ( mAB * mAB ) );
548
549 temp = ( mR * mR - mA * mA - mB * mB ) * ( mR * mR - mA * mA - mB * mB ) / 4.0 -
550 mA * mA * mB * mB;
551 if ( temp < 0 ) temp = 0;
552 double pR = sqrt( temp / ( mR * mR ) );
553
554 temp = ( mD * mD - mR * mR - mC * mC ) * ( mD * mD - mR * mR - mC * mC ) / 4.0 -
555 mR * mR * mC * mC;
556 if ( temp < 0 ) temp = 0;
557 double pD = sqrt( temp / ( mD * mD ) );
558
559 temp = ( mD * mD - mAB * mAB - mC * mC ) * ( mD * mD - mAB * mAB - mC * mC ) / 4.0 -
560 mAB * mAB * mC * mC;
561 if ( temp < 0 ) temp = 0;
562 double pDAB = sqrt( temp / ( mD * mD ) );
563
564 double fR = sqrt( 1.0 + 1.5 * 1.5 * pR * pR ) / sqrt( 1.0 + 1.5 * 1.5 * pAB * pAB );
565 double fD = sqrt( 1.0 + 5.0 * 5.0 * pD * pD ) / sqrt( 1.0 + 5.0 * 5.0 * pDAB * pDAB );
566 //-----------------------------------------------------------------
567
568 double pi, m2, m_pi2, ss, ppi2, p02, ppi, p0;
569 double d, hs, hm, dhdq, f, gamma_s, dr, di;
570
571 pi = 3.14159265358979;
572
573 m2 = m_r * m_r;
574 m_pi2 = m_pi * m_pi;
575 ss = sqrt( m2_ab );
576
577 ppi2 = ( m2_ab - 4. * m_pi2 ) / 4.;
578 p02 = ( m2 - 4. * m_pi2 ) / 4.;
579 if ( p02 < 0 ) p02 = 0;
580 if ( ppi2 < 0 ) ppi2 = 0;
581 p0 = sqrt( p02 );
582 ppi = sqrt( ppi2 );
583
584 d = 3. * m_pi2 / pi / p02 * log( ( m_r + 2. * p0 ) / 2. / m_pi ) + m_r / 2. / pi / p0 -
585 m_pi2 * m_r / pi / ( p0 * p0 * p0 );
586
587 hs = 2. * ppi / pi / ss * log( ( ss + 2. * ppi ) / 2. / m_pi );
588 hm = 2. * p0 / pi / m_r * log( ( m_r + 2. * p0 ) / 2. / m_pi );
589
590 dhdq = hm * ( 1. / 8. / p02 - 1. / 2. / m2 ) + 1. / 2. / pi / m2;
591
592 f = gamma_r * m_r * m_r / ( p0 * p0 * p0 ) *
593 ( ppi2 * ( hs - hm ) - p02 * ( m2_ab - m2 ) * dhdq );
594
595 gamma_s = gamma_r * m2 * ppi * ppi * ppi / ss / ( p0 * p0 * p0 );
596
597 dr = m2 - m2_ab + f;
598 di = gamma_s;
599
600 //----------------------------------------------------------------------------
601 num *= fR * fD * ( 1 + d * gamma_r / m_r );
602 //----------------------------------------------------------------------------
603 TComplex ampl = num / TComplex( dr, -di );
604 ampl = _ampl * TComplex( cos( _theta * pi180inv ), sin( _theta * pi180inv ) ) * ampl;
605
606 return ampl;
607
608} // end Babar_sakurai
609
610TComplex Dalitz::Babar_resAmp( double mAC, double mBC, double mAB, double mB, double mA,
611 double mC, double _ampl, double _theta, double _gamma,
612 double _bwm, int _spin ) {
613
614 double pi180inv = 3.1415926 / 180.;
615
616 TComplex ampl;
617
618 double mD = 1.86450;
619
620 double mR = _bwm;
621 double gammaR = _gamma;
622
623 double temp = ( mAB * mAB - mA * mA - mB * mB ) * ( mAB * mAB - mA * mA - mB * mB ) / 4.0 -
624 mA * mA * mB * mB;
625 if ( temp < 0 ) temp = 0;
626 double pAB = sqrt( temp / ( mAB * mAB ) );
627
628 temp = ( mR * mR - mA * mA - mB * mB ) * ( mR * mR - mA * mA - mB * mB ) / 4.0 -
629 mA * mA * mB * mB;
630 if ( temp < 0 ) temp = 0;
631 double pR = sqrt( temp / ( mR * mR ) );
632
633 temp = ( mD * mD - mR * mR - mC * mC ) * ( mD * mD - mR * mR - mC * mC ) / 4.0 -
634 mR * mR * mC * mC;
635 if ( temp < 0 ) temp = 0;
636 double pD = sqrt( temp / ( mD * mD ) );
637
638 temp = ( mD * mD - mAB * mAB - mC * mC ) * ( mD * mD - mAB * mAB - mC * mC ) / 4.0 -
639 mAB * mAB * mC * mC;
640 if ( temp < 0 ) temp = 0;
641 double pDAB = sqrt( temp / ( mD * mD ) );
642
643 double fR = 1;
644 double fD = 1;
645 int power = 0;
646 switch ( _spin )
647 {
648 case 0:
649 fR = 1.0;
650 fD = 1.0;
651 power = 1;
652 break;
653 case 1:
654 fR = sqrt( 1.0 + 1.5 * 1.5 * pR * pR ) / sqrt( 1.0 + 1.5 * 1.5 * pAB * pAB );
655 fD = sqrt( 1.0 + 5.0 * 5.0 * pD * pD ) / sqrt( 1.0 + 5.0 * 5.0 * pDAB * pDAB );
656 power = 3;
657 break;
658 case 2:
659 fR = sqrt( ( 9 + 3 * pow( ( 1.5 * pR ), 2 ) + pow( ( 1.5 * pR ), 4 ) ) /
660 ( 9 + 3 * pow( ( 1.5 * pAB ), 2 ) + pow( ( 1.5 * pAB ), 4 ) ) );
661 fD = sqrt( ( 9 + 3 * pow( ( 5.0 * pD ), 2 ) + pow( ( 5.0 * pD ), 4 ) ) /
662 ( 9 + 3 * pow( ( 5.0 * pDAB ), 2 ) + pow( ( 5.0 * pDAB ), 4 ) ) );
663 power = 5;
664 break;
665 }
666
667 double gammaAB = gammaR * pow( pAB / pR, power ) * ( mR / mAB ) * fR * fR;
668
669 switch ( _spin )
670 {
671 case 0:
672 ampl = _ampl * TComplex( cos( _theta * pi180inv ), sin( _theta * pi180inv ) ) * fR * fD /
673 ( mR * mR - mAB * mAB - TComplex( 0.0, mR * gammaAB ) );
674 break;
675 case 1:
676 ampl = _ampl * TComplex( cos( _theta * pi180inv ), sin( _theta * pi180inv ) ) * fR * fD *
677 ( mAC * mAC - mBC * mBC +
678 ( mD * mD - mC * mC ) * ( mB * mB - mA * mA ) / ( mAB * mAB ) ) /
679 ( mR * mR - mAB * mAB - TComplex( 0.0, mR * gammaAB ) );
680 break;
681 case 2:
682 ampl = _ampl * TComplex( cos( _theta * pi180inv ), sin( _theta * pi180inv ) ) * fR * fD /
683 ( mR * mR - mAB * mAB - TComplex( 0.0, mR * gammaAB ) ) *
684 ( pow( ( mBC * mBC - mAC * mAC +
685 ( mD * mD - mC * mC ) * ( mA * mA - mB * mB ) / ( mAB * mAB ) ),
686 2 ) -
687 ( 1.0 / 3.0 ) *
688 ( mAB * mAB - 2 * mD * mD - 2 * mC * mC +
689 pow( ( mD * mD - mC * mC ) / mAB, 2 ) ) *
690 ( mAB * mAB - 2 * mA * mA - 2 * mB * mB +
691 pow( ( mA * mA - mB * mB ) / mAB, 2 ) ) );
692 break;
693 }
694
695 return ampl;
696}
697
698TComplex Dalitz::Babar_Amplitude( double x, double y, double z ) {
699
700 // PRD 73, 112009(2006) Belle
701 // for D0 particle 1: K particle 2: pi- particle 3: pi+
702 // the right order is: (1,2), (1,3), (3,2)
703 double m_mass[4] = { 1.86450, 0.497648, 0.139570, 0.139570 }; // mass
704
705 TComplex D0( 0.0, 0.0 );
706
707 // x, y, z already squared, need to get back the mass!!!!
708 x = sqrt( x );
709 y = sqrt( y );
710 z = sqrt( z );
711
712 TComplex DK2piRes[17];
713 DK2piRes[0] = Babar_sakurai( x, y, z, 1.00, 0.0, 0.1464, 0.7758 ); // RHO(770)
714 DK2piRes[1] = Babar_resAmp( x, y, z, m_mass[3], m_mass[2], m_mass[1], 0.0391, 115.3, 0.00849,
715 0.78259, 1 ); // OMEGA
716 DK2piRes[2] = Babar_resAmp( x, y, z, m_mass[3], m_mass[2], m_mass[1], 0.482, -141.8, 0.044,
717 0.975, 0 ); // F_0(980)
718 DK2piRes[3] = Babar_resAmp( x, y, z, m_mass[3], m_mass[2], m_mass[1], 0.922, -21.3, 0.1851,
719 1.2754, 2 ); // F_2(1270)
720 DK2piRes[4] = Babar_resAmp( x, y, z, m_mass[3], m_mass[2], m_mass[1], 2.25, 113.2, 0.173,
721 1.434, 0 ); // F_0(1370) hep-ph 0009168
722 DK2piRes[5] = Babar_sakurai( x, y, z, 0.52, 38, 0.455, 1.406 ); // RHO(1450)
723 DK2piRes[6] = Babar_resAmp( x, y, z, m_mass[3], m_mass[2], m_mass[1], 1.36, -177.9, 0.383,
724 0.484, 0 ); // Sigma(600)
725 DK2piRes[7] = Babar_resAmp( x, y, z, m_mass[3], m_mass[2], m_mass[1], 0.340, 153.0, 0.088,
726 1.014, 0 ); // Sigma
727 DK2piRes[8] = Babar_resAmp( x, z, y, m_mass[3], m_mass[1], m_mass[2], 1.781, 131.0, 0.0508,
728 0.89166, 1 ); // K*(892)-
729 DK2piRes[9] = Babar_resAmp( x, z, y, m_mass[3], m_mass[1], m_mass[2], 0.52, 154, 0.232,
730 1.414, 1 ); // K*(1410)-
731 DK2piRes[10] = Babar_resAmp( x, z, y, m_mass[3], m_mass[1], m_mass[2], 2.45, -8.3, 0.294,
732 1.412, 0 ); // K*_0(1430)-
733 DK2piRes[11] = Babar_resAmp( x, z, y, m_mass[3], m_mass[1], m_mass[2], 1.05, -54.3, 0.0985,
734 1.4256, 2 ); // K*_2(1430)-
735 DK2piRes[12] = Babar_resAmp( x, z, y, m_mass[3], m_mass[1], m_mass[2], 0.89, -139, 0.322,
736 1.717, 1 ); // K*(1680)-
737 DK2piRes[13] = Babar_resAmp( y, z, x, m_mass[2], m_mass[1], m_mass[3], 0.180, -44.1, 0.0508,
738 0.89166, 1 ); // K*(892)+
739 DK2piRes[14] = Babar_resAmp( y, z, x, m_mass[2], m_mass[1], m_mass[3], 0.37, 18, 0.294,
740 1.412, 0 ); // K*_0(1430)+
741 DK2piRes[15] = Babar_resAmp( y, z, x, m_mass[2], m_mass[1], m_mass[3], 0.075, -104, 0.0985,
742 1.4256, 2 ); // K*_2(1430)+
743
744 double pi180inv = 3.1415926 / 180.;
745 DK2piRes[16] = TComplex( 3.53 * cos( 128 * pi180inv ), 3.53 * sin( 128 * pi180inv ) );
746
747 for ( int i = 0; i < 17; i++ ) { D0 += DK2piRes[i]; }
748
749 return D0;
750} // End Babar_amplitude
751
752int Dalitz::getBin( double mx, double my, double mz ) {
753
754 // For computing Dalitz variable kinemtaics
755 double m_mass_2[4] = { 1.86450 * 1.86450, 0.497648 * 0.497648, 0.139570 * 0.139570,
756 0.139570 * 0.139570 }; // mass square
757 double m_sum_m_2 = m_mass_2[0] + m_mass_2[1] + m_mass_2[2] + m_mass_2[3];
758
759 // Check if on DP
760 if ( !( Point_on_DP2( mx, my ) ) && !( Point_on_DP2( my, mx ) ) ) return -1;
761
762 // Over-ride user z, use computed z(x,y) instead
763 mz = m_sum_m_2 - mx - my;
764 if ( mz < 0 ) mz = 0;
765
766 // Determine phase
767 double thisPhase = Phase( mx, my, mz );
768
769 // Change this to bin that is found (-1 means not found)
770 int thisbin = -1;
771
772 // Determine the bin and increment
773 for ( int bin = 0; bin < N; bin++ )
774 {
775 if ( ( thisPhase >= ( bin - 0.5 ) * 2 * PI / N ) &&
776 ( thisPhase < ( bin + 0.5 ) * 2 * PI / N ) )
777 thisbin = bin;
778 }
779
780 return thisbin;
781
782} // end getBin
783
784// ------------ end Dalitz.cxx
m_pi
*********DOUBLE PRECISION m_pi
Definition BVR.h:11
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
TComplex
complex< double > TComplex
Definition Dalitz.h:14
arg
double arg(const EvtComplex &c)
Definition EvtComplex.hh:195
pi
double pi
Definition EvtPhokharaDef.hh:17
bin
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per bin
Definition FoamA.h:85
m_b
*********DOUBLE PRECISION m_pi INTEGER m_lenwt !max no of aux weights INTEGER m_phmax !maximum photon multiplicity ISR FSR *DOUBLE COMPLEX m_Pauli4 DOUBLE COMPLEX m_AmpBorn DOUBLE COMPLEX m_AmpBoxy DOUBLE COMPLEX m_AmpBorn1 DOUBLE COMPLEX m_AmpBorn2 DOUBLE COMPLEX m_AmpExpo2p DOUBLE COMPLEX m_Rmat DOUBLE COMPLEX m_BoxGZut !DOUBLE COMPLEX m_F1finPair2 !DOUBLE PRECISION m_Vcut DOUBLE PRECISION m_Alfinv DOUBLE PRECISION m_Lorin1 DOUBLE PRECISION m_Lorin2 DOUBLE PRECISION m_b
Definition GPS.h:30
sin
double sin(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:185
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
Dalitz::Amplitude
TComplex Amplitude(double x, double y, double z)
Definition Dalitz.cxx:30
Dalitz::Point_on_DP2
bool Point_on_DP2(double x, double y)
Definition Dalitz.cxx:193
Dalitz::resAmp
TComplex resAmp(double mAC, double mBC, double mAB, double mA, double mB, double mC, double _ampl, double _theta, double _gamma, double _bwm, int _spin)
Definition Dalitz.cxx:338
Dalitz::CLEO_resAmp
TComplex CLEO_resAmp(double mAC, double mBC, double mAB, double mA, double mB, double mC, double _ampl, double _theta, double _gamma, double _bwm, int _spin)
Definition Dalitz.cxx:241
Dalitz::CLEO_Amplitude
TComplex CLEO_Amplitude(double x, double y, double z)
Definition Dalitz.cxx:90
Dalitz::Dalitz
Dalitz()
Definition Dalitz.cxx:21
Dalitz::Phase
double Phase(double x, double y, double z, int Babar=1)
Definition Dalitz.cxx:116
Dalitz::Babar_resAmp
TComplex Babar_resAmp(double mAC, double mBC, double mAB, double mB, double mA, double mC, double _ampl, double _theta, double _gamma, double _bwm, int _spin)
Definition Dalitz.cxx:610
Dalitz::Babar_Amplitude
TComplex Babar_Amplitude(double x, double y, double z)
Definition Dalitz.cxx:698
Dalitz::sakurai
TComplex sakurai(double mkp, double mkm, double mpp, double _ampl, double _theta, double gamma_r, double m_r)
Definition Dalitz.cxx:448
Dalitz::getBin
int getBin(double mx, double my, double mz)
Definition Dalitz.cxx:752
Dalitz::f_980
TComplex f_980(double mPP, double mR, double _ampl, double _theta)
Definition Dalitz.cxx:426
Dalitz::Point_on_DP
bool Point_on_DP(double x, double y)
Definition Dalitz.cxx:143
Dalitz::Babar_sakurai
TComplex Babar_sakurai(double mkp, double mkm, double mpp, double _ampl, double _theta, double gamma_r, double m_r)
Definition Dalitz.cxx:511
m2
double double * m2
Definition qcdloop1.h:83