BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDToKSKSpi.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: EvtDToKSKSpi.cc
11// the necessary file: EvtDToKSKSpi.cc
12//
13// Description: D+ -> KS0 KS0 pi+,
14//
15//
16// Modification history:
17// Liaoyuan Dong Thu Sep 29 22:06:58 CST 2022 Module created
18//------------------------------------------------------------------------
19//
20#include "EvtDToKSKSpi.hh"
28#include <cmath>
29#include <iostream>
30#include <stdlib.h>
31using namespace std;
32
34
35void EvtDToKSKSpi::getName( std::string& model_name ) { model_name = "DToKSKSpi"; }
36
38
40 checkNArg( 0 );
41 checkNDaug( 3 );
43
44 mass[0] = 0.89167; // KpiSW
45 mass[1] = 0.89167; // K*892
46 width[0] = 0.05140; // KpiSW
47 width[1] = 0.05140; // K*892
48
49 rho[0] = 0.680648291453; // KpiSW
50 rho[1] = 1; // K*892
51 phi[0] = 1.837253259522; // KpiSW
52 phi[1] = 0; // K*892
53
54 spin[0] = 0; // KpiSW
55 spin[1] = 1; // K*892
56 modetype[0] = 13; // KpiSW
57 modetype[1] = 13; // K*892
58
59 /*
60 std::cout << "EvtDToKSKSpi (May 01, 2023) ==> Initialization" << std::endl;
61 for (int i=0; i<2; i++) {
62 cout << i << "rho,phi = " << rho[i] << ", "<< phi[i] << endl;
63 }
64 */
65 math_pi = 3.1415926;
66 GS1 = 0.636619783;
67 GS2 = 0.01860182466;
68 GS3 = 0.1591549458; // 1/(2*math_2pi)
69 GS4 = 0.00620060822; // mass_Pion2/math_pi
70
71 int G[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
72}
73
75 setProbMax( 361.0 ); // max 353.68 2022/11/08
76}
77
78void EvtDToKSKSpi::decay( EvtParticle* p ) { // phsp
79
80 //-------------------------------------------------begin-----------------------------------------
81 /*
82 //检查产生子正确后打开输出最大pdf
83 double maxprob = 0.0;
84 for(int ir=0;ir<=60000000;ir++){
85 p->initializePhaseSpace(getNDaug(),getDaugs());
86 EvtVector4R D1 = p->getDaug(0)->getP4();
87 EvtVector4R D2 = p->getDaug(1)->getP4();
88 EvtVector4R D3 = p->getDaug(2)->getP4();
89
90 double P1[4], P2[4], P3[4];
91 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
92 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
93 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
94
95 double value;
96 // value = calDalEva(P1, P2, P3);
97 int g0[2]={12,1};
98 int spin[2]={0,1};
99 int nstates=2;
100 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value);
101 if (value<0) continue;
102 if(value>maxprob) {
103 maxprob=value;
104 cout << "ir= " << ir << endl;
105 cout << "double P1[4] = {" << P1[0] <<","<< P1[1] <<","<< P1[2] <<","<< P1[3]
106 <<"};"<< endl; cout << "double P2[4] = {" << P2[0] <<","<< P2[1] <<","<< P2[2] <<","<< P2[3]
107 <<"};"<< endl; cout << "double P3[4] = {" << P3[0] <<","<< P3[1] <<","<< P3[2] <<","<< P3[3]
108 <<"};"<< endl; cout << "MAX====> " << maxprob << endl;
109 }
110 }
111 printf("MAXprob = %.10f\n",maxprob);
112 */
113 //---------------------------------------------endls-----------------------------------------------------
114
115 // 检验最大值关闭
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 );
123 P1[1] = D1.get( 1 );
124 P1[2] = D1.get( 2 );
125 P1[3] = D1.get( 3 );
126 P2[0] = D2.get( 0 );
127 P2[1] = D2.get( 1 );
128 P2[2] = D2.get( 2 );
129 P2[3] = D2.get( 3 );
130 P3[0] = D3.get( 0 );
131 P3[1] = D3.get( 1 );
132 P3[2] = D3.get( 2 );
133 P3[3] = D3.get( 3 );
134 // Check
135 // P1[0] = 0.817527;P1[1] = -0.0302783; P1[2] = 0.535352; P1[3] = -0.364866;
136 // P2[0] = 0.660657;P2[1] = 0.132324; P2[2] = -0.368085; P2[3] = 0.18912 ;
137 // P3[0] = 0.40678; P3[1] = 0.0991563 ; P3[2] = -0.281318; P3[3] = 0.238786;
138
139 double value;
140 int g0[2] = { 12, 1 };
141 int spin[2] = { 0, 1 };
142 int nstates = 2;
143 calEva( P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value );
144 setProb( value ); // 检验最大值关闭
145
146 return;
147}
148
149void EvtDToKSKSpi::Com_Multi( double a1[2], double a2[2], double res[2] ) {
150 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
151 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
152}
153void EvtDToKSKSpi::Com_Divide( double a1[2], double a2[2], double res[2] ) {
154 double tmp = a2[0] * a2[0] + a2[1] * a2[1];
155 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / tmp;
156 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / tmp;
157}
158//------------base---------------------------------
159double EvtDToKSKSpi::SCADot( double a1[4], double a2[4] ) {
160 double _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
161 return _cal;
162}
163double EvtDToKSKSpi::barrier( int l, double sa, double sb, double sc, double r, double mass ) {
164 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
165 if ( q < 0 ) q = 1e-16;
166 double z;
167 z = q * r * r;
168 double sa0;
169 sa0 = mass * mass;
170 double q0 = ( sa0 + sb - sc ) * ( sa0 + sb - sc ) / ( 4 * sa0 ) - sb;
171 if ( q0 < 0 ) q0 = 1e-16;
172 double z0 = q0 * r * r;
173 double F = 0.0;
174 if ( l == 0 ) F = 1;
175 if ( l == 1 ) F = sqrt( ( 1 + z0 ) / ( 1 + z ) );
176 if ( l == 2 ) F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) );
177 return F;
178}
179
180void EvtDToKSKSpi::calt1( double daug1[4], double daug2[4], double t1[4] ) {
181 double p, pq, tmp;
182 double pa[4], qa[4];
183 for ( int i = 0; i < 4; i++ )
184 {
185 pa[i] = daug1[i] + daug2[i];
186 qa[i] = daug1[i] - daug2[i];
187 }
188 p = SCADot( pa, pa );
189 pq = SCADot( pa, qa );
190 tmp = pq / p;
191 for ( int i = 0; i < 4; i++ ) { t1[i] = qa[i] - tmp * pa[i]; }
192}
193
194void EvtDToKSKSpi::calt2( double daug1[4], double daug2[4], double t2[4][4] ) {
195 double p, r;
196 double pa[4], t1[4];
197 calt1( daug1, daug2, t1 );
198 r = SCADot( t1, t1 ) / 3.0;
199 for ( int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
200 p = SCADot( pa, pa );
201 for ( int i = 0; i < 4; i++ )
202 {
203 for ( int j = 0; j < 4; j++ )
204 { t2[i][j] = t1[i] * t1[j] - r * ( G[i][j] - pa[i] * pa[j] / p ); }
205 }
206}
207
208double EvtDToKSKSpi::wid( double mass2, double mass, double sa, double sb, double sc,
209 double r2, int l ) {
210 double widm = 0.;
211 double m = sqrt( sa );
212 double tmp = sb - sc;
213 double tmp1 = sa + tmp;
214 double q = 0.25 * tmp1 * tmp1 / sa - sb;
215 if ( q < 0 ) q = 1e-16;
216 double tmp2 = mass2 + tmp;
217 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
218 if ( q0 < 0 ) q0 = 1e-16;
219 double z = q * r2;
220 double z0 = q0 * r2;
221 double t = q / q0;
222 if ( l == 0 ) { widm = sqrt( t ) * mass / m; }
223 else if ( l == 1 ) { widm = t * sqrt( t ) * mass / m * ( 1 + z0 ) / ( 1 + z ); }
224 else if ( l == 2 )
225 { widm = t * t * sqrt( t ) * mass / m * ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ); }
226 return widm;
227}
228double EvtDToKSKSpi::widl1( double mass2, double mass, double sa, double sb, double sc,
229 double r2 ) {
230 double widm = 0.;
231 double m = sqrt( sa );
232 double tmp = sb - sc;
233 double tmp1 = sa + tmp;
234 double q = 0.25 * tmp1 * tmp1 / sa - sb;
235 if ( q < 0 ) q = 1e-16;
236 double tmp2 = mass2 + tmp;
237 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
238 if ( q0 < 0 ) q0 = 1e-16;
239 double z = q * r2;
240 double z0 = q0 * r2;
241 double F = ( 1 + z0 ) / ( 1 + z );
242 double t = q / q0;
243 widm = t * sqrt( t ) * mass / m * F;
244 return widm;
245}
246void EvtDToKSKSpi::propagatorRBW( double mass2, double mass, double width, double sa,
247 double sb, double sc, double r2, int l, double prop[2] ) {
248 double a[2], b[2];
249 a[0] = 1;
250 a[1] = 0;
251 b[0] = mass2 - sa;
252 b[1] = -mass * width * wid( mass2, mass, sa, sb, sc, r2, l );
253 Com_Divide( a, b, prop );
254}
255
256void EvtDToKSKSpi::propagatorGS( double mass2, double mass, double width, double sa, double sb,
257 double sc, double r2, double prop[2] ) {
258 double a[2], b[2];
259 double tmp = sb - sc;
260 double tmp1 = sa + tmp;
261 double q2 = 0.25 * tmp1 * tmp1 / sa - sb;
262 if ( q2 < 0 ) q2 = 1e-16;
263
264 double tmp2 = mass2 + tmp;
265 double q02 = 0.25 * tmp2 * tmp2 / mass2 - sb;
266 if ( q02 < 0 ) q02 = 1e-16;
267
268 double q = sqrt( q2 );
269 double q0 = sqrt( q02 );
270 double m = sqrt( sa );
271 double q03 = q0 * q02;
272 double tmp3 = log( mass + 2 * q0 ) + 1.2926305904; // log(mpi0+mpip) = -1.2926305904;
273
274 double h = GS1 * q / m * ( log( m + 2 * q ) + 1.2926305904 );
275 double h0 = GS1 * q0 / mass * tmp3;
276 double dh = h0 * ( 0.125 / q02 - 0.5 / mass2 ) + GS3 / mass2;
277 double d = GS2 / q02 * tmp3 + GS3 * mass / q0 - GS4 * mass / q03;
278 double f = mass2 / q03 * ( q2 * ( h - h0 ) + ( mass2 - sa ) * q02 * dh );
279
280 a[0] = 1.0 + d * width / mass;
281 a[1] = 0.0;
282 b[0] = mass2 - sa + width * f;
283 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r2 );
284 Com_Divide( a, b, prop );
285}
286
287void EvtDToKSKSpi::KPiSLASS( double sa, double sb, double sc, double prop[2] ) {
288 const double m1430 = 1.441;
289 const double sa0 = 1.441 * 1.441; // m1430*m1430;
290 const double w1430 = 0.193;
291 const double Lass1 = 0.25 / sa0;
292 double tmp = sb - sc;
293 double tmp1 = sa0 + tmp;
294 double q0 = Lass1 * tmp1 * tmp1 - sb;
295 // if(q0<0) q0 = 1e-16;
296 double tmp2 = sa + tmp;
297 double qs = 0.25 * tmp2 * tmp2 / sa - sb;
298 double q = sqrt( qs );
299 double width = w1430 * q * m1430 / sqrt( sa * q0 );
300 double temp_R = atan( m1430 * width / ( sa0 - sa ) );
301 if ( temp_R < 0 ) temp_R += math_pi;
302 double deltaR = -109.7 * math_pi / 180 + temp_R;
303 double temp_F =
304 atan( 0.226 * q / ( 2.0 - 3.819 * qs ) ); // 2.0*0.113 = 0.226; 0.113*33.8 = 3.819
305 if ( temp_F < 0 ) temp_F += math_pi;
306 double deltaF = 0.1 * math_pi / 180 + temp_F;
307 double deltaS = deltaR + 2.0 * deltaF;
308 double t1 = 0.96 * sin( deltaF );
309 double t2 = sin( deltaR );
310 double CF[2], CS[2];
311 CF[0] = cos( deltaF );
312 CF[1] = sin( deltaF );
313 CS[0] = cos( deltaS );
314 CS[1] = sin( deltaS );
315 prop[0] = t1 * CF[0] + t2 * CS[0];
316 prop[1] = t1 * CF[1] + t2 * CS[1];
317}
318
319double EvtDToKSKSpi::DDalitz( double P1[4], double P2[4], double P3[4], int Ang,
320 double mass ) {
321 double pR[4], pD[4];
322 double temp_PDF, v_re;
323 temp_PDF = 0.0;
324 v_re = 0.0;
325 double B[2], s1, s2, s3, sR, sD;
326 for ( int i = 0; i < 4; i++ )
327 {
328 pR[i] = P1[i] + P2[i];
329 pD[i] = pR[i] + P3[i];
330 }
331 s1 = SCADot( P1, P1 );
332 s2 = SCADot( P2, P2 );
333 s3 = SCADot( P3, P3 );
334 sR = SCADot( pR, pR );
335 sD = SCADot( pD, pD );
336
337 int G[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
338
339 if ( Ang == 0 )
340 {
341 B[0] = 1;
342 B[1] = 1;
343 temp_PDF = 1;
344 }
345 if ( Ang == 1 )
346 {
347 B[0] = barrier( 1, sR, s1, s2, 3.0, mass );
348 B[1] = barrier( 1, sD, sR, s3, 5.0, 1.86966 );
349 double t1[4], T1[4];
350 calt1( P1, P2, t1 );
351 calt1( pR, P3, T1 );
352 temp_PDF = 0;
353 for ( int i = 0; i < 4; i++ ) { temp_PDF += t1[i] * T1[i] * G[i][i]; }
354 }
355 if ( Ang == 2 )
356 {
357 B[0] = barrier( 2, sR, s1, s2, 3.0, mass );
358 B[1] = barrier( 2, sD, sR, s3, 5.0, 1.86966 );
359 double t2[4][4], T2[4][4];
360 calt2( P1, P2, t2 );
361 calt2( pR, P3, T2 );
362 temp_PDF = 0;
363 for ( int i = 0; i < 4; i++ )
364 {
365 for ( int j = 0; j < 4; j++ ) { temp_PDF += t2[i][j] * T2[j][i] * G[i][i] * G[j][j]; }
366 }
367 }
368 v_re = temp_PDF * B[0] * B[1];
369 return v_re;
370}
371
372void EvtDToKSKSpi::calEva( double* Ks01, double* Ks02, double* Pic, double* mass1,
373 double* width1, double* amp, double* phase, int* g0, int* spin,
374 int* modetype, int nstates, double& Result ) {
375 double rRes = 9.0;
376 double P12[4], P23[4], P13[4];
377 double cof[2], amp_PDF[2], PDF[2];
378 double scpi, sks02, sks01;
379 double s12, s13, s23;
380 for ( int i = 0; i < 4; i++ )
381 {
382 P12[i] = Ks01[i] + Ks02[i];
383 P13[i] = Ks01[i] + Pic[i];
384 P23[i] = Ks02[i] + Pic[i];
385 }
386 sks01 = SCADot( Ks01, Ks01 );
387 sks02 = SCADot( Ks02, Ks02 );
388 scpi = SCADot( Pic, Pic );
389 s12 = SCADot( P12, P12 );
390 s13 = SCADot( P13, P13 );
391 s23 = SCADot( P23, P23 );
392 double pro[2], temp_PDF, amp_tmp[2], temp_PDF1, temp_PDF2, pro1[2], pro2[2];
393 double mass1sq;
394 amp_PDF[0] = 0;
395 amp_PDF[1] = 0;
396 PDF[0] = 0;
397 PDF[1] = 0;
398 amp_tmp[0] = 0;
399 amp_tmp[1] = 0;
400 for ( int i = 0; i < nstates; i++ )
401 {
402 amp_tmp[0] = 0;
403 amp_tmp[1] = 0;
404 mass1sq = mass1[i] * mass1[i];
405 cof[0] = amp[i] * cos( phase[i] );
406 cof[1] = amp[i] * sin( phase[i] );
407 temp_PDF = 0;
408 if ( modetype[i] == 12 )
409 {
410 temp_PDF = DDalitz( Ks01, Ks02, Pic, spin[i], mass1[i] );
411 if ( g0[i] == 1 )
412 propagatorRBW( mass1sq, mass1[i], width1[i], s12, sks01, sks02, rRes, spin[i], pro );
413 if ( g0[i] == 12 ) KPiSLASS( s12, sks01, sks02, pro );
414 if ( g0[i] == 0 )
415 {
416 pro[0] = 1;
417 pro[1] = 0;
418 }
419 amp_tmp[0] = temp_PDF * pro[0];
420 amp_tmp[1] = temp_PDF * pro[1];
421 }
422 if ( modetype[i] == 13 )
423 {
424 temp_PDF1 = DDalitz( Ks01, Pic, Ks02, spin[i], mass1[i] );
425 if ( g0[i] == 1 )
426 propagatorRBW( mass1sq, mass1[i], width1[i], s13, sks01, scpi, rRes, spin[i], pro1 );
427 if ( g0[i] == 12 ) KPiSLASS( s13, sks01, scpi, pro1 );
428 if ( g0[i] == 0 )
429 {
430 pro1[0] = 1;
431 pro1[1] = 0;
432 }
433 temp_PDF2 = DDalitz( Ks02, Pic, Ks01, spin[i], mass1[i] );
434 if ( g0[i] == 1 )
435 propagatorRBW( mass1sq, mass1[i], width1[i], s23, sks02, scpi, rRes, spin[i], pro2 );
436 if ( g0[i] == 12 ) KPiSLASS( s23, sks02, scpi, pro2 );
437 if ( g0[i] == 0 )
438 {
439 pro2[0] = 1;
440 pro2[1] = 0;
441 }
442 amp_tmp[0] = ( temp_PDF1 * pro1[0] + temp_PDF2 * pro2[0] );
443 amp_tmp[1] = ( temp_PDF1 * pro1[1] + temp_PDF2 * pro2[1] );
444 }
445 Com_Multi( amp_tmp, cof, amp_PDF );
446 PDF[0] += amp_PDF[0];
447 PDF[1] += amp_PDF[1];
448 }
449 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
450 if ( value <= 0 ) value = 1e-20;
451 Result = value;
452 // Check
453 // printf("%s %8.15f\n","AAAvalue = ",value);
454}
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 getName(std::string &name)
void decay(EvtParticle *p)
virtual ~EvtDToKSKSpi()
EvtDecayBase * clone()
void initProbMax()
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