00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "QCMCFilterAlg/Dalitz.h"
00014
00015 #include <complex>
00016
00017 using std::string;
00018
00019 double Dalitz::PI = 3.1415926;
00020
00021 Dalitz::Dalitz() {
00022 N = 8;
00023 }
00024
00025
00026 Dalitz::Dalitz(int binNum) {
00027 N = binNum;
00028 }
00029
00030
00031 TComplex Dalitz::Amplitude(double x, double y, double z) {
00032
00033
00034
00035
00036 double m_mass[4] = {1.86450, 0.497648, 0.139570, 0.139570};
00037
00038 TComplex D0(0.0,0.0);
00039
00040
00041 x = sqrt(x);
00042 y = sqrt(y);
00043 z = sqrt(z);
00044
00045
00046 TComplex DK2piRes[19];
00047
00048
00049 DK2piRes[0] = sakurai(x, y, z, 1.00, 0.0, 0.1503, 0.7758);
00050 DK2piRes[1] = resAmp(x, y, z, m_mass[3], m_mass[2], m_mass[1], 0.0314,110.8, 0.00849,0.78259,1);
00051 DK2piRes[2] = f_980(z, 0.980, 0.365, 201.9);
00052 DK2piRes[3] = resAmp(x, y, z, m_mass[3], m_mass[2], m_mass[1], 1.32, 348, 0.1851, 1.2754, 2);
00053 DK2piRes[4] = resAmp(x, y, z, m_mass[3], m_mass[2], m_mass[1], 1.44, 82, 0.173, 1.434, 0);
00054 DK2piRes[5] = sakurai(x, y, z, 0.66, 9, 0.400, 1.465);
00055 DK2piRes[6] = resAmp(x, y, z, m_mass[3], m_mass[2], m_mass[1], 1.43, 212, 0.454, 0.519, 0);
00056 DK2piRes[7] = resAmp(x, y, z, m_mass[3], m_mass[2], m_mass[1], 0.23, 237, 0.101, 1.050, 0);
00057 DK2piRes[8] = resAmp(x, z, y, m_mass[3], m_mass[1], m_mass[2], 1.644, 132.1, 0.0508, 0.89166,1);
00058 DK2piRes[9] = resAmp(x, z, y, m_mass[3], m_mass[1], m_mass[2], 0.61, 113, 0.232, 1.414, 1);
00059 DK2piRes[10] = resAmp(x, z, y, m_mass[3], m_mass[1], m_mass[2], 2.15, 353.6, 0.294, 1.412, 0);
00060 DK2piRes[11] = resAmp(x, z, y, m_mass[3], m_mass[1], m_mass[2], 0.88, 318.7, 0.0985, 1.4256, 2);
00061 DK2piRes[12] = resAmp(x, z, y, m_mass[3], m_mass[1], m_mass[2], 1.39, 103, 0.322, 1.717, 1);
00062 DK2piRes[13] = resAmp(y, z, x, m_mass[2], m_mass[1], m_mass[3], 0.144, 320.3, 0.0508, 0.89166,1);
00063 DK2piRes[14] = resAmp(y, z, x, m_mass[2], m_mass[1], m_mass[3], 0.45, 254, 0.232, 1.414, 1);
00064 DK2piRes[15] = resAmp(y, z, x, m_mass[2], m_mass[1], m_mass[3], 0.47, 88, 0.294, 1.412, 0);
00065 DK2piRes[16] = resAmp(y, z, x, m_mass[2], m_mass[1], m_mass[3], 0.25, 265, 0.0985, 1.4256, 2);
00066 DK2piRes[17] = resAmp(y, z, x, m_mass[2], m_mass[1], m_mass[3], 1.2, 118, 0.322, 1.717, 1);
00067
00068 double pi180inv = 3.1415926/180.;
00069 DK2piRes[18] = TComplex(3.0*cos(164*pi180inv),3.0*sin(164*pi180inv));
00070
00071 for(int i=0; i<19; i++){
00072 D0 += DK2piRes[i];
00073 }
00074
00075 return D0;
00076 }
00077
00078 TComplex Dalitz::CLEO_Amplitude(double x, double y, double z) {
00079
00080 double m_mass[4] = {1.86450, 0.497648, 0.139570, 0.139570};
00081
00082 TComplex D0(0.0,0.0);
00083
00084
00085 x = sqrt(x);
00086 y = sqrt(y);
00087 z = sqrt(z);
00088
00089
00090 TComplex DK2piRes[3];
00091
00092
00093 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);
00094 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);
00095 DK2piRes[2] = TComplex(1.0, 0.0) ;
00096
00097 for(int i=0; i<3; i++){
00098 D0 += DK2piRes[i];
00099 }
00100
00101 return D0;
00102 }
00103
00104 double Dalitz::Phase(double x, double y, double z, int Babar) {
00105
00106 TComplex D0(0,0);
00107 TComplex D0bar(0,0);
00108
00109 if (Babar == 1) {
00110 D0 = Babar_Amplitude(x, y, z);
00111 D0bar = Babar_Amplitude(y, x, z);
00112 } else{
00113 D0 = Amplitude(x, y, z);
00114 D0bar = Amplitude(y, x, z);
00115 }
00116
00117
00118 double deltaD = arg(D0) - arg(D0bar);
00119 if(x<y) deltaD = -deltaD;
00120
00121 if ( deltaD < -PI/N ) deltaD += 2*PI;
00122 if ( deltaD > (2*N-1)*PI/N ) deltaD -= 2*PI;
00123
00124 return deltaD;
00125
00126 }
00127
00128 bool Dalitz::Point_on_DP(double x, double y) {
00129
00130 double m_mass[4] = {1.86450, 0.497648, 0.139570, 0.139570};
00131 double m_mass2[4]= {1.86450*1.86450, 0.497648*0.497648,
00132 0.139570*0.139570, 0.139570*0.139570};
00133
00134 double m_XmaxDP = m_mass[0] - m_mass[3]; m_XmaxDP *= m_XmaxDP;
00135 double m_XminDP = m_mass[1] + m_mass[2]; m_XminDP *= m_XminDP;
00136
00137 if ( (x > m_XmaxDP) || (x < m_XminDP) ) return false;
00138
00139 double Low = 0;
00140 double Up = 0;
00141 double HInv_m12 = 0.5/sqrt(x);
00142 double E1 = HInv_m12*(x + m_mass2[1] - m_mass2[2]);
00143 double E3 = HInv_m12*(m_mass2[0] - m_mass2[3] - x);
00144 double E1_2 = E1*E1;
00145 double E3_2 = E3*E3;
00146
00147 if (E1 < m_mass[1]) { E1=m_mass[1]; E1_2=m_mass2[1]; }
00148 if (E3 < m_mass[3]) { E3=m_mass[3]; E3_2=m_mass2[3]; }
00149
00150 double temp = E1_2-m_mass2[1];
00151 if (temp < 0) temp = 0;
00152 double P1 = sqrt(temp);
00153 temp = E3_2 - m_mass2[3];
00154 if (temp<0) temp = 0;
00155 double P3 = sqrt(temp);
00156 double E13_2 = (E1+E3)*(E1+E3);
00157
00158
00159 Low = E13_2 - (P1+P3)*(P1+P3);
00160 Up = E13_2 - (P1-P3)*(P1-P3);
00161
00162 if ( (y > Up) || (y < Low) ) return false;
00163
00164 return true;
00165 }
00166
00167
00168 bool Dalitz::Point_on_DP2(double x,double y) {
00169
00170 double m_mass[4] = {1.86450, 0.497648, 0.139570, 0.139570};
00171 double m_mass2[4] = {1.86450*1.86450, 0.497648*0.497648,
00172 0.139570*0.139570, 0.139570*0.139570};
00173
00174 double m_XmaxDP = m_mass[0] - m_mass[3]; m_XmaxDP *= m_XmaxDP;
00175 double m_XminDP = m_mass[1] + m_mass[2]; m_XminDP *= m_XminDP;
00176
00177 if ( (x > m_XmaxDP) || (x < m_XminDP) ) return false;
00178
00179 double Low = 0;
00180 double Up = 0;
00181 double HInv_m12 = 0.5/sqrt(x);
00182 double E1 = HInv_m12*(x + m_mass2[1] - m_mass2[2]);
00183 double E3 = HInv_m12*(m_mass2[0] - m_mass2[3] - x);
00184 double E1_2 = E1*E1;
00185 double E3_2 = E3*E3;
00186
00187 if (E1 < m_mass[1]) { E1=m_mass[1]; E1_2=m_mass2[1]; }
00188 if (E3 < m_mass[3]) { E3=m_mass[3]; E3_2=m_mass2[3]; }
00189
00190 double temp = E1_2-m_mass2[1];
00191 if (temp < 0) temp = 0;
00192 double P1 = sqrt(temp);
00193 temp = E3_2-m_mass2[3];
00194 if (temp < 0) temp = 0;
00195 double P3 = sqrt(temp);
00196 double E13_2 = (E1+E3)*(E1+E3);
00197
00198 Low = E13_2 - (P1+P3)*(P1+P3);
00199 Up = E13_2 - (P1-P3)*(P1-P3);
00200
00201 if ( (y > (Up+0.05)) || (y < (Low-0.05)) ) return false;
00202
00203 return true;
00204 }
00205
00206
00207 TComplex Dalitz::CLEO_resAmp(double mAC, double mBC, double mAB,
00208 double mB , double mA , double mC ,
00209 double _ampl, double _theta, double _gamma, double _bwm, int _spin) {
00210
00211 double pi180inv = 3.1415926/180.;
00212
00213 TComplex ampl;
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244 double mD = 1.86484 ;
00245 double eA = ( mD*mD - mBC*mBC + mA*mA ) / (2.*mD) ;
00246 double eAB = ( mD*mD - mC*mC + mAB*mAB ) / (2.*mD) ;
00247
00248
00249 double pd = mD * eA ;
00250 double pq = mD * eAB ;
00251 double qd = mA*mA + 0.5 * ( mAB*mAB - mA*mA - mB*mB ) ;
00252 double mp2 = mD*mD ;
00253 double mq2 = mAB*mAB ;
00254 double md2 = mA*mA ;
00255 double cos_phi_0 = (pd*mq2-pq*qd)/sqrt((pq*pq-mq2*mp2)*(qd*qd-mq2*md2));
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267 switch (_spin) {
00268
00269 case 0 :
00270 ampl=(_ampl*TComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
00271 sqrt(_gamma/(2.*PI))*
00272 (1.0/(mAB-_bwm-TComplex(0.0,0.5*_gamma))));
00273 break;
00274
00275 case 1 :
00276 ampl=(_ampl*TComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
00277 sqrt(_gamma/(2.*PI))*
00278 (cos_phi_0/(mAB-_bwm-TComplex(0.0,0.5*_gamma))));
00279 break;
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293 default:
00294
00295 ampl = TComplex(0.0);
00296 break;
00297
00298 }
00299
00300 return ampl;
00301 }
00302
00303
00304 TComplex Dalitz::resAmp(double mAC, double mBC, double mAB,
00305 double mB , double mA , double mC ,
00306 double _ampl, double _theta, double _gamma, double _bwm, int _spin) {
00307
00308 double pi180inv = 3.1415926/180.;
00309
00310 TComplex ampl;
00311
00312 double mD = 1.86450;
00313
00314 double mR = _bwm;
00315 double gammaR = _gamma;
00316
00317 double temp = (mAB*mAB-mA*mA-mB*mB)*(mAB*mAB-mA*mA-mB*mB)/4.0- mA*mA*mB*mB;
00318 if (temp < 0) temp = 0;
00319 double pAB = sqrt( temp/(mAB*mAB) );
00320
00321 temp = (mR*mR-mA*mA-mB*mB)*(mR*mR-mA*mA-mB*mB)/4.0 - mA*mA*mB*mB;
00322 if (temp<0) temp = 0;
00323 double pR = sqrt( temp/(mR*mR));
00324
00325 temp = (mD*mD-mR*mR-mC*mC)*(mD*mD-mR*mR-mC*mC)/4.0 - mR*mR*mC*mC;
00326 if (temp < 0) temp = 0;
00327 double pD = sqrt( temp/(mD*mD) );
00328
00329 temp = (mD*mD-mAB*mAB-mC*mC)*(mD*mD-mAB*mAB-mC*mC)/4.0 - mAB*mAB*mC*mC;
00330 if (temp<0) temp = 0;
00331 double pDAB = sqrt( temp/(mD*mD));
00332
00333 double fR = 1;
00334 double fD = 1;
00335 int power = 0;
00336 switch (_spin) {
00337 case 0:
00338 fR = 1.0;
00339 fD = 1.0;
00340 power = 1;
00341 break;
00342 case 1:
00343 fR = sqrt(1.0+1.5*1.5*pR*pR)/sqrt(1.0+1.5*1.5*pAB*pAB);
00344 fD = sqrt(1.0+5.0*5.0*pD*pD)/sqrt(1.0+5.0*5.0*pDAB*pDAB);
00345 power = 3;
00346 break;
00347 case 2:
00348 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)) );
00349 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)) );
00350 power = 5;
00351 break;
00352 }
00353
00354 double gammaAB= gammaR*pow(pAB/pR,power)*(mR/mAB)*fR*fR;
00355
00356 switch (_spin) {
00357 case 0:
00358 ampl = _ampl*TComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
00359 fR*fD/(mR*mR-mAB*mAB-TComplex(0.0,mR*gammaAB));
00360 break;
00361 case 1:
00362 ampl = _ampl*TComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
00363 fR*fD*(mAC*mAC-mBC*mBC+(mD*mD-mC*mC)*(mB*mB-mA*mA)/(mR*mR))/
00364 (mR*mR-mAB*mAB-TComplex(0.0,mR*gammaAB));
00365 break;
00366 case 2:
00367 ampl = _ampl*TComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
00368 fR*fD/(mR*mR-mAB*mAB-TComplex(0.0,mR*gammaAB))*
00369 (pow((mBC*mBC-mAC*mAC+(mD*mD-mC*mC)*(mA*mA-mB*mB)/(mR*mR)),2)-
00370 (1.0/3.0)*(mAB*mAB-2*mD*mD-2*mC*mC+pow((mD*mD- mC*mC)/mR, 2))*
00371 (mAB*mAB-2*mA*mA-2*mB*mB+pow((mA*mA-mB*mB)/mR,2)));
00372 break;
00373 }
00374
00375 return ampl;
00376 }
00377
00378
00379 TComplex Dalitz::f_980(double mPP, double mR,
00380 double _ampl, double _theta ) {
00381
00382 double pi180inv = 3.1415926/180.;
00383 double mK = 0.493677;
00384 double mK0 = 0.497648;
00385 double mP = 0.13957;
00386
00387 double m2_PP = mPP*mPP;
00388 double gamma = 0.09*sqrt(m2_PP/4.-mP*mP);
00389 if( m2_PP/4.>mK*mK ) gamma = gamma + 0.02/2.*sqrt(m2_PP/4.-mK*mK);
00390 if( m2_PP/4.>mK0*mK0 ) gamma = gamma + 0.02/2.*sqrt(m2_PP/4.-mK0*mK0);
00391
00392
00393 TComplex ampl;
00394 ampl = _ampl*TComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
00395 1./TComplex(mR*mR-m2_PP, -mR*gamma);
00396
00397 return ampl;
00398
00399 }
00400
00401
00402 TComplex Dalitz::sakurai(double mkp, double mkm, double mpp, double _ampl, double _theta,
00403 double gamma_r, double m_r) {
00404
00405 double pi180inv = 3.1415926/180.;
00406 double m_pi = 0.139570;
00407 double m_k = 0.497648;
00408 double mD = 1.86450;
00409 double num, m_a, m_b, m_c, m2_ab, m2_ac, m2_bc;
00410 double m_ab, m_ac, m_bc, m2_a, m2_b, m2_c, m2_d;
00411 m_a = m_pi;
00412 m_b = m_pi;
00413 m_c = m_k;
00414 m_ab = mpp;
00415 m_ac = mkp;
00416 m_bc = mkm;
00417
00418 m2_ab = m_ab*m_ab;
00419 m2_ac = m_ac*m_ac;
00420 m2_bc = m_bc*m_bc;
00421 m2_a = m_a*m_a;
00422 m2_b = m_b*m_b;
00423 m2_c = m_c*m_c;
00424 m2_d = mD*mD;
00425
00426
00427 num = m2_ac-m2_bc+(m2_d-m2_c)*(m2_b-m2_a)/(m_r*m_r);
00428
00429 double pi, m2, m_pi2, ss, ppi2, p02, ppi, p0;
00430 double d, hs, hm, dhdq, f, gamma_s, dr, di;
00431
00432 pi = 3.14159265358979;
00433
00434 m2 = m_r*m_r;
00435 m_pi2 = m_pi*m_pi;
00436 ss = sqrt(m2_ab);
00437
00438 ppi2 = (m2_ab-4.*m_pi2)/4.;
00439 p02 = (m2-4.*m_pi2)/4.;
00440 p0 = sqrt(p02);
00441 ppi = sqrt(ppi2);
00442
00443 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);
00444
00445 hs = 2.*ppi/pi/ss*log((ss+2.*ppi)/2./m_pi);
00446 hm = 2.*p0/pi/m_r*log((m_r+2.*p0)/2./m_pi);
00447
00448 dhdq = hm*(1./8./p02 - 1./2./m2) + 1./2./pi/m2;
00449
00450 f = gamma_r*m_r*m_r/(p0*p0*p0)*(ppi2*(hs-hm) - p02*(m2_ab-m2)*dhdq);
00451
00452 gamma_s = gamma_r*m2*ppi*ppi*ppi/ss/(p0*p0*p0);
00453
00454 dr = m2-m2_ab+f;
00455 di = gamma_s;
00456
00457 TComplex ampl = num/TComplex(dr, -di);
00458 ampl = _ampl*TComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*ampl;
00459
00460 return ampl;
00461
00462 }
00463
00464 TComplex Dalitz::Babar_sakurai(double mkp, double mkm, double mpp, double _ampl, double _theta,
00465 double gamma_r, double m_r) {
00466
00467 double pi180inv = 3.1415926/180.;
00468 double m_pi = 0.139570;
00469 double m_k = 0.497648;
00470 double mD = 1.86450;
00471 double num, m_a, m_b, m_c, m2_ab, m2_ac, m2_bc;
00472 double m_ab, m_ac, m_bc, m2_a, m2_b, m2_c, m2_d;
00473 m_a = m_pi;
00474 m_b = m_pi;
00475 m_c = m_k;
00476 m_ab = mpp;
00477 m_ac = mkp;
00478 m_bc = mkm;
00479
00480 m2_ab = m_ab*m_ab;
00481 m2_ac = m_ac*m_ac;
00482 m2_bc = m_bc*m_bc;
00483 m2_a = m_a*m_a;
00484 m2_b = m_b*m_b;
00485 m2_c = m_c*m_c;
00486 m2_d = mD*mD;
00487
00488
00489 num=m2_ac-m2_bc+(m2_d-m2_c)*(m2_b-m2_a)/(m_r*m_r);
00490
00491
00492 double mAB = m_ab;
00493 double mA = m_a;
00494 double mB = m_b;
00495 double mC = m_c;
00496 double mR = m_r;
00497 double temp = (mAB*mAB-mA*mA-mB*mB)*(mAB*mAB-mA*mA-mB*mB)/4.0- mA*mA*mB*mB;
00498 if (temp < 0) temp = 0;
00499 double pAB = sqrt( temp/(mAB*mAB) );
00500
00501 temp = (mR*mR-mA*mA-mB*mB)*(mR*mR-mA*mA-mB*mB)/4.0 - mA*mA*mB*mB;
00502 if (temp < 0) temp = 0;
00503 double pR = sqrt( temp/(mR*mR));
00504
00505 temp = (mD*mD-mR*mR-mC*mC)*(mD*mD-mR*mR-mC*mC)/4.0 - mR*mR*mC*mC;
00506 if (temp < 0) temp = 0;
00507 double pD = sqrt( temp/(mD*mD) );
00508
00509 temp = (mD*mD-mAB*mAB-mC*mC)*(mD*mD-mAB*mAB-mC*mC)/4.0 - mAB*mAB*mC*mC;
00510 if (temp < 0) temp = 0;
00511 double pDAB = sqrt( temp/(mD*mD));
00512
00513 double fR = sqrt(1.0+1.5*1.5*pR*pR)/sqrt(1.0+1.5*1.5*pAB*pAB);
00514 double fD = sqrt(1.0+5.0*5.0*pD*pD)/sqrt(1.0+5.0*5.0*pDAB*pDAB);
00515
00516
00517 double pi,m2,m_pi2,ss,ppi2,p02,ppi,p0;
00518 double d,hs,hm,dhdq,f,gamma_s,dr,di;
00519
00520 pi = 3.14159265358979;
00521
00522 m2 = m_r*m_r;
00523 m_pi2 = m_pi*m_pi;
00524 ss = sqrt(m2_ab);
00525
00526 ppi2 = (m2_ab-4.*m_pi2)/4.;
00527 p02 = (m2-4.*m_pi2)/4.;
00528 if (p02 < 0) p02 = 0;
00529 if (ppi2 < 0) ppi2 = 0;
00530 p0 = sqrt(p02);
00531 ppi = sqrt(ppi2);
00532
00533 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);
00534
00535 hs = 2.*ppi/pi/ss*log((ss+2.*ppi)/2./m_pi);
00536 hm = 2.*p0/pi/m_r*log((m_r+2.*p0)/2./m_pi);
00537
00538 dhdq = hm*(1./8./p02 - 1./2./m2) + 1./2./pi/m2;
00539
00540 f = gamma_r*m_r*m_r/(p0*p0*p0)*(ppi2*(hs-hm) - p02*(m2_ab-m2)*dhdq);
00541
00542 gamma_s = gamma_r*m2*ppi*ppi*ppi/ss/(p0*p0*p0);
00543
00544 dr = m2-m2_ab+f;
00545 di = gamma_s;
00546
00547
00548 num *= fR*fD*(1+d*gamma_r/m_r);
00549
00550 TComplex ampl = num/TComplex(dr, -di);
00551 ampl = _ampl*TComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*ampl;
00552
00553 return ampl;
00554
00555 }
00556
00557
00558 TComplex Dalitz::Babar_resAmp(double mAC, double mBC, double mAB,
00559 double mB , double mA , double mC ,
00560 double _ampl, double _theta, double _gamma, double _bwm, int _spin) {
00561
00562 double pi180inv = 3.1415926/180.;
00563
00564 TComplex ampl;
00565
00566 double mD = 1.86450;
00567
00568 double mR = _bwm;
00569 double gammaR = _gamma;
00570
00571 double temp = (mAB*mAB-mA*mA-mB*mB)*(mAB*mAB-mA*mA-mB*mB)/4.0- mA*mA*mB*mB;
00572 if (temp < 0) temp = 0;
00573 double pAB = sqrt( temp/(mAB*mAB) );
00574
00575 temp = (mR*mR-mA*mA-mB*mB)*(mR*mR-mA*mA-mB*mB)/4.0 - mA*mA*mB*mB;
00576 if (temp < 0) temp = 0;
00577 double pR = sqrt( temp/(mR*mR));
00578
00579 temp = (mD*mD-mR*mR-mC*mC)*(mD*mD-mR*mR-mC*mC)/4.0 - mR*mR*mC*mC;
00580 if (temp < 0) temp = 0;
00581 double pD = sqrt( temp/(mD*mD) );
00582
00583 temp = (mD*mD-mAB*mAB-mC*mC)*(mD*mD-mAB*mAB-mC*mC)/4.0 - mAB*mAB*mC*mC;
00584 if (temp < 0) temp = 0;
00585 double pDAB = sqrt( temp/(mD*mD));
00586
00587 double fR = 1;
00588 double fD = 1;
00589 int power = 0;
00590 switch (_spin) {
00591 case 0:
00592 fR = 1.0;
00593 fD = 1.0;
00594 power = 1;
00595 break;
00596 case 1:
00597 fR = sqrt(1.0+1.5*1.5*pR*pR)/sqrt(1.0+1.5*1.5*pAB*pAB);
00598 fD = sqrt(1.0+5.0*5.0*pD*pD)/sqrt(1.0+5.0*5.0*pDAB*pDAB);
00599 power = 3;
00600 break;
00601 case 2:
00602 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)) );
00603 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)) );
00604 power = 5;
00605 break;
00606 }
00607
00608 double gammaAB = gammaR*pow(pAB/pR,power)*(mR/mAB)*fR*fR;
00609
00610 switch (_spin) {
00611 case 0:
00612 ampl = _ampl*TComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
00613 fR*fD/(mR*mR-mAB*mAB-TComplex(0.0,mR*gammaAB));
00614 break;
00615 case 1:
00616 ampl = _ampl*TComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
00617 fR*fD*(mAC*mAC-mBC*mBC+(mD*mD-mC*mC)*(mB*mB-mA*mA)/(mAB*mAB))/
00618 (mR*mR-mAB*mAB-TComplex(0.0,mR*gammaAB));
00619 break;
00620 case 2:
00621 ampl = _ampl*TComplex(cos(_theta*pi180inv),sin(_theta*pi180inv))*
00622 fR*fD/(mR*mR-mAB*mAB-TComplex(0.0,mR*gammaAB))*
00623 (pow((mBC*mBC-mAC*mAC+(mD*mD-mC*mC)*(mA*mA-mB*mB)/(mAB*mAB)),2)-
00624 (1.0/3.0)*(mAB*mAB-2*mD*mD-2*mC*mC+pow((mD*mD- mC*mC)/mAB, 2))*
00625 (mAB*mAB-2*mA*mA-2*mB*mB+pow((mA*mA-mB*mB)/mAB,2)));
00626 break;
00627 }
00628
00629 return ampl;
00630
00631 }
00632
00633 TComplex Dalitz::Babar_Amplitude(double x, double y, double z) {
00634
00635
00636
00637
00638 double m_mass[4] = { 1.86450, 0.497648, 0.139570, 0.139570};
00639
00640 TComplex D0(0.0,0.0);
00641
00642
00643 x = sqrt(x);
00644 y = sqrt(y);
00645 z = sqrt(z);
00646
00647 TComplex DK2piRes[17];
00648 DK2piRes[0] = Babar_sakurai(x, y, z, 1.00, 0.0, 0.1464, 0.7758);
00649 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);
00650 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);
00651 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);
00652 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);
00653 DK2piRes[5] = Babar_sakurai(x, y, z, 0.52, 38, 0.455, 1.406);
00654 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);
00655 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);
00656 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);
00657 DK2piRes[9] = Babar_resAmp(x, z, y, m_mass[3], m_mass[1], m_mass[2], 0.52, 154, 0.232, 1.414, 1);
00658 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);
00659 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);
00660 DK2piRes[12] = Babar_resAmp(x, z, y, m_mass[3], m_mass[1], m_mass[2], 0.89, -139, 0.322, 1.717, 1);
00661 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);
00662 DK2piRes[14] = Babar_resAmp(y, z, x, m_mass[2], m_mass[1], m_mass[3], 0.37, 18, 0.294, 1.412, 0);
00663 DK2piRes[15] = Babar_resAmp(y, z, x, m_mass[2], m_mass[1], m_mass[3], 0.075, -104, 0.0985, 1.4256, 2);
00664
00665 double pi180inv = 3.1415926/180.;
00666 DK2piRes[16] = TComplex(3.53*cos(128*pi180inv),3.53*sin(128*pi180inv));
00667
00668 for(int i=0; i<17; i++){
00669 D0 += DK2piRes[i];
00670 }
00671
00672 return D0;
00673 }
00674
00675 int Dalitz::getBin(double mx, double my, double mz) {
00676
00677
00678 double m_mass_2[4]={ 1.86450*1.86450, 0.497648*0.497648,
00679 0.139570*0.139570, 0.139570*0.139570};
00680 double m_sum_m_2 = m_mass_2[0] + m_mass_2[1] + m_mass_2[2] + m_mass_2[3];
00681
00682
00683 if( !(Point_on_DP2(mx, my)) && !(Point_on_DP2(my, mx)) ) return -1;
00684
00685
00686 mz = m_sum_m_2 - mx - my;
00687 if (mz < 0) mz = 0;
00688
00689
00690 double thisPhase = Phase(mx, my, mz);
00691
00692
00693 int thisbin = -1;
00694
00695
00696 for (int bin = 0; bin < N; bin++) {
00697 if((thisPhase >= (bin-0.5)*2*PI/N) && (thisPhase < (bin+0.5)*2*PI/N)) thisbin = bin;
00698 }
00699
00700 return thisbin;
00701
00702 }
00703
00704