BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDToKKpi.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2// Environment:
3// This software is part of models developed at BES collaboration
4// based on the EvtGen framework. If you use all or part
5// of it, please give an appropriate acknowledgement.
6//
7// Copyright Information: See EvtGen/BesCopyright
8// Copyright (A) 2006 Ping Rong-Gang @IHEP
9//
10// Module: EvtDToKKpi.cc
11// the necessary file: Evt.hh
12//
13// Description: D+ -> K- K+ pi+,
14//
15//
16// Modification history:
17// Liaoyuan Dong Fri Dec 02 23:35:56 CST 2022 Module created
18//------------------------------------------------------------------------
19
20#include "EvtDToKKpi.hh"
21#include "../EvtGenBase/EvtConst.hh"
22#include "../EvtGenBase/EvtDecayTable.hh"
23#include "../EvtGenBase/EvtGenKine.hh"
24#include "../EvtGenBase/EvtPDL.hh"
25#include "../EvtGenBase/EvtParticle.hh"
26#include "../EvtGenBase/EvtPatches.hh"
27#include "../EvtGenBase/EvtReport.hh"
28#include "../EvtGenBase/EvtResonance.hh"
29#include "../EvtGenBase/EvtResonance2.hh"
30#include <fstream>
31#include <stdlib.h>
32#include <string>
33using std::endl;
34
35EvtDToKKpi::~EvtDToKKpi() {}
36
37void EvtDToKKpi::getName( std::string& model_name ) { model_name = "DToKKpi"; }
38
39EvtDecayBase* EvtDToKKpi::clone() { return new EvtDToKKpi; }
40
41void EvtDToKKpi::init() {
42 // check that there are 0 arguments
43 checkNArg( 0 );
44 checkNDaug( 3 );
45 checkSpinParent( EvtSpinType::SCALAR );
46 checkSpinDaughter( 0, EvtSpinType::SCALAR );
47 checkSpinDaughter( 1, EvtSpinType::SCALAR );
48 checkSpinDaughter( 2, EvtSpinType::SCALAR );
49
50 mass[0] = 0.89581; // Ksta0
51 mass[1] = 1.019461; // phi1020
52 mass[2] = 0.965; // f0(980)
53 mass[3] = 1.474; // a0(1450)
54 mass[4] = 1.68; // phi1680
55 mass[5] = 1.4324; // K21430
56 mass[6] = 0.89581; // KPi Swave
57 mass[7] = 1.3183; // a21320
58 width[0] = 0.0474; // Ksta0
59 width[1] = 0.004266; // phi1020
60 width[2] = 0.0400; // f0(980)
61 width[3] = 0.265; // a01450
62 width[4] = 0.15; // phi1680
63 width[5] = 0.109; // K21430
64 width[6] = 0.0474; // KPi Swave
65 width[7] = 0.107; // a01320
66
67 rho[0] = 1; // K*892
68 phi[0] = 0; //
69
70 rho[1] = 1.121155232058; // phi1020
71 phi[1] = 3.322105212534; //
72
73 rho[2] = 1.214264427529; // f0980
74 phi[2] = -5.321965906572;
75
76 rho[3] = 1.086952630367; // a01450
77 phi[3] = -4.213398459878;
78
79 rho[4] = 1.121071578335; // phi1680
80 phi[4] = -1.055595448689; //
81
82 rho[5] = -1.473195648049; // K21430
83 phi[5] = 5.841826108667;
84
85 rho[6] = 5.396864431453; // KPISW
86 phi[6] = -0.360227442474;
87
88 rho[7] = 0.436480711396; // a21320
89 phi[7] = 0.778423040673;
90
91 spin[0] = 1; // K*892
92 spin[1] = 1; // phi1020
93 spin[2] = 0; // f0980
94 spin[3] = 0; // a01450
95 spin[4] = 1; // phi1680
96 spin[5] = 2; // K21430
97 spin[6] = 0; // KPISW
98 spin[7] = 2; // a21320
99
100 modetype[0] = 13; // K*892
101 modetype[1] = 12; // phi1020
102 modetype[2] = 12; // f0980
103 modetype[3] = 12; // a01450
104 modetype[4] = 12; // phi1680
105 modetype[5] = 13; // K21430 error
106 modetype[6] = 13; // KPiSW
107 modetype[7] = 12; // a21320 error
108
109 /*
110 std::cout << "EvtDToKKpi (May 01, 2023) ==> Initialization" << std::endl;
111 for (int i=0; i<8; i++) {
112 cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
113 }
114 */
115
116 mDp = 1.86966;
117 mK0 = 0.497611;
118 mKa = 0.49368;
119 mPi = 0.13957;
120 mK02 = 0.237616707;
121 mPi2 = 0.01947978;
122 mass_EtaP = 0.95778;
123 mass_Kaon = 0.49368;
124
125 math_pi = 3.1415926;
126 mass_Pion = 0.13957;
127 mass_Pion2 = 0.0194797849;
128 mass_2Pion = 0.27914;
129 math_2pi = 6.2831852;
130 rD2 = 25.0; // 5*5
131 rRes2 = 9.0; // 3*3
132 // g1 = 0.5468;
133 g2 = 0.23;
134 GS1 = 0.636619783;
135 GS2 = 0.01860182466;
136 GS3 = 0.1591549458;
137 GS4 = 0.00620060822;
138
139 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
140 for ( int i = 0; i < 4; i++ )
141 {
142 for ( int j = 0; j < 4; j++ ) { G[i][j] = GG[i][j]; }
143 }
144}
145
146void EvtDToKKpi::initProbMax() {
147 setProbMax( 11700.0 ); // MAXprob = 11690.0512849074
148}
149
150void EvtDToKKpi::decay( EvtParticle* p ) {
151 //==================================================get
152 // max============================================
153 /*
154 double maxprob = 0.0;
155 for(int ir=0;ir<=60000000;ir++){
156 p->initializePhaseSpace(getNDaug(),getDaugs());
157 EvtVector4R D1 = p->getDaug(0)->getP4();
158 EvtVector4R D2 = p->getDaug(1)->getP4();
159 EvtVector4R D3 = p->getDaug(2)->getP4();
160
161 double P1[4], P2[4], P3[4];
162 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
163 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
164 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
165
166 double value;
167 int nstates=8;
168 int g0[8]={1,1,4,1,1,1,3,1};
169 double r0[8]={3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0};
170 double r1[8]={5.0,5.0,5.0,5.0,5.0,5.0,5.0,5.0};
171
172 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value,r0,r1);
173
174 if (value<0) continue;
175 if(value>maxprob) {
176 maxprob=value;
177 cout << "ir= " << ir << endl;
178 cout << "MAX====> " << maxprob << endl;
179 }
180 }
181 printf("MAXprob = %.10f\n",maxprob);
182 */
183 //==================================================get
184 // max============================================
185
186 p->initializePhaseSpace( getNDaug(), getDaugs() );
187 EvtVector4R D1 = p->getDaug( 0 )->getP4();
188 EvtVector4R D2 = p->getDaug( 1 )->getP4();
189 EvtVector4R D3 = p->getDaug( 2 )->getP4();
190
191 double P1[4], P2[4], P3[4];
192 P1[0] = D1.get( 0 );
193 P1[1] = D1.get( 1 );
194 P1[2] = D1.get( 2 );
195 P1[3] = D1.get( 3 );
196 P2[0] = D2.get( 0 );
197 P2[1] = D2.get( 1 );
198 P2[2] = D2.get( 2 );
199 P2[3] = D2.get( 3 );
200 P3[0] = D3.get( 0 );
201 P3[1] = D3.get( 1 );
202 P3[2] = D3.get( 2 );
203 P3[3] = D3.get( 3 );
204
205 // Check
206 // P1[0] = 0.694086; P1[1] = 0.343301; P1[2] = 0.263795; P1[3] = -0.224934;
207 // P2[0] = 0.719665; P2[1] = -0.165392; P2[2] = -0.495843; P2[3] = -0.0314011;
208 // P3[0] = 0.474568; P3[1] = 0.0149532 ; P3[2] = 0.405429; P3[3] = 0.202826;
209
210 double value;
211 int nstates = 8;
212 int g0[8] = { 1, 1, 4, 1, 1, 1, 3, 1 };
213 double r0[8] = { 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0 };
214 double r1[8] = { 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0 };
215
216 calEva( P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value, r0, r1 );
217 setProb( value );
218 return;
219}
220
221void EvtDToKKpi::Com_Multi( double a1[2], double a2[2], double res[2] ) {
222 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
223 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
224}
225void EvtDToKKpi::Com_Divide( double a1[2], double a2[2], double res[2] ) {
226 double tmp = a2[0] * a2[0] + a2[1] * a2[1];
227 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / tmp;
228 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / tmp;
229}
230//------------base---------------------------------
231double EvtDToKKpi::SCADot( double a1[4], double a2[4] ) {
232 double _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
233 return _cal;
234}
235double EvtDToKKpi::barrier( int l, double sa, double sb, double sc, double r, double mass ) {
236 double q = fabs( ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb );
237 double z;
238 z = q * r * r;
239 double sa0;
240 sa0 = mass * mass;
241 double q0 = fabs( ( sa0 + sb - sc ) * ( sa0 + sb - sc ) / ( 4 * sa0 ) - sb );
242 double z0 = q0 * r * r;
243 double F = 0.0;
244 if ( l == 0 ) F = 1;
245 if ( l == 1 ) F = sqrt( ( 1 + z0 ) / ( 1 + z ) );
246 if ( l == 2 ) F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) );
247 return F;
248}
249//------------------spin-------------------------------------------
250void EvtDToKKpi::calt1( double daug1[4], double daug2[4], double t1[4] ) {
251 double p, pq, tmp;
252 double pa[4], qa[4];
253 for ( int i = 0; i < 4; i++ )
254 {
255 pa[i] = daug1[i] + daug2[i];
256 qa[i] = daug1[i] - daug2[i];
257 }
258 p = SCADot( pa, pa );
259 pq = SCADot( pa, qa );
260 tmp = pq / p;
261 for ( int i = 0; i < 4; i++ ) { t1[i] = qa[i] - tmp * pa[i]; }
262}
263void EvtDToKKpi::calt2( double daug1[4], double daug2[4], double t2[4][4] ) {
264 double p, r;
265 double pa[4], t1[4];
266 calt1( daug1, daug2, t1 );
267 r = SCADot( t1, t1 ) / 3.0;
268 for ( int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
269 p = SCADot( pa, pa );
270 for ( int i = 0; i < 4; i++ )
271 {
272 for ( int j = 0; j < 4; j++ )
273 { t2[i][j] = t1[i] * t1[j] - r * ( G[i][j] - pa[i] * pa[j] / p ); }
274 }
275}
276//-------------------prop--------------------------------------------
277void EvtDToKKpi::propagator( double mass2, double mass, double width, double sx,
278 double prop[2] ) {
279 double a[2], b[2];
280 a[0] = 1;
281 a[1] = 0;
282 b[0] = mass2 - sx;
283 b[1] = -mass * width;
284 Com_Divide( a, b, prop );
285}
286double EvtDToKKpi::wid( double mass2, double mass, double sa, double sb, double sc, double r2,
287 int l ) {
288 double widm = 0.;
289 double m = sqrt( sa );
290 double tmp = sb - sc;
291 double tmp1 = sa + tmp;
292 double q = fabs( 0.25 * tmp1 * tmp1 / sa - sb );
293 double tmp2 = mass2 + tmp;
294 double q0 = fabs( 0.25 * tmp2 * tmp2 / mass2 - sb );
295 double z = q * r2;
296 double z0 = q0 * r2;
297 double t = q / q0;
298 if ( l == 0 ) { widm = sqrt( t ) * mass / m; }
299 else if ( l == 1 ) { widm = t * sqrt( t ) * mass / m * ( 1 + z0 ) / ( 1 + z ); }
300 else if ( l == 2 )
301 { widm = t * t * sqrt( t ) * mass / m * ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ); }
302 return widm;
303}
304double EvtDToKKpi::widl1( double mass2, double mass, double sa, double sb, double sc,
305 double r2 ) {
306 double widm = 0.;
307 double m = sqrt( sa );
308 double tmp = sb - sc;
309 double tmp1 = sa + tmp;
310 double q = fabs( 0.25 * tmp1 * tmp1 / sa - sb );
311 double tmp2 = mass2 + tmp;
312 double q0 = fabs( 0.25 * tmp2 * tmp2 / mass2 - sb );
313 double z = q * r2;
314 double z0 = q0 * r2;
315 double F = ( 1 + z0 ) / ( 1 + z );
316 double t = q / q0;
317 widm = t * sqrt( t ) * mass / m * F;
318 return widm;
319}
320void EvtDToKKpi::propagatorRBW( double mass, double width, double sa, double sb, double sc,
321 double r2, int l, double prop[2] ) {
322
323 double a[2], b[2];
324 double mass2 = mass * mass;
325 a[0] = 1;
326 a[1] = 0;
327 b[0] = mass2 - sa;
328 b[1] = -mass * width * wid( mass2, mass, sa, sb, sc, r2, l );
329 Com_Divide( a, b, prop );
330}
331
332void EvtDToKKpi::propagatorFlatte( double mass, double width, double sa, double prop[2] ) {
333 double q2_Pi, q2_Ka;
334 double rhoPi[2], rhoKa[2];
335
336 q2_Pi = 0.25 * sa - mPi * mPi;
337 q2_Ka = 0.25 * sa - mKa * mKa;
338
339 if ( q2_Pi > 0 )
340 {
341 rhoPi[0] = 2.0 * sqrt( q2_Pi / sa );
342 rhoPi[1] = 0.0;
343 }
344 if ( q2_Pi <= 0 )
345 {
346 rhoPi[0] = 0.0;
347 rhoPi[1] = 2.0 * sqrt( -q2_Pi / sa );
348 }
349
350 if ( q2_Ka > 0 )
351 {
352 rhoKa[0] = 2.0 * sqrt( q2_Ka / sa );
353 rhoKa[1] = 0.0;
354 }
355 if ( q2_Ka <= 0 )
356 {
357 rhoKa[0] = 0.0;
358 rhoKa[1] = 2.0 * sqrt( -q2_Ka / sa );
359 }
360
361 double a[2], b[2];
362 a[0] = 1;
363 a[1] = 0;
364 b[0] = mass * mass - sa + 0.165 * rhoPi[1] + 0.69465 * rhoKa[1];
365 b[1] = -( 0.165 * rhoPi[0] + 0.69465 * rhoKa[0] );
366 Com_Divide( a, b, prop );
367}
368
369void EvtDToKKpi::rhoab( double sa, double sb, double sc, double res[2] ) {
370 double tmp = sa + sb - sc;
371 double q = 0.25 * tmp * tmp / sa - sb;
372 if ( q >= 0 )
373 {
374 res[0] = 2.0 * sqrt( q / sa );
375 res[1] = 0.0;
376 }
377 else
378 {
379 res[0] = 0.0;
380 res[1] = 2.0 * sqrt( -q / sa );
381 }
382}
383void EvtDToKKpi::rho4Pi( double sa, double res[2] ) {
384 double temp = 1.0 - 0.3116765584 / sa; // 0.3116765584=0.13957*0.13957*16
385 if ( temp >= 0 )
386 {
387 res[0] = sqrt( temp ) / ( 1.0 + exp( 9.8 - 3.5 * sa ) );
388 res[1] = 0.0;
389 }
390 else
391 {
392 res[0] = 0.0;
393 res[1] = sqrt( -temp ) / ( 1.0 + exp( 9.8 - 3.5 * sa ) ); //?????????????????????
394 }
395}
396
397void EvtDToKKpi::propagatorsigma500( double sa, double sb, double sc, double prop[2] ) {
398 double f = 0.5843 + 1.6663 * sa;
399 const double M = 0.9264; // M(f0500)
400 const double mass2 = 0.85821696; // M*M
401 const double mpi2d2 = 0.00973989245; // mass_Pion2/2 = 0.0194797849/2
402 double g1 = f * ( sa - mpi2d2 ) / ( mass2 - mpi2d2 ) * exp( ( mass2 - sa ) / 1.082 );
403 double rho1s[2], rho1M[2], rho2s[2], rho2M[2], rho1[2], rho2[2];
404 rhoab( sa, sb, sc, rho1s );
405 rhoab( mass2, sb, sc, rho1M );
406 rho4Pi( sa, rho2s );
407 rho4Pi( mass2, rho2M );
408 Com_Divide( rho1s, rho1M, rho1 );
409 Com_Divide( rho2s, rho2M, rho2 );
410 double a[2], b[2];
411 a[0] = 1.0;
412 a[1] = 0.0;
413 b[0] = mass2 - sa + M * ( g1 * rho1[1] + 0.0024 * rho2[1] );
414 b[1] = -M * ( g1 * rho1[0] + 0.0024 * rho2[0] );
415 Com_Divide( a, b, prop );
416}
417void EvtDToKKpi::Flatte_rhoab( double sa, double sb, double rho[2] ) {
418 double q = 1.0 - ( 4 * sb / sa );
419
420 if ( q > 0 )
421 {
422 rho[0] = sqrt( q );
423 rho[1] = 0;
424 }
425 else if ( q < 0 )
426 {
427 rho[0] = 0;
428 rho[1] = sqrt( -q );
429 }
430}
431
432void EvtDToKKpi::propagator980( double mass, double sx, double* sb, double prop[2] ) {
433
434 double gpipi1 = 2.0 / 3.0 * 0.165;
435 double gpipi2 = 1.0 / 3.0 * 0.165;
436 double gKK1 = 0.5 * 0.69465;
437 double gKK2 = 0.5 * 0.69465;
438
439 double unit[2] = { 1.0 };
440 double ci[2] = { 0, 1 };
441 double rho1[2];
442 Flatte_rhoab( sx, sb[0], rho1 );
443 double rho2[2];
444 Flatte_rhoab( sx, sb[1], rho2 );
445 double rho3[2];
446 Flatte_rhoab( sx, sb[2], rho3 );
447 double rho4[2];
448 Flatte_rhoab( sx, sb[3], rho4 );
449
450 double tmp1[2] = { gpipi1, 0 };
451 double tmp11[2];
452 double tmp2[2] = { gpipi2, 0 };
453 double tmp22[2];
454 double tmp3[2] = { gKK1, 0 };
455 double tmp33[2];
456 double tmp4[2] = { gKK2, 0 };
457 double tmp44[2];
458
459 Com_Multi( tmp1, rho1, tmp11 );
460 Com_Multi( tmp2, rho2, tmp22 );
461 Com_Multi( tmp3, rho3, tmp33 );
462 Com_Multi( tmp4, rho4, tmp44 );
463
464 double tmp5[2] = { tmp11[0] + tmp22[0] + tmp33[0] + tmp44[0],
465 tmp11[1] + tmp22[1] + tmp33[1] + tmp44[1] };
466 double tmp51[2];
467 Com_Multi( tmp5, ci, tmp51 );
468 double tmp6[2] = { mass * mass - sx - tmp51[0], -1.0 * tmp51[1] };
469
470 Com_Divide( unit, tmp6, prop );
471}
472
473void EvtDToKKpi::KPiSLASS( double sa, double sb, double sc, double prop[2] ) {
474 const double m1430 = 1.441;
475 const double sa0 = 2.076481; // m1430*m1430;
476 const double w1430 = 0.193;
477 const double Lass1 = 0.25 / sa0;
478 double tmp = sb - sc;
479 double tmp1 = sa0 + tmp;
480 double q0 = fabs( Lass1 * tmp1 * tmp1 - sb );
481 double tmp2 = sa + tmp;
482 double qs = 0.25 * tmp2 * tmp2 / sa - sb;
483 double q = sqrt( qs );
484 double width = w1430 * q * m1430 / sqrt( sa * q0 );
485 double temp_R = atan( m1430 * width / ( sa0 - sa ) );
486 if ( temp_R < 0 ) temp_R += math_pi;
487 double deltaR = -109.7 * math_pi / 180.0 + temp_R;
488 double temp_F =
489 atan( 0.226 * q / ( 2.0 - 3.8194 * qs ) ); // 2.0*0.113 = 0.226; -33.8*0.113 = -3.8194
490 if ( temp_F < 0 ) temp_F += math_pi;
491 double deltaF = 0.1 * math_pi / 180.0 + temp_F;
492 double deltaS = deltaR + 2.0 * deltaF;
493 double t1 = 0.96 * sin( deltaF );
494 double t2 = sin( deltaR );
495 double CF[2], CS[2];
496 CF[0] = cos( deltaF );
497 CF[1] = sin( deltaF );
498 CS[0] = cos( deltaS );
499 CS[1] = sin( deltaS );
500 prop[0] = t1 * CF[0] + t2 * CS[0];
501 prop[1] = t1 * CF[1] + t2 * CS[1];
502}
503//------------GS---used by rho----------------------------
504void EvtDToKKpi::propagatorGS( double mass2, double mass, double width, double sa, double sb,
505 double sc, double r2, double prop[2] ) {
506 double a[2], b[2];
507 double tmp = sb - sc;
508 double tmp1 = sa + tmp;
509 double q2 = fabs( 0.25 * tmp1 * tmp1 / sa - sb );
510
511 double tmp2 = mass2 + tmp;
512 double q02 = fabs( 0.25 * tmp2 * tmp2 / mass2 - sb );
513
514 double q = sqrt( q2 );
515 double q0 = sqrt( q02 );
516 double m = sqrt( sa );
517 double q03 = q0 * q02;
518 double tmp3 = log( mass + 2 * q0 ) + 1.2926305904; // log(mpi0+mpip) = -1.2926305904;
519
520 double h = GS1 * q / m * ( log( m + 2 * q ) + 1.2926305904 );
521 double h0 = GS1 * q0 / mass * tmp3;
522 double dh = h0 * ( 0.125 / q02 - 0.5 / mass2 ) + GS3 / mass2;
523 double d = GS2 / q02 * tmp3 + GS3 * mass / q0 - GS4 * mass / q03;
524 double f = mass2 / q03 * ( q2 * ( h - h0 ) + ( mass2 - sa ) * q02 * dh );
525
526 a[0] = 1.0 + d * width / mass;
527 a[1] = 0.0;
528 b[0] = mass2 - sa + width * f;
529 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r2 );
530 Com_Divide( a, b, prop );
531}
532
533double EvtDToKKpi::DDalitz( double P1[4], double P2[4], double P3[4], int Ang, double mass ) {
534 double pR[4], pD[4];
535 double temp_PDF, v_re;
536 temp_PDF = 0.0;
537 v_re = 0.0;
538 double B[2], s1, s2, s3, sR, sD;
539 for ( int i = 0; i < 4; i++ )
540 {
541 pR[i] = P1[i] + P2[i];
542 pD[i] = pR[i] + P3[i];
543 }
544 s1 = SCADot( P1, P1 );
545 s2 = SCADot( P2, P2 );
546 s3 = SCADot( P3, P3 );
547 sR = SCADot( pR, pR );
548 sD = SCADot( pD, pD );
549 int G[4][4];
550 for ( int i = 0; i != 4; i++ )
551 {
552 for ( int j = 0; j != 4; j++ )
553 {
554 if ( i == j )
555 {
556 if ( i == 0 ) G[i][j] = 1;
557 else G[i][j] = -1;
558 }
559 else G[i][j] = 0;
560 }
561 }
562 if ( Ang == 0 )
563 {
564 B[0] = 1;
565 B[1] = 1;
566 temp_PDF = 1;
567 }
568 if ( Ang == 1 )
569 {
570 B[0] = barrier( 1, sR, s1, s2, 3.0, mass );
571 B[1] = barrier( 1, sD, sR, s3, 5.0, mDp );
572 double t1[4], T1[4];
573 calt1( P1, P2, t1 );
574 calt1( pR, P3, T1 );
575 temp_PDF = 0;
576 for ( int i = 0; i < 4; i++ ) { temp_PDF += t1[i] * T1[i] * G[i][i]; }
577 }
578 if ( Ang == 2 )
579 {
580 B[0] = barrier( 2, sR, s1, s2, 3.0, mass );
581 B[1] = barrier( 2, sD, sR, s3, 5.0, mDp );
582 double t2[4][4], T2[4][4];
583 calt2( P1, P2, t2 );
584 calt2( pR, P3, T2 );
585 temp_PDF = 0;
586 for ( int i = 0; i < 4; i++ )
587 {
588 for ( int j = 0; j < 4; j++ ) { temp_PDF += t2[i][j] * T2[j][i] * G[i][i] * G[j][j]; }
589 }
590 }
591 v_re = temp_PDF * B[0] * B[1];
592 return v_re;
593}
594
595void EvtDToKKpi::calEva( double* K01, double* K02, double* Pi, double* mass1, double* width1,
596 double* amp, double* phase, int* g0, int* spin, int* modetype,
597 int nstates, double& Result, double* r0, double* r1 ) {
598
599 double P12[4], P23[4], P13[4];
600 double cof[2], amp_PDF[2], PDF[2];
601 double Amp_KPiS[2];
602 double s1, s2, s3;
603 double s12, s13, s23;
604 for ( int i = 0; i < 4; i++ )
605 {
606 P12[i] = K01[i] + K02[i];
607 P13[i] = K01[i] + Pi[i];
608 P23[i] = K02[i] + Pi[i];
609 }
610 s1 = SCADot( K01, K01 );
611 s2 = SCADot( K02, K02 );
612 s3 = SCADot( Pi, Pi );
613
614 s12 = SCADot( P12, P12 );
615 s13 = SCADot( P13, P13 );
616 s23 = SCADot( P23, P23 );
617
618 double pro[2], temp_PDF, amp_tmp[2];
619 double mass1sq;
620 amp_PDF[0] = 0;
621 amp_PDF[1] = 0;
622 PDF[0] = 0;
623 PDF[1] = 0;
624 amp_tmp[0] = 0;
625 amp_tmp[1] = 0;
626 for ( int i = 0; i < nstates; i++ )
627 {
628 amp_tmp[0] = 0;
629 amp_tmp[1] = 0;
630 mass1sq = mass1[i] * mass1[i];
631 cof[0] = amp[i] * cos( phase[i] );
632 cof[1] = amp[i] * sin( phase[i] );
633 temp_PDF = 0;
634
635 if ( modetype[i] == 12 )
636 {
637 temp_PDF = DDalitz( K01, K02, Pi, spin[i], mass1[i] );
638 if ( g0[i] == 1 ) propagatorRBW( mass1[i], width1[i], s12, s1, s2, rRes2, spin[i], pro );
639 if ( g0[i] == 4 ) propagatorFlatte( mass1[i], width1[i], s12, pro ); // Only for f0(980)
640 amp_tmp[0] = temp_PDF * pro[0];
641 amp_tmp[1] = temp_PDF * pro[1];
642 }
643
644 if ( modetype[i] == 13 )
645 {
646 temp_PDF = DDalitz( K01, Pi, K02, spin[i], mass1[i] );
647 if ( g0[i] == 1 ) propagatorRBW( mass1[i], width1[i], s13, s1, s3, rRes2, spin[i], pro );
648 if ( g0[i] == 3 ) KPiSLASS( s13, s1, s3, pro );
649
650 amp_tmp[0] = temp_PDF * pro[0];
651 amp_tmp[1] = temp_PDF * pro[1];
652 }
653
654 Com_Multi( amp_tmp, cof, amp_PDF );
655 PDF[0] += amp_PDF[0];
656 PDF[1] += amp_PDF[1];
657 }
658 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
659 if ( value <= 0 ) value = 1e-20;
660 Result = value;
661 // printf("%s %8.15f\n","AAABBBvalue = ",value);
662}
mass
double mass
Definition CosmicGenerator.cxx:128
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
g1
TF1 * g1
Definition DataBase/tau_mode.c:58
exp
EvtComplex exp(const EvtComplex &c)
Definition EvtComplex.hh:209
unit
*******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 saves r n generator level $ !Flag for chat level in !Latex Output unit
Definition FoamA.h:90
sin
double sin(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:185
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
q
****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
CS
TCrossPart * CS
Definition Mcgpj.cxx:51
EvtDToKKpi::getName
void getName(std::string &name)
Definition EvtDToKKpi.cc:37
EvtDToKKpi::initProbMax
void initProbMax()
Definition EvtDToKKpi.cc:146
EvtDToKKpi::clone
EvtDecayBase * clone()
Definition EvtDToKKpi.cc:39
EvtDToKKpi::decay
void decay(EvtParticle *p)
Definition EvtDToKKpi.cc:150
EvtDToKKpi::init
void init()
Definition EvtDToKKpi.cc:41
EvtDToKKpi::~EvtDToKKpi
virtual ~EvtDToKKpi()
Definition EvtDToKKpi.cc:35
EvtDecayBase::checkSpinDaughter
void checkSpinDaughter(int d1, EvtSpinType::spintype sp)
Definition EvtDecayBase.cc:492
EvtDecayBase::checkSpinParent
void checkSpinParent(EvtSpinType::spintype sp)
Definition EvtDecayBase.cc:479
EvtDecayBase::setProbMax
void setProbMax(double prbmx)
Definition EvtDecayBase.cc:261
EvtDecayBase::checkNDaug
void checkNDaug(int d1, int d2=-1)
Definition EvtDecayBase.cc:465
EvtDecayBase::getDaugs
EvtId * getDaugs()
Definition EvtDecayBase.hh:68
EvtDecayBase::checkNArg
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
Definition EvtDecayBase.cc:449
EvtDecayProb::setProb
void setProb(double prob)
Definition EvtDecayProb.hh:40
EvtParticle::getP4
const EvtVector4R & getP4() const
Definition EvtParticle.cc:108
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition EvtParticle.cc:80
EvtParticle::initializePhaseSpace
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
Definition EvtParticle.cc:1042
EvtVector4R::get
double get(int i) const
Definition EvtVector4R.hh:171
t
int t()
Definition t.c:1