BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDToKSpieta.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: EvtDToKSpieta.cc
11// the necessary file: EvtDToKSpieta.hh
12//
13// Description: D+ -> KS0 pi+ eta, see BAM-614
14//
15// Modification history:
16//
17// Liaoyuan Dong Nov. 15, 2022 Module created
18//------------------------------------------------------------------------
19#include "EvtDToKSpieta.hh"
20#include "../EvtGenBase/EvtComplex.hh"
21#include "../EvtGenBase/EvtDecayTable.hh"
22#include "../EvtGenBase/EvtGenKine.hh"
23#include "../EvtGenBase/EvtPDL.hh"
24#include "../EvtGenBase/EvtParticle.hh"
25#include "../EvtGenBase/EvtPatches.hh"
26#include "../EvtGenBase/EvtReport.hh"
27#include <cmath>
28#include <iostream>
29#include <stdlib.h>
30using namespace std;
31
32EvtDToKSpieta::~EvtDToKSpieta() {}
33
34void EvtDToKSpieta::getName( std::string& model_name ) { model_name = "DToKSpieta"; }
35
36EvtDecayBase* EvtDToKSpieta::clone() { return new EvtDToKSpieta; }
37
38void EvtDToKSpieta::init() {
39 checkNArg( 0 );
40 checkNDaug( 3 );
41 checkSpinParent( EvtSpinType::SCALAR );
42 ma0 = 1.004;
43 Ga0 = 0.0732;
44 mKn_1430 = 1.423;
45 GKn_1430 = 0.1815;
46 mK1270 = 1.272;
47 mK1400 = 1.403;
48 GK1270 = 0.09;
49 GK1400 = 0.174;
50 phi_omega = -0.02;
51
52 rho_omega = 0.00294;
53
54 rho[0] = 1; // a0980-pieta
55 phi[0] = 0;
56
57 rho[1] = 0.30586; // K1430-kseta
58 phi[1] = 2.5805;
59
60 spin[0] = 0;
61 spin[1] = 0;
62
63 modetype[0] = 23;
64 modetype[1] = 13;
65
66 /*
67 cout << "Initializing EvtDToKSpieta" << endl;
68 for (int i=0; i<2; i++) {
69 cout << i << " rho,phi = " << rho[i] << ", "<< phi[i] << endl;
70 }
71 */
72
73 mass_Ks = 0.4977;
74 mass_Eta = 0.547862;
75 mD = 1.86965;
76 rD = 5;
77 metap = 0.95778;
78 mkstr = 0.89594;
79 mk0 = 0.497611;
80 mass_Kaon = 0.493677;
81 mass_Pion = 0.13957;
82 mass_Pi0 = 0.1349768;
83 math_pi = 3.1415926;
84
85 pi = 3.1415926;
86 mpi = 0.13957;
87 g1 = 0.5468;
88 g2 = 0.23;
89
90 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
91 int EE[4][4][4][4] = {
92 { { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
93 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 1 }, { 0, 0, -1, 0 } },
94 { { 0, 0, 0, 0 }, { 0, 0, 0, -1 }, { 0, 0, 0, 0 }, { 0, 1, 0, 0 } },
95 { { 0, 0, 0, 0 }, { 0, 0, 1, 0 }, { 0, -1, 0, 0 }, { 0, 0, 0, 0 } } },
96 { { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, -1 }, { 0, 0, 1, 0 } },
97 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
98 { { 0, 0, 0, 1 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { -1, 0, 0, 0 } },
99 { { 0, 0, -1, 0 }, { 0, 0, 0, 0 }, { 1, 0, 0, 0 }, { 0, 0, 0, 0 } } },
100 { { { 0, 0, 0, 0 }, { 0, 0, 0, 1 }, { 0, 0, 0, 0 }, { 0, -1, 0, 0 } },
101 { { 0, 0, 0, -1 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 1, 0, 0, 0 } },
102 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
103 { { 0, 1, 0, 0 }, { -1, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } } },
104 { { { 0, 0, 0, 0 }, { 0, 0, -1, 0 }, { 0, 1, 0, 0 }, { 0, 0, 0, 0 } },
105 { { 0, 0, 1, 0 }, { 0, 0, 0, 0 }, { -1, 0, 0, 0 }, { 0, 0, 0, 0 } },
106 { { 0, -1, 0, 0 }, { 1, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
107 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } } } };
108 for ( int i = 0; i < 4; i++ )
109 {
110 for ( int j = 0; j < 4; j++ )
111 {
112 G[i][j] = GG[i][j];
113 for ( int k = 0; k < 4; k++ )
114 {
115 for ( int l = 0; l < 4; l++ ) { E[i][j][k][l] = EE[i][j][k][l]; }
116 }
117 }
118 }
119}
120
121void EvtDToKSpieta::initProbMax() {
122 setProbMax( 21.0 ); // setProbMax(1680.0);
123}
124
125void EvtDToKSpieta::decay( EvtParticle* p ) {
126 //-----------for max value------------------
127 /* double maxprob = 0.0;
128 for(int ir=0;ir<=60000000;ir++){
129 p->initializePhaseSpace(getNDaug(),getDaugs());
130 EvtVector4R ks0 = p->getDaug(0)->getP4();
131 EvtVector4R pic = p->getDaug(1)->getP4();
132 EvtVector4R pi0 = p->getDaug(2)->getP4();
133
134 double Ks0[4],Pic[4],Pi0[4];
135 Ks0[0] = ks0.get(0); Pic[0] = pic.get(0); Pi0[0] =
136 pi0.get(0); Ks0[1] = ks0.get(1); Pic[1] = pic.get(1); Pi0[1] = pi0.get(1); Ks0[2] =
137 ks0.get(2); Pic[2] = pic.get(2); Pi0[2] = pi0.get(2); Ks0[3] = ks0.get(3); Pic[3] =
138 pic.get(3); Pi0[3] = pi0.get(3); double value; calPDF(Ks0, Pic, Pi0, value);
139 if(value>maxprob) {
140 maxprob=value;
141 std::cout << "Max PDF = " << ir << " prob= " << value <<
142 std::endl;
143 }
144 }
145 std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
146 return;*/
147 //-----------------------------------------------
148 p->initializePhaseSpace( getNDaug(), getDaugs() );
149 EvtVector4R ks0 = p->getDaug( 0 )->getP4();
150 EvtVector4R pic = p->getDaug( 1 )->getP4();
151 EvtVector4R pi0 = p->getDaug( 2 )->getP4();
152
153 double Ks0[4], Pic[4], Pi0[4];
154 Ks0[0] = ks0.get( 0 );
155 Pic[0] = pic.get( 0 );
156 Pi0[0] = pi0.get( 0 );
157 Ks0[1] = ks0.get( 1 );
158 Pic[1] = pic.get( 1 );
159 Pi0[1] = pi0.get( 1 );
160 Ks0[2] = ks0.get( 2 );
161 Pic[2] = pic.get( 2 );
162 Pi0[2] = pi0.get( 2 );
163 Ks0[3] = ks0.get( 3 );
164 Pic[3] = pic.get( 3 );
165 Pi0[3] = pi0.get( 3 );
166
167 double value;
168 calPDF( Ks0, Pic, Pi0, value );
169 setProb( value );
170 return;
171}
172
173void EvtDToKSpieta::calPDF( double Ks0[], double Pic[], double Pi0[], double& Result ) {
174 double cof[2], amp_tmp1[2], amp_tmp2[2], amp_tmp[2], amp_PDF[2], PDF[2];
175 double flag[3], mass_R[2], width_R[2];
176 double pro[2];
177 double sa[3], sb[3], sc[3], B[3];
178 double t1D[4], t1V[4], t1A[4];
179 double pS[4], pV[4], pD[4], pA[4];
180 double pro1[2], pro2[2], proKPi_S[2];
181
182 double rD2 = 25.0;
183 double rRes2 = 9.0;
184 double rRes = 9.0;
185 double mass1[2] = { ma0, mKn_1430 };
186 double width1[2] = { Ga0, GKn_1430 };
187 double g0[2] = { 3, 1 };
188 double temp_PDF = 0;
189 // PDF[0]=0;
190 // PDF[1]=0;
191 double P12[4], P23[4], P13[4];
192 double scpi, snpi, sks0;
193 double s12, s13, s23;
194 for ( int i = 0; i < 4; i++ )
195 {
196 P12[i] = Pic[i] + Ks0[i];
197 P13[i] = Pi0[i] + Ks0[i];
198 P23[i] = Pic[i] + Pi0[i];
199 }
200 scpi = SCADot( Pic, Pic );
201 snpi = SCADot( Pi0, Pi0 );
202 sks0 = SCADot( Ks0, Ks0 );
203 s12 = SCADot( P12, P12 );
204 s13 = SCADot( P13, P13 );
205 s23 = SCADot( P23, P23 );
206 double mass1sq;
207 double Amp_KPiS[2];
208 amp_PDF[0] = 0;
209 amp_PDF[1] = 0;
210 PDF[0] = 0;
211 PDF[1] = 0;
212 // amp_tmp[0] = 0;
213 // amp_tmp[1] = 0;
214
215 for ( int i = 0; i < 2; i++ )
216 {
217 amp_tmp[0] = 0;
218 amp_tmp[1] = 0;
219 mass1sq = mass1[i] * mass1[i];
220 cof[0] = rho[i] * cos( phi[i] );
221 cof[1] = rho[i] * sin( phi[i] );
222 temp_PDF = 0;
223 if ( modetype[i] == 23 )
224 {
225 temp_PDF = DDalitz( Pic, Pi0, Ks0, spin[i], mass1[i] );
226 if ( g0[i] == 1 )
227 propagatorRBW( mass1sq, mass1[i], width1[i], s23, scpi, snpi, rRes, spin[i], pro );
228 if ( g0[i] == 2 ) KPiSLASS( s23, scpi, snpi, pro );
229 if ( g0[i] == 3 ) propagatorFlatte( mass1[i], width1[i], s23, scpi, snpi, 0, pro );
230 if ( g0[i] == 0 )
231 {
232 pro[0] = 1;
233 pro[1] = 0;
234 }
235 amp_tmp[0] = temp_PDF * pro[0];
236 amp_tmp[1] = temp_PDF * pro[1];
237 }
238 if ( modetype[i] == 12 )
239 {
240 temp_PDF = DDalitz( Ks0, Pic, Pi0, spin[i], mass1[i] );
241 if ( g0[i] == 1 )
242 propagatorRBW( mass1sq, mass1[i], width1[i], s12, sks0, scpi, rRes, spin[i], pro );
243 if ( g0[i] == 2 )
244 {
245 propagatorFlatte( mass1[i], width1[i], s12, sks0, scpi, 1, pro );
246 pro[0] = pro[0] * ( 0.01 + 0.990 * 0.990 ) / ( 0.01 + s12 );
247 pro[1] = pro[1] * ( 0.01 + 0.990 * 0.990 ) / ( 0.01 + s12 );
248 }
249 if ( g0[i] == 3 )
250 propagatorFlatte( mass1[i], width1[i], s12, sks0, scpi, 1, pro ); // Only for a0(980)
251
252 if ( g0[i] == 0 )
253 {
254 pro[0] = 1;
255 pro[1] = 0;
256 }
257 amp_tmp[0] = temp_PDF * pro[0];
258 amp_tmp[1] = temp_PDF * pro[1];
259 }
260 if ( modetype[i] == 13 )
261 {
262 temp_PDF = DDalitz( Ks0, Pi0, Pic, spin[i], mass1[i] );
263 if ( g0[i] == 1 )
264 propagatorRBW( mass1sq, mass1[i], width1[i], s13, sks0, snpi, rRes, spin[i], pro );
265 if ( g0[i] == 2 )
266 { // K*1430 Flatte
267 double skm2[2] = { sks0, mass_Ks * mass_Ks };
268 double spi2[2] = { snpi, mass_Eta * mass_Eta };
269 propagatorKstr1430( mass1[i], s13, skm2, spi2, pro );
270 };
271 if ( g0[i] == 3 )
272 propagatorFlatte( mass1[i], width1[i], s13, sks0, snpi, 1, pro ); // Only for a0(980)
273 if ( g0[i] == 0 )
274 {
275 pro[0] = 1;
276 pro[1] = 0;
277 }
278 amp_tmp[0] = temp_PDF * pro[0];
279 amp_tmp[1] = temp_PDF * pro[1];
280 }
281 if ( modetype[i] == 132 )
282 {
283 KPiSLASS( s12, sks0, scpi, Amp_KPiS );
284 amp_tmp[0] = Amp_KPiS[0];
285 amp_tmp[1] = Amp_KPiS[1];
286 }
287 // cout<<amp_tmp[0]<<" "<<amp_tmp[1]<<endl;
288 Com_Multi( amp_tmp, cof, amp_PDF );
289 PDF[0] += amp_PDF[0];
290 PDF[1] += amp_PDF[1];
291 }
292
293 // printf("PDF[0] : %.10f \n",PDF[0]);
294 // printf("PDF[1] : %.10f \n",PDF[1]);
295 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
296 // printf("Value : %.10f \n",value);
297
298 Result = value;
299}
300
301void EvtDToKSpieta::Com_Multi( double a1[2], double a2[2], double res[2] ) {
302 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
303 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
304}
305void EvtDToKSpieta::Com_Divide( double a1[2], double a2[2], double res[2] ) {
306 double tmp = a2[0] * a2[0] + a2[1] * a2[1];
307 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / tmp;
308 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / tmp;
309}
310//------------base---------------------------------
311double EvtDToKSpieta::SCADot( double a1[4], double a2[4] ) {
312 double _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
313 return _cal;
314}
315
316double EvtDToKSpieta::barrier( int l, double sa, double sb, double sc, double r,
317 double mass ) {
318 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
319 if ( q < 0 ) q = 1e-16;
320 double z;
321 z = q * r * r;
322 double sa0;
323 sa0 = mass * mass;
324 double q0 = ( sa0 + sb - sc ) * ( sa0 + sb - sc ) / ( 4 * sa0 ) - sb;
325 if ( q0 < 0 ) q0 = 1e-16;
326 double z0 = q0 * r * r;
327 double F = 0.0;
328 if ( l == 0 ) F = 1;
329 if ( l == 1 ) F = sqrt( ( 1 + z0 ) / ( 1 + z ) );
330 if ( l == 2 ) F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) );
331 return F;
332}
333
334double EvtDToKSpieta::Barrier( int l, double sa, double sb, double sc, double r2 ) {
335 double F;
336 double tmp = sa + sb - sc;
337 double q = 0.25 * tmp * tmp / sa - sb;
338 if ( q < 0 ) q = 1e-16;
339 double z = q * r2;
340 if ( l == 1 ) { F = sqrt( 2.0 * z / ( 1.0 + z ) ); }
341 else if ( l == 2 )
342 {
343 double z2 = z * z;
344 F = sqrt( 13.0 * z2 / ( 9.0 + 3.0 * z + z2 ) );
345 }
346 else { F = 1.0; }
347 return F;
348}
349
350//------------------spin-------------------------------------------
351void EvtDToKSpieta::calt1( double daug1[4], double daug2[4], double t1[4] ) {
352 double p, pq, tmp;
353 double pa[4], qa[4];
354 for ( int i = 0; i < 4; i++ )
355 {
356 pa[i] = daug1[i] + daug2[i];
357 qa[i] = daug1[i] - daug2[i];
358 }
359 p = SCADot( pa, pa );
360 pq = SCADot( pa, qa );
361 tmp = pq / p;
362 for ( int i = 0; i < 4; i++ ) { t1[i] = qa[i] - tmp * pa[i]; }
363}
364void EvtDToKSpieta::calt2( double daug1[4], double daug2[4], double t2[4][4] ) {
365 double p, r;
366 double pa[4], t1[4];
367 calt1( daug1, daug2, t1 );
368 r = SCADot( t1, t1 ) / 3.0;
369 for ( int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
370 p = SCADot( pa, pa );
371 for ( int i = 0; i < 4; i++ )
372 {
373 for ( int j = 0; j < 4; j++ )
374 { t2[i][j] = t1[i] * t1[j] - r * ( G[i][j] - pa[i] * pa[j] / p ); }
375 }
376}
377//-------------------prop--------------------------------------------
378void EvtDToKSpieta::propagator( double mass2, double mass, double width, double sx,
379 double prop[2] ) {
380 double a[2], b[2];
381 a[0] = 1;
382 a[1] = 0;
383 b[0] = mass2 - sx;
384 b[1] = -mass * width;
385 Com_Divide( a, b, prop );
386}
387double EvtDToKSpieta::wid( double mass2, double mass, double sa, double sb, double sc,
388 double r2, int l ) {
389 double widm = 0.;
390 double m = sqrt( sa );
391 double tmp = sb - sc;
392 double tmp1 = sa + tmp;
393 double q = 0.25 * tmp1 * tmp1 / sa - sb;
394 if ( q < 0 ) q = 1e-16;
395 double tmp2 = mass2 + tmp;
396 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
397 if ( q0 < 0 ) q0 = 1e-16;
398 double z = q * r2;
399 double z0 = q0 * r2;
400 double t = q / q0;
401 if ( l == 0 ) { widm = sqrt( t ) * mass / m; }
402 else if ( l == 1 ) { widm = t * sqrt( t ) * mass / m * ( 1 + z0 ) / ( 1 + z ); }
403 else if ( l == 2 )
404 { widm = t * t * sqrt( t ) * mass / m * ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ); }
405 return widm;
406}
407double EvtDToKSpieta::widl1( double mass2, double mass, double sa, double sb, double sc,
408 double r2 ) {
409 double widm = 0.;
410 double m = sqrt( sa );
411 double tmp = sb - sc;
412 double tmp1 = sa + tmp;
413 double q = 0.25 * tmp1 * tmp1 / sa - sb;
414 if ( q < 0 ) q = 1e-16;
415 double tmp2 = mass2 + tmp;
416 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
417 if ( q0 < 0 ) q0 = 1e-16;
418 double z = q * r2;
419 double z0 = q0 * r2;
420 double F = ( 1 + z0 ) / ( 1 + z );
421 double t = q / q0;
422 widm = t * sqrt( t ) * mass / m * F;
423 return widm;
424}
425void EvtDToKSpieta::propagatorRBW( double mass2, double mass, double width, double sa,
426 double sb, double sc, double r2, int l, double prop[2] ) {
427 double a[2], b[2];
428 a[0] = 1;
429 a[1] = 0;
430 double q = 0.25 * ( sa + sb - sc ) * ( sa + sb - sc ) / sa - sb;
431
432 b[0] = mass2 - sa;
433 b[1] = -mass * width * wid( mass2, mass, sa, sb, sc, r2, l );
434 Com_Divide( a, b, prop );
435}
436
437void EvtDToKSpieta::propagatorFlatte( double mass, double width, double sa, double sb,
438 double sc, int r, double prop[2] ) {
439 double q, qKsK, qetapi;
440 // double qKsK,qetapi;
441 double rhoab[2], rhoKsK[2];
442 q = 0.25 * ( sa + sb - sc ) * ( sa + sb - sc ) / sa - sb;
443 qetapi = 0.25 * ( sa + 0.547862 * 0.547862 - 0.13957039 * 0.13957039 ) *
444 ( sa + 0.547862 * 0.547862 - 0.13957039 * 0.13957039 ) / sa -
445 0.547862 * 0.547862;
446 if ( r == 0 ) qKsK = 0.25 * sa - 0.49368 * 0.49368;
447 if ( r == 1 ) qKsK = 0.25 * ( sa + sb - sc ) * ( sa + sb - sc ) / sa - sb;
448 if ( qetapi > 0 )
449 {
450 rhoab[0] = 2 * sqrt( qetapi / sa );
451 rhoab[1] = 0;
452 }
453 if ( qetapi < 0 )
454 {
455 rhoab[0] = 0;
456 rhoab[1] = 2 * sqrt( -qetapi / sa );
457 }
458 if ( qKsK > 0 )
459 {
460 rhoKsK[0] = 2 * sqrt( qKsK / sa );
461 rhoKsK[1] = 0;
462 }
463 if ( qKsK < 0 )
464 {
465 rhoKsK[0] = 0;
466 rhoKsK[1] = 2 * sqrt( -qKsK / sa );
467 }
468 double a[2], b[2];
469 a[0] = 1;
470 a[1] = 0;
471 b[0] = mass * mass - sa + 0.341 * rhoab[1] + 0.892 * 0.341 * rhoKsK[1];
472 b[1] = -( 0.341 * rhoab[0] + 0.892 * 0.341 * rhoKsK[0] );
473 Com_Divide( a, b, prop );
474}
475
476void EvtDToKSpieta::PiPiSWASS( double sa, double sb, double sc, double prop[2] ) {
477 double tmp = sb - sc;
478 double tmp2 = sa + tmp;
479 double qs = 0.25 * tmp2 * tmp2 / sa - sb;
480 double q = sqrt( qs );
481 double a0 = -0.11 / mass_Pion;
482 prop[0] = 1 / ( 1 + a0 * a0 * q * q );
483 prop[1] = a0 * q / ( 1 + a0 * a0 * q * q );
484}
485void EvtDToKSpieta::KPiSLASS( double sa, double sb, double sc, double prop[2] ) {
486 const double m1430 = 1.441;
487 const double sa0 = 1.441 * 1.441; // m1430*m1430;
488 const double w1430 = 0.193;
489 const double Lass1 = 0.25 / sa0;
490 double tmp = sb - sc;
491 double tmp1 = sa0 + tmp;
492 double q0 = Lass1 * tmp1 * tmp1 - sb;
493 if ( q0 < 0 ) q0 = 1e-16;
494 double tmp2 = sa + tmp;
495 double qs = 0.25 * tmp2 * tmp2 / sa - sb;
496 double q = sqrt( qs );
497 double width = w1430 * q * m1430 / sqrt( sa * q0 );
498 double temp_R = atan( m1430 * width / ( sa0 - sa ) );
499 if ( temp_R < 0 ) temp_R += math_pi;
500 double deltaR = -1.915 + temp_R; // fiR=-109.7
501 double temp_F =
502 atan( 0.226 * q / ( 2.0 - 3.819 * qs ) ); // 2.0*0.113 = 0.226; 0.113*33.8 = 3.819
503 if ( temp_F < 0 ) temp_F += math_pi;
504 double deltaF = 0.002 + temp_F; // fiF=0.1
505 double deltaS = deltaR + 2.0 * deltaF;
506 double t1 = 0.96 * sin( deltaF );
507 double t2 = sin( deltaR );
508 double CF[2], CS[2];
509 CF[0] = cos( deltaF );
510 CF[1] = sin( deltaF );
511 CS[0] = cos( deltaS );
512 CS[1] = sin( deltaS );
513 prop[0] = t1 * CF[0] + t2 * CS[0];
514 prop[1] = t1 * CF[1] + t2 * CS[1];
515}
516
517void EvtDToKSpieta::Flatte_rhoab( double sa, double sb, double sc, double rho[2] ) {
518 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
519 if ( q > 0 )
520 {
521 rho[0] = 2 * sqrt( q / sa );
522 rho[1] = 0;
523 }
524 else if ( q < 0 )
525 {
526 rho[0] = 0;
527 rho[1] = 2 * sqrt( -q / sa );
528 }
529}
530
531void EvtDToKSpieta::propagatorKstr1430( double mass, double sx, double* sb, double* sc,
532 double prop[2] ) // K*1430 Flatte
533{
534 double unit[2] = { 1.0 };
535 double ci[2] = { 0, 1 };
536 double rho1[2];
537 Flatte_rhoab( sx, sb[0], sc[0], rho1 );
538 double rho2[2];
539 Flatte_rhoab( sx, sb[1], sc[1], rho2 );
540 double gKPi_Kstr1430 = 0.2990, gEtaPK_Kstr1430 = 0.0529;
541 double tmp1[2] = { gKPi_Kstr1430, 0 };
542 double tmp11[2];
543 double tmp2[2] = { gEtaPK_Kstr1430, 0 };
544 double tmp22[2];
545 Com_Multi( tmp1, rho1, tmp11 );
546 Com_Multi( tmp2, rho2, tmp22 );
547 double tmp3[2] = { tmp11[0] + tmp22[0], tmp11[1] + tmp22[1] };
548 double tmp31[2];
549 Com_Multi( tmp3, ci, tmp31 );
550 double tmp4[2] = { mass * mass - sx - tmp31[0], -1.0 * tmp31[1] };
551 Com_Divide( unit, tmp4, prop );
552}
553double EvtDToKSpieta::DDalitz( double P1[4], double P2[4], double P3[4], int Ang,
554 double mass ) {
555 double pR[4], pD[4];
556 double temp_PDF, v_re;
557 temp_PDF = 0.0;
558 v_re = 0.0;
559 double B[2], s1, s2, s3, sR, sD;
560 for ( int i = 0; i < 4; i++ )
561 {
562 pR[i] = P1[i] + P2[i];
563 pD[i] = pR[i] + P3[i];
564 }
565 s1 = SCADot( P1, P1 );
566 s2 = SCADot( P2, P2 );
567 s3 = SCADot( P3, P3 );
568 sR = SCADot( pR, pR );
569 sD = SCADot( pD, pD );
570 int G[4][4];
571 for ( int i = 0; i != 4; i++ )
572 {
573 for ( int j = 0; j != 4; j++ )
574 {
575 if ( i == j )
576 {
577 if ( i == 0 ) G[i][j] = 1;
578 else G[i][j] = -1;
579 }
580 else G[i][j] = 0;
581 }
582 }
583 if ( Ang == 0 )
584 {
585 B[0] = 1;
586 B[1] = 1;
587 temp_PDF = 1;
588 }
589 if ( Ang == 1 )
590 {
591 B[0] = barrier( 1, sR, s1, s2, 3.0, mass );
592 B[1] = barrier( 1, sD, sR, s3, 5.0, 1.869 );
593 // B[0] = Barrier(1,sR,s1,s2,9.0);
594 // B[1] = Barrier(1,sD,sR,s3,25.0);
595 double t1[4], T1[4];
596 calt1( P1, P2, t1 );
597 calt1( pR, P3, T1 );
598 temp_PDF = 0;
599 for ( int i = 0; i < 4; i++ ) { temp_PDF += t1[i] * T1[i] * G[i][i]; }
600 }
601 if ( Ang == 2 )
602 {
603 B[0] = barrier( 2, sR, s1, s2, 3.0, mass );
604 B[1] = barrier( 2, sD, sR, s3, 5.0, 1.869 );
605 // B[0] = Barrier(2,sR,s1,s2,9.0);
606 // B[1] = Barrier(2,sD,sR,s3,25.0);
607 double t2[4][4], T2[4][4];
608 calt2( P1, P2, t2 );
609 calt2( pR, P3, T2 );
610 temp_PDF = 0;
611 for ( int i = 0; i < 4; i++ )
612 {
613 for ( int j = 0; j < 4; j++ ) { temp_PDF += t2[i][j] * T2[j][i] * G[i][i] * G[j][j]; }
614 }
615 }
616 v_re = temp_PDF * B[0] * B[1];
617 return v_re;
618}
mass
double mass
Definition CosmicGenerator.cxx:128
a0
character *LEPTONflag integer iresonances real zeta5 real a0
Definition Ekhara/src/ekhara/common.h:4
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
EvtDToKSpieta::~EvtDToKSpieta
virtual ~EvtDToKSpieta()
Definition EvtDToKSpieta.cc:32
EvtDToKSpieta::getName
void getName(std::string &name)
Definition EvtDToKSpieta.cc:34
EvtDToKSpieta::decay
void decay(EvtParticle *p)
Definition EvtDToKSpieta.cc:125
EvtDToKSpieta::clone
EvtDecayBase * clone()
Definition EvtDToKSpieta.cc:36
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