BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDToPiPi0Etap.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: EvtD0ToEtapipi.cc
11//
12// Description: Routine to handle three-body decays of D0/D0_bar or D+/D-
13//
14// Modification history:
15//
16// Liaoyuan Dong Jul 3 23:56:42 2023 Module created
17// Xian Kui Oct 5 20:45:00 2023 DToPiPi0Etap
18//------------------------------------------------------------------------
19
20#include "EvtDToPiPi0Etap.hh"
21#include "../EvtGenBase/EvtComplex.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/EvtVector4R.hh"
29#include <stdlib.h>
30using namespace std;
31
32EvtDToPiPi0Etap::~EvtDToPiPi0Etap() {}
33
34void EvtDToPiPi0Etap::getName( std::string& model_name ) { model_name = "DToPiPi0Etap"; }
35
36EvtDecayBase* EvtDToPiPi0Etap::clone() { return new EvtDToPiPi0Etap; }
37
38void EvtDToPiPi0Etap::init() {
39
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 modetype[0] = 12;
48 modetype[1] = 12;
49 rho[0] = 1.0; // rho+
50 phi[0] = 0.0;
51 rho[1] = 8.3996e+00; // (Pi+Pi0) D-wave
52 phi[1] = 4.8088e+00;
53
54 mass[0] = 0.77511;
55 width[0] = 0.1491;
56
57 mass[1] = 0.27454719;
58 width[1] = 0.1491; // TBD
59
60 massDp = 1.86966;
61 massD0 = 1.86484;
62
63 massEtap = 0.95778;
64 massPip = 0.13957;
65 massPi0 = 0.13498;
66
67 mathPi = 3.1415926;
68
69 // GS
70 GS1 = 0.6366197830;
71 GS2 = 0.0186018247;
72 GS3 = 0.1591549431; // 1/(2*mathPi)
73 GS4 = 0.0062006081; // massPip*massPip/mathPi
74
75 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
76
77 for ( int i = 0; i < 4; i++ )
78 {
79 for ( int j = 0; j < 4; j++ ) { G[i][j] = GG[i][j]; }
80 }
81}
82
83void EvtDToPiPi0Etap::initProbMax() {
84 setProbMax( 48.9 ); // The Max Prob = 48.895731256754906724
85}
86
87void EvtDToPiPi0Etap::decay( EvtParticle* p ) {
88
89 // This piece of code could in principle be used to calculate maximum
90 // probablity on fly. But as it uses high number of points and model
91 // deals with single final state, we keep hardcoded number for now rather
92 // than adapting code to work here.
93
94 /*
95 double maxprob = 0.0;
96 for(int ir=0;ir<=60000000;ir++){
97 p->initializePhaseSpace(getNDaug(),getDaugs());
98 EvtVector4R D1 = p->getDaug(0)->getP4();
99 EvtVector4R D2 = p->getDaug(1)->getP4();
100 EvtVector4R D3 = p->getDaug(2)->getP4();
101
102 double P1[4], P2[4], P3[4];
103 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
104 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
105 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
106
107 double value;
108 int g0[2]={1,0};
109 int spin[2]={1,2};
110 int nstates=2;
111
112 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates,
113 value); if(value>maxprob) { maxprob=value; cout << "ir = " << ir << " maxProb= " << value
114 << endl;
115 }
116 }
117 cout << "The maxProb = " << maxprob << endl;
118 */
119
120 p->initializePhaseSpace( getNDaug(), getDaugs() );
121 EvtVector4R D1 = p->getDaug( 0 )->getP4();
122 EvtVector4R D2 = p->getDaug( 1 )->getP4();
123 EvtVector4R D3 = p->getDaug( 2 )->getP4();
124
125 double P1[4], P2[4], P3[4];
126
127 P1[0] = D1.get( 0 );
128 P1[1] = D1.get( 1 );
129 P1[2] = D1.get( 2 );
130 P1[3] = D1.get( 3 );
131 P2[0] = D2.get( 0 );
132 P2[1] = D2.get( 1 );
133 P2[2] = D2.get( 2 );
134 P2[3] = D2.get( 3 );
135 P3[0] = D3.get( 0 );
136 P3[1] = D3.get( 1 );
137 P3[2] = D3.get( 2 );
138 P3[3] = D3.get( 3 );
139
140 double value;
141 int g0[2] = { 1, 0 };
142 int spin[2] = { 1, 2 };
143 int nstates = 2;
144 calEva( P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value );
145
146 setProb( value );
147 return;
148}
149
150/* *
151 * Base
152 * */
153
154void EvtDToPiPi0Etap::Com_Multi( double a1[2], double a2[2], double res[2] ) {
155 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
156 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
157}
158
159void EvtDToPiPi0Etap::Com_Divide( double a1[2], double a2[2], double res[2] ) {
160 double tmp = a2[0] * a2[0] + a2[1] * a2[1];
161 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / tmp;
162 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / tmp;
163}
164
165double EvtDToPiPi0Etap::SCADot( double a1[4], double a2[4] ) {
166 double _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
167 return _cal;
168}
169
170double EvtDToPiPi0Etap::barrier( int l, double sa, double sb, double sc, double r,
171 double mass ) {
172 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
173 if ( q < 0 ) q = -q;
174 double z;
175 z = q * r * r;
176 double sa0;
177 sa0 = mass * mass;
178 double q0 = ( sa0 + sb - sc ) * ( sa0 + sb - sc ) / ( 4 * sa0 ) - sb;
179 if ( q0 < 0 ) q0 = -1 * q0;
180 double z0 = q0 * r * r;
181 double F = 0.0;
182 if ( l == 0 ) F = 1;
183 if ( l == 1 ) F = sqrt( ( 1 + z0 ) / ( 1 + z ) );
184 if ( l == 2 ) F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) );
185 return F;
186}
187
188/* *
189 * Spin
190 * */
191
192void EvtDToPiPi0Etap::calt1( double daug1[4], double daug2[4], double t1[4] ) {
193 double p, pq, tmp;
194 double pa[4], qa[4];
195 for ( int i = 0; i < 4; i++ )
196 {
197 pa[i] = daug1[i] + daug2[i];
198 qa[i] = daug1[i] - daug2[i];
199 }
200 p = SCADot( pa, pa );
201 pq = SCADot( pa, qa );
202 tmp = pq / p;
203 for ( int i = 0; i < 4; i++ ) { t1[i] = qa[i] - tmp * pa[i]; }
204}
205
206void EvtDToPiPi0Etap::calt2( double daug1[4], double daug2[4], double t2[4][4] ) {
207 double p, r;
208 double pa[4], t1[4];
209 calt1( daug1, daug2, t1 );
210 r = SCADot( t1, t1 ) / 3.0;
211 for ( int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
212 p = SCADot( pa, pa );
213 for ( int i = 0; i < 4; i++ )
214 {
215 for ( int j = 0; j < 4; j++ )
216 { t2[i][j] = t1[i] * t1[j] - r * ( G[i][j] - pa[i] * pa[j] / p ); }
217 }
218}
219
220/* *
221 * Propagator
222 * */
223
224double EvtDToPiPi0Etap::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 = -q;
232 double tmp2 = mass2 + tmp;
233 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
234 if ( q0 < 0 ) q0 = -q0;
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}
244
245double EvtDToPiPi0Etap::widl1( double mass2, double mass, double sa, double sb, double sc,
246 double r2 ) {
247 double widm = 0.;
248 double m = sqrt( sa );
249 double tmp = sb - sc;
250 double tmp1 = sa + tmp;
251 double q = 0.25 * tmp1 * tmp1 / sa - sb;
252 if ( q < 0 ) q = -q;
253 double tmp2 = mass2 + tmp;
254 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
255 if ( q0 < 0 ) q0 = -q0;
256 double z = q * r2;
257 double z0 = q0 * r2;
258 double F = ( 1 + z0 ) / ( 1 + z );
259 double t = q / q0;
260 widm = t * sqrt( t ) * mass / m * F;
261 return widm;
262}
263
264void EvtDToPiPi0Etap::propagatorRBW( double mass2, double mass, double width, double sa,
265 double sb, double sc, double r2, int l, double prop[2] ) {
266 double a[2], b[2];
267 a[0] = 1;
268 a[1] = 0;
269
270 b[0] = mass2 - sa;
271 b[1] = -mass * width * wid( mass2, mass, sa, sb, sc, r2, l );
272 Com_Divide( a, b, prop );
273}
274
275void EvtDToPiPi0Etap::propagatorGS( double mass2, double mass, double width, double sa,
276 double sb, double sc, double r2, double prop[2] ) {
277 double a[2], b[2];
278 double tmp = sb - sc;
279 double tmp1 = sa + tmp;
280 double q2 = 0.25 * tmp1 * tmp1 / sa - sb;
281 if ( q2 < 0 ) q2 = 1e-16;
282
283 double tmp2 = mass2 + tmp;
284 double q02 = 0.25 * tmp2 * tmp2 / mass2 - sb;
285 if ( q02 < 0 ) q02 = 1e-16;
286
287 double q = sqrt( q2 );
288 double q0 = sqrt( q02 );
289 double m = sqrt( sa );
290 double q03 = q0 * q02;
291 double tmp3 = log( mass + 2 * q0 ) + 1.292621885; // ln(massPip + massPi0) = -1.292621885;
292
293 double h = GS1 * q / m * ( log( m + 2 * q ) + 1.292621885 );
294 double h0 = GS1 * q0 / mass * tmp3;
295 double dh = h0 * ( 0.125 / q02 - 0.5 / mass2 ) + GS3 / mass2;
296 double d = GS2 / q02 * tmp3 + GS3 * mass / q0 - GS4 * mass / q03;
297 double f = mass2 / q03 * ( q2 * ( h - h0 ) + ( mass2 - sa ) * q02 * dh );
298
299 a[0] = 1.0 + d * width / mass;
300 a[1] = 0.0;
301 b[0] = mass2 - sa + width * f;
302 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r2 );
303 Com_Divide( a, b, prop );
304}
305
306double EvtDToPiPi0Etap::DDalitz( double P1[4], double P2[4], double P3[4], int Ang,
307 double mass ) {
308 double pR[4], pD[4];
309 double temp_PDF, v_re;
310 temp_PDF = 0.0;
311 v_re = 0.0;
312 double B[2], s1, s2, s3, sR, sD;
313 for ( int i = 0; i < 4; i++ )
314 {
315 pR[i] = P1[i] + P2[i];
316 pD[i] = pR[i] + P3[i];
317 }
318 s1 = SCADot( P1, P1 );
319 s2 = SCADot( P2, P2 );
320 s3 = SCADot( P3, P3 );
321 sR = SCADot( pR, pR );
322 sD = SCADot( pD, pD );
323 int G[4][4];
324 for ( int i = 0; i != 4; i++ )
325 {
326 for ( int j = 0; j != 4; j++ )
327 {
328 if ( i == j )
329 {
330 if ( i == 0 ) G[i][j] = 1;
331 else G[i][j] = -1;
332 }
333 else G[i][j] = 0;
334 }
335 }
336
337 if ( Ang == 0 )
338 {
339 B[0] = 1;
340 B[1] = 1;
341 temp_PDF = 1;
342 }
343
344 if ( Ang == 1 )
345 {
346 B[0] = barrier( 1, sR, s1, s2, 3.0, mass );
347 B[1] = barrier( 1, sD, sR, s3, 5.0, massDp );
348 double t1[4], T1[4];
349 calt1( P1, P2, t1 );
350 calt1( pR, P3, T1 );
351 temp_PDF = 0;
352 for ( int i = 0; i < 4; i++ ) { temp_PDF += t1[i] * T1[i] * G[i][i]; }
353 }
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, massDp );
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
369 v_re = temp_PDF * B[0] * B[1];
370 return v_re;
371}
372
373/* *
374 * Calculation
375 * */
376
377void EvtDToPiPi0Etap::calEva( double* Pip, double* Pi0, double* Etap, double* mass,
378 double* width, double* amp, double* phase, int* g0, int* spin,
379 int* modetype, int nstates, double& Result ) {
380 double P12[4], P23[4], P13[4];
381 double coEff[2], tempPDF[2], ultPDF[2];
382 double rRes = 9;
383
384 double sPi0, sEtap, sPip;
385 double s12, s13, s23;
386
387 for ( int i = 0; i < 4; i++ )
388 {
389 P12[i] = Pip[i] + Pi0[i];
390 P13[i] = Pip[i] + Etap[i];
391 P23[i] = Pi0[i] + Etap[i];
392 }
393
394 sPip = SCADot( Pip, Pip );
395 sPi0 = SCADot( Pi0, Pi0 );
396 sEtap = SCADot( Etap, Etap );
397
398 s12 = SCADot( P12, P12 );
399 s13 = SCADot( P13, P13 );
400 s23 = SCADot( P23, P23 );
401
402 double prop[2], partAmp, tempAmp[2];
403 double massSq;
404
405 ultPDF[0] = 0;
406 ultPDF[1] = 0;
407
408 tempPDF[0] = 0;
409 tempPDF[1] = 0;
410
411 for ( int i = 0; i < nstates; i++ )
412 {
413 tempAmp[0] = 0;
414 tempAmp[1] = 0;
415 partAmp = 0;
416
417 massSq = mass[i] * mass[i];
418 coEff[0] = amp[i] * cos( phase[i] );
419 coEff[1] = amp[i] * sin( phase[i] );
420
421 if ( modetype[i] == 12 )
422 {
423 partAmp = DDalitz( Pip, Pi0, Etap, spin[i], mass[i] );
424 if ( g0[i] == 1 ) propagatorGS( massSq, mass[i], width[i], s12, sPip, sPi0, rRes, prop );
425 if ( g0[i] == 0 )
426 {
427 prop[0] = 1;
428 prop[1] = 0;
429 }
430 tempAmp[0] = partAmp * prop[0];
431 tempAmp[1] = partAmp * prop[1];
432 }
433
434 if ( modetype[i] == 13 )
435 {
436 partAmp = DDalitz( Pip, Etap, Pi0, spin[i], mass[i] );
437 if ( g0[i] == 1 )
438 propagatorRBW( massSq, mass[i], width[i], s13, sPip, sEtap, rRes, spin[i], prop );
439 if ( g0[i] == 0 )
440 {
441 prop[0] = 1;
442 prop[1] = 0;
443 }
444 tempAmp[0] = partAmp * prop[0];
445 tempAmp[1] = partAmp * prop[1];
446 }
447
448 if ( modetype[i] == 23 )
449 {
450 partAmp = DDalitz( Pi0, Etap, Pip, spin[i], mass[i] );
451 if ( g0[i] == 1 )
452 propagatorRBW( massSq, mass[i], width[i], s23, sPi0, sEtap, rRes, spin[i], prop );
453 if ( g0[i] == 0 )
454 {
455 prop[0] = 1;
456 prop[1] = 0;
457 }
458 tempAmp[0] = partAmp * prop[0];
459 tempAmp[1] = partAmp * prop[1];
460 }
461
462 if ( modetype[i] == 1111 )
463 {
464 tempAmp[0] = 1;
465 tempAmp[1] = 0;
466 }
467
468 Com_Multi( tempAmp, coEff, tempPDF );
469
470 ultPDF[0] += tempPDF[0];
471 ultPDF[1] += tempPDF[1];
472 }
473
474 double value = ultPDF[0] * ultPDF[0] + ultPDF[1] * ultPDF[1];
475
476 if ( value <= 0 ) value = 1e-20;
477 Result = value;
478}
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
EvtDToPiPi0Etap::~EvtDToPiPi0Etap
virtual ~EvtDToPiPi0Etap()
Definition EvtDToPiPi0Etap.cc:32
EvtDToPiPi0Etap::decay
void decay(EvtParticle *p)
Definition EvtDToPiPi0Etap.cc:87
EvtDToPiPi0Etap::clone
EvtDecayBase * clone()
Definition EvtDToPiPi0Etap.cc:36
EvtDToPiPi0Etap::getName
void getName(std::string &name)
Definition EvtDToPiPi0Etap.cc:34
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