BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtD0ToKSpi0eta.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: EvtD0ToKSpi0eta.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 "EvtD0ToKSpi0eta.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
35EvtD0ToKSpi0eta::~EvtD0ToKSpi0eta() {}
36
37void EvtD0ToKSpi0eta::getName( std::string& model_name ) { model_name = "D0ToKSpi0eta"; }
38
39EvtDecayBase* EvtD0ToKSpi0eta::clone() { return new EvtD0ToKSpi0eta; }
40
41void EvtD0ToKSpi0eta::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.98; // a0980
51 mass[1] = 0.89255; // K892
52 mass[2] = 1.414; // K1410
53 mass[3] = 1.414; // KPISwave
54 mass[4] = 1.414; // non-Dwave
55
56 width[0] = 0.756; // a0980
57 width[1] = 0.0473; // k892
58 width[2] = 0.232; // k1410
59 width[3] = 0.232; // kpiswave
60 width[4] = 0.232; // non-Dwave
61
62 rho[0] = 1; // a0980
63 phi[0] = 0; //
64
65 rho[1] = 0.20146; // k892
66 phi[1] = 2.9328; //
67
68 rho[2] = 0.81611; // k1410
69 phi[2] = -2.0637;
70
71 rho[3] = 0.49421; // kpiswave
72 phi[3] = -4.4523;
73
74 rho[4] = 5.4953; // non-Dwave
75 phi[4] = 1.7266; //
76
77 spin[0] = 0; // a0980
78 spin[1] = 1; // k892
79 spin[2] = 1; // k1410
80 spin[3] = 0; // kpiSwave
81 spin[4] = 2; // nonD
82
83 modetype[0] = 23; // a0980
84 modetype[1] = 12; // k892
85 modetype[2] = 12; // k1410
86 modetype[3] = 12; // kpis
87 modetype[4] = 13; // non D
88
89 /*
90 std::cout << "EvtD0ToKSpi0eta ==> Initialization" << std::endl;
91 for (int i=0; i<5; i++) {
92 cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
93 }
94 */
95
96 mD0 = 1.864;
97 mK0 = 0.497611;
98 mKa = 0.49368;
99 mPi = 0.13957;
100 mK02 = 0.237616707;
101 mPi2 = 0.01947978;
102 mass_EtaP = 0.95778;
103 mass_Kaon = 0.49368;
104
105 math_pi = 3.1415926;
106 mass_Pion = 0.13957;
107 mass_Pion2 = 0.0194797849;
108 mass_2Pion = 0.27914;
109 math_2pi = 6.2831852;
110 rD2 = 25.0; // 5*5
111 rRes2 = 9.0; // 3*3
112 // g1 = 0.5468;
113 g2 = 0.23;
114 GS1 = 0.636619783;
115 GS2 = 0.01860182466;
116 GS3 = 0.1591549458;
117 GS4 = 0.00620060822;
118
119 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
120 for ( int i = 0; i < 4; i++ )
121 {
122 for ( int j = 0; j < 4; j++ ) { G[i][j] = GG[i][j]; }
123 }
124}
125
126void EvtD0ToKSpi0eta::initProbMax() {
127 setProbMax( 625.0 ); // MAXprob = 621.9269871585
128}
129
130void EvtD0ToKSpi0eta::decay( EvtParticle* p ) {
131 //==================================================get
132 // max============================================
133 /*
134 double maxprob = 0.0;
135 for(int ir=0;ir<=60000000;ir++){
136 p->initializePhaseSpace(getNDaug(),getDaugs());
137 EvtVector4R D1 = p->getDaug(0)->getP4();
138 EvtVector4R D2 = p->getDaug(1)->getP4();
139 EvtVector4R D3 = p->getDaug(2)->getP4();
140
141 double P1[4], P2[4], P3[4];
142 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
143 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
144 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
145
146 double value;
147 int nstates=8;
148 int g0[8]={1,1,4,1,1,1,3,1};
149 double r0[8]={3.0,3.0,3.0,3.0,3.0,3.0,3.0,3.0};
150 double r1[8]={5.0,5.0,5.0,5.0,5.0,5.0,5.0,5.0};
151
152 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value,r0,r1);
153
154 if (value<0) continue;
155 if(value>maxprob) {
156 maxprob=value;
157 cout << "ir= " << ir << endl;
158 cout << "MAX====> " << maxprob << endl;
159 }
160 }
161 printf("MAXprob = %.10f\n",maxprob);
162 */
163 //==================================================get
164 // max============================================
165
166 p->initializePhaseSpace( getNDaug(), getDaugs() );
167 EvtVector4R D1 = p->getDaug( 0 )->getP4();
168 EvtVector4R D2 = p->getDaug( 1 )->getP4();
169 EvtVector4R D3 = p->getDaug( 2 )->getP4();
170
171 double P1[4], P2[4], P3[4];
172 P1[0] = D1.get( 0 );
173 P1[1] = D1.get( 1 );
174 P1[2] = D1.get( 2 );
175 P1[3] = D1.get( 3 );
176 P2[0] = D2.get( 0 );
177 P2[1] = D2.get( 1 );
178 P2[2] = D2.get( 2 );
179 P2[3] = D2.get( 3 );
180 P3[0] = D3.get( 0 );
181 P3[1] = D3.get( 1 );
182 P3[2] = D3.get( 2 );
183 P3[3] = D3.get( 3 );
184
185 // Check
186 // P1[0] = 0.734098420544862; P2[0] = 0.352732515142708; P3[0] =
187 // 0.799420828049644;
188 // P1[1] = 0.409477641275596; P2[1] = 0.078790398165110; P3[1] = -0.431210425541509;
189 // P1[2] = 0.301451605847801; P2[2] = 0.308201450489482; P3[2] = -0.378552010532733;
190 // P1[3] =-0.180693640379551; P2[3] =-0.070704784074596; P3[3] = 0.09759914740263;
191
192 double value;
193 int nstates = 5;
194 int g0[5] = { 3, 1, 1, 2, 0 };
195 double r0[5] = { 3.0, 3.0, 3.0, 3.0, 3.0 };
196 double r1[5] = { 5.0, 5.0, 5.0, 5.0, 5.0 };
197
198 calEva( P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value, r0, r1 );
199 setProb( value );
200 return;
201
202 //////
203}
204
205void EvtD0ToKSpi0eta::Com_Multi( double a1[2], double a2[2], double res[2] ) {
206 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
207 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
208}
209void EvtD0ToKSpi0eta::Com_Divide( double a1[2], double a2[2], double res[2] ) {
210 double tmp = a2[0] * a2[0] + a2[1] * a2[1];
211 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / tmp;
212 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / tmp;
213}
214//------------base---------------------------------
215double EvtD0ToKSpi0eta::SCADot( double a1[4], double a2[4] ) {
216 double _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
217 return _cal;
218}
219double EvtD0ToKSpi0eta::barrier( int l, double sa, double sb, double sc, double r,
220 double mass ) {
221 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
222 if ( q < 0 ) q = 1e-16;
223 double z;
224 z = q * r * r;
225 double sa0;
226 sa0 = mass * mass;
227 double q0 = ( sa0 + sb - sc ) * ( sa0 + sb - sc ) / ( 4 * sa0 ) - sb;
228 if ( q0 < 0 ) q0 = 1e-16;
229 double z0 = q0 * r * r;
230 double F = 0.0;
231 if ( l == 0 ) F = 1;
232 if ( l == 1 ) F = sqrt( ( 1 + z0 ) / ( 1 + z ) );
233 if ( l == 2 ) F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) );
234 return F;
235}
236//------------------spin-------------------------------------------
237void EvtD0ToKSpi0eta::calt1( double daug1[4], double daug2[4], double t1[4] ) {
238 double p, pq, tmp;
239 double pa[4], qa[4];
240 for ( int i = 0; i < 4; i++ )
241 {
242 pa[i] = daug1[i] + daug2[i];
243 qa[i] = daug1[i] - daug2[i];
244 }
245 p = SCADot( pa, pa );
246 pq = SCADot( pa, qa );
247 tmp = pq / p;
248 for ( int i = 0; i < 4; i++ ) { t1[i] = qa[i] - tmp * pa[i]; }
249}
250void EvtD0ToKSpi0eta::calt2( double daug1[4], double daug2[4], double t2[4][4] ) {
251 double p, r;
252 double pa[4], t1[4];
253 calt1( daug1, daug2, t1 );
254 r = SCADot( t1, t1 ) / 3.0;
255 for ( int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
256 p = SCADot( pa, pa );
257 for ( int i = 0; i < 4; i++ )
258 {
259 for ( int j = 0; j < 4; j++ )
260 { t2[i][j] = t1[i] * t1[j] - r * ( G[i][j] - pa[i] * pa[j] / p ); }
261 }
262}
263//-------------------prop--------------------------------------------
264void EvtD0ToKSpi0eta::propagator( double mass2, double mass, double width, double sx,
265 double prop[2] ) {
266 double a[2], b[2];
267 a[0] = 1;
268 a[1] = 0;
269 b[0] = mass2 - sx;
270 b[1] = -mass * width;
271 Com_Divide( a, b, prop );
272}
273double EvtD0ToKSpi0eta::wid( double mass2, double mass, double sa, double sb, double sc,
274 double r2, int l ) {
275 double widm = 0.;
276 double m = sqrt( sa );
277 double tmp = sb - sc;
278 double tmp1 = sa + tmp;
279 double q = 0.25 * tmp1 * tmp1 / sa - sb;
280 if ( q < 0 ) q = 1e-16;
281 double tmp2 = mass2 + tmp;
282 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
283 if ( q0 < 0 ) q0 = 1e-16;
284 double z = q * r2;
285 double z0 = q0 * r2;
286 double t = q / q0;
287 if ( l == 0 ) { widm = sqrt( t ) * mass / m; }
288 else if ( l == 1 ) { widm = t * sqrt( t ) * mass / m * ( 1 + z0 ) / ( 1 + z ); }
289 else if ( l == 2 )
290 { widm = t * t * sqrt( t ) * mass / m * ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ); }
291 return widm;
292}
293double EvtD0ToKSpi0eta::widl1( double mass2, double mass, double sa, double sb, double sc,
294 double r2 ) {
295 double widm = 0.;
296 double m = sqrt( sa );
297 double tmp = sb - sc;
298 double tmp1 = sa + tmp;
299 double q = 0.25 * tmp1 * tmp1 / sa - sb;
300 if ( q < 0 ) q = 1e-16;
301 double tmp2 = mass2 + tmp;
302 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
303 if ( q0 < 0 ) q0 = 1e-16;
304 double z = q * r2;
305 double z0 = q0 * r2;
306 double F = ( 1 + z0 ) / ( 1 + z );
307 double t = q / q0;
308 widm = t * sqrt( t ) * mass / m * F;
309 return widm;
310}
311void EvtD0ToKSpi0eta::propagatorRBW( double mass, double width, double sa, double sb,
312 double sc, double r2, int l, double prop[2] ) {
313
314 double a[2], b[2];
315 double mass2 = mass * mass;
316 a[0] = 1;
317 a[1] = 0;
318 b[0] = mass2 - sa;
319 b[1] = -mass * width * wid( mass2, mass, sa, sb, sc, r2, l );
320 Com_Divide( a, b, prop );
321}
322
323void EvtD0ToKSpi0eta::KPiSLASS( double sa, double sb, double sc, double prop[2] ) {
324 const double m1430 = 1.441;
325 const double sa0 = 1.441 * 1.441; // m1430*m1430;
326 const double w1430 = 0.193;
327 const double Lass1 = 0.25 / sa0;
328 double tmp = sb - sc;
329 double tmp1 = sa0 + tmp;
330 double q0 = Lass1 * tmp1 * tmp1 - sb;
331 if ( q0 < 0 ) q0 = 1e-16;
332 double tmp2 = sa + tmp;
333 double qs = 0.25 * tmp2 * tmp2 / sa - sb;
334 double q = sqrt( qs );
335 double width = w1430 * q * m1430 / sqrt( sa * q0 );
336 double temp_R = atan( m1430 * width / ( sa0 - sa ) );
337 if ( temp_R < 0 ) temp_R += math_pi;
338 double deltaR = -1.915 + temp_R; // fiR=-109.7
339 double temp_F =
340 atan( 0.226 * q / ( 2.0 - 3.819 * qs ) ); // 2.0*0.113 = 0.226; 0.113*33.8 = 3.819
341 if ( temp_F < 0 ) temp_F += math_pi;
342 double deltaF = 0.002 + temp_F; // fiF=0.1
343 double deltaS = deltaR + 2.0 * deltaF;
344 double t1 = 0.96 * sin( deltaF );
345 double t2 = sin( deltaR );
346 double CF[2], CS[2];
347 CF[0] = cos( deltaF );
348 CF[1] = sin( deltaF );
349 CS[0] = cos( deltaS );
350 CS[1] = sin( deltaS );
351 prop[0] = t1 * CF[0] + t2 * CS[0];
352 prop[1] = t1 * CF[1] + t2 * CS[1];
353}
354//------------GS---used by rho----------------------------
355void EvtD0ToKSpi0eta::propagatorGS( double mass2, double mass, double width, double sa,
356 double sb, double sc, double r2, double prop[2] ) {
357 double a[2], b[2];
358 double tmp = sb - sc;
359 double tmp1 = sa + tmp;
360 double q2 = 0.25 * tmp1 * tmp1 / sa - sb;
361 if ( q2 < 0 ) q2 = 1e-16;
362
363 double tmp2 = mass2 + tmp;
364 double q02 = 0.25 * tmp2 * tmp2 / mass2 - sb;
365 if ( q02 < 0 ) q02 = 1e-16;
366
367 double q = sqrt( q2 );
368 double q0 = sqrt( q02 );
369 double m = sqrt( sa );
370 double q03 = q0 * q02;
371 // double tmp3 = log(mass+2*q0)+1.2760418309; // log(mass_2Pion) = 1.2760418309;
372 double tmp3 = log( mass + 2 * q0 ) + 1.2926305904; // log(mpi0+mpip) = -1.2926305904;
373
374 double h = GS1 * q / m * ( log( m + 2 * q ) + 1.2926305904 );
375 double h0 = GS1 * q0 / mass * tmp3;
376 double dh = h0 * ( 0.125 / q02 - 0.5 / mass2 ) + GS3 / mass2;
377 double d = GS2 / q02 * tmp3 + GS3 * mass / q0 - GS4 * mass / q03;
378 double f = mass2 / q03 * ( q2 * ( h - h0 ) + ( mass2 - sa ) * q02 * dh );
379
380 a[0] = 1.0 + d * width / mass;
381 a[1] = 0.0;
382 b[0] = mass2 - sa + width * f;
383 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r2 );
384 Com_Divide( a, b, prop );
385}
386void EvtD0ToKSpi0eta::propagatorFlatte( double mass, double width, double sa, double sb,
387 double sc, int r, double prop[2] ) {
388 double qKK, qetapi, qKSKS;
389 double rhoab[2], rhoKsK[2];
390 // q = 0.25*(sa+sb-sc)*(sa+sb-sc)/sa-sb;
391 qetapi = 0.25 * ( sa + 0.547862 * 0.547862 - 0.1349768 * 0.1349768 ) *
392 ( sa + 0.547862 * 0.547862 - 0.1349768 * 0.1349768 ) / sa -
393 0.547862 * 0.547862; // eta pi0
394 qKK = 1 - ( 2 * 0.49368 / sa ) * ( 2 * 0.49368 / sa ); // K+K-
395 qKSKS = 1 - ( 2 * 0.497614 / sa ) * ( 2 * 0.497614 / sa ); // KS0KS0
396 // eta pi0
397 if ( qetapi > 0 )
398 {
399 rhoab[0] = 2 * sqrt( qetapi / sa );
400 rhoab[1] = 0;
401 }
402 if ( qetapi < 0 )
403 {
404 rhoab[0] = 0;
405 rhoab[1] = 2 * sqrt( -qetapi / sa );
406 }
407 // KK&&KSKS
408 if ( qKK > 0 && qKSKS > 0 )
409 {
410 rhoKsK[0] = 0.5 * ( sqrt( qKK ) + sqrt( qKSKS ) );
411 rhoKsK[1] = 0;
412 }
413 if ( qKK < 0 && qKSKS < 0 )
414 {
415 rhoKsK[0] = 0;
416 rhoKsK[1] = 0.5 * ( sqrt( -qKK ) + sqrt( -qKSKS ) );
417 }
418 if ( qKK < 0 && qKSKS > 0 )
419 {
420 rhoKsK[0] = 0.5 * sqrt( qKSKS );
421 rhoKsK[1] = 0.5 * sqrt( -qKK );
422 }
423 if ( qKK > 0 && qKSKS < 0 )
424 {
425 rhoKsK[0] = 0.5 * sqrt( qKK );
426 rhoKsK[1] = 0.5 * sqrt( -qKSKS );
427 }
428 double a[2], b[2];
429 a[0] = 1;
430 a[1] = 0;
431 b[0] = mass * mass - sa + 0.341 * rhoab[1] +
432 0.892 * 0.341 * rhoKsK[1]; // （考虑两项:eta pi 和kk）
433 b[1] = -( 0.341 * rhoab[0] + 0.892 * 0.341 * rhoKsK[0] );
434 Com_Divide( a, b, prop );
435}
436
437double EvtD0ToKSpi0eta::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 {
457 for ( int j = 0; j != 4; j++ )
458 {
459 if ( i == j )
460 {
461 if ( i == 0 ) G[i][j] = 1;
462 else G[i][j] = -1;
463 }
464 else G[i][j] = 0;
465 }
466 }
467 if ( Ang == 0 )
468 {
469 B[0] = 1;
470 B[1] = 1;
471 temp_PDF = 1;
472 }
473 if ( Ang == 1 )
474 {
475 B[0] = barrier( 1, sR, s1, s2, 3.0, mass );
476 B[1] = barrier( 1, sD, sR, s3, 5.0, 1.864 );
477 double t1[4], T1[4];
478 calt1( P1, P2, t1 );
479 calt1( pR, P3, T1 );
480 temp_PDF = 0;
481 for ( int i = 0; i < 4; i++ ) { temp_PDF += t1[i] * T1[i] * G[i][i]; }
482 }
483 if ( Ang == 2 )
484 {
485 B[0] = barrier( 2, sR, s1, s2, 3.0, mass );
486 B[1] = barrier( 2, sD, sR, s3, 5.0, 1.864 );
487 double t2[4][4], T2[4][4];
488 calt2( P1, P2, t2 );
489 calt2( pR, P3, T2 );
490 temp_PDF = 0;
491 for ( int i = 0; i < 4; i++ )
492 {
493 for ( int j = 0; j < 4; j++ ) { temp_PDF += t2[i][j] * T2[j][i] * G[i][i] * G[j][j]; }
494 }
495 }
496 v_re = temp_PDF * B[0] * B[1];
497 return v_re;
498}
499
500void EvtD0ToKSpi0eta::calEva( double* KS0, double* Pi0, double* Eta, double* mass1,
501 double* width1, double* amp, double* phase, int* g0, int* spin,
502 int* modetype, int nstates, double& Result, double* r0,
503 double* r1 ) {
504 double rRes = 9.0;
505 double P12[4], P23[4], P13[4];
506 double cof[2], amp_PDF[2], PDF[2];
507 double seta, spi0, sks0;
508 double s12, s13, s23;
509 for ( int i = 0; i < 4; i++ )
510 {
511 P12[i] = KS0[i] + Pi0[i];
512 P13[i] = KS0[i] + Eta[i];
513 P23[i] = Pi0[i] + Eta[i];
514 }
515 sks0 = SCADot( KS0, KS0 );
516 spi0 = SCADot( Pi0, Pi0 );
517 seta = SCADot( Eta, Eta );
518 s12 = SCADot( P12, P12 );
519 s13 = SCADot( P13, P13 );
520 s23 = SCADot( P23, P23 );
521 double pro[2], temp_PDF, amp_tmp[2];
522 double mass1sq;
523 double Amp_KPiS[2];
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 if ( modetype[i] == 23 )
539 {
540 temp_PDF = DDalitz( Pi0, Eta, KS0, spin[i], mass1[i] ); // 计算An里面的Sn
541 // if(g0[i]==1) propagatorGS(mass1sq,mass1[i],width1[i],s23,spi0,seta,9.0,pro);
542 if ( g0[i] == 1 )
543 propagatorRBW( mass1[i], width1[i], s23, spi0, seta, rRes, spin[i],
544 pro ); // 利用RBW函数计算了传播子对应An中的Pn
545 if ( g0[i] == 2 )
546 KPiSLASS( s23, spi0, seta, pro ); // 利用K pi s波（不同于传播子的参数化的一种方式），
547 if ( g0[i] == 3 )
548 propagatorFlatte( mass1[i], width1[i], s23, spi0, seta, 0, pro ); // spi0:pi0 seta:eta
549 if ( g0[i] == 0 )
550 {
551 pro[0] = 1;
552 pro[1] = 0;
553 }
554 amp_tmp[0] = temp_PDF * pro[0];
555 amp_tmp[1] = temp_PDF * pro[1];
556 }
557
558 if ( modetype[i] == 12 )
559 {
560 temp_PDF = DDalitz( KS0, Pi0, Eta, spin[i], mass1[i] );
561 if ( g0[i] == 1 )
562 propagatorRBW( mass1[i], width1[i], s12, sks0, spi0, rRes, spin[i], pro );
563 if ( g0[i] == 2 ) KPiSLASS( s12, sks0, spi0, pro );
564
565 if ( g0[i] == 0 )
566 {
567 pro[0] = 1; // 指的是pro的实部为1
568 pro[1] = 0; // 虚部为0------------------------pro0==1&&pro1==0即为非共振态
569 }
570 amp_tmp[0] = temp_PDF * pro[0];
571 amp_tmp[1] = temp_PDF * pro[1];
572 }
573
574 if ( modetype[i] == 13 )
575 {
576 // temp_PDF = 0;
577 temp_PDF = DDalitz( KS0, Eta, Pi0, spin[i], mass1[i] );
578 if ( g0[i] == 1 )
579 propagatorRBW( mass1[i], width1[i], s13, sks0, seta, rRes, spin[i], pro );
580 // if(g0[i]==2) KEtaSLASS(s13,sks0,seta,pro);
581
582 if ( g0[i] == 0 )
583 {
584 pro[0] = 1;
585 pro[1] = 0;
586 }
587 // amp_tmp[0] = amp_tmp[0] - temp_PDF*pro[0];
588 // amp_tmp[1] = amp_tmp[1] - temp_PDF*pro[1];
589 amp_tmp[0] = temp_PDF * pro[0];
590 amp_tmp[1] = temp_PDF * pro[1];
591 }
592 Com_Multi( amp_tmp, cof, amp_PDF );
593 PDF[0] += amp_PDF[0];
594 PDF[1] += amp_PDF[1];
595 }
596 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
597 if ( value <= 0 ) value = 1e-20;
598 Result = value;
599 // Check
600 // printf("%s %8.14f\n","AAABBBvalue = ",value);
601}
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
EvtD0ToKSpi0eta::~EvtD0ToKSpi0eta
virtual ~EvtD0ToKSpi0eta()
Definition EvtD0ToKSpi0eta.cc:35
EvtD0ToKSpi0eta::getName
void getName(std::string &name)
Definition EvtD0ToKSpi0eta.cc:37
EvtD0ToKSpi0eta::clone
EvtDecayBase * clone()
Definition EvtD0ToKSpi0eta.cc:39
EvtD0ToKSpi0eta::decay
void decay(EvtParticle *p)
Definition EvtD0ToKSpi0eta.cc:130
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