BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDsToKSKpi0.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: EvtDsToKSKpi0.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 "EvtDsToKSKpi0.hh"
18#include "../EvtGenBase/EvtConst.hh"
19#include "../EvtGenBase/EvtDecayTable.hh"
20#include "../EvtGenBase/EvtGenKine.hh"
21#include "../EvtGenBase/EvtPDL.hh"
22#include "../EvtGenBase/EvtParticle.hh"
23#include "../EvtGenBase/EvtPatches.hh"
24#include "../EvtGenBase/EvtReport.hh"
25#include "../EvtGenBase/EvtResonance.hh"
26#include "../EvtGenBase/EvtResonance2.hh"
27#include <fstream>
28#include <stdlib.h>
29#include <string>
30using std::endl;
31
32EvtDsToKSKpi0::~EvtDsToKSKpi0() {}
33
34void EvtDsToKSKpi0::getName( std::string& model_name ) { model_name = "DsToKSKpi0"; }
35
36EvtDecayBase* EvtDsToKSKpi0::clone() { return new EvtDsToKSKpi0; }
37
38void EvtDsToKSKpi0::init() {
39 // check that there are 0 arguments
40 checkNArg( 0 );
41 checkNDaug( 3 );
42 checkSpinParent( EvtSpinType::SCALAR );
43 checkSpinDaughter( 0, EvtSpinType::SCALAR );
44 checkSpinDaughter( 1, EvtSpinType::SCALAR );
45 checkSpinDaughter( 2, EvtSpinType::SCALAR );
46
47 phi[0] = 0;
48 rho[0] = 0;
49 phi[1] = 0;
50 rho[1] = 1;
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[5] = 0;
58 rho[5] = 0;
59 phi[0] = 5.2793e+00;
60 phi[2] = -6.4411e+00;
61 phi[3] = 2.2286e-01;
62 phi[5] = -2.5771e+00;
63
64 rho[0] = 8.5648e-01;
65 rho[2] = 6.7257e-01;
66 rho[3] = 1.8923e+00;
67 rho[5] = 2.3549e+00;
68
69 modetype[0] = 12;
70 modetype[1] = 13;
71 modetype[2] = 23;
72 modetype[3] = 13;
73 modetype[4] = 23;
74 modetype[5] = 12;
75
76 // cout << "DsToKSKpi0 :" << endl;
77 // for (int i=0; i<6; i++) {
78 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
79 // }
80
81 width[0] = 0.0732;
82 width[1] = 0.0473;
83 width[2] = 0.0503;
84 width[3] = 0.232;
85 width[4] = 0.1767;
86 width[5] = 0.0993;
87
88 mass[0] = 1.0007;
89 mass[1] = 0.89555;
90 mass[2] = 0.89176;
91 mass[3] = 1.414;
92 mass[4] = 1.600;
93 mass[5] = 1.8154;
94
95 mD = 1.86486;
96 mDs = 1.9683;
97 rRes = 9.0;
98 rD = 5.0;
99 metap = 0.95778;
100 mkstr = 0.89594;
101 mk0 = 0.497614;
102 mass_Kaon = 0.49368;
103 mass_Pion = 0.13957;
104 mass_Pi0 = 0.1349766;
105 math_pi = 3.1415926;
106
107 GS1 = 0.636619783;
108 GS2 = 0.01860182466;
109 GS3 = 0.1591549458; // 1/(2*math_2pi)
110 GS4 = 0.00620060822; // mass_Pion2/math_pi
111
112 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
113 for ( int i = 0; i < 4; i++ )
114 {
115 for ( int j = 0; j < 4; j++ ) { G[i][j] = GG[i][j]; }
116 }
117}
118
119void EvtDsToKSKpi0::initProbMax() { setProbMax( 1871.0 ); }
120
121void EvtDsToKSKpi0::decay( EvtParticle* p ) {
122 /*
123 double maxprob = 0.0;
124 for(int ir=0;ir<=60000000;ir++){
125 p->initializePhaseSpace(getNDaug(),getDaugs());
126 EvtVector4R D1 = p->getDaug(0)->getP4();
127 EvtVector4R D2 = p->getDaug(1)->getP4();
128 EvtVector4R D3 = p->getDaug(2)->getP4();
129
130 double P1[4], P2[4], P3[4];
131 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
132 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
133 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
134
135 double value;
136 //value = calDalEva(P1, P2, P3);
137 int g0[6]={1,1,1,1,1,1};
138 int spin[6]={0,1,1,1,1,0};
139
140 int nstates=6;
141 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value);
142 if (value<0) continue;
143 if(value>maxprob) {
144 maxprob=value;
145 cout << "ir= " << ir << endl;
146 cout << "double P1[4] = {" << P1[0] <<","<< P1[1] <<","<< P1[2] <<","<< P1[3]
147 <<"};"<<" "<< sqrt(P1[0]*P1[0]-P1[1]*P1[1]-P1[2]*P1[2]-P1[3]*P1[3]) <<endl; cout <<
148 "double P2[4] = {" << P2[0] <<","<< P2[1]
149 <<","<< P2[2] <<","<< P2[3] <<"};"<< " "<<
150 sqrt(P2[0]*P2[0]-P2[1]*P2[1]-P2[2]*P2[2]-P2[3]*P2[3]) <<endl; cout << "double P3[4] = {"
151 << P3[0] <<","<< P3[1] <<","<< P3[2] <<","<< P3[3] <<"};"<< " "<<
152 sqrt(P3[0]*P3[0]-P3[1]*P3[1]-P3[2]*P3[2]-P3[3]*P3[3]) <<endl; cout << "MAX====> " <<
153 maxprob << endl;
154
155
156 }
157
158 }
159 printf("MAXprob = %.10f\n",maxprob);
160 */
161
162 p->initializePhaseSpace( getNDaug(), getDaugs() );
163 EvtVector4R D1 = p->getDaug( 0 )->getP4();
164 EvtVector4R D2 = p->getDaug( 1 )->getP4();
165 EvtVector4R D3 = p->getDaug( 2 )->getP4();
166
167 double P1[4], P2[4], P3[4];
168 P1[0] = D1.get( 0 );
169 P1[1] = D1.get( 1 );
170 P1[2] = D1.get( 2 );
171 P1[3] = D1.get( 3 );
172 P2[0] = D2.get( 0 );
173 P2[1] = D2.get( 1 );
174 P2[2] = D2.get( 2 );
175 P2[3] = D2.get( 3 );
176 P3[0] = D3.get( 0 );
177 P3[1] = D3.get( 1 );
178 P3[2] = D3.get( 2 );
179 P3[3] = D3.get( 3 );
180
181 // P1[0] = 0.619505234038183; P1[1] = -0.091552038742807; P1[2] = -0.326369642554372; P1[3]
182 // = 0.145835944211028; P2[0] = 0.537869938150145; P2[1] = 0.060139034763020; P2[2] =
183 // -0.082936309719630; P2[3] = -0.187328468377079; P3[0] = 0.875173341522316; P3[1] =
184 // 0.176856301894967; P3[2] = 0.816309543635941; P3[3] = -0.223764071089403;
185 // P1[0] = 0.958942320949657; P1[1] =-0.130886193017869 ; P1[2] = 0.458368295934040; P1[3] =
186 // -0.666871795529042; P2[0] = 0.840895041827643; P2[1] = -0.279686661949455; P2[2] =
187 // -0.173310031669695; P2[3] = 0.595924907258905; P3[0] = 0.205192165542905 ; P3[1] =
188 // 0.148385551896119; P3[2] = -0.008476135905980; P3[3] = 0.042367724031334; P1[0] = ; P1[1]
189 // = ; P1[2] = ; P1[3] = ;
190 // P2[0] = ; P2[1] = ; P2[2] = ; P2[3] = ;
191 // P3[0] = ; P3[1] = ; P3[2] = ; P3[3] = ;
192 // P1[0] =0.61748841695829326248 ; P1[1] =-0.16431886625773378663 ; P1[2] =
193 // -0.015984626326897594106; P1[3] = -0.32623606810397931532; P2[0] =0.61165109530060246534 ;
194 // P2[1] = -0.16009556658150744801; P2[2] = -0.016518497076036600668; P2[3] =
195 // -0.32325478223947284873; P3[0] =0.73916048774110421071 ; P3[1] = 0.32441443283924120689;
196 // P3[2] = 0.032503123402934194774; P3[3] = 0.64949085034345221956;
197
198 double value;
199 int g0[6] = { 1, 1, 1, 1, 1, 1 };
200 int spin[6] = { 0, 1, 1, 1, 1, 0 };
201 int nstates = 6;
202 calEva( P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value );
203
204 setProb( value );
205
206 return;
207}
208
209void EvtDsToKSKpi0::Com_Multi( double a1[2], double a2[2], double res[2] ) {
210 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
211 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
212}
213void EvtDsToKSKpi0::Com_Divide( double a1[2], double a2[2], double res[2] ) {
214 double tmp = a2[0] * a2[0] + a2[1] * a2[1];
215 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / tmp;
216 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / tmp;
217}
218//------------base---------------------------------
219double EvtDsToKSKpi0::SCADot( double a1[4], double a2[4] ) {
220 double _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
221 return _cal;
222}
223
224double EvtDsToKSKpi0::barrier( int l, double sa, double sb, double sc, double r,
225 double mass ) {
226 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
227 if ( q < 0 ) q = 1e-16;
228 double z;
229 z = q * r * r;
230 double sa0;
231 sa0 = mass * mass;
232 double q0 = ( sa0 + sb - sc ) * ( sa0 + sb - sc ) / ( 4 * sa0 ) - sb;
233 if ( q0 < 0 ) q0 = 1e-16;
234 double z0 = q0 * r * r;
235 double F = 0.0;
236 if ( l == 0 ) F = 1;
237 if ( l == 1 ) F = sqrt( ( 1 + z0 ) / ( 1 + z ) );
238 if ( l == 2 ) F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) );
239 return F;
240}
241/*
242double EvtDsToKSKpi0::Barrier(int l, double sa, double sb, double sc, double r2)
243{
244 double F;
245 double tmp = sa+sb-sc;
246 double q = 0.25*tmp*tmp/sa-sb;
247 if (q < 0) q = 1e-16;
248 double z = q*r2;
249 if (l==1) {
250 F = sqrt(2.0*z/(1.0+z));
251 }
252 else if (l==2) {
253 double z2 = z*z;
254 F = sqrt(13.0*z2/(9.0+3.0*z+z2));
255 } else {
256 F = 1.0;
257 }
258 return F;
259}
260*/
261void EvtDsToKSKpi0::calt1( double daug1[4], double daug2[4], double t1[4] ) {
262 double p, pq, tmp;
263 double pa[4], qa[4];
264 for ( int i = 0; i < 4; i++ )
265 {
266 pa[i] = daug1[i] + daug2[i];
267 qa[i] = daug1[i] - daug2[i];
268 }
269 p = SCADot( pa, pa );
270 pq = SCADot( pa, qa );
271 tmp = pq / p;
272 for ( int i = 0; i < 4; i++ ) { t1[i] = qa[i] - tmp * pa[i]; }
273}
274void EvtDsToKSKpi0::calt2( double daug1[4], double daug2[4], double t2[4][4] ) {
275 double p, r;
276 double pa[4], t1[4];
277 calt1( daug1, daug2, t1 );
278 r = SCADot( t1, t1 ) / 3.0;
279 for ( int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
280 p = SCADot( pa, pa );
281 for ( int i = 0; i < 4; i++ )
282 {
283 for ( int j = 0; j < 4; j++ )
284 { t2[i][j] = t1[i] * t1[j] - r * ( G[i][j] - pa[i] * pa[j] / p ); }
285 }
286}
287//-------------------prop--------------------------------------------
288
289double EvtDsToKSKpi0::wid( double mass2, double mass, double sa, double sb, double sc,
290 double r2, int l ) {
291 double widm = 0.;
292 double m = sqrt( sa );
293 double tmp = sb - sc;
294 double tmp1 = sa + tmp;
295 double q = 0.25 * tmp1 * tmp1 / sa - sb;
296 if ( q < 0 ) q = 1e-16;
297 double tmp2 = mass2 + tmp;
298 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
299 if ( q0 < 0 ) q0 = 1e-16;
300 double z = q * r2;
301 double z0 = q0 * r2;
302 double t = q / q0;
303 if ( l == 0 ) { widm = sqrt( t ) * mass / m; }
304 else if ( l == 1 ) { widm = t * sqrt( t ) * mass / m * ( 1 + z0 ) / ( 1 + z ); }
305 else if ( l == 2 )
306 { widm = t * t * sqrt( t ) * mass / m * ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ); }
307 return widm;
308}
309double EvtDsToKSKpi0::widl1( double mass2, double mass, double sa, double sb, double sc,
310 double r2 ) {
311 double widm = 0.;
312 double m = sqrt( sa );
313 double tmp = sb - sc;
314 double tmp1 = sa + tmp;
315 double q = 0.25 * tmp1 * tmp1 / sa - sb;
316 if ( q < 0 ) q = 1e-16;
317 double tmp2 = mass2 + tmp;
318 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
319 if ( q0 < 0 ) q0 = 1e-16;
320 double z = q * r2;
321 double z0 = q0 * r2;
322 double F = ( 1 + z0 ) / ( 1 + z );
323 double t = q / q0;
324 widm = t * sqrt( t ) * mass / m * F;
325 return widm;
326}
327void EvtDsToKSKpi0::propagatorRBW( double mass2, double mass, double width, double sa,
328 double sb, double sc, double r2, int l, double prop[2] ) {
329 double a[2], b[2];
330 a[0] = 1;
331 a[1] = 0;
332 b[0] = mass2 - sa;
333 b[1] = -mass * width * wid( mass2, mass, sa, sb, sc, r2, l );
334 Com_Divide( a, b, prop );
335}
336
337void EvtDsToKSKpi0::propagatorGS( double mass2, double mass, double width, double sa,
338 double sb, double sc, double r2, double prop[2] ) {
339 double a[2], b[2];
340 double tmp = sb - sc;
341 double tmp1 = sa + tmp;
342 double q2 = 0.25 * tmp1 * tmp1 / sa - sb;
343 if ( q2 < 0 ) q2 = 1e-16;
344
345 double tmp2 = mass2 + tmp;
346 double q02 = 0.25 * tmp2 * tmp2 / mass2 - sb;
347 if ( q02 < 0 ) q02 = 1e-16;
348
349 double q = sqrt( q2 );
350 double q0 = sqrt( q02 );
351 double m = sqrt( sa );
352 double q03 = q0 * q02;
353 double tmp3 = log( mass + 2 * q0 ) + 1.2926305904; // log(mpi0+mpip) = -1.2926305904;
354
355 double h = GS1 * q / m * ( log( m + 2 * q ) + 1.2926305904 );
356 double h0 = GS1 * q0 / mass * tmp3;
357 double dh = h0 * ( 0.125 / q02 - 0.5 / mass2 ) + GS3 / mass2;
358 double d = GS2 / q02 * tmp3 + GS3 * mass / q0 - GS4 * mass / q03;
359 double f = mass2 / q03 * ( q2 * ( h - h0 ) + ( mass2 - sa ) * q02 * dh );
360
361 a[0] = 1.0 + d * width / mass;
362 a[1] = 0.0;
363 b[0] = mass2 - sa + width * f;
364 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r2 );
365 Com_Divide( a, b, prop );
366}
367
368void EvtDsToKSKpi0::propagatorFlatte( double mass, double width, double sa, double sb,
369 double sc, int r, double prop[2] ) {
370 double q, qKsK, qetapi;
371 double rhoab[2], rhoKsK[2];
372 q = 0.25 * ( sa + sb - sc ) * ( sa + sb - sc ) / sa - sb;
373 qetapi = 0.25 * ( sa + 0.547862 - 0.13957039 ) * ( sa + 0.547862 - 0.13957039 ) / sa -
374 0.547862 * 0.547862;
375 if ( r == 0 ) qKsK = 0.25 * sa - 0.49368 * 0.49368;
376 if ( r == 1 ) qKsK = 0.25 * ( sa + sb - sc ) * ( sa + sb - sc ) / sa - sb;
377 if ( qetapi > 0 )
378 {
379 rhoab[0] = 2 * sqrt( qetapi / sa );
380 rhoab[1] = 0;
381 }
382 if ( qetapi < 0 )
383 {
384 rhoab[0] = 0;
385 rhoab[1] = 2 * sqrt( -qetapi / sa );
386 }
387 if ( qKsK > 0 )
388 {
389 rhoKsK[0] = 2 * sqrt( qKsK / sa );
390 rhoKsK[1] = 0;
391 }
392 if ( qKsK < 0 )
393 {
394 rhoKsK[0] = 0;
395 rhoKsK[1] = 2 * sqrt( -qKsK / sa );
396 }
397 double a[2], b[2];
398 a[0] = 1;
399 a[1] = 0;
400 b[0] = mass * mass - sa + 0.341 * rhoab[1] + 0.892 * 0.341 * rhoKsK[1];
401 b[1] = -( 0.341 * rhoab[0] + 0.892 * 0.341 * rhoKsK[0] );
402 Com_Divide( a, b, prop );
403}
404
405void EvtDsToKSKpi0::KPiSLASS( double sa, double sb, double sc, double prop[2] ) {
406 const double m1430 = 1.441;
407 const double sa0 = 1.441 * 1.441; // m1430*m1430;
408 const double w1430 = 0.193;
409 const double Lass1 = 0.25 / sa0;
410 double tmp = sb - sc;
411 double tmp1 = sa0 + tmp;
412 double q0 = Lass1 * tmp1 * tmp1 - sb;
413 if ( q0 < 0 ) q0 = 1e-16;
414 double tmp2 = sa + tmp;
415 double qs = 0.25 * tmp2 * tmp2 / sa - sb;
416 double q = sqrt( qs );
417 double width = w1430 * q * m1430 / sqrt( sa * q0 );
418 double temp_R = atan( m1430 * width / ( sa0 - sa ) );
419 if ( temp_R < 0 ) temp_R += math_pi;
420 double deltaR = -109.7 + temp_R;
421 double temp_F =
422 atan( 0.226 * q / ( 2.0 - 3.819 * qs ) ); // 2.0*0.113 = 0.226; 0.113*33.8 = 3.819
423 if ( temp_F < 0 ) temp_F += math_pi;
424 double deltaF = 0.1 + temp_F;
425 double deltaS = deltaR + 2.0 * deltaF;
426 double t1 = 0.96 * sin( deltaF );
427 double t2 = sin( deltaR );
428 double CF[2], CS[2];
429 CF[0] = cos( deltaF );
430 CF[1] = sin( deltaF );
431 CS[0] = cos( deltaS );
432 CS[1] = sin( deltaS );
433 prop[0] = t1 * CF[0] + t2 * CS[0];
434 prop[1] = t1 * CF[1] + t2 * CS[1];
435}
436
437double EvtDsToKSKpi0::DDalitz( double P1[4], double P2[4], double P3[4], int Ang,
438 double mass ) {
439 double pR[4], pD[4];
440 double temp_PDF, v_re;
441 temp_PDF = 0.0;
442 v_re = 0.0;
443 double B[2], s1, s2, s3, sR, sD;
444 for ( int i = 0; i < 4; i++ )
445 {
446 pR[i] = P1[i] + P2[i];
447 pD[i] = pR[i] + P3[i];
448 }
449 s1 = SCADot( P1, P1 );
450 s2 = SCADot( P2, P2 );
451 s3 = SCADot( P3, P3 );
452 sR = SCADot( pR, pR );
453 sD = SCADot( pD, pD );
454 // int G[4][4];
455 // for(int i=0; i!=4; i++){
456 // for(int j=0; j!=4; j++){
457 // if(i==j){
458 // if(i==0) G[i][j] = 1;
459 // else G[i][j] = -1;
460 // }
461 // else G[i][j] = 0;
462 // }
463 // }
464 if ( Ang == 0 )
465 {
466 B[0] = 1;
467 B[1] = 1;
468 temp_PDF = 1;
469 }
470 if ( Ang == 1 )
471 {
472 B[0] = barrier( 1, sR, s1, s2, 3.0, mass );
473 B[1] = barrier( 1, sD, sR, s3, 5.0, 1.9683 );
474 // B[0] = Barrier(1,sR,s1,s2,9.0);
475 // B[1] = Barrier(1,sD,sR,s3,25.0);
476 double t1[4], T1[4];
477 calt1( P1, P2, t1 );
478 calt1( pR, P3, T1 );
479 temp_PDF = 0;
480 for ( int i = 0; i < 4; i++ ) { temp_PDF += t1[i] * T1[i] * G[i][i]; }
481 }
482 if ( Ang == 2 )
483 {
484 B[0] = barrier( 2, sR, s1, s2, 3.0, mass );
485 B[1] = barrier( 2, sD, sR, s3, 5.0, 1.9683 );
486 // B[0] = Barrier(2,sR,s1,s2,9.0);
487 // B[1] = Barrier(2,sD,sR,s3,25.0);
488 double t2[4][4], T2[4][4];
489 calt2( P1, P2, t2 );
490 calt2( pR, P3, T2 );
491 temp_PDF = 0;
492 for ( int i = 0; i < 4; i++ )
493 {
494 for ( int j = 0; j < 4; j++ ) { temp_PDF += t2[i][j] * T2[j][i] * G[i][i] * G[j][j]; }
495 }
496 }
497 v_re = temp_PDF * B[0] * B[1];
498 return v_re;
499}
500
501void EvtDsToKSKpi0::calEva( double* Ks0, double* Kc, double* Pi0, double* mass1,
502 double* width1, double* amp, double* phase, int* g0, int* spin,
503 int* modetype, int nstates, double& Result ) {
504 // double Ks0[4], Pic[4], Pi0[4];
505 // double rRes = 9.0;
506 double P12[4], P23[4], P13[4];
507 double cof[2], amp_PDF[2], PDF[2];
508 double snpi, sck, sks0;
509 double s12, s13, s23;
510 for ( int i = 0; i < 4; i++ )
511 {
512 P12[i] = Kc[i] + Ks0[i];
513 P13[i] = Pi0[i] + Ks0[i];
514 P23[i] = Kc[i] + Pi0[i];
515 }
516 sck = SCADot( Kc, Kc );
517 snpi = SCADot( Pi0, Pi0 );
518 sks0 = SCADot( Ks0, Ks0 );
519 s12 = SCADot( P12, P12 );
520 s13 = SCADot( P13, P13 );
521 s23 = SCADot( P23, P23 );
522 double pro[2], temp_PDF, amp_tmp[2], temp_PDF1, temp_PDF2, pro1[2], pro2[2];
523 double mass1sq;
524 amp_PDF[0] = 0;
525 amp_PDF[1] = 0;
526 PDF[0] = 0;
527 PDF[1] = 0;
528 amp_tmp[0] = 0;
529 amp_tmp[1] = 0;
530 for ( int i = 0; i < nstates; i++ )
531 {
532 amp_tmp[0] = 0;
533 amp_tmp[1] = 0;
534 mass1sq = mass1[i] * mass1[i];
535 cof[0] = amp[i] * cos( phase[i] );
536 cof[1] = amp[i] * sin( phase[i] );
537 temp_PDF = 0;
538
539 if ( modetype[i] == 23 )
540 {
541 temp_PDF = DDalitz( Kc, Pi0, Ks0, spin[i], mass1[i] );
542 if ( g0[i] == 1 )
543 propagatorRBW( mass1sq, mass1[i], width1[i], s23, sck, snpi, rRes, spin[i], pro );
544 if ( g0[i] == 2 ) KPiSLASS( s23, sck, snpi, pro );
545 if ( g0[i] == 0 )
546 {
547 pro[0] = 1;
548 pro[1] = 0;
549 }
550 amp_tmp[0] = temp_PDF * pro[0];
551 amp_tmp[1] = temp_PDF * pro[1];
552 }
553 if ( modetype[i] == 12 )
554 {
555 temp_PDF = DDalitz( Ks0, Kc, Pi0, spin[i], mass1[i] );
556 if ( g0[i] == 1 )
557 propagatorRBW( mass1sq, mass1[i], width1[i], s12, sks0, sck, rRes, spin[i], pro );
558 if ( g0[i] == 2 ) KPiSLASS( s12, sks0, sck, pro );
559 if ( g0[i] == 3 )
560 propagatorFlatte( mass1[i], width1[i], s12, sks0, sck, 1, pro ); // Only for a0(980)
561 if ( g0[i] == 0 )
562 {
563 pro[0] = 1;
564 pro[1] = 0;
565 }
566 amp_tmp[0] = temp_PDF * pro[0];
567 amp_tmp[1] = temp_PDF * pro[1];
568 }
569 if ( modetype[i] == 13 )
570 {
571 temp_PDF = DDalitz( Ks0, Pi0, Kc, spin[i], mass1[i] );
572 if ( g0[i] == 1 )
573 propagatorRBW( mass1sq, mass1[i], width1[i], s13, sks0, snpi, rRes, spin[i], pro );
574 if ( g0[i] == 2 ) KPiSLASS( s13, sks0, snpi, pro );
575 if ( g0[i] == 0 )
576 {
577 pro[0] = 1;
578 pro[1] = 0;
579 }
580 amp_tmp[0] = temp_PDF * pro[0];
581 amp_tmp[1] = temp_PDF * pro[1];
582 }
583 Com_Multi( amp_tmp, cof, amp_PDF );
584 PDF[0] += amp_PDF[0];
585 PDF[1] += amp_PDF[1];
586 }
587 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
588 // cout<<"value = "<<value<<endl;
589 if ( value <= 0 ) value = 1e-20;
590 Result = value;
591 // cout<<"value = "<<value<<" Result = "<<Result<<endl;
592}
mass
double mass
Definition CosmicGenerator.cxx:128
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
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
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
EvtDsToKSKpi0::decay
void decay(EvtParticle *p)
Definition EvtDsToKSKpi0.cc:121
EvtDsToKSKpi0::clone
EvtDecayBase * clone()
Definition EvtDsToKSKpi0.cc:36
EvtDsToKSKpi0::~EvtDsToKSKpi0
virtual ~EvtDsToKSKpi0()
Definition EvtDsToKSKpi0.cc:32
EvtDsToKSKpi0::getName
void getName(std::string &name)
Definition EvtDsToKSKpi0.cc:34
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