BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDsToEta3pi.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: EvtDsToEta3pi.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 "EvtDsToEta3pi.hh"
18#include "../EvtGenBase/EvtConst.hh"
19#include "../EvtGenBase/EvtDecayTable.hh"
20#include "../EvtGenBase/EvtFlatte.hh"
21#include "../EvtGenBase/EvtGenKine.hh"
22#include "../EvtGenBase/EvtPDL.hh"
23#include "../EvtGenBase/EvtParticle.hh"
24#include "../EvtGenBase/EvtPatches.hh"
25#include "../EvtGenBase/EvtReport.hh"
26#include "../EvtGenBase/EvtResonance.hh"
27#include "../EvtGenBase/EvtResonance2.hh"
28#include <fstream>
29#include <stdlib.h>
30#include <string>
31using std::endl;
32
33EvtDsToEta3pi::~EvtDsToEta3pi() {}
34
35void EvtDsToEta3pi::getName( std::string& model_name ) { model_name = "DsToEta3pi"; }
36
37EvtDecayBase* EvtDsToEta3pi::clone() { return new EvtDsToEta3pi; }
38
39void EvtDsToEta3pi::init() {
40
41 // check that there are 0 arguments
42 checkNArg( 0 );
43 checkNDaug( 4 );
44
45 checkSpinParent( EvtSpinType::SCALAR );
46 checkSpinDaughter( 0, EvtSpinType::SCALAR );
47 checkSpinDaughter( 1, EvtSpinType::SCALAR );
48 checkSpinDaughter( 2, EvtSpinType::SCALAR );
49 checkSpinDaughter( 3, EvtSpinType::SCALAR );
50
51 int mode = 0;
52 phi[0] = 2.4600e+00;
53 modetype[0] = 1;
54 phi[1] = 0.0000e+00;
55 modetype[1] = 2;
56 phi[2] = -1.5394e+00;
57 modetype[2] = 3;
58 phi[3] = -6.0185e+00;
59 modetype[3] = 4;
60 phi[4] = -6.0185e+00;
61 modetype[4] = 5;
62 phi[5] = -6.1365e+00;
63 modetype[5] = 6;
64 phi[6] = -6.1365e+00;
65 modetype[6] = 7;
66 phi[7] = 1.3514e+00;
67 modetype[7] = 8;
68 phi[8] = 2.4527e+00;
69 modetype[8] = 9;
70 phi[9] = -2.0956e+00;
71 modetype[9] = 10;
72 phi[10] = -2.0956e+00;
73 modetype[10] = 11;
74
75 rho[0] = -1.5823e+00;
76 rho[1] = 10.000e+00;
77 rho[2] = 1.3750e+01;
78 rho[3] = -3.7257e-01;
79 rho[4] = -3.7257e-01;
80 rho[5] = -3.6293e+00;
81 rho[6] = -3.6293e+00;
82 rho[7] = -8.9184e+00;
83 rho[8] = -1.4292e+01;
84 rho[9] = 1.2518e+00;
85 rho[10] = 1.2518e+00;
86
87 mass1[0] = 9.9000e-01;
88 mass1[1] = 7.7526e-01;
89 mass1[2] = 5.2600e-01;
90 mass1[3] = 9.9000e-01;
91 mass1[4] = 9.9000e-01;
92 mass1[5] = 9.9000e-01;
93 mass1[6] = 9.9000e-01;
94 mass1[7] = 9.6500e-01;
95 mass1[8] = 5.2600e-01;
96 mass1[9] = 9.9000e-01;
97 mass1[10] = 9.9000e-01;
98
99 mass2[0] = 7.7526e-01;
100 mass2[1] = 1.3462e+00;
101 mass2[2] = 1.3462e+00;
102 mass2[3] = 1.4050e+00;
103 mass2[4] = 1.4050e+00;
104 mass2[5] = 1.8000e+00;
105 mass2[6] = 1.8000e+00;
106 mass2[7] = 1.8000e+00;
107 mass2[8] = 1.8000e+00;
108 mass2[9] = 1.4502e+00;
109 mass2[10] = 1.4502e+00;
110
111 width1[0] = 7.5000e-02;
112 width1[1] = 1.4910e-01;
113 width1[2] = 5.0000e-01;
114 width1[3] = 7.5000e-02;
115 width1[4] = 7.5000e-02;
116 width1[5] = 7.5000e-02;
117 width1[6] = 7.5000e-02;
118 width1[7] = 7.5000e-02;
119 width1[8] = 5.0000e-01;
120 width1[9] = 7.5000e-02;
121 width1[10] = 7.5000e-02;
122
123 width2[0] = 1.4910e-01;
124 width2[1] = 3.6900e-01;
125 width2[2] = 3.6900e-01;
126 width2[3] = 5.1000e-02;
127 width2[4] = 5.1000e-02;
128 width2[5] = 3.8600e-01;
129 width2[6] = 3.8600e-01;
130 width2[7] = 3.8600e-01;
131 width2[8] = 3.8600e-01;
132 width2[9] = 5.4000e-02;
133 width2[10] = 5.4000e-02;
134
135 // cout << "DsToEta3pi :"<< endl;
136 // for (int i=0; i<11; i++) {
137 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << " m1= " << mass1[i] << " m2= "
138 // << mass2[i] << " w1= "
139 // << width1[i] << " w2= " << width2[i] << endl;
140 // }
141
142 rRes = 3.0;
143 rRes2 = 9.0;
144 rD = 5.0;
145 rD2 = 25.0;
146 mass_Ks = 0.497611;
147 mass_Kaon = 0.49368;
148 mass_Pion = 0.13957;
149 mass_Pi0 = 0.1349766;
150 meta = 0.547862;
151 GS1 = 0.636619783;
152 GS2 = 0.01860182466;
153 GS3 = 0.1591549458; // 1/(2*math_2pi)
154 GS4 = 0.00620060822;
155 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
156 int EE[4][4][4][4] = {
157 { { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
158 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 1 }, { 0, 0, -1, 0 } },
159 { { 0, 0, 0, 0 }, { 0, 0, 0, -1 }, { 0, 0, 0, 0 }, { 0, 1, 0, 0 } },
160 { { 0, 0, 0, 0 }, { 0, 0, 1, 0 }, { 0, -1, 0, 0 }, { 0, 0, 0, 0 } } },
161 { { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, -1 }, { 0, 0, 1, 0 } },
162 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
163 { { 0, 0, 0, 1 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { -1, 0, 0, 0 } },
164 { { 0, 0, -1, 0 }, { 0, 0, 0, 0 }, { 1, 0, 0, 0 }, { 0, 0, 0, 0 } } },
165 { { { 0, 0, 0, 0 }, { 0, 0, 0, 1 }, { 0, 0, 0, 0 }, { 0, -1, 0, 0 } },
166 { { 0, 0, 0, -1 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 1, 0, 0, 0 } },
167 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
168 { { 0, 1, 0, 0 }, { -1, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } } },
169 { { { 0, 0, 0, 0 }, { 0, 0, -1, 0 }, { 0, 1, 0, 0 }, { 0, 0, 0, 0 } },
170 { { 0, 0, 1, 0 }, { 0, 0, 0, 0 }, { -1, 0, 0, 0 }, { 0, 0, 0, 0 } },
171 { { 0, -1, 0, 0 }, { 1, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
172 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } } } };
173 for ( int i = 0; i < 4; i++ )
174 {
175 for ( int j = 0; j < 4; j++ )
176 {
177 G[i][j] = GG[i][j];
178 for ( int k = 0; k < 4; k++ )
179 {
180 for ( int l = 0; l < 4; l++ ) { E[i][j][k][l] = EE[i][j][k][l]; }
181 }
182 }
183 }
184}
185
186void EvtDsToEta3pi::initProbMax() { setProbMax( 108000.0 ); }
187
188void EvtDsToEta3pi::decay( EvtParticle* p ) {
189
190 // double maxprob = 0.0;
191 // for(int ir=0;ir<=60000000;ir++){
192 // p->initializePhaseSpace(getNDaug(),getDaugs());
193 // EvtVector4R _pip1 = p->getDaug(0)->getP4();
194 // EvtVector4R _pip2 = p->getDaug(1)->getP4();
195 // EvtVector4R _pim = p->getDaug(2)->getP4();
196 // EvtVector4R _eta = p->getDaug(3)->getP4();
197 //
198 // double _Pip1[4],_Pip2[4],_Pim[4],_Eta[4];
199 // _Pip1[1]=_pip1.get(1);_Pip1[2]=_pip1.get(2);_Pip1[3]=_pip1.get(3);
200 // _Pip2[1]=_pip2.get(1);_Pip2[2]=_pip2.get(2);_Pip2[3]=_pip2.get(3);
201 // _Pim[1] = _pim.get(1);_Pim[2] = _pim.get(2);_Pim[3] = _pim.get(3);
202 // _Eta[1] = _eta.get(1);_Eta[2] = _eta.get(2);_Eta[3] = _eta.get(3);
203 //
204 // double _prob;
205 // int g0[11]={1,1,1,1,1,1,1,1,1,1,1};
206 // int g1[11]={1,1,1,1,1,1,1,1,1,1,1};
207 // int g2[11]={0,0,0,0,0,0,0,0,0,0,0};
208 // int nstates=11;
209 // calEvaMy(_Pip1, _Pip2, _Pim, _Eta, mass1, mass2, width1, width2, rho, phi, g0,g1,g2,
210 // modetype, nstates, _prob); if(_prob>maxprob) {
211 // maxprob=_prob;
212 // printf("nrun = %d,maxprob = %.10f,prob = %.10f\n\n",ir,maxprob,_prob);
213 // }
214 // }
215
216 p->initializePhaseSpace( getNDaug(), getDaugs() );
217 EvtVector4R pip1 = p->getDaug( 0 )->getP4();
218 EvtVector4R pip2 = p->getDaug( 1 )->getP4();
219 EvtVector4R pim = p->getDaug( 2 )->getP4();
220 EvtVector4R eta = p->getDaug( 3 )->getP4();
221
222 double Pip1[4], Pip2[4], Pim[4], Eta[4];
223 Pip1[1] = pip1.get( 1 );
224 Pip1[2] = pip1.get( 2 );
225 Pip1[3] = pip1.get( 3 );
226 Pip2[1] = pip2.get( 1 );
227 Pip2[2] = pip2.get( 2 );
228 Pip2[3] = pip2.get( 3 );
229 Pim[1] = pim.get( 1 );
230 Pim[2] = pim.get( 2 );
231 Pim[3] = pim.get( 3 );
232 Eta[1] = eta.get( 1 );
233 Eta[2] = eta.get( 2 );
234 Eta[3] = eta.get( 3 );
235
236 double prob;
237 int g0[11] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 };
238 int g1[11] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 };
239 int g2[11] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
240 int nstates = 11;
241
242 calEvaMy( Pip1, Pip2, Pim, Eta, mass1, mass2, width1, width2, rho, phi, g0, g1, g2, modetype,
243 nstates, prob );
244 setProb( prob );
245
246 return;
247}
248
249void EvtDsToEta3pi::calEvaMy( double* Pip1, double* Pip2, double* Pim, double* Eta,
250 double* mass1, double* mass2, double* width1, double* width2,
251 double* amp, double* phase, int* g0, int* g1, int* g2,
252 int* modetype, int nstates, double& Result ) {
253 double prho1[4], prho2[4], pa11[4], pa12[4], pa01[4], pa02[4], pD[4], pf11[4], pf12[4],
254 pa0m[4], peta12951[4], peta12952[4], psigma1[4], psigma2[4];
255 Pip1[0] = sqrt( mass_Pion * mass_Pion + Pip1[1] * Pip1[1] + Pip1[2] * Pip1[2] +
256 Pip1[3] * Pip1[3] );
257 Pip2[0] = sqrt( mass_Pion * mass_Pion + Pip2[1] * Pip2[1] + Pip2[2] * Pip2[2] +
258 Pip2[3] * Pip2[3] );
259 Pim[0] = sqrt( mass_Pion * mass_Pion + Pim[1] * Pim[1] + Pim[2] * Pim[2] + Pim[3] * Pim[3] );
260 Eta[0] = sqrt( meta * meta + Eta[1] * Eta[1] + Eta[2] * Eta[2] + Eta[3] * Eta[3] );
261
262 for ( int i = 0; i != 4; i++ )
263 {
264 prho1[i] = Pim[i] + Pip1[i];
265 prho2[i] = Pim[i] + Pip2[i];
266 pa11[i] = prho1[i] + Pip2[i];
267 pa12[i] = prho2[i] + Pip1[i];
268 psigma1[i] = Pim[i] + Pip1[i];
269 psigma2[i] = Pim[i] + Pip2[i];
270 pa01[i] = Eta[i] + Pip1[i];
271 pa02[i] = Eta[i] + Pip2[i];
272 pa0m[i] = Eta[i] + Pim[i];
273 pf11[i] = Pim[i] + Pip1[i] + Eta[i];
274 pf12[i] = Pim[i] + Pip2[i] + Eta[i];
275 pD[i] = pa11[i] + Eta[i];
276 peta12951[i] = Pim[i] + Pip1[i] + Eta[i];
277 peta12952[i] = Pim[i] + Pip2[i] + Eta[i];
278 }
279 double cof[2], amp_tmp1[2], amp_tmp2[2], amp_tmp[2], amp_PDF[2], PDF[2];
280 double amp_Omg, phs_Omg;
281 amp_PDF[0] = 0;
282 amp_PDF[1] = 0;
283 PDF[0] = 0;
284 PDF[1] = 0;
285 amp_PDF[0] = 0;
286 amp_PDF[1] = 0;
287 PDF[0] = 0;
288 PDF[1] = 0;
289 double spion1, spion2, spim, seta, sa11, sa12, srho1, srho2, sa01, sa02, sD, sf11, sf12,
290 sa0m, seta12951, seta12952, ssigma1, ssigma2;
291 double t1rho1[4], t2rho1[4][4], t1rho2[4], t2rho2[4][4], t1a01[4], t1a02[4], t2a11[4][4],
292 t2a12[4][4], t1d01[4], t1d02[4], t1f11[4], t2f11[4][4], t1f12[4], t2f12[4][4], t1a11[4],
293 t1a12[4], t1f11_2[4], t1f12_2[4];
294 double temp_PDF, tmp1, tmp2, temp[2], tmp3, tmp4, temp_PDF1;
295 double pro[2], pro0[2], pro1[2], pro2[2], pro3[2], pro4[2];
296 double t1D[4], t2D[4][4], B[3], t1Tm[4], A[2], Bc[3];
297
298 double mass1sq, mass2sq;
299
300 seta = SCADot( Eta, Eta );
301 spion1 = SCADot( Pip1, Pip1 );
302 spion2 = SCADot( Pip2, Pip2 );
303 spim = SCADot( Pim, Pim );
304
305 srho1 = SCADot( prho1, prho1 );
306 srho2 = SCADot( prho2, prho2 );
307 sa11 = SCADot( pa11, pa11 );
308 sa12 = SCADot( pa12, pa12 );
309
310 sa01 = SCADot( pa01, pa01 );
311 sa02 = SCADot( pa02, pa02 );
312 sa0m = SCADot( pa0m, pa0m );
313 sf11 = SCADot( pf11, pf11 );
314 sf12 = SCADot( pf12, pf12 );
315 sD = SCADot( pD, pD );
316 seta12951 = SCADot( peta12951, peta12951 );
317 seta12952 = SCADot( peta12952, peta12952 );
318 ssigma1 = SCADot( psigma1, psigma1 );
319 ssigma2 = SCADot( psigma2, psigma2 );
320 double seta2[2] = { mass_Ks * mass_Ks, seta };
321 double spion12[2] = { mass_Kaon * mass_Kaon, spion1 };
322 double spion22[2] = { mass_Kaon * mass_Kaon, spion2 };
323 double spim2[2] = { mass_Kaon * mass_Kaon, spim };
324
325 calt1( Pip1, Pim, t1rho1 );
326 calt2( Pip1, Pim, t2rho1 );
327 calt1( Pip2, Pim, t1rho2 );
328 calt2( Pip2, Pim, t2rho2 );
329
330 calt1( Pip1, Eta, t1a01 );
331 calt1( Pip2, Eta, t1a02 );
332
333 calt1( pa01, prho2, t1d01 );
334 calt1( pa02, prho1, t1d02 );
335
336 calt1( pa0m, Pip1, t1f11 );
337 calt1( pa0m, Pip2, t1f12 );
338
339 calt1( psigma1, Pip2, t1a11 );
340 calt1( psigma2, Pip1, t1a12 );
341
342 calt1( pa01, Pim, t1f11_2 );
343 calt1( pa02, Pim, t1f12_2 );
344
345 calt2( prho1, Pip2, t2a11 );
346 calt2( prho2, Pip1, t2a12 );
347
348 calt2( pa01, Pip1, t2f11 );
349 calt2( pa02, Pip2, t2f12 );
350
351 for ( int i = 0; i < nstates; i++ )
352 {
353 int d = 0;
354 cof[0] = amp[i] * cos( phase[i] );
355 cof[1] = amp[i] * sin( phase[i] );
356 mass1sq = mass1[i] * mass1[i];
357 mass2sq = mass2[i] * mass2[i];
358 temp_PDF = 0;
359 temp_PDF1 = 0;
360 amp_tmp[0] = 0;
361 amp_tmp[1] = 0;
362 amp_tmp1[0] = 0;
363 amp_tmp1[1] = 0;
364 amp_tmp2[0] = 0;
365 amp_tmp2[1] = 0;
366 tmp1 = 0;
367 tmp2 = 0;
368 tmp3 = 0;
369 tmp4 = 0;
370
371 if ( modetype[i] == 1 )
372 {
373 if ( g0[i] == 1 ) propagatora0980( mass1[i], sa01, seta2, spion12, pro0 );
374 if ( g0[i] == 0 )
375 {
376 pro0[0] = 1;
377 pro0[1] = 0;
378 }
379 if ( g1[i] == 1 )
380 propagatorGS( mass2sq, mass2[i], width2[i], srho2, spion2, spim, rRes2, pro1 );
381 if ( g1[i] == 0 )
382 {
383 pro1[0] = 1;
384 pro1[1] = 0;
385 }
386 Com_Multi( pro0, pro1, pro );
387
388 for ( int a = 0; a < 4; a++ ) { temp_PDF += G[a][a] * t1d01[a] * t1rho2[a]; }
389
390 B[1] = Barrier( 1, srho2, spion2, spim, rRes2 );
391 B[2] = Barrier( 1, sD, sa01, srho2, rD2 );
392 tmp1 = B[1] * B[2] * temp_PDF;
393 amp_tmp1[0] = tmp1 * pro[0];
394 amp_tmp1[1] = tmp1 * pro[1];
395
396 temp_PDF = 0;
397 if ( g0[i] == 1 ) propagatora0980( mass1[i], sa02, seta2, spion22, pro0 );
398 if ( g0[i] == 0 )
399 {
400 pro0[0] = 1;
401 pro0[1] = 0;
402 }
403 if ( g1[i] == 1 )
404 propagatorGS( mass2sq, mass2[i], width2[i], srho1, spion1, spim, rRes2, pro1 );
405 if ( g1[i] == 0 )
406 {
407 pro1[0] = 1;
408 pro1[1] = 0;
409 }
410 Com_Multi( pro0, pro1, pro );
411
412 for ( int a = 0; a < 4; a++ ) { temp_PDF += G[a][a] * t1d02[a] * t1rho1[a]; }
413
414 B[1] = Barrier( 1, srho1, spion1, spim, rRes2 );
415 B[2] = Barrier( 1, sD, sa02, srho1, rD2 );
416 tmp2 = B[1] * B[2] * temp_PDF;
417 amp_tmp2[0] = tmp2 * pro[0];
418 amp_tmp2[1] = tmp2 * pro[1];
419 }
420 if ( modetype[i] == 2 )
421 {
422 if ( g0[i] == 1 )
423 propagatorGS( mass1sq, mass1[i], width1[i], srho1, spion1, spim, rRes2, pro0 );
424 if ( g0[i] == 0 )
425 {
426 pro0[0] = 1;
427 pro0[1] = 0;
428 }
429 if ( g1[i] == 1 )
430 propagatorRBW( mass2sq, mass2[i], width2[i], sa11, srho1, spion2, rRes2, g2[i], pro1 );
431 if ( g1[i] == 0 )
432 {
433 pro1[0] = 1;
434 pro1[1] = 0;
435 }
436 Com_Multi( pro0, pro1, pro );
437 calt1( pa11, Eta, t1D );
438 B[0] = Barrier( 1, srho1, spion1, spim, rRes2 );
439 B[2] = Barrier( 1, sD, sa11, seta, rD2 );
440 if ( g2[i] == 0 )
441 {
442 for ( int a = 0; a < 4; a++ )
443 {
444 for ( int j = 0; j < 4; j++ )
445 {
446 temp_PDF += t1D[a] * ( G[a][j] - pa11[a] * pa11[j] / sa11 ) * t1rho1[j] * G[a][a] *
447 G[j][j];
448 }
449 }
450 tmp1 = B[0] * B[2] * temp_PDF;
451 }
452 // if(g2[i]==2)
453 // {
454 // for(int a=0; a<4; a++)
455 // {
456 // for(int j=0; j<4; j++)
457 // {
458 // temp_PDF += t1D[a]*t2a11[a][j]*t1rho1[j]*G[a][a]*G[j][j];
459 // }
460 // }
461 // B[1] = Barrier(2,sa11,srho1,spion2,rRes2);
462 // tmp1 = B[0]*B[1]*B[2]*temp_PDF;
463 // }
464 amp_tmp1[0] = tmp1 * pro[0];
465 amp_tmp1[1] = tmp1 * pro[1];
466
467 temp_PDF = 0;
468 if ( g0[i] == 1 )
469 propagatorGS( mass1sq, mass1[i], width1[i], srho2, spion2, spim, rRes2, pro0 );
470 if ( g0[i] == 0 )
471 {
472 pro0[0] = 1;
473 pro0[1] = 0;
474 }
475 if ( g1[i] == 1 )
476 propagatorRBW( mass2sq, mass2[i], width2[i], sa12, srho2, spion1, rRes2, g2[i], pro1 );
477 if ( g1[i] == 0 )
478 {
479 pro1[0] = 1;
480 pro1[1] = 0;
481 }
482 Com_Multi( pro0, pro1, pro );
483 calt1( pa12, Eta, t1D );
484 B[0] = Barrier( 1, srho2, spion2, spim, rRes2 );
485 B[2] = Barrier( 1, sD, sa12, seta, rD2 );
486 if ( g2[i] == 0 )
487 {
488 for ( int a = 0; a < 4; a++ )
489 {
490 for ( int j = 0; j < 4; j++ )
491 {
492 temp_PDF += t1D[a] * ( G[a][j] - pa12[a] * pa12[j] / sa12 ) * t1rho2[j] * G[a][a] *
493 G[j][j];
494 }
495 }
496 tmp2 = B[0] * B[2] * temp_PDF;
497 }
498 // if(g2[i]==2)
499 // {
500 // for(int a=0; a<4; a++)
501 // {
502 // for(int j=0; j<4; j++)
503 // {
504 // temp_PDF += t1D[a]*t2a12[a][j]*t1rho2[j]*G[a][a]*G[j][j];
505 // }
506 // }
507 // B[1] = Barrier(2,sa12,srho2,spion1,rRes2);
508 // tmp2 = B[0]*B[1]*B[2]*temp_PDF;
509 // }
510 amp_tmp2[0] = tmp2 * pro[0];
511 amp_tmp2[1] = tmp2 * pro[1];
512 }
513 if ( modetype[i] == 3 )
514 {
515 if ( g0[i] == 1 ) propagatorsigma500( ssigma1, spion1, spim, pro0 );
516 if ( g0[i] == 0 )
517 {
518 pro0[0] = 1;
519 pro0[1] = 0;
520 }
521 if ( g1[i] == 1 )
522 propagatorRBW( mass2sq, mass2[i], width2[i], sa11, ssigma1, spion2, rRes2, 1, pro1 );
523 if ( g1[i] == 0 )
524 {
525 pro1[0] = 1;
526 pro1[1] = 0;
527 }
528 // printf("pro0 = %.10f, %.10f\n", pro0[0], pro0[1]);
529 // printf("pro1 = %.10f, %.10f\n", pro1[0], pro1[1]);
530 Com_Multi( pro0, pro1, pro );
531 calt1( pa11, Eta, t1D );
532 B[2] = Barrier( 1, sD, sa11, seta, rD2 );
533 B[1] = Barrier( 1, sa11, ssigma1, spion2, rRes2 );
534 for ( int a = 0; a < 4; a++ ) { temp_PDF += t1D[a] * t1a11[a] * G[a][a]; }
535 tmp1 = B[1] * B[2] * temp_PDF;
536 amp_tmp1[0] = tmp1 * pro[0];
537 amp_tmp1[1] = tmp1 * pro[1];
538 temp_PDF = 0;
539 if ( g0[i] == 1 ) propagatorsigma500( ssigma2, spion2, spim, pro0 );
540 if ( g0[i] == 0 )
541 {
542 pro0[0] = 1;
543 pro0[1] = 0;
544 }
545 if ( g1[i] == 1 )
546 propagatorRBW( mass2sq, mass2[i], width2[i], sa12, ssigma2, spion1, rRes2, 1, pro1 );
547 if ( g1[i] == 0 )
548 {
549 pro1[0] = 1;
550 pro1[1] = 0;
551 }
552 // printf("pro0 = %f, %f\n", pro0[0], pro0[1]);
553 // printf("pro1 = %f, %f\n", pro1[0], pro1[1]);
554 Com_Multi( pro0, pro1, pro );
555 calt1( pa12, Eta, t1D );
556 B[2] = Barrier( 1, sD, sa12, seta, rD2 );
557 B[1] = Barrier( 1, sa12, ssigma2, spion1, rRes2 );
558 for ( int a = 0; a < 4; a++ ) { temp_PDF += t1D[a] * t1a12[a] * G[a][a]; }
559 tmp2 = B[1] * B[2] * temp_PDF;
560 amp_tmp2[0] = tmp2 * pro[0];
561 amp_tmp2[1] = tmp2 * pro[1];
562 }
563 if ( modetype[i] == 4 )
564 {
565 if ( g0[i] == 1 ) propagatora0980( mass1[i], sa0m, spim2, seta2, pro0 );
566 if ( g0[i] == 0 )
567 {
568 pro0[0] = 1;
569 pro0[1] = 0;
570 pro2[0] = 1;
571 pro2[1] = 0;
572 }
573 if ( g1[i] == 1 )
574 propagatorRBW( mass2sq, mass2[i], width2[i], seta12951, sa0m, spion1, rRes2, 0, pro1 );
575 if ( g1[i] == 0 )
576 {
577 pro1[0] = 1;
578 pro1[1] = 0;
579 pro3[0] = 1;
580 pro3[1] = 0;
581 }
582 Com_Multi( pro0, pro1, pro );
583 calt1( pf11, Pip2, t1D );
584 tmp1 = 1;
585 amp_tmp1[0] = tmp1 * pro[0];
586 amp_tmp1[1] = tmp1 * pro[1];
587
588 temp_PDF = 0;
589 if ( g0[i] == 1 ) propagatora0980( mass1[i], sa0m, spim2, seta2, pro0 );
590 if ( g0[i] == 0 )
591 {
592 pro0[0] = 1;
593 pro0[1] = 0;
594 pro2[0] = 1;
595 pro2[1] = 0;
596 }
597 if ( g1[i] == 1 )
598 propagatorRBW( mass2sq, mass2[i], width2[i], seta12952, sa0m, spion2, rRes2, 0, pro1 );
599 if ( g1[i] == 0 )
600 {
601 pro1[0] = 1;
602 pro1[1] = 0;
603 pro3[0] = 1;
604 pro3[1] = 0;
605 }
606 Com_Multi( pro0, pro1, pro );
607 calt1( pf12, Pip1, t1D );
608 tmp2 = 1;
609 amp_tmp2[0] = tmp2 * pro[0];
610 amp_tmp2[1] = tmp2 * pro[1];
611 }
612 if ( modetype[i] == 5 )
613 {
614 if ( g0[i] == 1 ) propagatora0980( mass1[i], sa01, spion12, seta2, pro2 );
615 if ( g0[i] == 0 )
616 {
617 pro0[0] = 1;
618 pro0[1] = 0;
619 pro2[0] = 1;
620 pro2[1] = 0;
621 }
622 if ( g1[i] == 1 )
623 propagatorRBW( mass2sq, mass2[i], width2[i], seta12951, sa01, spim, rRes2, 0, pro3 );
624 if ( g1[i] == 0 )
625 {
626 pro1[0] = 1;
627 pro1[1] = 0;
628 pro3[0] = 1;
629 pro3[1] = 0;
630 }
631 Com_Multi( pro2, pro3, pro4 );
632 calt1( pf11, Pip2, t1D );
633 tmp3 = 1;
634 amp_tmp1[0] = tmp3 * pro4[0];
635 amp_tmp1[1] = tmp3 * pro4[1];
636
637 temp_PDF1 = 0;
638 if ( g0[i] == 1 ) propagatora0980( mass1[i], sa02, spion22, seta2, pro2 );
639 if ( g0[i] == 0 )
640 {
641 pro0[0] = 1;
642 pro0[1] = 0;
643 pro2[0] = 1;
644 pro2[1] = 0;
645 }
646 if ( g1[i] == 1 )
647 propagatorRBW( mass2sq, mass2[i], width2[i], seta12952, sa02, spim, rRes2, 0, pro3 );
648 if ( g1[i] == 0 )
649 {
650 pro1[0] = 1;
651 pro1[1] = 0;
652 pro3[0] = 1;
653 pro3[1] = 0;
654 }
655 Com_Multi( pro2, pro3, pro4 );
656 calt1( pf12, Pip1, t1D );
657 tmp4 = 1;
658 amp_tmp2[0] = tmp4 * pro4[0];
659 amp_tmp2[1] = tmp4 * pro4[1];
660 }
661 if ( modetype[i] == 6 )
662 {
663 if ( g0[i] == 1 ) propagatora0980( mass1[i], sa0m, spim2, seta2, pro0 );
664 if ( g0[i] == 0 )
665 {
666 pro0[0] = 1;
667 pro0[1] = 0;
668 pro2[0] = 1;
669 pro2[1] = 0;
670 }
671 // if(g1[i]==0)
672 // propagatorRBW(mass2sq,mass2[i],width2[i],seta12951,sa0m,spion1,rRes2,0,pro1);
673 if ( g1[i] == 1 )
674 {
675 pro1[0] = 1;
676 pro1[1] = 0;
677 pro3[0] = 1;
678 pro3[1] = 0;
679 }
680 Com_Multi( pro0, pro1, pro );
681 calt1( pf11, Pip2, t1D );
682 tmp1 = 1;
683 amp_tmp1[0] = tmp1 * pro[0];
684 amp_tmp1[1] = tmp1 * pro[1];
685 temp_PDF = 0;
686 if ( g0[i] == 1 ) propagatora0980( mass1[i], sa0m, spim2, seta2, pro0 );
687 if ( g0[i] == 0 )
688 {
689 pro0[0] = 1;
690 pro0[1] = 0;
691 pro2[0] = 1;
692 pro2[1] = 0;
693 }
694 // if(g1[i]==0)
695 // propagatorRBW(mass2sq,mass2[i],width2[i],seta12952,sa0m,spion2,rRes2,0,pro1);
696 if ( g1[i] == 1 )
697 {
698 pro1[0] = 1;
699 pro1[1] = 0;
700 pro3[0] = 1;
701 pro3[1] = 0;
702 }
703 Com_Multi( pro0, pro1, pro );
704 calt1( pf12, Pip1, t1D );
705 tmp2 = 1;
706 amp_tmp2[0] = tmp2 * pro[0];
707 amp_tmp2[1] = tmp2 * pro[1];
708 }
709 if ( modetype[i] == 7 )
710 {
711
712 if ( g0[i] == 1 ) propagatora0980( mass1[i], sa01, spion12, seta2, pro2 );
713 if ( g0[i] == 0 )
714 {
715 pro0[0] = 1;
716 pro0[1] = 0;
717 pro2[0] = 1;
718 pro2[1] = 0;
719 }
720 // if(g1[i]==0)
721 // propagatorRBW(mass2sq,mass2[i],width2[i],seta12951,sa01,spim,rRes2,0,pro3);
722 if ( g1[i] == 1 )
723 {
724 pro1[0] = 1;
725 pro1[1] = 0;
726 pro3[0] = 1;
727 pro3[1] = 0;
728 }
729 Com_Multi( pro2, pro3, pro4 );
730 calt1( pf11, Pip2, t1D );
731 tmp3 = 1;
732 amp_tmp1[0] = tmp3 * pro4[0];
733 amp_tmp1[1] = tmp3 * pro4[1];
734 temp_PDF1 = 0;
735 if ( g0[i] == 1 ) propagatora0980( mass1[i], sa02, spion22, seta2, pro2 );
736 if ( g0[i] == 0 )
737 {
738 pro0[0] = 1;
739 pro0[1] = 0;
740 pro2[0] = 1;
741 pro2[1] = 0;
742 }
743 // if(g1[i]==0)
744 // propagatorRBW(mass2sq,mass2[i],width2[i],seta12952,sa02,spim,rRes2,0,pro3);
745 if ( g1[i] == 1 )
746 {
747 pro1[0] = 1;
748 pro1[1] = 0;
749 pro3[0] = 1;
750 pro3[1] = 0;
751 }
752 Com_Multi( pro2, pro3, pro4 );
753 calt1( pf12, Pip1, t1D );
754 tmp4 = 1;
755 amp_tmp2[0] = tmp4 * pro4[0];
756 amp_tmp2[1] = tmp4 * pro4[1];
757 }
758 if ( modetype[i] == 8 )
759 {
760 if ( g0[i] == 1 ) propagator980( mass1[i], srho1, spion12, spim2, pro0 );
761 if ( g0[i] == 0 )
762 {
763 pro0[0] = 1;
764 pro0[1] = 0;
765 }
766 // if(g1[i]==0)
767 // propagatorRBW(mass2sq,mass2[i],width2[i],seta12951,srho1,seta,rRes2,0,pro1);
768 if ( g1[i] == 1 )
769 {
770 pro1[0] = 1;
771 pro1[1] = 0;
772 }
773 Com_Multi( pro0, pro1, pro );
774 // printf("pro0 = %f, %f\n", pro0[0], pro0[1]);
775 // printf("pro1 = %f, %f\n", pro1[0], pro1[1]);
776 // printf("pro = %f, %f\n", pro[0], pro[1]);
777 tmp1 = 1;
778 amp_tmp1[0] = tmp1 * pro[0];
779 amp_tmp1[1] = tmp1 * pro[1];
780 temp_PDF = 0;
781 if ( g0[i] == 1 ) propagator980( mass1[i], srho2, spion22, spim2, pro0 );
782 if ( g0[i] == 0 )
783 {
784 pro0[0] = 1;
785 pro0[1] = 0;
786 }
787 // if(g1[i]==0)
788 // propagatorRBW(mass2sq,mass2[i],width2[i],seta12952,srho2,seta,rRes2,0,pro1);
789 if ( g1[i] == 1 )
790 {
791 pro1[0] = 1;
792 pro1[1] = 0;
793 }
794 Com_Multi( pro0, pro1, pro );
795 // printf("pro0 = %f, %f\n", pro0[0], pro0[1]);
796 // printf("pro1 = %f, %f\n", pro1[0], pro1[1]);
797 // printf("pro = %f, %f\n", pro[0], pro[1]);
798 calt1( pf12, Pip1, t1D );
799 tmp2 = 1;
800 amp_tmp2[0] = tmp2 * pro[0];
801 amp_tmp2[1] = tmp2 * pro[1];
802 }
803 if ( modetype[i] == 9 )
804 {
805 if ( g0[i] == 1 ) propagatorsigma500( ssigma1, spion1, spim, pro0 );
806 if ( g0[i] == 0 )
807 {
808 pro0[0] = 1;
809 pro0[1] = 0;
810 }
811 // if(g1[i]==1)
812 // propagatorRBW(mass2sq,mass2[i],width2[i],seta12951,ssigma1,seta,rRes2,0,pro1);
813 if ( g1[i] == 1 )
814 {
815 pro1[0] = 1;
816 pro1[1] = 0;
817 }
818 Com_Multi( pro0, pro1, pro );
819 calt1( pf11, Pip2, t1D );
820 tmp1 = 1;
821 amp_tmp1[0] = tmp1 * pro[0];
822 amp_tmp1[1] = tmp1 * pro[1];
823 temp_PDF = 0;
824 if ( g0[i] == 1 ) propagatorsigma500( ssigma2, spion2, spim, pro0 );
825 if ( g0[i] == 0 )
826 {
827 pro0[0] = 1;
828 pro0[1] = 0;
829 }
830 // if(g1[i]==1)
831 // propagatorRBW(mass2sq,mass2[i],width2[i],seta12952,ssigma2,seta,rRes2,0,pro1);
832 if ( g1[i] == 1 )
833 {
834 pro1[0] = 1;
835 pro1[1] = 0;
836 }
837 Com_Multi( pro0, pro1, pro );
838 calt1( pf12, Pip1, t1D );
839 tmp2 = 1;
840 amp_tmp2[0] = tmp2 * pro[0];
841 amp_tmp2[1] = tmp2 * pro[1];
842 }
843
844 if ( modetype[i] == 10 )
845 {
846 if ( g0[i] == 1 ) propagatora0980( mass1[i], sa0m, seta2, spim2, pro0 );
847 if ( g0[i] == 0 )
848 {
849 pro0[0] = 1;
850 pro0[1] = 0;
851 pro2[0] = 1;
852 pro2[1] = 0;
853 }
854 if ( g1[i] == 1 )
855 propagatorRBW( mass2sq, mass2[i], width2[i], sf11, sa0m, spion1, rRes2, 1, pro1 );
856 if ( g1[i] == 0 )
857 {
858 pro1[0] = 1;
859 pro1[1] = 0;
860 pro3[0] = 1;
861 pro3[1] = 0;
862 }
863 Com_Multi( pro0, pro1, pro );
864 calt1( pf11, Pip2, t1D );
865 B[1] = Barrier( 1, sf11, sa0m, spion1, rRes2 );
866 B[2] = Barrier( 1, sD, sf11, spion2, rD2 );
867 for ( int a = 0; a < 4; a++ ) { temp_PDF += t1D[a] * t1f11[a] * G[a][a]; }
868 tmp1 = B[1] * B[2] * temp_PDF;
869 amp_tmp1[0] = tmp1 * pro[0];
870 amp_tmp1[1] = tmp1 * pro[1];
871
872 temp_PDF = 0;
873 if ( g0[i] == 1 ) propagatora0980( mass1[i], sa0m, seta2, spim2, pro0 );
874 if ( g0[i] == 0 )
875 {
876 pro0[0] = 1;
877 pro0[1] = 0;
878 pro2[0] = 1;
879 pro2[1] = 0;
880 }
881 if ( g1[i] == 1 )
882 propagatorRBW( mass2sq, mass2[i], width2[i], sf12, sa0m, spion2, rRes2, 1, pro1 );
883 if ( g1[i] == 0 )
884 {
885 pro1[0] = 1;
886 pro1[1] = 0;
887 pro3[0] = 1;
888 pro3[1] = 0;
889 }
890 Com_Multi( pro0, pro1, pro );
891 calt1( pf12, Pip1, t1D );
892 B[1] = Barrier( 1, sf12, sa0m, spion2, rRes2 );
893 B[2] = Barrier( 1, sD, sf12, spion1, rD2 );
894 for ( int a = 0; a < 4; a++ ) { temp_PDF += t1D[a] * t1f12[a] * G[a][a]; }
895 tmp2 = B[1] * B[2] * temp_PDF;
896 amp_tmp2[0] = tmp2 * pro[0];
897 amp_tmp2[1] = tmp2 * pro[1];
898 }
899
900 if ( modetype[i] == 11 )
901 {
902 if ( g0[i] == 1 ) propagatora0980( mass1[i], sa01, seta2, spion12, pro2 );
903 if ( g0[i] == 0 )
904 {
905 pro0[0] = 1;
906 pro0[1] = 0;
907 pro2[0] = 1;
908 pro2[1] = 0;
909 }
910 if ( g1[i] == 1 )
911 propagatorRBW( mass2sq, mass2[i], width2[i], sf11, sa01, spim, rRes2, 1, pro3 );
912 if ( g1[i] == 0 )
913 {
914 pro1[0] = 1;
915 pro1[1] = 0;
916 pro3[0] = 1;
917 pro3[1] = 0;
918 }
919 Com_Multi( pro2, pro3, pro4 );
920 calt1( pf11, Pip2, t1D );
921 Bc[1] = Barrier( 1, sf11, sa01, spim, rRes2 );
922 Bc[2] = Barrier( 1, sD, sf11, spion2, rD2 );
923 for ( int a = 0; a < 4; a++ ) { temp_PDF1 += t1D[a] * t1f11_2[a] * G[a][a]; }
924 tmp3 = Bc[1] * Bc[2] * temp_PDF1;
925 amp_tmp1[0] = tmp3 * pro4[0];
926 amp_tmp1[1] = tmp3 * pro4[1];
927 // printf("amp_tmp1 = %f , %f\n", amp_tmp1[0], amp_tmp1[1]);
928
929 temp_PDF1 = 0;
930 if ( g0[i] == 1 ) propagatora0980( mass1[i], sa02, seta2, spion22, pro2 );
931 if ( g0[i] == 0 )
932 {
933 pro0[0] = 1;
934 pro0[1] = 0;
935 pro2[0] = 1;
936 pro2[1] = 0;
937 }
938 if ( g1[i] == 1 )
939 propagatorRBW( mass2sq, mass2[i], width2[i], sf12, sa02, spim, rRes2, 1, pro3 );
940 if ( g1[i] == 0 )
941 {
942 pro1[0] = 1;
943 pro1[1] = 0;
944 pro3[0] = 1;
945 pro3[1] = 0;
946 }
947 Com_Multi( pro2, pro3, pro4 );
948 calt1( pf12, Pip1, t1D );
949 Bc[1] = Barrier( 1, sf12, sa02, spim, rRes2 );
950 Bc[2] = Barrier( 1, sD, sf12, spion1, rD2 );
951 for ( int a = 0; a < 4; a++ ) { temp_PDF1 += t1D[a] * t1f12_2[a] * G[a][a]; }
952 tmp4 = Bc[1] * Bc[2] * temp_PDF1;
953 amp_tmp2[0] = tmp4 * pro4[0];
954 amp_tmp2[1] = tmp4 * pro4[1];
955 // printf("amp_tmp2 = %f , %f\n", amp_tmp2[0], amp_tmp2[1]);
956 }
957 amp_tmp[0] = amp_tmp1[0] + amp_tmp2[0];
958 amp_tmp[1] = amp_tmp1[1] + amp_tmp2[1];
959 Com_Multi( amp_tmp, cof, amp_PDF );
960 PDF[0] += amp_PDF[0];
961 PDF[1] += amp_PDF[1];
962 }
963
964 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
965 if ( value <= 0 ) { value = 1e-20; }
966 Result = value;
967}
968
969void EvtDsToEta3pi::Com_Multi( double a1[2], double a2[2], double res[2] ) {
970 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
971 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
972}
973void EvtDsToEta3pi::Com_Divide( double a1[2], double a2[2], double res[2] ) {
974 double tmp = a2[0] * a2[0] + a2[1] * a2[1];
975 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / tmp;
976 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / tmp;
977}
978double EvtDsToEta3pi::SCADot( double a1[4], double a2[4] ) {
979 double _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
980 return _cal;
981}
982double EvtDsToEta3pi::Barrier( int l, double sa, double sb, double sc, double r2 ) {
983 double F;
984 double tmp = sa + sb - sc;
985 double q = 0.25 * tmp * tmp / sa - sb;
986 if ( q < 0 ) q = 1e-16;
987 double z = q * r2;
988 if ( l == 1 ) { F = sqrt( 2.0 * z / ( 1.0 + z ) ); }
989 else if ( l == 2 )
990 {
991 double z2 = z * z;
992 F = sqrt( 13.0 * z2 / ( 9.0 + 3.0 * z + z2 ) );
993 }
994 else { F = 1.0; }
995 return F;
996}
997void EvtDsToEta3pi::calt1( double daug1[4], double daug2[4], double t1[4] ) {
998 double p, pq, tmp;
999 double pa[4], qa[4];
1000 for ( int i = 0; i < 4; i++ )
1001 {
1002 pa[i] = daug1[i] + daug2[i];
1003 qa[i] = daug1[i] - daug2[i];
1004 }
1005 p = SCADot( pa, pa );
1006 pq = SCADot( pa, qa );
1007 tmp = pq / p;
1008 for ( int i = 0; i < 4; i++ ) { t1[i] = qa[i] - tmp * pa[i]; }
1009}
1010
1011void EvtDsToEta3pi::calt2( double daug1[4], double daug2[4], double t2[4][4] ) {
1012 double p, r;
1013 double pa[4], t1[4];
1014 calt1( daug1, daug2, t1 );
1015 r = SCADot( t1, t1 );
1016 for ( int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
1017 p = SCADot( pa, pa );
1018 for ( int i = 0; i < 4; i++ )
1019 {
1020 for ( int j = 0; j < 4; j++ )
1021 { t2[i][j] = t1[i] * t1[j] - 1.0 / 3 * r * ( G[i][j] - pa[i] * pa[j] / p ); }
1022 }
1023}
1024
1025void EvtDsToEta3pi::propagator( double mass2, double mass, double width, double sx,
1026 double prop[2] ) {
1027 double a[2], b[2];
1028 a[0] = 1;
1029 a[1] = 0;
1030 b[0] = mass2 - sx;
1031 b[1] = -mass * width;
1032 Com_Divide( a, b, prop );
1033}
1034double EvtDsToEta3pi::wid( double mass2, double mass, double sa, double sb, double sc,
1035 double r2, int l ) {
1036 double widm = 0.;
1037 double m = sqrt( sa );
1038 double tmp = sb - sc;
1039 double tmp1 = sa + tmp;
1040 double q = 0.25 * tmp1 * tmp1 / sa - sb;
1041 if ( q < 0 ) q = 1e-16;
1042 double tmp2 = mass2 + tmp;
1043 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
1044 if ( q0 < 0 ) q0 = 1e-16;
1045 double z = q * r2;
1046 double z0 = q0 * r2;
1047 double t = q / q0;
1048 if ( l == 0 ) { widm = sqrt( t ) * mass / m; }
1049 else if ( l == 1 ) { widm = t * sqrt( t ) * mass / m * ( 1 + z0 ) / ( 1 + z ); }
1050 else if ( l == 2 )
1051 { widm = t * t * sqrt( t ) * mass / m * ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ); }
1052 return widm;
1053}
1054double EvtDsToEta3pi::widl1( double mass2, double mass, double sa, double sb, double sc,
1055 double r2 ) {
1056 double widm = 0.;
1057 double m = sqrt( sa );
1058 double tmp = sb - sc;
1059 double tmp1 = sa + tmp;
1060 double q = 0.25 * tmp1 * tmp1 / sa - sb;
1061 if ( q < 0 ) q = 1e-16;
1062 double tmp2 = mass2 + tmp;
1063 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
1064 if ( q0 < 0 ) q0 = 1e-16;
1065 double z = q * r2;
1066 double z0 = q0 * r2;
1067 double F = ( 1 + z0 ) / ( 1 + z );
1068 double t = q / q0;
1069 widm = t * sqrt( t ) * mass / m * F;
1070 return widm;
1071}
1072void EvtDsToEta3pi::propagatorRBW( double mass2, double mass, double width, double sa,
1073 double sb, double sc, double r2, int l, double prop[2] ) {
1074 double a[2], b[2];
1075 a[0] = 1;
1076 a[1] = 0;
1077 b[0] = mass2 - sa;
1078 b[1] = -mass * width * wid( mass2, mass, sa, sb, sc, r2, l );
1079 Com_Divide( a, b, prop );
1080}
1081
1082void EvtDsToEta3pi::propagatorGS( double mass2, double mass, double width, double sa,
1083 double sb, double sc, double r2, double prop[2] ) {
1084 double a[2], b[2];
1085 double tmp = sb - sc;
1086 double tmp1 = sa + tmp;
1087 double q2 = 0.25 * tmp1 * tmp1 / sa - sb;
1088 if ( q2 < 0 ) q2 = 1e-16;
1089
1090 double tmp2 = mass2 + tmp;
1091 double q02 = 0.25 * tmp2 * tmp2 / mass2 - sb;
1092 if ( q02 < 0 ) q02 = 1e-16;
1093
1094 double q = sqrt( q2 );
1095 double q0 = sqrt( q02 );
1096 double m = sqrt( sa );
1097 double q03 = q0 * q02;
1098 double tmp3 = log( mass + 2 * q0 ) + 1.2760418309; // log(mass_2Pion) = 1.2760418309;
1099
1100 double h = GS1 * q / m * ( log( m + 2 * q ) + 1.2760418309 );
1101 double h0 = GS1 * q0 / mass * tmp3;
1102 double dh = h0 * ( 0.125 / q02 - 0.5 / mass2 ) + GS3 / mass2;
1103 double d = GS2 / q02 * tmp3 + GS3 * mass / q0 - GS4 * mass / q03;
1104 double f = mass2 / q03 * ( q2 * ( h - h0 ) + ( mass2 - sa ) * q02 * dh );
1105
1106 a[0] = 1.0 + d * width / mass;
1107 a[1] = 0.0;
1108 b[0] = mass2 - sa + width * f;
1109 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r2 );
1110 Com_Divide( a, b, prop );
1111}
1112void EvtDsToEta3pi::rhoab( double sa, double sb, double sc, double res[2] ) {
1113 double tmp = sa + sb - sc;
1114 double q = 0.25 * tmp * tmp / sa - sb;
1115 if ( q >= 0 )
1116 {
1117 res[0] = 2.0 * sqrt( q / sa );
1118 res[1] = 0.0;
1119 }
1120 else
1121 {
1122 res[0] = 0.0;
1123 res[1] = 2.0 * sqrt( -q / sa );
1124 }
1125}
1126void EvtDsToEta3pi::rho4Pi( double sa, double res[2] ) {
1127 double temp = 1.0 - 0.3116765584 / sa; // 0.3116765584=0.13957*0.13957*16
1128 if ( temp >= 0 )
1129 {
1130 res[0] = sqrt( temp ) / ( 1.0 + exp( 9.8 - 3.5 * sa ) );
1131 res[1] = 0.0;
1132 }
1133 else
1134 {
1135 res[0] = 0.0;
1136 res[1] = sqrt( -temp ) / ( 1.0 + exp( 9.8 - 3.5 * sa ) );
1137 }
1138}
1139void EvtDsToEta3pi::propagatorsigma500( double sa, double sb, double sc, double prop[2] ) {
1140 double f = 0.5843 + 1.6663 * sa;
1141 const double M = 0.9264;
1142 const double mass2 = 0.85821696; // M*M
1143 const double mpi2d2 = 0.00973989245;
1144 double g1 = f * ( sa - mpi2d2 ) / ( mass2 - mpi2d2 ) * exp( ( mass2 - sa ) / 1.082 );
1145 double rho1s[2], rho1M[2], rho2s[2], rho2M[2], rho1[2], rho2[2];
1146 rhoab( sa, sb, sc, rho1s );
1147 rhoab( mass2, sb, sc, rho1M );
1148 rho4Pi( sa, rho2s );
1149 rho4Pi( mass2, rho2M );
1150 Com_Divide( rho1s, rho1M, rho1 );
1151 Com_Divide( rho2s, rho2M, rho2 );
1152 double a[2], b[2];
1153 a[0] = 1.0;
1154 a[1] = 0.0;
1155 b[0] = mass2 - sa + M * ( g1 * rho1[1] + 0.0024 * rho2[1] );
1156 b[1] = -M * ( g1 * rho1[0] + 0.0024 * rho2[0] );
1157 Com_Divide( a, b, prop );
1158}
1159void EvtDsToEta3pi::propagatora0_980( double mass2, double mass, double width, double sa,
1160 double sb, double sc, int charge, double prop[2] ) {
1161 double q, qKsK;
1162 double rho[2], rhoKsK[2];
1163 rhoab( sa, sb, sc, rho );
1164 if ( charge == 0 ) { qKsK = 0.25 * sa - 0.2437199424; }
1165 else if ( charge == 1 )
1166 {
1167 double tmp1 = sa + 3.899750596e-03;
1168 qKsK = 0.25 * tmp1 * tmp1 / sa - 0.247619692996;
1169 }
1170 if ( qKsK > 0 )
1171 {
1172 rhoKsK[0] = 2 * sqrt( qKsK / sa );
1173 rhoKsK[1] = 0;
1174 }
1175 else if ( qKsK < 0 )
1176 {
1177 rhoKsK[0] = 0;
1178 rhoKsK[1] = 2 * sqrt( -qKsK / sa );
1179 }
1180 double a[2], b[2];
1181 a[0] = 1;
1182 a[1] = 0;
1183 b[0] = mass2 - sa + 0.341 * rho[1] + 0.304172 * rhoKsK[1]; // 0.304172=0.892*0.341
1184 b[1] = -0.341 * rho[0] - 0.304172 * rhoKsK[0];
1185 Com_Divide( a, b, prop );
1186}
1187void EvtDsToEta3pi::Flatte_rhoab( double sa, double sb, double sc, double rho[2] ) {
1188 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
1189 double b[2];
1190 if ( q > 0 )
1191 {
1192 rho[0] = 2 * sqrt( q / sa );
1193 rho[1] = 0;
1194 }
1195 else if ( q < 0 )
1196 {
1197 rho[0] = 0;
1198 rho[1] = 2 * sqrt( -q / sa );
1199 }
1200}
1201void EvtDsToEta3pi::propagator980( double mass, double sx, double* sb, double* sc,
1202 double prop[2] ) {
1203 double unit[2] = { 1.0 };
1204 double ci[2] = { 0, 1 };
1205 double rho1[2];
1206 Flatte_rhoab( sx, sb[0], sc[0], rho1 );
1207 double rho2[2];
1208 Flatte_rhoab( sx, sb[1], sc[1], rho2 );
1209
1210 double gK_f980 = 0.69466, gPi_f980 = 0.165;
1211 double tmp1[2] = { gK_f980, 0 };
1212 double tmp11[2];
1213 double tmp2[2] = { gPi_f980, 0 };
1214 double tmp22[2];
1215 Com_Multi( tmp1, rho1, tmp11 );
1216 Com_Multi( tmp2, rho2, tmp22 );
1217 double tmp3[2] = { tmp11[0] + tmp22[0], tmp11[1] + tmp22[1] };
1218 double tmp31[2];
1219 Com_Multi( tmp3, ci, tmp31 );
1220 double tmp4[2] = { mass * mass - sx - tmp31[0], -1.0 * tmp31[1] };
1221 Com_Divide( unit, tmp4, prop );
1222}
1223
1224void EvtDsToEta3pi::propagatora0980( double mass, double sx, double* sb, double* sc,
1225 double prop[2] ) {
1226 double unit[2] = { 1.0 };
1227 double ci[2] = { 0, 1 };
1228 double rho1[2];
1229 Flatte_rhoab( sx, sb[0], sc[0], rho1 );
1230 double rho2[2];
1231 Flatte_rhoab( sx, sb[1], sc[1], rho2 );
1232
1233 double gKK_a980 = 0.892 * 0.341, gPiEta_a980 = 0.341;
1234 double tmp1[2] = { gKK_a980, 0 };
1235 double tmp11[2];
1236 double tmp2[2] = { gPiEta_a980, 0 };
1237 double tmp22[2];
1238 Com_Multi( tmp1, rho1, tmp11 );
1239 Com_Multi( tmp2, rho2, tmp22 );
1240 double tmp3[2] = { tmp11[0] + tmp22[0], tmp11[1] + tmp22[1] };
1241 double tmp31[2];
1242 Com_Multi( tmp3, ci, tmp31 );
1243 double tmp4[2] = { mass * mass - sx - tmp31[0], -1.0 * tmp31[1] };
1244 Com_Divide( unit, tmp4, prop );
1245}
mass
double mass
Definition CosmicGenerator.cxx:128
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
g1
TF1 * g1
Definition DataBase/tau_mode.c:58
exp
EvtComplex exp(const EvtComplex &c)
Definition EvtComplex.hh:209
mass_Pion
const double mass_Pion
Definition EvtD0ToKSpipipi0.cc:1561
meta
double meta
Definition EvtPhokharaDef.hh:22
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
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
EvtDsToEta3pi::~EvtDsToEta3pi
virtual ~EvtDsToEta3pi()
Definition EvtDsToEta3pi.cc:33
EvtDsToEta3pi::clone
EvtDecayBase * clone()
Definition EvtDsToEta3pi.cc:37
EvtDsToEta3pi::getName
void getName(std::string &name)
Definition EvtDsToEta3pi.cc:35
EvtDsToEta3pi::decay
void decay(EvtParticle *p)
Definition EvtDsToEta3pi.cc:188
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