BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtD0ToKpiEtap.cc
Go to the documentation of this file.
1// Environment:
2// This software is part of the EvtGen package developed jointly
3// for the BaBar and CLEO collaborations. If you use all or part
4// of it, please give an appropriate acknowledgement.
5//
6// Copyright Information: See EvtGen/COPYRIGHT
7// Copyright (C) 1998 Caltech, UCSB
8//
9// Module: EvtD0ToKpiEtap.cc
10//
11// Description: Routine to handle three-body decays of D0/D0_bar or D+/D-
12//
13// Modification history:
14//
15// Liaoyuan Dong Jul 3 23:56:42 2023 Module created
16// Mar 16 16:57:48 2024 Module updated
17//------------------------------------------------------------------------
18//
19#include "EvtD0ToKpiEtap.hh"
28#include <stdlib.h>
29using namespace std;
30
32
33void EvtD0ToKpiEtap::getName( std::string& model_name ) { model_name = "D0ToKpiEtap"; }
34
37 checkNArg( 0 );
38 checkNDaug( 3 );
43
44 phi[0] = 0.0;
45 rho[0] = 1.0;
46 phi[1] = -3.3097e+00;
47 rho[1] = 1.6085e-01; // 892
48 phi[2] = 3.2189e+00;
49 rho[2] = 7.4439e+00; // 1680
50 modetype[0] = 132;
51 modetype[1] = 12;
52 modetype[2] = 12;
53
54 mass[2] = 1.4273;
55 mass[1] = 0.89555;
56 mass[0] = 1.414;
57 width[2] = 0.100;
58 width[1] = 0.0473;
59 width[0] = 0.232;
60 const double mk0 = 0.497614;
61 const double mass_Kaon = 0.49368;
62 const double mass_Pion = 0.13957;
63 const double mass_Pi0 = 0.1349766;
64 const double meta = 0.547862;
65 const double metap = 0.95778;
66 mpi = 0.13957;
67 mD = 1.86486;
68 mDs = 1.9683;
69 rRes = 9.0;
70 rD = 5.0;
71 math_pi = 3.1415926;
72 GS1 = 0.636619783;
73 GS2 = 0.01860182466;
74 GS3 = 0.1591549458; // 1/(2*math_2pi)
75 GS4 = 0.00620060822; // mass_Pion2/math_pi
76
77 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
78 for ( int i = 0; i < 4; i++ )
79 {
80 for ( int j = 0; j < 4; j++ ) { G[i][j] = GG[i][j]; }
81 }
82}
84
86
87 // This piece of code could in principle be used to calculate maximum
88 // probablity on fly. But as it uses high number of points and model
89 // deals with single final state, we keep hardcoded number for now rather
90 // than adapting code to work here.
91 /*
92 double maxprob = 0.0;
93 for(int ir=0;ir<=60000000;ir++){
94 p->initializePhaseSpace(getNDaug(),getDaugs());
95 EvtVector4R D1 = p->getDaug(0)->getP4();
96 EvtVector4R D2 = p->getDaug(1)->getP4();
97 EvtVector4R D3 = p->getDaug(2)->getP4();
98
99 double P1[4], P2[4], P3[4];
100 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
101 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
102 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
103
104 double value;
105 int g0[3]={0,1,0};
106 int spin[3]={0,1,2};
107 /////////////////////
108 int nstates=3;
109 //
110 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value);
111
112 // value = calEva(P1, P2, P3);
113 if(value>maxprob) {
114 maxprob=value;
115 cout << "ir = " << ir << " maxProb= " << value << endl;
116 }
117 }
118 cout << "maxProb = " << maxprob << endl;
119
120 */
121 //////////////////////////////////////////////////////////
122
124 EvtVector4R D1 = p->getDaug( 0 )->getP4();
125 EvtVector4R D2 = p->getDaug( 1 )->getP4();
126 EvtVector4R D3 = p->getDaug( 2 )->getP4();
127
128 double P1[4], P2[4], P3[4];
129 P1[0] = D1.get( 0 );
130 P1[1] = D1.get( 1 );
131 P1[2] = D1.get( 2 );
132 P1[3] = D1.get( 3 );
133 P2[0] = D2.get( 0 );
134 P2[1] = D2.get( 1 );
135 P2[2] = D2.get( 2 );
136 P2[3] = D2.get( 3 );
137 P3[0] = D3.get( 0 );
138 P3[1] = D3.get( 1 );
139 P3[2] = D3.get( 2 );
140 P3[3] = D3.get( 3 );
141
142 double value;
143 int g0[3] = { 0, 1, 0 };
144 int spin[3] = { 0, 1, 2 };
145 /////////////////////
146 int nstates = 3;
147 calEva( P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value );
148
149 setProb( value );
150 return;
151}
152
153void EvtD0ToKpiEtap::Com_Multi( double a1[2], double a2[2], double res[2] ) {
154 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
155 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
156}
157void EvtD0ToKpiEtap::Com_Divide( double a1[2], double a2[2], double res[2] ) {
158 double tmp = a2[0] * a2[0] + a2[1] * a2[1];
159 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / tmp;
160 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / tmp;
161}
162//------------base---------------------------------
163double EvtD0ToKpiEtap::SCADot( double a1[4], double a2[4] ) {
164 double _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
165 return _cal;
166}
167double EvtD0ToKpiEtap::barrier( int l, double sa, double sb, double sc, double r,
168 double mass ) {
169 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
170 if ( q < 0 ) q = -q;
171 double z;
172 z = q * r * r;
173 double sa0;
174 sa0 = mass * mass;
175 double q0 = ( sa0 + sb - sc ) * ( sa0 + sb - sc ) / ( 4 * sa0 ) - sb;
176 if ( q0 < 0 ) q0 = -1 * q0;
177 double z0 = q0 * r * r;
178 double F = 0.0;
179 if ( l == 0 ) F = 1;
180 if ( l == 1 ) F = sqrt( ( 1 + z0 ) / ( 1 + z ) );
181 if ( l == 2 ) F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) );
182 return F;
183}
184void EvtD0ToKpiEtap::calt1( double daug1[4], double daug2[4], double t1[4] ) {
185 double p, pq, tmp;
186 double pa[4], qa[4];
187 for ( int i = 0; i < 4; i++ )
188 {
189 pa[i] = daug1[i] + daug2[i];
190 qa[i] = daug1[i] - daug2[i];
191 }
192 p = SCADot( pa, pa );
193 pq = SCADot( pa, qa );
194 tmp = pq / p;
195 for ( int i = 0; i < 4; i++ ) { t1[i] = qa[i] - tmp * pa[i]; }
196}
197void EvtD0ToKpiEtap::calt2( double daug1[4], double daug2[4], double t2[4][4] ) {
198 double p, r;
199 double pa[4], t1[4];
200 calt1( daug1, daug2, t1 );
201 r = SCADot( t1, t1 ) / 3.0;
202 for ( int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
203 p = SCADot( pa, pa );
204 for ( int i = 0; i < 4; i++ )
205 {
206 for ( int j = 0; j < 4; j++ )
207 { t2[i][j] = t1[i] * t1[j] - r * ( G[i][j] - pa[i] * pa[j] / p ); }
208 }
209}
210//-------------------prop--------------------------------------------
211double EvtD0ToKpiEtap::wid( double mass2, double mass, double sa, double sb, double sc,
212 double r2, int l ) {
213 double widm = 0.;
214 double m = sqrt( sa );
215 double tmp = sb - sc;
216 double tmp1 = sa + tmp;
217 double q = 0.25 * tmp1 * tmp1 / sa - sb;
218 if ( q < 0 ) q = -q;
219 double tmp2 = mass2 + tmp;
220 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
221 if ( q0 < 0 ) q0 = -q0;
222 double z = q * r2;
223 double z0 = q0 * r2;
224 double t = q / q0;
225 if ( l == 0 ) { widm = sqrt( t ) * mass / m; }
226 else if ( l == 1 ) { widm = t * sqrt( t ) * mass / m * ( 1 + z0 ) / ( 1 + z ); }
227 else if ( l == 2 )
228 { widm = t * t * sqrt( t ) * mass / m * ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ); }
229 return widm;
230}
231double EvtD0ToKpiEtap::widl1( double mass2, double mass, double sa, double sb, double sc,
232 double r2 ) {
233 double widm = 0.;
234 double m = sqrt( sa );
235 double tmp = sb - sc;
236 double tmp1 = sa + tmp;
237 double q = 0.25 * tmp1 * tmp1 / sa - sb;
238 if ( q < 0 ) q = -q;
239 double tmp2 = mass2 + tmp;
240 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
241 if ( q0 < 0 ) q0 = -q0;
242 double z = q * r2;
243 double z0 = q0 * r2;
244 double F = ( 1 + z0 ) / ( 1 + z );
245 double t = q / q0;
246 widm = t * sqrt( t ) * mass / m * F;
247 return widm;
248}
249void EvtD0ToKpiEtap::propagatorRBW( double mass2, double mass, double width, double sa,
250 double sb, double sc, double r2, int l, double prop[2] ) {
251 double a[2], b[2];
252 a[0] = 1;
253 a[1] = 0;
254 double q = 0.25 * ( sa + sb - sc ) * ( sa + sb - sc ) / sa - sb;
255
256 b[0] = mass2 - sa;
257 b[1] = -mass * width * wid( mass2, mass, sa, sb, sc, r2, l );
258 Com_Divide( a, b, prop );
259}
260
261void EvtD0ToKpiEtap::propagatorFlatte( double mass, double width, double sa, double sb,
262 double sc, int r, double prop[2] ) {
263 double q, qKsK, qetapi;
264 // double qKsK,qetapi;
265 double rhoab[2], rhoKsK[2];
266 q = 0.25 * ( sa + sb - sc ) * ( sa + sb - sc ) / sa - sb;
267 qetapi = 0.25 * ( sa + 0.547862 - 0.13957039 ) * ( sa + 0.547862 - 0.13957039 ) / sa -
268 0.547862 * 0.547862;
269 if ( r == 0 ) qKsK = 0.25 * sa - 0.49368 * 0.49368;
270 if ( r == 1 ) qKsK = 0.25 * ( sa + sb - sc ) * ( sa + sb - sc ) / sa - sb;
271 if ( qetapi > 0 )
272 {
273 rhoab[0] = 2 * sqrt( qetapi / sa );
274 rhoab[1] = 0;
275 }
276 if ( qetapi < 0 )
277 {
278 rhoab[0] = 0;
279 rhoab[1] = 2 * sqrt( -qetapi / sa );
280 }
281 if ( qKsK > 0 )
282 {
283 rhoKsK[0] = 2 * sqrt( qKsK / sa );
284 rhoKsK[1] = 0;
285 }
286 if ( qKsK < 0 )
287 {
288 rhoKsK[0] = 0;
289 rhoKsK[1] = 2 * sqrt( -qKsK / sa );
290 }
291 double a[2], b[2];
292 a[0] = 1;
293 a[1] = 0;
294 b[0] = mass * mass - sa + 0.341 * rhoab[1] + 0.892 * 0.341 * rhoKsK[1];
295 b[1] = -( 0.341 * rhoab[0] + 0.892 * 0.341 * rhoKsK[0] );
296 Com_Divide( a, b, prop );
297}
298
299void EvtD0ToKpiEtap::propagatorGS( double mass2, double mass, double width, double sa,
300 double sb, double sc, double r2, double prop[2] ) {
301 double a[2], b[2];
302 double tmp = sb - sc;
303 double tmp1 = sa + tmp;
304 double q2 = 0.25 * tmp1 * tmp1 / sa - sb;
305 if ( q2 < 0 ) q2 = 1e-16;
306
307 double tmp2 = mass2 + tmp;
308 double q02 = 0.25 * tmp2 * tmp2 / mass2 - sb;
309 if ( q02 < 0 ) q02 = 1e-16;
310
311 double q = sqrt( q2 );
312 double q0 = sqrt( q02 );
313 double m = sqrt( sa );
314 double q03 = q0 * q02;
315 // double tmp3 = log(mass+2*q0)+1.2760418309; // log(mass_2Pion) = 1.2760418309;
316 double tmp3 = log( mass + 2 * q0 ) + 1.2926305904; // log(mpi0+mpip) = -1.2926305904;
317
318 double h = GS1 * q / m * ( log( m + 2 * q ) + 1.2926305904 );
319 double h0 = GS1 * q0 / mass * tmp3;
320 double dh = h0 * ( 0.125 / q02 - 0.5 / mass2 ) + GS3 / mass2;
321 double d = GS2 / q02 * tmp3 + GS3 * mass / q0 - GS4 * mass / q03;
322 double f = mass2 / q03 * ( q2 * ( h - h0 ) + ( mass2 - sa ) * q02 * dh );
323
324 a[0] = 1.0 + d * width / mass;
325 a[1] = 0.0;
326 b[0] = mass2 - sa + width * f;
327 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r2 );
328 Com_Divide( a, b, prop );
329}
330
331void EvtD0ToKpiEtap::KPiSLASS( double sa, double sb, double sc, double prop[2] ) {
332 const double m1430 = 1.441;
333 const double sa0 = 1.441 * 1.441; // m1430*m1430;
334 const double w1430 = 0.193;
335 const double Lass1 = 0.25 / sa0;
336 double tmp = sb - sc;
337 double tmp1 = sa0 + tmp;
338 double q0 = Lass1 * tmp1 * tmp1 - sb;
339 if ( q0 < 0 ) q0 = -1 * q0;
340 double tmp2 = sa + tmp;
341 double qs = 0.25 * tmp2 * tmp2 / sa - sb;
342 double q = sqrt( qs );
343 double width = w1430 * q * m1430 / sqrt( sa * q0 );
344 double temp_R = atan( m1430 * width / ( sa0 - sa ) );
345 if ( temp_R < 0 ) temp_R += math_pi;
346 double deltaR = -1.915 + temp_R; // fiR=-109.7
347 double temp_F =
348 atan( 0.226 * q / ( 2.0 - 3.819 * qs ) ); // 2.0*0.113 = 0.226; 0.113*33.8 = 3.819
349 if ( temp_F < 0 ) temp_F += math_pi;
350 double deltaF = 0.002 + temp_F; // fiF=0.1
351 double deltaS = deltaR + 2.0 * deltaF;
352 double t1 = 0.96 * sin( deltaF );
353 double t2 = sin( deltaR );
354 double CF[2], CS[2];
355 CF[0] = cos( deltaF );
356 CF[1] = sin( deltaF );
357 CS[0] = cos( deltaS );
358 CS[1] = sin( deltaS );
359 prop[0] = t1 * CF[0] + t2 * CS[0];
360 prop[1] = t1 * CF[1] + t2 * CS[1];
361}
362double EvtD0ToKpiEtap::DDalitz( double P1[4], double P2[4], double P3[4], int Ang,
363 double mass ) {
364 double pR[4], pD[4];
365 double temp_PDF, v_re;
366 temp_PDF = 0.0;
367 v_re = 0.0;
368 double B[2], s1, s2, s3, sR, sD;
369 for ( int i = 0; i < 4; i++ )
370 {
371 pR[i] = P1[i] + P2[i];
372 pD[i] = pR[i] + P3[i];
373 }
374 s1 = SCADot( P1, P1 );
375 s2 = SCADot( P2, P2 );
376 s3 = SCADot( P3, P3 );
377 sR = SCADot( pR, pR );
378 sD = SCADot( pD, pD );
379 if ( Ang == 0 )
380 {
381 B[0] = 1;
382 B[1] = 1;
383 temp_PDF = 1;
384 }
385 if ( Ang == 1 )
386 {
387 B[0] = barrier( 1, sR, s1, s2, 3.0, mass );
388 B[1] = barrier( 1, sD, sR, s3, 5.0, 1.86484 );
389 // B[0] = Barrier(1,sR,s1,s2,9.0);
390 // B[1] = Barrier(1,sD,sR,s3,25.0);
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.86484 );
401 // B[0] = Barrier(2,sR,s1,s2,9.0);
402 // B[1] = Barrier(2,sD,sR,s3,25.0);
403 double t2[4][4], T2[4][4];
404 calt2( P1, P2, t2 );
405 calt2( pR, P3, T2 );
406 temp_PDF = 0;
407 for ( int i = 0; i < 4; i++ )
408 {
409 for ( int j = 0; j < 4; j++ ) { temp_PDF += t2[i][j] * T2[j][i] * G[i][i] * G[j][j]; }
410 }
411 }
412 v_re = temp_PDF * B[0] * B[1];
413 return v_re;
414}
415
416void EvtD0ToKpiEtap::calEva( double* Ks0, double* Kc, double* Pi0, double* mass1,
417 double* width1, double* amp, double* phase, int* g0, int* spin,
418 int* modetype, int nstates, double& Result ) {
419 double P12[4], P23[4], P13[4];
420 double cof[2], amp_PDF[2], PDF[2];
421 double snpi, sck, sks0;
422 double s12, s13, s23;
423 for ( int i = 0; i < 4; i++ )
424 {
425 P12[i] = Kc[i] + Ks0[i];
426 P13[i] = Pi0[i] + Ks0[i];
427 P23[i] = Kc[i] + Pi0[i];
428 }
429 sck = SCADot( Kc, Kc );
430 snpi = SCADot( Pi0, Pi0 );
431 sks0 = SCADot( Ks0, Ks0 );
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
452 if ( modetype[i] == 12 )
453 {
454 temp_PDF = DDalitz( Ks0, Kc, Pi0, spin[i], mass1[i] );
455 if ( g0[i] == 1 )
456 propagatorRBW( mass1sq, mass1[i], width1[i], s12, sks0, sck, rRes, spin[i], pro );
457 // if(g0[i]==2) KPiSLASS(s12,sks0,scpi,pro);
458 if ( g0[i] == 2 )
459 {
460 propagatorFlatte( mass1[i], width1[i], s12, sks0, sck, 1, pro );
461 pro[0] = pro[0] * ( 0.01 + 0.990 * 0.990 ) / ( 0.01 + s12 );
462 pro[1] = pro[1] * ( 0.01 + 0.990 * 0.990 ) / ( 0.01 + s12 );
463 // pro[0]=pro[0]*(0.6*0.6)/(0.6*0.6)+(s12-0.990*0.990)*(s12-0.990*0.990);
464 // pro[1]=pro[1]*(0.6*0.6)/(0.6*0.6)+(s12-0.990*0.990)*(s12-0.990*0.990);
465 }
466 if ( g0[i] == 3 )
467 propagatorFlatte( mass1[i], width1[i], s12, sks0, sck, 1, pro ); // Only for a0(980)
468
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
478 if ( modetype[i] == 132 )
479 {
480 KPiSLASS( s12, sks0, sck, pro );
481 amp_tmp[0] = pro[0];
482 amp_tmp[1] = pro[1];
483 }
484 Com_Multi( amp_tmp, cof, amp_PDF );
485 PDF[0] += amp_PDF[0];
486 PDF[1] += amp_PDF[1];
487 }
488 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
489 if ( value <= 0 ) value = 1e-20;
490 Result = value;
491}
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 decay(EvtParticle *p)
EvtDecayBase * clone()
virtual ~EvtD0ToKpiEtap()
void getName(std::string &name)
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)
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