BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDsToKSpipi0.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: EvtDsToKSpipi0.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 "EvtDsToKSpipi0.hh"
27#include <fstream>
28#include <stdlib.h>
29#include <string>
30using std::endl;
31
33
34void EvtDsToKSpipi0::getName( std::string& model_name ) { model_name = "DsToKSpipi0"; }
35
37
39 // check that there are 0 arguments
40 checkNArg( 0 );
41 checkNDaug( 3 );
46
47 phi[0] = 1.9740e-01;
48 phi[1] = 3.1956e+00;
49 phi[3] = 1.9523e-01;
50 phi[4] = 2.1602e+00;
51
52 rho[0] = 3.0648e-01;
53 rho[1] = 4.0708e-01;
54 rho[3] = 7.5484e-01;
55 rho[4] = 2.6511e+00;
56
57 phi[2] = 0;
58 rho[2] = 1;
59 modetype[0] = 12;
60 modetype[1] = 13;
61 modetype[2] = 23;
62 modetype[3] = 131;
63 modetype[4] = 231;
64
65 // cout << "DsToKSpipi0 :" << endl;
66 // for (int i=0; i<5; i++) {
67 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
68 // }
69
70 width[0] = 5.0800e-02;
71 width[1] = 4.7300e-02;
72 width[2] = 1.5020e-01;
73 width[3] = 2.3200e-01;
74 width[4] = 4.0000e-01;
75 mass[0] = 8.9166e-01;
76 mass[1] = 8.9555e-01;
77 mass[2] = 7.6650e-01;
78 mass[3] = 1.4140e+00;
79 mass[4] = 1.4650e+00;
80
81 mD = 1.86486;
82 mDs = 1.9683;
83 rRes = 9.0;
84 rD = 5.0;
85 metap = 0.95778;
86 mkstr = 0.89594;
87 mk0 = 0.497614;
88 mass_Kaon = 0.49368;
89 mass_Pion = 0.13957;
90 mass_Pi0 = 0.1349766;
91 math_pi = 3.1415926;
92 // mrho = 0.77549;
93 // Grho = 0.1491;
94 ma0 = 0.99;
95 Ga0 = 0.0756;
96 meta = 0.547862;
97
98 GS1 = 0.636619783;
99 GS2 = 0.01860182466;
100 GS3 = 0.1591549458; // 1/(2*math_2pi)
101 GS4 = 0.00620060822; // mass_Pion2/math_pi
102
103 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
104 for ( int i = 0; i < 4; i++ )
105 {
106 for ( int j = 0; j < 4; j++ ) { G[i][j] = GG[i][j]; }
107 }
108}
109
111
113 /*
114 double maxprob = 0.0;
115 for(int ir=0;ir<=6000000;ir++){
116 p->initializePhaseSpace(getNDaug(),getDaugs());
117 EvtVector4R D1 = p->getDaug(0)->getP4();
118 EvtVector4R D2 = p->getDaug(1)->getP4();
119 EvtVector4R D3 = p->getDaug(2)->getP4();
120
121 double P1[4], P2[4], P3[4];
122 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
123 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
124 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
125
126 double value;
127 //value = calDalEva(P1, P2, P3);
128 int g0[5]={1,1,1,1,1};
129 int nstates=5;
130 calEva(P1, P2, P3, mass, width, rho, phi, g0, modetype, nstates, value);
131 if (value<0) continue;
132 if(value>maxprob) {
133 maxprob=value;
134 cout << "ir= " << ir << endl;
135 cout << "double P1[4] = {" << P1[0] <<","<< P1[1] <<","<< P1[2] <<","<< P1[3]
136 <<"};"<< endl; cout << "double P2[4] = {" << P2[0] <<","<< P2[1] <<","<< P2[2] <<","<<
137 P2[3] <<"};"<< endl; cout << "double P3[4] = {" << P3[0] <<","<< P3[1] <<","<< P3[2]
138 <<","<< P3[3] <<"};"<< endl; cout << "MAX====> " << maxprob << endl;
139 }
140 }
141 printf("MAXprob = %.10f\n",maxprob);
142 */
144 EvtVector4R D1 = p->getDaug( 0 )->getP4();
145 EvtVector4R D2 = p->getDaug( 1 )->getP4();
146 EvtVector4R D3 = p->getDaug( 2 )->getP4();
147
148 double P1[4], P2[4], P3[4];
149 P1[0] = D1.get( 0 );
150 P1[1] = D1.get( 1 );
151 P1[2] = D1.get( 2 );
152 P1[3] = D1.get( 3 );
153 P2[0] = D2.get( 0 );
154 P2[1] = D2.get( 1 );
155 P2[2] = D2.get( 2 );
156 P2[3] = D2.get( 3 );
157 P3[0] = D3.get( 0 );
158 P3[1] = D3.get( 1 );
159 P3[2] = D3.get( 2 );
160 P3[3] = D3.get( 3 );
161 // P1[0] = 0.840899; P1[1] = -0.624276; P1[2] = 0.0113601 ; P1[3] = -0.263898;
162 // P2[0] = 0.89303 ; P2[1] = 0.860432 ; P2[2] = -0.0491258; P2[3] = 0.187791 ;
163 // P3[0] = 0.288373; P3[1] = 0.102614 ; P3[2] = -0.132757 ; P3[3] = -0.191797;
164 // P1[0] = 0.817635; P1[1] = 0.140538 ; P1[2] = -0.501451 ; P1[3] = 0.386916;
165 // P2[0] = 0.952772; P2[1] = -0.164752; P2[2] = 0.437717 ; P2[3] = -0.818264;
166 // P3[0] = 0.260117; P3[1] = 0.213518 ; P3[2] = -0.0606403; P3[3] = -0.0132164;
167
168 double value;
169 int g0[5] = { 1, 1, 1, 1, 1 };
170 int nstates = 5;
171 calEva( P1, P2, P3, mass, width, rho, phi, g0, modetype, nstates, value );
172
173 // value = calEva(P1, P2, P3);
174 setProb( value );
175 return;
176}
177
178void EvtDsToKSpipi0::Com_Multi( double a1[2], double a2[2], double res[2] ) {
179 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
180 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
181}
182void EvtDsToKSpipi0::Com_Divide( double a1[2], double a2[2], double res[2] ) {
183 double tmp = a2[0] * a2[0] + a2[1] * a2[1];
184 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / tmp;
185 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / tmp;
186}
187//------------base---------------------------------
188double EvtDsToKSpipi0::SCADot( double a1[4], double a2[4] ) {
189 double _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
190 return _cal;
191}
192double EvtDsToKSpipi0::barrier( int l, double sa, double sb, double sc, double r,
193 double mass ) {
194 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
195 if ( q < 0 ) q = 1e-16;
196 double z;
197 z = q * r * r;
198 double sa0;
199 sa0 = mass * mass;
200 double q0 = ( sa0 + sb - sc ) * ( sa0 + sb - sc ) / ( 4 * sa0 ) - sb;
201 if ( q0 < 0 ) q0 = 1e-16;
202 double z0 = q0 * r * r;
203 double F = 0.0;
204 if ( l == 0 ) F = 1;
205 if ( l == 1 ) F = sqrt( ( 1 + z0 ) / ( 1 + z ) );
206 if ( l == 2 ) F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) );
207 return F;
208}
209
210void EvtDsToKSpipi0::calt1( double daug1[4], double daug2[4], double t1[4] ) {
211 double p, pq, tmp;
212 double pa[4], qa[4];
213 for ( int i = 0; i < 4; i++ )
214 {
215 pa[i] = daug1[i] + daug2[i];
216 qa[i] = daug1[i] - daug2[i];
217 }
218 p = SCADot( pa, pa );
219 pq = SCADot( pa, qa );
220 tmp = pq / p;
221 for ( int i = 0; i < 4; i++ ) { t1[i] = qa[i] - tmp * pa[i]; }
222}
223void EvtDsToKSpipi0::calt2( double daug1[4], double daug2[4], double t2[4][4] ) {
224 double p, r;
225 double pa[4], t1[4];
226 calt1( daug1, daug2, t1 );
227 r = SCADot( t1, t1 ) / 3.0;
228 for ( int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
229 p = SCADot( pa, pa );
230 for ( int i = 0; i < 4; i++ )
231 {
232 for ( int j = 0; j < 4; j++ )
233 { t2[i][j] = t1[i] * t1[j] - r * ( G[i][j] - pa[i] * pa[j] / p ); }
234 }
235}
236//-------------------prop--------------------------------------------
237
238double EvtDsToKSpipi0::wid( double mass2, double mass, double sa, double sb, double sc,
239 double r2, int l ) {
240 double widm = 0.;
241 double m = sqrt( sa );
242 double tmp = sb - sc;
243 double tmp1 = sa + tmp;
244 double q = 0.25 * tmp1 * tmp1 / sa - sb;
245 if ( q < 0 ) q = 1e-16;
246 double tmp2 = mass2 + tmp;
247 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
248 if ( q0 < 0 ) q0 = 1e-16;
249 double z = q * r2;
250 double z0 = q0 * r2;
251 double t = q / q0;
252 if ( l == 0 ) { widm = sqrt( t ) * mass / m; }
253 else if ( l == 1 ) { widm = t * sqrt( t ) * mass / m * ( 1 + z0 ) / ( 1 + z ); }
254 else if ( l == 2 )
255 { widm = t * t * sqrt( t ) * mass / m * ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ); }
256 return widm;
257}
258double EvtDsToKSpipi0::widl1( double mass2, double mass, double sa, double sb, double sc,
259 double r2 ) {
260 double widm = 0.;
261 double m = sqrt( sa );
262 double tmp = sb - sc;
263 double tmp1 = sa + tmp;
264 double q = 0.25 * tmp1 * tmp1 / sa - sb;
265 if ( q < 0 ) q = 1e-16;
266 double tmp2 = mass2 + tmp;
267 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
268 if ( q0 < 0 ) q0 = 1e-16;
269 double z = q * r2;
270 double z0 = q0 * r2;
271 double F = ( 1 + z0 ) / ( 1 + z );
272 double t = q / q0;
273 widm = t * sqrt( t ) * mass / m * F;
274 return widm;
275}
276void EvtDsToKSpipi0::propagatorRBW( double mass2, double mass, double width, double sa,
277 double sb, double sc, double r2, int l, double prop[2] ) {
278 double a[2], b[2];
279 a[0] = 1;
280 a[1] = 0;
281 b[0] = mass2 - sa;
282 b[1] = -mass * width * wid( mass2, mass, sa, sb, sc, r2, l );
283 Com_Divide( a, b, prop );
284}
285
286void EvtDsToKSpipi0::propagatorGS( double mass2, double mass, double width, double sa,
287 double sb, double sc, double r2, double prop[2] ) {
288 double a[2], b[2];
289 double tmp = sb - sc;
290 double tmp1 = sa + tmp;
291 double q2 = 0.25 * tmp1 * tmp1 / sa - sb;
292 if ( q2 < 0 ) q2 = 1e-16;
293
294 double tmp2 = mass2 + tmp;
295 double q02 = 0.25 * tmp2 * tmp2 / mass2 - sb;
296 if ( q02 < 0 ) q02 = 1e-16;
297
298 double q = sqrt( q2 );
299 double q0 = sqrt( q02 );
300 double m = sqrt( sa );
301 double q03 = q0 * q02;
302 double tmp3 = log( mass + 2 * q0 ) + 1.2926305904; // log(mpi0+mpip) = -1.2926305904;
303
304 double h = GS1 * q / m * ( log( m + 2 * q ) + 1.2926305904 );
305 double h0 = GS1 * q0 / mass * tmp3;
306 double dh = h0 * ( 0.125 / q02 - 0.5 / mass2 ) + GS3 / mass2;
307 double d = GS2 / q02 * tmp3 + GS3 * mass / q0 - GS4 * mass / q03;
308 double f = mass2 / q03 * ( q2 * ( h - h0 ) + ( mass2 - sa ) * q02 * dh );
309
310 a[0] = 1.0 + d * width / mass;
311 a[1] = 0.0;
312 b[0] = mass2 - sa + width * f;
313 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r2 );
314 Com_Divide( a, b, prop );
315}
316
317void EvtDsToKSpipi0::KPiSLASS( double sa, double sb, double sc, double prop[2] ) {
318 const double m1430 = 1.441;
319 const double sa0 = 1.441 * 1.441; // m1430*m1430;
320 const double w1430 = 0.193;
321 const double Lass1 = 0.25 / sa0;
322 double tmp = sb - sc;
323 double tmp1 = sa0 + tmp;
324 double q0 = Lass1 * tmp1 * tmp1 - sb;
325 if ( q0 < 0 ) q0 = 1e-16;
326 double tmp2 = sa + tmp;
327 double qs = 0.25 * tmp2 * tmp2 / sa - sb;
328 double q = sqrt( qs );
329 double width = w1430 * q * m1430 / sqrt( sa * q0 );
330 double temp_R = atan( m1430 * width / ( sa0 - sa ) );
331 if ( temp_R < 0 ) temp_R += math_pi;
332 double deltaR = -109.7 + temp_R;
333 double temp_F =
334 atan( 0.226 * q / ( 2.0 - 3.819 * qs ) ); // 2.0*0.113 = 0.226; 0.113*33.8 = 3.819
335 if ( temp_F < 0 ) temp_F += math_pi;
336 double deltaF = 0.1 + temp_F;
337 double deltaS = deltaR + 2.0 * deltaF;
338 double t1 = 0.96 * sin( deltaF );
339 double t2 = sin( deltaR );
340 double CF[2], CS[2];
341 CF[0] = cos( deltaF );
342 CF[1] = sin( deltaF );
343 CS[0] = cos( deltaS );
344 CS[1] = sin( deltaS );
345 prop[0] = t1 * CF[0] + t2 * CS[0];
346 prop[1] = t1 * CF[1] + t2 * CS[1];
347}
348
349double EvtDsToKSpipi0::DDalitz( double P1[4], double P2[4], double P3[4], int Ang,
350 double mass ) {
351 double pR[4], pD[4];
352 double temp_PDF, v_re;
353 temp_PDF = 0.0;
354 v_re = 0.0;
355 double B[2], s1, s2, s3, sR, sD;
356 for ( int i = 0; i < 4; i++ )
357 {
358 pR[i] = P1[i] + P2[i];
359 pD[i] = pR[i] + P3[i];
360 }
361 s1 = SCADot( P1, P1 );
362 s2 = SCADot( P2, P2 );
363 s3 = SCADot( P3, P3 );
364 sR = SCADot( pR, pR );
365 sD = SCADot( pD, pD );
366 // int G[4][4];
367 // for(int i=0; i!=4; i++){
368 // for(int j=0; j!=4; j++){
369 // if(i==j){
370 // if(i==0) G[i][j] = 1;
371 // else G[i][j] = -1;
372 // }
373 // else G[i][j] = 0;
374 // }
375 // }
376 if ( Ang == 0 )
377 {
378 B[0] = 1;
379 B[1] = 1;
380 temp_PDF = 1;
381 }
382 if ( Ang == 1 )
383 {
384 B[0] = barrier( 1, sR, s1, s2, 3.0, mass );
385 B[1] = barrier( 1, sD, sR, s3, 5.0, 1.9683 );
386 double t1[4], T1[4];
387 calt1( P1, P2, t1 );
388 calt1( pR, P3, T1 );
389 temp_PDF = 0;
390 for ( int i = 0; i < 4; i++ ) { temp_PDF += t1[i] * T1[i] * G[i][i]; }
391 }
392 if ( Ang == 2 )
393 {
394 B[0] = barrier( 2, sR, s1, s2, 3.0, mass );
395 B[1] = barrier( 2, sD, sR, s3, 5.0, 1.9683 );
396 double t2[4][4], T2[4][4];
397 calt2( P1, P2, t2 );
398 calt2( pR, P3, T2 );
399 temp_PDF = 0;
400 for ( int i = 0; i < 4; i++ )
401 {
402 for ( int j = 0; j < 4; j++ ) { temp_PDF += t2[i][j] * T2[j][i] * G[i][i] * G[j][j]; }
403 }
404 }
405 v_re = temp_PDF * B[0] * B[1];
406 return v_re;
407}
408
409void EvtDsToKSpipi0::calEva( double* Ks0, double* Pic, double* Pi0, double* mass1,
410 double* width1, double* amp, double* phase, int* g0,
411 int* modetype, int nstates, double& Result ) {
412 // double Ks0[4], Pic[4], Pi0[4];
413 // double rRes = 9.0;
414 double P12[4], P23[4], P13[4];
415 double cof[2], amp_PDF[2], PDF[2];
416 double scpi, snpi, sks0;
417 double s12, s13, s23;
418 for ( int i = 0; i < 4; i++ )
419 {
420 P12[i] = Pic[i] + Ks0[i];
421 P13[i] = Pi0[i] + Ks0[i];
422 P23[i] = Pic[i] + Pi0[i];
423 }
424 scpi = SCADot( Pic, Pic );
425 snpi = SCADot( Pi0, Pi0 );
426 sks0 = SCADot( Ks0, Ks0 );
427 s12 = SCADot( P12, P12 );
428 s13 = SCADot( P13, P13 );
429 s23 = SCADot( P23, P23 );
430 double pro[2], temp_PDF, amp_tmp[2];
431 double mass1sq;
432 amp_PDF[0] = 0;
433 amp_PDF[1] = 0;
434 PDF[0] = 0;
435 PDF[1] = 0;
436 amp_tmp[0] = 0;
437 amp_tmp[1] = 0;
438 for ( int i = 0; i < nstates; i++ )
439 {
440 amp_tmp[0] = 0;
441 amp_tmp[1] = 0;
442 mass1sq = mass1[i] * mass1[i];
443 cof[0] = amp[i] * cos( phase[i] );
444 cof[1] = amp[i] * sin( phase[i] );
445 temp_PDF = 0;
446 if ( modetype[i] == 23 || modetype[i] == 231 )
447 {
448 temp_PDF = DDalitz( Pic, Pi0, Ks0, 1, mass1[i] );
449 if ( g0[i] == 1 )
450 propagatorGS( mass1sq, mass1[i], width1[i], s23, scpi, snpi, 9.0, pro );
451 if ( g0[i] == 0 )
452 {
453 pro[0] = 1;
454 pro[1] = 0;
455 }
456 amp_tmp[0] = temp_PDF * pro[0];
457 amp_tmp[1] = temp_PDF * pro[1];
458 }
459 if ( modetype[i] == 12 )
460 {
461 temp_PDF = DDalitz( Ks0, Pic, Pi0, 1, mass1[i] );
462 if ( g0[i] == 1 )
463 propagatorRBW( mass1sq, mass1[i], width1[i], s12, sks0, scpi, rRes, 1, pro );
464 if ( g0[i] == 12 ) KPiSLASS( s12, sks0, scpi, pro );
465 if ( g0[i] == 0 )
466 {
467 pro[0] = 1;
468 pro[1] = 0;
469 }
470 amp_tmp[0] = temp_PDF * pro[0];
471 amp_tmp[1] = temp_PDF * pro[1];
472 }
473 if ( modetype[i] == 13 || modetype[i] == 131 )
474 {
475 temp_PDF = DDalitz( Ks0, Pi0, Pic, 1, mass1[i] );
476 if ( g0[i] == 1 )
477 propagatorRBW( mass1sq, mass1[i], width1[i], s13, sks0, snpi, rRes, 1, pro );
478 if ( g0[i] == 12 ) KPiSLASS( s13, sks0, snpi, pro );
479 if ( g0[i] == 0 )
480 {
481 pro[0] = 1;
482 pro[1] = 0;
483 }
484 amp_tmp[0] = temp_PDF * pro[0];
485 amp_tmp[1] = temp_PDF * pro[1];
486 }
487 Com_Multi( amp_tmp, cof, amp_PDF );
488 PDF[0] += amp_PDF[0];
489 PDF[1] += amp_PDF[1];
490 }
491 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
492 if ( value <= 0 ) value = 1e-20;
493 Result = value;
494 // cout<<"value = "<<value<<" Result = "<<Result<<endl;
495}
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)
EvtDecayBase * clone()
virtual ~EvtDsToKSpipi0()
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