BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDsToKSKSpi.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2// Environment:
3// This software is part of the EvtGen package developed jointly
4// for the BaBar and CLEO collaborations. If you use all or part
5// of it, please give an appropriate acknowledgement.
6//
7// Copyright Information: See EvtGen/COPYRIGHT
8// Copyright (C) 1998 Caltech, UCSB
9//
10// Module: EvtDsToKSKSpi.cc
11//
12// Description: Routine to handle three-body decays of Ds+/Ds-
13//
14// Modification history:
15// Liaoyuan Dong Sun Jun 02 01:23:25 2024 Module created
16//------------------------------------------------------------------------
17#include "EvtDsToKSKSpi.hh"
27#include <fstream>
28#include <stdlib.h>
29#include <string>
30using std::endl;
31
33
34void EvtDsToKSKSpi::getName( std::string& model_name ) { model_name = "DsToKSKSpi"; }
35
37
39 // check that there are 0 arguments
40 checkNArg( 0 );
41 checkNDaug( 3 );
46
47 phi[0] = 0;
48 rho[0] = 1;
49 phi[1] = 0;
50 rho[1] = 0;
51 phi[2] = 0;
52 rho[2] = 0;
53 phi[3] = 0;
54 rho[3] = 0;
55 phi[4] = 0;
56 rho[4] = 0;
57 phi[0] = 0.0;
58 phi[1] = 2.3351;
59
60 rho[0] = 1.0;
61 rho[1] = 3.3891;
62
63 modetype[0] = 13;
64 modetype[1] = 12;
65 modetype[2] = 12;
66 modetype[3] = 13;
67 modetype[4] = 13;
68
69 // cout << "DsToKSKSpi :" << endl;
70 // for (int i=0; i<5; i++) {
71 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
72 // }
73
74 width[0] = 5.0800e-02;
75 width[1] = 1.4041e-01; // 1.2300e-01
76 width[2] = 2.6500e-01;
77 width[3] = 2.7000e-01;
78 width[4] = 5.0800e-02;
79
80 mass[0] = 8.9166e-01;
81 mass[1] = 1.7228e+00;
82 mass[2] = 1.3500e+00;
83 mass[3] = 1.4250e+00;
84 mass[4] = 8.9166e-01;
85
86 mD = 1.86486;
87 mDs = 1.9683;
88 rRes = 9.0;
89 rD = 5.0;
90 metap = 0.95778;
91 mkstr = 0.89594;
92 mk0 = 0.497614;
93 mass_Kaon = 0.49368;
94 mass_Pion = 0.13957;
95 mass_Pi0 = 0.1349766;
96 math_pi = 3.1415926;
97 // mrho = 0.77549;
98 // Grho = 0.1491;
99 ma0 = 0.99;
100 Ga0 = 0.0756;
101 meta = 0.547862;
102
103 GS1 = 0.636619783;
104 GS2 = 0.01860182466;
105 GS3 = 0.1591549458; // 1/(2*math_2pi)
106 GS4 = 0.00620060822; // mass_Pion2/math_pi
107
108 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
109 for ( int i = 0; i < 4; i++ )
110 {
111 for ( int j = 0; j < 4; j++ ) { G[i][j] = GG[i][j]; }
112 }
113}
114
116
118 /*
119 double maxprob = 0.0;
120 for(int ir=0;ir<=60000000;ir++){
121 p->initializePhaseSpace(getNDaug(),getDaugs());
122 EvtVector4R D1 = p->getDaug(0)->getP4();
123 EvtVector4R D2 = p->getDaug(1)->getP4();
124 EvtVector4R D3 = p->getDaug(2)->getP4();
125
126 double P1[4], P2[4], P3[4];
127 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
128 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
129 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
130
131 double value;
132 // value = calDalEva(P1, P2, P3);
133 int g0[5]={1,1,1,1,12};
134 int spin[5]={1,0,0,0,0};
135 int nstates=5;
136 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value);
137 if (value<0) continue;
138 if(value>maxprob) {
139 maxprob=value;
140 cout << "ir= " << ir << endl;
141 cout << "double P1[4] = {" << P1[0] <<","<< P1[1] <<","<< P1[2] <<","<< P1[3]
142 <<"};"<< endl; cout << "double P2[4] = {" << P2[0] <<","<< P2[1] <<","<< P2[2] <<","<< P2[3]
143 <<"};"<< endl; cout << "double P3[4] = {" << P3[0] <<","<< P3[1] <<","<< P3[2] <<","<< P3[3]
144 <<"};"<< endl; cout << "MAX====> " << maxprob << endl;
145 }
146 }
147 printf("MAXprob = %.10f\n",maxprob);
148 */
149
151 EvtVector4R D1 = p->getDaug( 0 )->getP4();
152 EvtVector4R D2 = p->getDaug( 1 )->getP4();
153 EvtVector4R D3 = p->getDaug( 2 )->getP4();
154
155 double P1[4], P2[4], P3[4];
156 P1[0] = D1.get( 0 );
157 P1[1] = D1.get( 1 );
158 P1[2] = D1.get( 2 );
159 P1[3] = D1.get( 3 );
160 P2[0] = D2.get( 0 );
161 P2[1] = D2.get( 1 );
162 P2[2] = D2.get( 2 );
163 P2[3] = D2.get( 3 );
164 P3[0] = D3.get( 0 );
165 P3[1] = D3.get( 1 );
166 P3[2] = D3.get( 2 );
167 P3[3] = D3.get( 3 );
168 // P1[0] = 0.770687885731853; P1[1] = 0.085867043334332; P1[2] =-0.550183164280213; P1[3]
169 // =-0.190434925445587; P2[0] = 0.840768850984605; P2[1] =-0.240025188001374; P2[2] =
170 // 0.631078772181167; P2[3] = 0.058310035302289; P3[0] = 0.393776666168180; P3[1]
171 // =-0.159532664043002; P3[2] =-0.300922058149019; P3[3] = 0.139912371490140;
172
173 double value;
174 int g0[5] = { 1, 1, 1, 1, 12 };
175 int spin[5] = { 1, 0, 0, 0, 0 };
176 int nstates = 5;
177 calEva( P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value );
178
179 setProb( value );
180 return;
181}
182
183void EvtDsToKSKSpi::Com_Multi( double a1[2], double a2[2], double res[2] ) {
184 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
185 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
186}
187void EvtDsToKSKSpi::Com_Divide( double a1[2], double a2[2], double res[2] ) {
188 double tmp = a2[0] * a2[0] + a2[1] * a2[1];
189 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / tmp;
190 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / tmp;
191}
192//------------base---------------------------------
193double EvtDsToKSKSpi::SCADot( double a1[4], double a2[4] ) {
194 double _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
195 return _cal;
196}
197double EvtDsToKSKSpi::barrier( int l, double sa, double sb, double sc, double r,
198 double mass ) {
199 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
200 if ( q < 0 ) q = 1e-16;
201 double z;
202 z = q * r * r;
203 double sa0;
204 sa0 = mass * mass;
205 double q0 = ( sa0 + sb - sc ) * ( sa0 + sb - sc ) / ( 4 * sa0 ) - sb;
206 if ( q0 < 0 ) q0 = 1e-16;
207 double z0 = q0 * r * r;
208 double F = 0.0;
209 if ( l == 0 ) F = 1;
210 if ( l == 1 ) F = sqrt( ( 1 + z0 ) / ( 1 + z ) );
211 if ( l == 2 ) F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) );
212 return F;
213}
214
215void EvtDsToKSKSpi::calt1( double daug1[4], double daug2[4], double t1[4] ) {
216 double p, pq, tmp;
217 double pa[4], qa[4];
218 for ( int i = 0; i < 4; i++ )
219 {
220 pa[i] = daug1[i] + daug2[i];
221 qa[i] = daug1[i] - daug2[i];
222 }
223 p = SCADot( pa, pa );
224 pq = SCADot( pa, qa );
225 tmp = pq / p;
226 for ( int i = 0; i < 4; i++ ) { t1[i] = qa[i] - tmp * pa[i]; }
227}
228void EvtDsToKSKSpi::calt2( double daug1[4], double daug2[4], double t2[4][4] ) {
229 double p, r;
230 double pa[4], t1[4];
231 calt1( daug1, daug2, t1 );
232 r = SCADot( t1, t1 ) / 3.0;
233 for ( int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
234 p = SCADot( pa, pa );
235 for ( int i = 0; i < 4; i++ )
236 {
237 for ( int j = 0; j < 4; j++ )
238 { t2[i][j] = t1[i] * t1[j] - r * ( G[i][j] - pa[i] * pa[j] / p ); }
239 }
240}
241//-------------------prop--------------------------------------------
242
243double EvtDsToKSKSpi::wid( double mass2, double mass, double sa, double sb, double sc,
244 double r2, int l ) {
245 double widm = 0.;
246 double m = sqrt( sa );
247 double tmp = sb - sc;
248 double tmp1 = sa + tmp;
249 double q = 0.25 * tmp1 * tmp1 / sa - sb;
250 if ( q < 0 ) q = 1e-16;
251 double tmp2 = mass2 + tmp;
252 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
253 if ( q0 < 0 ) q0 = 1e-16;
254 double z = q * r2;
255 double z0 = q0 * r2;
256 double t = q / q0;
257 if ( l == 0 ) { widm = sqrt( t ) * mass / m; }
258 else if ( l == 1 ) { widm = t * sqrt( t ) * mass / m * ( 1 + z0 ) / ( 1 + z ); }
259 else if ( l == 2 )
260 { widm = t * t * sqrt( t ) * mass / m * ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ); }
261 return widm;
262}
263double EvtDsToKSKSpi::widl1( double mass2, double mass, double sa, double sb, double sc,
264 double r2 ) {
265 double widm = 0.;
266 double m = sqrt( sa );
267 double tmp = sb - sc;
268 double tmp1 = sa + tmp;
269 double q = 0.25 * tmp1 * tmp1 / sa - sb;
270 if ( q < 0 ) q = 1e-16;
271 double tmp2 = mass2 + tmp;
272 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
273 if ( q0 < 0 ) q0 = 1e-16;
274 double z = q * r2;
275 double z0 = q0 * r2;
276 double F = ( 1 + z0 ) / ( 1 + z );
277 double t = q / q0;
278 widm = t * sqrt( t ) * mass / m * F;
279 return widm;
280}
281void EvtDsToKSKSpi::propagatorRBW( double mass2, double mass, double width, double sa,
282 double sb, double sc, double r2, int l, double prop[2] ) {
283 double a[2], b[2];
284 a[0] = 1;
285 a[1] = 0;
286 b[0] = mass2 - sa;
287 b[1] = -mass * width * wid( mass2, mass, sa, sb, sc, r2, l );
288 Com_Divide( a, b, prop );
289}
290
291void EvtDsToKSKSpi::propagatorGS( double mass2, double mass, double width, double sa,
292 double sb, double sc, double r2, double prop[2] ) {
293 double a[2], b[2];
294 double tmp = sb - sc;
295 double tmp1 = sa + tmp;
296 double q2 = 0.25 * tmp1 * tmp1 / sa - sb;
297 if ( q2 < 0 ) q2 = 1e-16;
298
299 double tmp2 = mass2 + tmp;
300 double q02 = 0.25 * tmp2 * tmp2 / mass2 - sb;
301 if ( q02 < 0 ) q02 = 1e-16;
302
303 double q = sqrt( q2 );
304 double q0 = sqrt( q02 );
305 double m = sqrt( sa );
306 double q03 = q0 * q02;
307 double tmp3 = log( mass + 2 * q0 ) + 1.2926305904; // log(mpi0+mpip) = -1.2926305904;
308
309 double h = GS1 * q / m * ( log( m + 2 * q ) + 1.2926305904 );
310 double h0 = GS1 * q0 / mass * tmp3;
311 double dh = h0 * ( 0.125 / q02 - 0.5 / mass2 ) + GS3 / mass2;
312 double d = GS2 / q02 * tmp3 + GS3 * mass / q0 - GS4 * mass / q03;
313 double f = mass2 / q03 * ( q2 * ( h - h0 ) + ( mass2 - sa ) * q02 * dh );
314
315 a[0] = 1.0 + d * width / mass;
316 a[1] = 0.0;
317 b[0] = mass2 - sa + width * f;
318 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r2 );
319 Com_Divide( a, b, prop );
320}
321
322void EvtDsToKSKSpi::KPiSLASS( double sa, double sb, double sc, double prop[2] ) {
323 const double m1430 = 1.441;
324 const double sa0 = 1.441 * 1.441; // m1430*m1430;
325 const double w1430 = 0.193;
326 const double Lass1 = 0.25 / sa0;
327 double tmp = sb - sc;
328 double tmp1 = sa0 + tmp;
329 double q0 = Lass1 * tmp1 * tmp1 - sb;
330 if ( q0 < 0 ) q0 = 1e-16;
331 double tmp2 = sa + tmp;
332 double qs = 0.25 * tmp2 * tmp2 / sa - sb;
333 double q = sqrt( qs );
334 double width = w1430 * q * m1430 / sqrt( sa * q0 );
335 double temp_R = atan( m1430 * width / ( sa0 - sa ) );
336 if ( temp_R < 0 ) temp_R += math_pi;
337 double deltaR = -109.7 + temp_R;
338 double temp_F =
339 atan( 0.226 * q / ( 2.0 - 3.819 * qs ) ); // 2.0*0.113 = 0.226; 0.113*33.8 = 3.819
340 if ( temp_F < 0 ) temp_F += math_pi;
341 double deltaF = 0.1 + temp_F;
342 double deltaS = deltaR + 2.0 * deltaF;
343 double t1 = 0.96 * sin( deltaF );
344 double t2 = sin( deltaR );
345 double CF[2], CS[2];
346 CF[0] = cos( deltaF );
347 CF[1] = sin( deltaF );
348 CS[0] = cos( deltaS );
349 CS[1] = sin( deltaS );
350 prop[0] = t1 * CF[0] + t2 * CS[0];
351 prop[1] = t1 * CF[1] + t2 * CS[1];
352}
353
354double EvtDsToKSKSpi::DDalitz( double P1[4], double P2[4], double P3[4], int Ang,
355 double mass ) {
356 double pR[4], pD[4];
357 double temp_PDF, v_re;
358 temp_PDF = 0.0;
359 v_re = 0.0;
360 double B[2], s1, s2, s3, sR, sD;
361 for ( int i = 0; i < 4; i++ )
362 {
363 pR[i] = P1[i] + P2[i];
364 pD[i] = pR[i] + P3[i];
365 }
366 s1 = SCADot( P1, P1 );
367 s2 = SCADot( P2, P2 );
368 s3 = SCADot( P3, P3 );
369 sR = SCADot( pR, pR );
370 sD = SCADot( pD, pD );
371 // int G[4][4];
372 // for(int i=0; i!=4; i++){
373 // for(int j=0; j!=4; j++){
374 // if(i==j){
375 // if(i==0) G[i][j] = 1;
376 // else G[i][j] = -1;
377 // }
378 // else G[i][j] = 0;
379 // }
380 // }
381 if ( Ang == 0 )
382 {
383 B[0] = 1;
384 B[1] = 1;
385 temp_PDF = 1;
386 }
387 if ( Ang == 1 )
388 {
389 B[0] = barrier( 1, sR, s1, s2, 3.0, mass );
390 B[1] = barrier( 1, sD, sR, s3, 5.0, 1.9683 );
391 double t1[4], T1[4];
392 calt1( P1, P2, t1 );
393 calt1( pR, P3, T1 );
394 temp_PDF = 0;
395 for ( int i = 0; i < 4; i++ ) { temp_PDF += t1[i] * T1[i] * G[i][i]; }
396 }
397 if ( Ang == 2 )
398 {
399 B[0] = barrier( 2, sR, s1, s2, 3.0, mass );
400 B[1] = barrier( 2, sD, sR, s3, 5.0, 1.9683 );
401 double t2[4][4], T2[4][4];
402 calt2( P1, P2, t2 );
403 calt2( pR, P3, T2 );
404 temp_PDF = 0;
405 for ( int i = 0; i < 4; i++ )
406 {
407 for ( int j = 0; j < 4; j++ ) { temp_PDF += t2[i][j] * T2[j][i] * G[i][i] * G[j][j]; }
408 }
409 }
410 v_re = temp_PDF * B[0] * B[1];
411 return v_re;
412}
413
414void EvtDsToKSKSpi::calEva( double* Ks01, double* Ks02, double* Pic, double* mass1,
415 double* width1, double* amp, double* phase, int* g0, int* spin,
416 int* modetype, int nstates, double& Result ) {
417 // double Ks0[4], Pic[4], Pi0[4];
418 // double rRes = 9.0;
419 double P12[4], P23[4], P13[4];
420 double cof[2], amp_PDF[2], PDF[2];
421 double scpi, sks02, sks01;
422 double s12, s13, s23;
423 for ( int i = 0; i < 4; i++ )
424 {
425 P12[i] = Ks01[i] + Ks02[i];
426 P13[i] = Ks01[i] + Pic[i];
427 P23[i] = Ks02[i] + Pic[i];
428 }
429 scpi = SCADot( Pic, Pic );
430 sks01 = SCADot( Ks01, Ks01 );
431 sks02 = SCADot( Ks02, Ks02 );
432 s12 = SCADot( P12, P12 );
433 s13 = SCADot( P13, P13 );
434 s23 = SCADot( P23, P23 );
435 double pro[2], temp_PDF, amp_tmp[2], temp_PDF1, temp_PDF2, pro1[2], pro2[2];
436 double mass1sq;
437 amp_PDF[0] = 0;
438 amp_PDF[1] = 0;
439 PDF[0] = 0;
440 PDF[1] = 0;
441 amp_tmp[0] = 0;
442 amp_tmp[1] = 0;
443 for ( int i = 0; i < nstates; i++ )
444 {
445 amp_tmp[0] = 0;
446 amp_tmp[1] = 0;
447 mass1sq = mass1[i] * mass1[i];
448 cof[0] = amp[i] * cos( phase[i] );
449 cof[1] = amp[i] * sin( phase[i] );
450 temp_PDF = 0;
451 /* if(modetype[i] == 23){
452 temp_PDF = DDalitz( Ks02, Pi0, Ks0, 1, mass1[i]);
453 if(g0[i]==1) propagatorGS(mass1sq,mass1[i],width1[i],s23,scpi,snpi,9.0,pro);
454 if(g0[i]==0){
455 pro[0] = 1;
456 pro[1] = 0;
457 }
458 amp_tmp[0] = temp_PDF*pro[0];
459 amp_tmp[1] = temp_PDF*pro[1];
460
461 }
462 */
463 if ( modetype[i] == 12 )
464 {
465 temp_PDF = DDalitz( Ks01, Ks02, Pic, spin[i], mass1[i] );
466 if ( g0[i] == 1 )
467 propagatorRBW( mass1sq, mass1[i], width1[i], s12, sks01, sks02, rRes, spin[i], pro );
468 if ( g0[i] == 12 ) KPiSLASS( s12, sks01, sks02, pro );
469 if ( g0[i] == 0 )
470 {
471 pro[0] = 1;
472 pro[1] = 0;
473 }
474 amp_tmp[0] = temp_PDF * pro[0];
475 amp_tmp[1] = temp_PDF * pro[1];
476 }
477 if ( modetype[i] == 13 )
478 {
479 temp_PDF1 = DDalitz( Ks01, Pic, Ks02, spin[i], mass1[i] );
480 if ( g0[i] == 1 )
481 propagatorRBW( mass1sq, mass1[i], width1[i], s13, sks01, scpi, rRes, spin[i], pro1 );
482 if ( g0[i] == 12 ) KPiSLASS( s13, sks01, scpi, pro1 );
483 if ( g0[i] == 0 )
484 {
485 pro1[0] = 1;
486 pro1[1] = 0;
487 }
488 temp_PDF2 = DDalitz( Ks02, Pic, Ks01, spin[i], mass1[i] );
489 if ( g0[i] == 1 )
490 propagatorRBW( mass1sq, mass1[i], width1[i], s23, sks02, scpi, rRes, spin[i], pro2 );
491 if ( g0[i] == 12 ) KPiSLASS( s23, sks02, scpi, pro2 );
492 if ( g0[i] == 0 )
493 {
494 pro2[0] = 1;
495 pro2[1] = 0;
496 }
497 amp_tmp[0] = ( temp_PDF1 * pro1[0] + temp_PDF2 * pro2[0] );
498 amp_tmp[1] = ( temp_PDF1 * pro1[1] + temp_PDF2 * pro2[1] );
499 }
500 Com_Multi( amp_tmp, cof, amp_PDF );
501 PDF[0] += amp_PDF[0];
502 PDF[1] += amp_PDF[1];
503 }
504 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
505 // cout<<"value = "<<value<<endl;
506 if ( value <= 0 ) value = 1e-20;
507 Result = value;
508 // cout<<"value = "<<value<<" Result = "<<Result<<endl;
509}
double mass
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
****INTEGER imax DOUBLE PRECISION m_pi *DOUBLE PRECISION m_amfin DOUBLE PRECISION m_Chfin DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_sinw2 DOUBLE PRECISION m_GFermi DOUBLE PRECISION m_MfinMin DOUBLE PRECISION m_ta2 INTEGER m_out INTEGER m_KeyFSR INTEGER m_KeyQCD *COMMON c_Semalib $ !copy of input $ !CMS energy $ !beam mass $ !final mass $ !beam charge $ !final charge $ !smallest final mass $ !Z mass $ !Z width $ !EW mixing angle $ !Gmu Fermi $ alphaQED at q
Definition KKsem.h:33
TCrossPart * CS
Definition Mcgpj.cxx:51
void checkSpinDaughter(int d1, EvtSpinType::spintype sp)
void checkSpinParent(EvtSpinType::spintype sp)
void setProbMax(double prbmx)
void checkNDaug(int d1, int d2=-1)
EvtId * getDaugs()
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
void setProb(double prob)
void decay(EvtParticle *p)
virtual ~EvtDsToKSKSpi()
EvtDecayBase * clone()
void getName(std::string &name)
const EvtVector4R & getP4() const
EvtParticle * getDaug(int i)
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
double get(int i) const
int t()
Definition t.c:1