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