BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDsToEtappipi0.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: EvtDsToEtappipi0.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 "EvtDsToEtappipi0.hh"
27#include <fstream>
28#include <stdlib.h>
29#include <string>
30using std::endl;
31
33
34void EvtDsToEtappipi0::getName( std::string& model_name ) { model_name = "DsToEtappipi0"; }
35
37
39 // check that there are 0 arguments
40 checkNArg( 0 );
41 checkNDaug( 3 );
46
47 phi[0] = 0;
48 rho[0] = 1;
49 phi[1] = 0;
50 rho[1] = 0;
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
60 phi[0] = 0.0;
61 rho[0] = 1.0;
62
63 modetype[0] = 23;
64
65 // cout << "DsToEtappipi0 :" << endl;
66 // for (int i=0; i<6; i++) {
67 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
68 // }
69
70 width[0] = 0.1491;
71 mass[0] = 0.77526;
72
73 mD = 1.86486;
74 mDs = 1.9683;
75 rRes = 9.0;
76 rD = 5.0;
77 metap = 0.95778;
78 mkstr = 0.89594;
79 mk0 = 0.497614;
80 mass_Kaon = 0.49368;
81 mass_Pion = 0.13957;
82 mass_Pi0 = 0.1349766;
83 math_pi = 3.1415926;
84 ma0 = 0.99;
85 Ga0 = 0.0756;
86 meta = 0.547862;
87
88 GS1 = 0.636619783;
89 GS2 = 0.01860182466;
90 GS3 = 0.1591549458; // 1/(2*math_2pi)
91 GS4 = 0.00620060822; // mass_Pion2/math_pi
92
93 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
94 for ( int i = 0; i < 4; i++ )
95 {
96 for ( int j = 0; j < 4; j++ ) { G[i][j] = GG[i][j]; }
97 }
98}
99
101
103 /*
104 double maxprob = 0.0;
105 for(int ir=0;ir<=40000000;ir++){
106 p->initializePhaseSpace(getNDaug(),getDaugs());
107 EvtVector4R D1 = p->getDaug(0)->getP4();
108 EvtVector4R D2 = p->getDaug(1)->getP4();
109 EvtVector4R D3 = p->getDaug(2)->getP4();
110
111 double P1[4], P2[4], P3[4];
112 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
113 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
114 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
115
116 double value;
117 //value = calDalEva(P1, P2, P3);
118 int g0[6]={1};
119 int spin[6]={1};
120 int nstates=1;
121 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value);
122 if (value<0) continue;
123 if(value>maxprob) {
124 maxprob=value;
125 cout << "ir= " << ir << endl;
126 cout << "double P1[4] = {" << P1[0] <<","<< P1[1] <<","<< P1[2] <<","<< P1[3] <<"};"<< endl;
127 cout << "double P2[4] = {" << P2[0] <<","<< P2[1] <<","<< P2[2] <<","<< P2[3] <<"};"<< endl;
128 cout << "double P3[4] = {" << P3[0] <<","<< P3[1] <<","<< P3[2] <<","<< P3[3] <<"};"<< endl;
129 cout << "MAX====> " << maxprob << endl;
130 }
131 }
132 printf("MAXprob = %.10f\n",maxprob);
133 */
135 EvtVector4R D1 = p->getDaug( 0 )->getP4();
136 EvtVector4R D2 = p->getDaug( 1 )->getP4();
137 EvtVector4R D3 = p->getDaug( 2 )->getP4();
138
139 double P1[4], P2[4], P3[4];
140 P1[0] = D1.get( 0 );
141 P1[1] = D1.get( 1 );
142 P1[2] = D1.get( 2 );
143 P1[3] = D1.get( 3 );
144 P2[0] = D2.get( 0 );
145 P2[1] = D2.get( 1 );
146 P2[2] = D2.get( 2 );
147 P2[3] = D2.get( 3 );
148 P3[0] = D3.get( 0 );
149 P3[1] = D3.get( 1 );
150 P3[2] = D3.get( 2 );
151 P3[3] = D3.get( 3 );
152
153 double value;
154 int g0[6] = { 1 };
155 int spin[6] = { 1 };
156 int nstates = 1;
157
158 calEva( P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value );
159
160 setProb( value );
161 return;
162}
163
164void EvtDsToEtappipi0::Com_Multi( double a1[2], double a2[2], double res[2] ) {
165 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
166 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
167}
168void EvtDsToEtappipi0::Com_Divide( double a1[2], double a2[2], double res[2] ) {
169 double tmp = a2[0] * a2[0] + a2[1] * a2[1];
170 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / tmp;
171 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / tmp;
172}
173//------------base---------------------------------
174double EvtDsToEtappipi0::SCADot( double a1[4], double a2[4] ) {
175 double _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
176 return _cal;
177}
178double EvtDsToEtappipi0::barrier( int l, double sa, double sb, double sc, double r,
179 double mass ) {
180 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
181 if ( q < 0 ) q = 1e-16;
182 double z;
183 z = q * r * r;
184 double sa0;
185 sa0 = mass * mass;
186 double q0 = ( sa0 + sb - sc ) * ( sa0 + sb - sc ) / ( 4 * sa0 ) - sb;
187 if ( q0 < 0 ) q0 = 1e-16;
188 double z0 = q0 * r * r;
189 double F = 0.0;
190 if ( l == 0 ) F = 1;
191 if ( l == 1 ) F = sqrt( ( 1 + z0 ) / ( 1 + z ) );
192 if ( l == 2 ) F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) );
193 return F;
194}
195
196void EvtDsToEtappipi0::calt1( double daug1[4], double daug2[4], double t1[4] ) {
197 double p, pq, tmp;
198 double pa[4], qa[4];
199 for ( int i = 0; i < 4; i++ )
200 {
201 pa[i] = daug1[i] + daug2[i];
202 qa[i] = daug1[i] - daug2[i];
203 }
204 p = SCADot( pa, pa );
205 pq = SCADot( pa, qa );
206 tmp = pq / p;
207 for ( int i = 0; i < 4; i++ ) { t1[i] = qa[i] - tmp * pa[i]; }
208}
209void EvtDsToEtappipi0::calt2( double daug1[4], double daug2[4], double t2[4][4] ) {
210 double p, r;
211 double pa[4], t1[4];
212 calt1( daug1, daug2, t1 );
213 r = SCADot( t1, t1 ) / 3.0;
214 for ( int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
215 p = SCADot( pa, pa );
216 for ( int i = 0; i < 4; i++ )
217 {
218 for ( int j = 0; j < 4; j++ )
219 { t2[i][j] = t1[i] * t1[j] - r * ( G[i][j] - pa[i] * pa[j] / p ); }
220 }
221}
222//-------------------prop--------------------------------------------
223
224double EvtDsToEtappipi0::wid( double mass2, double mass, double sa, double sb, double sc,
225 double r2, int l ) {
226 double widm = 0.;
227 double m = sqrt( sa );
228 double tmp = sb - sc;
229 double tmp1 = sa + tmp;
230 double q = 0.25 * tmp1 * tmp1 / sa - sb;
231 if ( q < 0 ) q = 1e-16;
232 double tmp2 = mass2 + tmp;
233 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
234 if ( q0 < 0 ) q0 = 1e-16;
235 double z = q * r2;
236 double z0 = q0 * r2;
237 double t = q / q0;
238 if ( l == 0 ) { widm = sqrt( t ) * mass / m; }
239 else if ( l == 1 ) { widm = t * sqrt( t ) * mass / m * ( 1 + z0 ) / ( 1 + z ); }
240 else if ( l == 2 )
241 { widm = t * t * sqrt( t ) * mass / m * ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ); }
242 return widm;
243}
244double EvtDsToEtappipi0::widl1( double mass2, double mass, double sa, double sb, double sc,
245 double r2 ) {
246 double widm = 0.;
247 double m = sqrt( sa );
248 double tmp = sb - sc;
249 double tmp1 = sa + tmp;
250 double q = 0.25 * tmp1 * tmp1 / sa - sb;
251 if ( q < 0 ) q = 1e-16;
252 double tmp2 = mass2 + tmp;
253 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
254 if ( q0 < 0 ) q0 = 1e-16;
255 double z = q * r2;
256 double z0 = q0 * r2;
257 double F = ( 1 + z0 ) / ( 1 + z );
258 double t = q / q0;
259 widm = t * sqrt( t ) * mass / m * F;
260 return widm;
261}
262void EvtDsToEtappipi0::propagatorRBW( double mass2, double mass, double width, double sa,
263 double sb, double sc, double r2, int l,
264 double prop[2] ) {
265 double a[2], b[2];
266 a[0] = 1;
267 a[1] = 0;
268 b[0] = mass2 - sa;
269 b[1] = -mass * width * wid( mass2, mass, sa, sb, sc, r2, l );
270 Com_Divide( a, b, prop );
271}
272
273void EvtDsToEtappipi0::propagatorGS( double mass2, double mass, double width, double sa,
274 double sb, double sc, double r2, double prop[2] ) {
275 double a[2], b[2];
276 double tmp = sb - sc;
277 double tmp1 = sa + tmp;
278 double q2 = 0.25 * tmp1 * tmp1 / sa - sb;
279 if ( q2 < 0 ) q2 = 1e-16;
280
281 double tmp2 = mass2 + tmp;
282 double q02 = 0.25 * tmp2 * tmp2 / mass2 - sb;
283 if ( q02 < 0 ) q02 = 1e-16;
284
285 double q = sqrt( q2 );
286 double q0 = sqrt( q02 );
287 double m = sqrt( sa );
288 double q03 = q0 * q02;
289 double tmp3 = log( mass + 2 * q0 ) + 1.2926305904; // log(mpi0+mpip) = -1.2926305904;
290
291 double h = GS1 * q / m * ( log( m + 2 * q ) + 1.2926305904 );
292 double h0 = GS1 * q0 / mass * tmp3;
293 double dh = h0 * ( 0.125 / q02 - 0.5 / mass2 ) + GS3 / mass2;
294 double d = GS2 / q02 * tmp3 + GS3 * mass / q0 - GS4 * mass / q03;
295 double f = mass2 / q03 * ( q2 * ( h - h0 ) + ( mass2 - sa ) * q02 * dh );
296
297 a[0] = 1.0 + d * width / mass;
298 a[1] = 0.0;
299 b[0] = mass2 - sa + width * f;
300 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r2 );
301 Com_Divide( a, b, prop );
302}
303
304void EvtDsToEtappipi0::PiPiSWASS( double sa, double sb, double sc, double prop[2] ) {
305 double tmp = sb - sc;
306 double tmp2 = sa + tmp;
307 double qs = 0.25 * tmp2 * tmp2 / sa - sb;
308 double q = sqrt( qs );
309 double a0 = -0.11 / mass_Pion;
310 prop[0] = 1 / ( 1 + a0 * a0 * q * q );
311 prop[1] = a0 * q / ( 1 + a0 * a0 * q * q );
312}
313
314void EvtDsToEtappipi0::Flatte_rhoab( double sa, double sb, double rho[2] ) {
315 double q = 1.0 - ( 4 * sb / sa );
316
317 if ( q > 0 )
318 {
319 rho[0] = sqrt( q );
320 rho[1] = 0;
321 }
322 else if ( q < 0 )
323 {
324 rho[0] = 0;
325 rho[1] = sqrt( -q );
326 }
327}
328
329void EvtDsToEtappipi0::propagator980( double mass, double sx, double* sb, double prop[2] ) {
330
331 double gpipi1 = 2.0 / 3.0 * 0.165;
332 double gpipi2 = 1.0 / 3.0 * 0.165;
333 double gKK1 = 0.5 * 0.69466;
334 double gKK2 = 0.5 * 0.69466;
335
336 double unit[2] = { 1.0 };
337 double ci[2] = { 0, 1 };
338 double rho1[2];
339 Flatte_rhoab( sx, sb[0], rho1 );
340 double rho2[2];
341 Flatte_rhoab( sx, sb[1], rho2 );
342 double rho3[2];
343 Flatte_rhoab( sx, sb[2], rho3 );
344 double rho4[2];
345 Flatte_rhoab( sx, sb[3], rho4 );
346
347 double tmp1[2] = { gpipi1, 0 };
348 double tmp11[2];
349 double tmp2[2] = { gpipi2, 0 };
350 double tmp22[2];
351 double tmp3[2] = { gKK1, 0 };
352 double tmp33[2];
353 double tmp4[2] = { gKK2, 0 };
354 double tmp44[2];
355
356 Com_Multi( tmp1, rho1, tmp11 );
357 Com_Multi( tmp2, rho2, tmp22 );
358 Com_Multi( tmp3, rho3, tmp33 );
359 Com_Multi( tmp4, rho4, tmp44 );
360
361 double tmp5[2] = { tmp11[0] + tmp22[0] + tmp33[0] + tmp44[0],
362 tmp11[1] + tmp22[1] + tmp33[1] + tmp44[1] };
363 double tmp51[2];
364 Com_Multi( tmp5, ci, tmp51 );
365 double tmp6[2] = { mass * mass - sx - tmp51[0], -1.0 * tmp51[1] };
366
367 Com_Divide( unit, tmp6, prop );
368}
369
370double EvtDsToEtappipi0::DDalitz( double P1[4], double P2[4], double P3[4], int Ang,
371 double mass ) {
372 double pR[4], pD[4];
373 double temp_PDF, v_re;
374 temp_PDF = 0.0;
375 v_re = 0.0;
376 double B[2], s1, s2, s3, sR, sD;
377 for ( int i = 0; i < 4; i++ )
378 {
379 pR[i] = P1[i] + P2[i];
380 pD[i] = pR[i] + P3[i];
381 }
382 s1 = SCADot( P1, P1 );
383 s2 = SCADot( P2, P2 );
384 s3 = SCADot( P3, P3 );
385 sR = SCADot( pR, pR );
386 sD = SCADot( pD, pD );
387 // int G[4][4];
388 // for(int i=0; i!=4; i++){
389 // for(int j=0; j!=4; j++){
390 // if(i==j){
391 // if(i==0) G[i][j] = 1;
392 // else G[i][j] = -1;
393 // }
394 // else G[i][j] = 0;
395 // }
396 // }
397 if ( Ang == 0 )
398 {
399 B[0] = 1;
400 B[1] = 1;
401 temp_PDF = 1;
402 }
403 if ( Ang == 1 )
404 {
405 B[0] = barrier( 1, sR, s1, s2, 3.0, mass );
406 B[1] = barrier( 1, sD, sR, s3, 5.0, 1.9683 );
407 double t1[4], T1[4];
408 calt1( P1, P2, t1 );
409 calt1( pR, P3, T1 );
410 temp_PDF = 0;
411 for ( int i = 0; i < 4; i++ ) { temp_PDF += t1[i] * T1[i] * G[i][i]; }
412 }
413 if ( Ang == 2 )
414 {
415 B[0] = barrier( 2, sR, s1, s2, 3.0, mass );
416 B[1] = barrier( 2, sD, sR, s3, 5.0, 1.9683 );
417 double t2[4][4], T2[4][4];
418 calt2( P1, P2, t2 );
419 calt2( pR, P3, T2 );
420 temp_PDF = 0;
421 for ( int i = 0; i < 4; i++ )
422 {
423 for ( int j = 0; j < 4; j++ ) { temp_PDF += t2[i][j] * T2[j][i] * G[i][i] * G[j][j]; }
424 }
425 }
426 v_re = temp_PDF * B[0] * B[1];
427 return v_re;
428}
429
430void EvtDsToEtappipi0::calEva( double* EtaP, double* Pi, double* Pi0, double* mass1,
431 double* width1, double* amp, double* phase, int* g0, int* spin,
432 int* modetype, int nstates, double& Result ) {
433 double P12[4], P23[4], P13[4];
434 double cof[2], amp_PDF[2], PDF[2];
435 double spi, spi0, setap;
436 double s12, s13, s23;
437 for ( int i = 0; i < 4; i++ )
438 {
439 P23[i] = Pi[i] + Pi0[i];
440 P12[i] = EtaP[i] + Pi[i];
441 P13[i] = EtaP[i] + Pi0[i];
442 }
443 setap = SCADot( EtaP, EtaP );
444 spi = SCADot( Pi, Pi );
445 spi0 = SCADot( Pi0, Pi0 );
446 s12 = SCADot( P12, P12 );
447 s13 = SCADot( P13, P13 );
448 s23 = SCADot( P23, P23 );
449
450 double pro[2], temp_PDF, amp_tmp[2];
451 double mass1sq;
452 amp_PDF[0] = 0;
453 amp_PDF[1] = 0;
454 PDF[0] = 0;
455 PDF[1] = 0;
456 amp_tmp[0] = 0;
457 amp_tmp[1] = 0;
458 for ( int i = 0; i < nstates; i++ )
459 {
460 amp_tmp[0] = 0;
461 amp_tmp[1] = 0;
462 mass1sq = mass1[i] * mass1[i];
463 cof[0] = amp[i] * cos( phase[i] );
464 cof[1] = amp[i] * sin( phase[i] );
465 temp_PDF = 0;
466
467 if ( modetype[i] == 23 )
468 {
469 temp_PDF = DDalitz( Pi, Pi0, EtaP, spin[i], mass1[i] );
470 if ( g0[i] == 1 ) propagatorGS( mass1sq, mass1[i], width1[i], s23, spi, spi0, 9.0, pro );
471 if ( g0[i] == 12 ) PiPiSWASS( s23, spi, spi0, pro );
472 if ( g0[i] == 0 )
473 {
474 pro[0] = 1;
475 pro[1] = 0;
476 }
477 // printf("cof1 = %.15f \n",cof[1]);
478
479 amp_tmp[0] = temp_PDF * pro[0];
480 amp_tmp[1] = temp_PDF * pro[1];
481 }
482 Com_Multi( amp_tmp, cof, amp_PDF );
483 PDF[0] += amp_PDF[0];
484 PDF[1] += amp_PDF[1];
485 }
486 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
487 // printf("value = %.15f\n",value);
488 Result = value;
489}
double mass
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
character *LEPTONflag integer iresonances real zeta5 real a0
*******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
****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
void checkSpinDaughter(int d1, EvtSpinType::spintype sp)
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)
EvtDecayBase * clone()
virtual ~EvtDsToEtappipi0()
void decay(EvtParticle *p)
void getName(std::string &name)
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