BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtD0ToKpipi0pi0.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: EvtD0ToKpipi0pi0.cc
11// the necessary file: EvtD0ToKpipi0pi0.hh
12//
13// Description: D0 -> K- pi+ pi0 pi0,
14// see PHYSICAL REVIEW D 99, 092008 (2019)
15//
16// Modification history:
17//
18// Liaoyuan Dong Jan.15, 2020 Module created
19//
20//------------------------------------------------------------------------
21#include "EvtD0ToKpipi0pi0.hh"
22#include "../EvtGenBase/EvtComplex.hh"
23#include "../EvtGenBase/EvtDecayTable.hh"
24#include "../EvtGenBase/EvtGenKine.hh"
25#include "../EvtGenBase/EvtPDL.hh"
26#include "../EvtGenBase/EvtParticle.hh"
27#include "../EvtGenBase/EvtPatches.hh"
28#include "../EvtGenBase/EvtReport.hh"
29#include <stdlib.h>
30
31EvtD0ToKpipi0pi0::~EvtD0ToKpipi0pi0() {}
32
33void EvtD0ToKpipi0pi0::getName( std::string& model_name ) { model_name = "D0ToKpipi0pi0"; }
34
35EvtDecayBase* EvtD0ToKpipi0pi0::clone() { return new EvtD0ToKpipi0pi0; }
36
37void EvtD0ToKpipi0pi0::init() {
38 checkNArg( 0 );
39 checkNDaug( 4 );
40 checkSpinParent( EvtSpinType::SCALAR );
41 /*
42 checkSpinDaughter(0,EvtSpinType::SCALAR);
43 checkSpinDaughter(1,EvtSpinType::SCALAR);
44 checkSpinDaughter(3,EvtSpinType::SCALAR);
45 checkSpinDaughter(4,EvtSpinType::SCALAR);
46 */
47 // std::cout << "Initializing EvtD0ToKpipi0pi0" << std::endl;
48
49 mod[0] = 1;
50 rho[0] = 2.02;
51 phi[0] = -0.75;
52 mod[1] = 1;
53 rho[1] = 1.66;
54 phi[1] = -2.90;
55 mod[2] = 0;
56 rho[2] = 0;
57 phi[2] = 0;
58 mod[3] = 0;
59 rho[3] = 0;
60 phi[3] = 0;
61 mod[4] = 0;
62 rho[4] = 0;
63 phi[4] = 0;
64 mod[5] = 0;
65 rho[5] = 0;
66 phi[5] = 0;
67 mod[6] = 0;
68 rho[6] = 0;
69 phi[6] = 0;
70 mod[7] = 1;
71 rho[7] = 1;
72 phi[7] = 0;
73 mod[8] = 1;
74 rho[8] = 0.842;
75 phi[8] = -2.05;
76 mod[9] = 1;
77 rho[9] = 0.0218;
78 phi[9] = 1.84;
79 mod[10] = 0;
80 rho[10] = 0;
81 phi[10] = 0;
82 mod[11] = 0;
83 rho[11] = 0;
84 phi[11] = 0;
85 mod[12] = 0;
86 rho[12] = 0;
87 phi[12] = 0;
88 mod[13] = 1;
89 rho[13] = 0.0336;
90 phi[13] = -1.55;
91 mod[14] = 1;
92 rho[14] = 0.109;
93 phi[14] = -1.35;
94 mod[15] = 1;
95 rho[15] = 0.196;
96 phi[15] = -2.07;
97 mod[16] = 0;
98 rho[16] = 0;
99 phi[16] = 0;
100 mod[17] = 0;
101 rho[17] = 0;
102 phi[17] = 0;
103 mod[18] = 0;
104 rho[18] = 0;
105 phi[18] = 0;
106 mod[19] = 1;
107 rho[19] = 0.363;
108 phi[19] = 1.93;
109 mod[20] = 0;
110 rho[20] = 0;
111 phi[20] = 0;
112 mod[21] = 0;
113 rho[21] = 0;
114 phi[21] = 0;
115 mod[22] = 0;
116 rho[22] = 0;
117 phi[22] = 0;
118 mod[23] = 1;
119 rho[23] = 0.555;
120 phi[23] = 0.44;
121 mod[24] = 1;
122 rho[24] = 0.526;
123 phi[24] = -1.84;
124 mod[25] = 0;
125 rho[25] = 0;
126 phi[25] = 0;
127 mod[26] = 1;
128 rho[26] = 1;
129 phi[26] = 0.64;
130 mod[27] = 0;
131 rho[27] = 0;
132 phi[27] = 0;
133 mod[28] = 0;
134 rho[28] = 0;
135 phi[28] = 0;
136 mod[29] = 1;
137 rho[29] = 3.34;
138 phi[29] = -0.02;
139 mod[30] = 0;
140 rho[30] = 0;
141 phi[30] = 0;
142 mod[31] = 0;
143 rho[31] = 0;
144 phi[31] = 0;
145 mod[32] = 0;
146 rho[32] = 0;
147 phi[32] = 0;
148 mod[33] = 0;
149 rho[33] = 0;
150 phi[33] = 0;
151 mod[34] = 1;
152 rho[34] = 1.76;
153 phi[34] = -2.39;
154 mod[35] = 1;
155 rho[35] = 0.175;
156 phi[35] = 1.59;
157 mod[36] = 1;
158 rho[36] = 0.397;
159 phi[36] = 1.45;
160 mod[37] = 0;
161 rho[37] = 0;
162 phi[37] = 0;
163 mod[38] = 0;
164 rho[38] = 0;
165 phi[38] = 0;
166 mod[39] = 0;
167 rho[39] = 0;
168 phi[39] = 0;
169 mod[40] = 0;
170 rho[40] = 0;
171 phi[40] = 0;
172 mod[41] = 1;
173 rho[41] = 1.02;
174 phi[41] = 0.52;
175 mod[42] = 0;
176 rho[42] = 0;
177 phi[42] = 0;
178 mod[43] = 0;
179 rho[43] = 0;
180 phi[43] = 0;
181 mod[44] = 0;
182 rho[44] = 0;
183 phi[44] = 0;
184 mod[45] = 0;
185 rho[45] = 0;
186 phi[45] = 0;
187 mod[46] = 1;
188 rho[46] = 0.146;
189 phi[46] = 1.24;
190 mod[47] = 1;
191 rho[47] = 0.0978;
192 phi[47] = -2.89;
193 mod[48] = 1;
194 rho[48] = 0.233;
195 phi[48] = 2.41;
196 mod[49] = 0;
197 rho[49] = 0;
198 phi[49] = 0;
199 mod[50] = 1;
200 rho[50] = 0.424;
201 phi[50] = -0.94;
202 mod[51] = 1;
203 rho[51] = 1.03;
204 phi[51] = -1.93;
205 mod[52] = 0;
206 rho[52] = 0;
207 phi[52] = 0;
208 mod[53] = 0;
209 rho[53] = 0;
210 phi[53] = 0;
211 mod[54] = 1;
212 rho[54] = 0.474;
213 phi[54] = -1.17;
214 mod[55] = 0;
215 rho[55] = 0;
216 phi[55] = 0;
217 mod[56] = 0;
218 rho[56] = 0;
219 phi[56] = 0;
220 mod[57] = 0;
221 rho[57] = 0;
222 phi[57] = 0;
223 mod[58] = 0;
224 rho[58] = 0;
225 phi[58] = 0;
226 mod[59] = 0;
227 rho[59] = 0;
228 phi[59] = 0;
229 mod[60] = 0;
230 rho[60] = 0;
231 phi[60] = 0;
232 mod[61] = 1;
233 rho[61] = 6.74;
234 phi[61] = -1.74;
235 mod[62] = 0;
236 rho[62] = 0;
237 phi[62] = 0;
238 mod[63] = 0;
239 rho[63] = 0;
240 phi[63] = 0;
241 mod[64] = 0;
242 rho[64] = 0;
243 phi[64] = 0;
244 mod[65] = 0;
245 rho[65] = 0;
246 phi[65] = 0;
247 mod[66] = 1;
248 rho[66] = 1.54;
249 phi[66] = -2.93;
250 mod[67] = 1;
251 rho[67] = 1.36;
252 phi[67] = 2.23;
253
254 mass[0] = 1.23;
255 width[0] = 0.50204;
256 mass[1] = 1.2723;
257 width[1] = 0.09;
258 mass[2] = 0.89166;
259 width[2] = 0.0508;
260 mass[3] = 0.89581;
261 width[3] = 0.0474;
262 mass[4] = 0.77511;
263 width[4] = 0.1491;
264
265 // for (int i=0; i<68; i++) {
266 // cout << i << " rho,phi = " << rho[i] << ", "<< phi[i] << endl;
267 // }
268 mD = 1.86486;
269 rRes = 3.0;
270 rD = 5.0;
271 metap = 0.95778;
272 mk0 = 0.497614;
273 mass_Kaon = 0.49368;
274 mass_Pion = 0.13957;
275 math_pi = 3.1415926;
276 mass_Pi0 = 0.1349766;
277 mkstrm = 0.89594;
278 mkstr0 = 0.89594;
279
280 pi = 3.1415926;
281 mpi = 0.13957;
282 g1 = 0.5468;
283 g2 = 0.23;
284
285 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
286 int EE[4][4][4][4] = {
287 { { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
288 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 1 }, { 0, 0, -1, 0 } },
289 { { 0, 0, 0, 0 }, { 0, 0, 0, -1 }, { 0, 0, 0, 0 }, { 0, 1, 0, 0 } },
290 { { 0, 0, 0, 0 }, { 0, 0, 1, 0 }, { 0, -1, 0, 0 }, { 0, 0, 0, 0 } } },
291 { { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, -1 }, { 0, 0, 1, 0 } },
292 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
293 { { 0, 0, 0, 1 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { -1, 0, 0, 0 } },
294 { { 0, 0, -1, 0 }, { 0, 0, 0, 0 }, { 1, 0, 0, 0 }, { 0, 0, 0, 0 } } },
295 { { { 0, 0, 0, 0 }, { 0, 0, 0, 1 }, { 0, 0, 0, 0 }, { 0, -1, 0, 0 } },
296 { { 0, 0, 0, -1 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 1, 0, 0, 0 } },
297 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
298 { { 0, 1, 0, 0 }, { -1, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } } },
299 { { { 0, 0, 0, 0 }, { 0, 0, -1, 0 }, { 0, 1, 0, 0 }, { 0, 0, 0, 0 } },
300 { { 0, 0, 1, 0 }, { 0, 0, 0, 0 }, { -1, 0, 0, 0 }, { 0, 0, 0, 0 } },
301 { { 0, -1, 0, 0 }, { 1, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
302 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } } } };
303 for ( int i = 0; i < 4; i++ )
304 {
305 for ( int j = 0; j < 4; j++ )
306 {
307 G[i][j] = GG[i][j];
308 for ( int k = 0; k < 4; k++ )
309 {
310 for ( int l = 0; l < 4; l++ ) { E[i][j][k][l] = EE[i][j][k][l]; }
311 }
312 }
313 }
314}
315
316void EvtD0ToKpipi0pi0::initProbMax() { setProbMax( 8060.0 ); }
317
318void EvtD0ToKpipi0pi0::decay( EvtParticle* p ) {
319 /*
320 double maxprob = 0.0;
321 for(int ir=0;ir<=60000000;ir++){
322 p->initializePhaseSpace(getNDaug(),getDaugs());
323 EvtVector4R Km0 = p->getDaug(0)->getP4();
324 EvtVector4R pi1 = p->getDaug(1)->getP4();
325 EvtVector4R pi2 = p->getDaug(2)->getP4();
326 EvtVector4R pi3 = p->getDaug(3)->getP4();
327 double Km[4],Pip[4],Pi01[4],Pi02[4];
328 Km[0] = Km0.get(0); Pip[0] = pi1.get(0); Pi01[0] = pi2.get(0); Pi02[0] = pi3.get(0);
329 Km[1] = Km0.get(1); Pip[1] = pi1.get(1); Pi01[1] = pi2.get(1); Pi02[1] = pi3.get(1);
330 Km[2] = Km0.get(2); Pip[2] = pi1.get(2); Pi01[2] = pi2.get(2); Pi02[2] = pi3.get(2);
331 Km[3] = Km0.get(3); Pip[3] = pi1.get(3); Pi01[3] = pi2.get(3); Pi02[3] = pi3.get(3);
332 double Prob = calPDF(Km, Pip, Pi01, Pi02);
333 if(Prob>maxprob) {
334 maxprob=Prob;
335 std::cout << "Max PDF = " << ir << " prob= " << Prob << std::endl;
336 }
337 }
338 std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
339 */
340 p->initializePhaseSpace( getNDaug(), getDaugs() );
341 EvtVector4R Km0 = p->getDaug( 0 )->getP4();
342 EvtVector4R pi1 = p->getDaug( 1 )->getP4();
343 EvtVector4R pi2 = p->getDaug( 2 )->getP4();
344 EvtVector4R pi3 = p->getDaug( 3 )->getP4();
345
346 double Km[4], Pip[4], Pi01[4], Pi02[4];
347 Km[0] = Km0.get( 0 );
348 Pip[0] = pi1.get( 0 );
349 Pi01[0] = pi2.get( 0 );
350 Pi02[0] = pi3.get( 0 );
351 Km[1] = Km0.get( 1 );
352 Pip[1] = pi1.get( 1 );
353 Pi01[1] = pi2.get( 1 );
354 Pi02[1] = pi3.get( 1 );
355 Km[2] = Km0.get( 2 );
356 Pip[2] = pi1.get( 2 );
357 Pi01[2] = pi2.get( 2 );
358 Pi02[2] = pi3.get( 2 );
359 Km[3] = Km0.get( 3 );
360 Pip[3] = pi1.get( 3 );
361 Pi01[3] = pi2.get( 3 );
362 Pi02[3] = pi3.get( 3 );
363 double prob = calPDF( Km, Pip, Pi01, Pi02 );
364 setProb( prob );
365 return;
366}
367
368double EvtD0ToKpipi0pi0::calPDF( double Km[], double Pip[], double Pi01[], double Pi02[] ) {
369 Km[0] = sqrt( mass_Kaon * mass_Kaon + Km[1] * Km[1] + Km[2] * Km[2] + Km[3] * Km[3] );
370 Pip[0] = sqrt( mass_Pion * mass_Pion + Pip[1] * Pip[1] + Pip[2] * Pip[2] + Pip[3] * Pip[3] );
371 Pi01[0] =
372 sqrt( mass_Pi0 * mass_Pi0 + Pi01[1] * Pi01[1] + Pi01[2] * Pi01[2] + Pi01[3] * Pi01[3] );
373 Pi02[0] =
374 sqrt( mass_Pi0 * mass_Pi0 + Pi02[1] * Pi02[1] + Pi02[2] * Pi02[2] + Pi02[3] * Pi02[3] );
375
376 EvtComplex PDF[68];
377 int g[3];
378 //-------PHSP---------
379 PDF[0] = PHSP( Km, Pip ) + PHSP( Km, Pip );
380 PDF[1] = PHSP( Km, Pi01 ) + PHSP( Km, Pi02 );
381 //-------D2PP,P2VP------
382 // PDF[2] = D2PP_P2VP(Pi01,Pi02,Km,Pip,0) + D2PP_P2VP(Pi02,Pi01,Km,Pip,0);
383 // PDF[3] = D2PP_P2VP(Pi01,Pip,Km,Pi02,10) + D2PP_P2VP(Pi02,Pip,Km,Pi01,10);
384 // PDF[4] = D2PP_P2VP(Pip,Pi01,Km,Pi02,20) + D2PP_P2VP(Pip,Pi02,Km,Pi01,20);
385 // PDF[5] = D2PP_P2VP(Pi01,Km,Pip,Pi02,1) + D2PP_P2VP(Pi02,Km,Pip,Pi01,1);
386 // PDF[6] = D2PP_P2VP(Km,Pi01,Pip,Pi02,11) + D2PP_P2VP(Km,Pi02,Pip,Pi01,11);
387 //----------D2AP,A2VP--------------
388 g[0] = 1;
389 g[1] = 1;
390 g[2] = 0;
391 PDF[7] = D2AP_A2VP( Km, Pi01, Pip, Pi02, g, 0 ) + D2AP_A2VP( Km, Pi02, Pip, Pi01, g, 0 );
392 g[2] = 2;
393 PDF[8] = D2AP_A2VP( Km, Pi01, Pip, Pi02, g, 0 ) + D2AP_A2VP( Km, Pi02, Pip, Pi01, g, 0 );
394 g[2] = 0;
395 PDF[9] = D2AP_A2VP( Pip, Pi01, Km, Pi02, g, 1 ) + D2AP_A2VP( Pip, Pi02, Km, Pi01, g, 1 );
396 // g[2] = 2;
397 // PDF[10] = D2AP_A2VP(Pip,Pi01,Km,Pi02,g,1) + D2AP_A2VP(Pip,Pi02,Km,Pi01,g,1);
398 // g[2] = 0;
399 // PDF[11] = D2AP_A2VP(Pi01,Pip,Km,Pi02,g,21) + D2AP_A2VP(Pi02,Pip,Km,Pi01,g,21);
400 // g[2] = 2;
401 // PDF[12] = D2AP_A2VP(Pi01,Pip,Km,Pi02,g,21) + D2AP_A2VP(Pi02,Pip,Km,Pi01,g,21);
402 g[2] = 0;
403 PDF[13] = D2AP_A2VP( Pi01, Pi02, Km, Pip, g, 31 ) + D2AP_A2VP( Pi02, Pi01, Km, Pip, g, 31 );
404 g[2] = 2;
405 PDF[14] = D2AP_A2VP( Pi01, Pi02, Km, Pip, g, 31 ) + D2AP_A2VP( Pi02, Pi01, Km, Pip, g, 31 );
406 g[2] = 0;
407 PDF[15] = D2AP_A2VP( Pi01, Km, Pip, Pi02, g, 3 ) + D2AP_A2VP( Pi02, Km, Pip, Pi01, g, 3 );
408 // g[2] = 2;
409 // PDF[16] = D2AP_A2VP(Pi01,Km,Pip,Pi02,g,3) + D2AP_A2VP(Pi02,Km,Pip,Pi01,g,3);
410 g[0] = 1;
411 g[1] = 0;
412 g[2] = 0;
413 // PDF[17] = D2AP_A2VP(Km,Pi01,Pip,Pi02,g,0) + D2AP_A2VP(Km,Pi02,Pip,Pi01,g,0);
414 // g[2] = 2;
415 // PDF[18] = D2AP_A2VP(Km,Pi01,Pip,Pi02,g,0) + D2AP_A2VP(Km,Pi02,Pip,Pi01,g,0);
416 g[2] = 0;
417 PDF[19] = D2AP_A2VP( Pip, Pi01, Km, Pi02, g, 1 ) + D2AP_A2VP( Pip, Pi02, Km, Pi01, g, 1 );
418 // g[2] = 2;
419 // PDF[20] = D2AP_A2VP(Pip,Pi01,Km,Pi02,g,1) + D2AP_A2VP(Pip,Pi02,Km,Pi01,g,1);
420 // g[2] = 0;
421 // PDF[21] = D2AP_A2VP(Pi01,Pip,Km,Pi02,g,21) + D2AP_A2VP(Pi02,Pip,Km,Pi01,g,21);
422 // g[2] = 2;
423 // PDF[22] = D2AP_A2VP(Pi01,Pip,Km,Pi02,g,21) + D2AP_A2VP(Pi02,Pip,Km,Pi01,g,21);
424 g[2] = 0;
425 PDF[23] = D2AP_A2VP( Pi01, Pi02, Km, Pip, g, 31 ) + D2AP_A2VP( Pi02, Pi01, Km, Pip, g, 31 );
426 g[2] = 2;
427 PDF[24] = D2AP_A2VP( Pi01, Pi02, Km, Pip, g, 31 ) + D2AP_A2VP( Pi02, Pi01, Km, Pip, g, 31 );
428 // g[2] = 0;
429 // PDF[25] = D2AP_A2VP(Pi01,Km,Pip,Pi02,g,3) + D2AP_A2VP(Pi02,Km,Pip,Pi01,g,3);
430 g[2] = 2;
431 PDF[26] = D2AP_A2VP( Pi01, Km, Pip, Pi02, g, 3 ) + D2AP_A2VP( Pi02, Km, Pip, Pi01, g, 3 );
432 //--------D2AP,A2SP-----------------------------------
433 // PDF[27] = D2AP_A2SP(Km,Pi01,Pip,Pi02,0) + D2AP_A2SP(Km,Pi02,Pip,Pi01,0);
434 // PDF[28] = D2AP_A2SP(Km,Pip,Pi01,Pi02,10) + D2AP_A2SP(Km,Pip,Pi02,Pi01,10);
435 PDF[29] = D2AP_A2SP( Pi01, Pi02, Km, Pip, 1 ) + D2AP_A2SP( Pi02, Pi01, Km, Pip, 1 );
436 // PDF[30] = D2AP_A2SP(Pi01,Pip,Km,Pi02,11) + D2AP_A2SP(Pi02,Pip,Km,Pi01,11);
437 // PDF[31] = D2AP_A2SP(Pip,Pi01,Km,Pi02,21) + D2AP_A2SP(Pip,Pi02,Km,Pi01,21);
438 // PDF[32] = D2AP_A2SP(Pi01,Km,Pip,Pi02,31) + D2AP_A2SP(Pi02,Km,Pip,Pi01,31);
439 // PDF[33] = D2AP_A2SP(Pip,Km,Pi01,Pi02,41) + D2AP_A2SP(Pip,Km,Pi02,Pi01,41);
440 //--------D2VS----------------------------
441 PDF[34] = D2VS( Pip, Pi01, Km, Pi02, 1, 0 ) + D2VS( Pip, Pi02, Km, Pi01, 1, 0 );
442 PDF[35] = D2VS( Km, Pi01, Pip, Pi02, 1, 1 ) + D2VS( Km, Pi02, Pip, Pi01, 1, 1 );
443 PDF[36] = D2VS( Km, Pip, Pi01, Pi02, 1, 11 ) + D2VS( Km, Pip, Pi02, Pi01, 1, 11 );
444 // PDF[37] = D2VS(Pi01,Pi02,Km,Pip,0,10) + D2VS(Pi02,Pi01,Km,Pip,0,10);
445 // PDF[38] = D2VS(Pip,Pi01,Km,Pi02,0,0) + D2VS(Pip,Pi02,Km,Pi01,0,0);
446 // PDF[39] = D2VS(Km,Pip,Pi01,Pi02,0,11) + D2VS(Km,Pip,Pi02,Pi01,0,11);
447 // PDF[40] = D2VS(Km,Pi01,Pip,Pi02,0,1) + D2VS(Km,Pi02,Pip,Pi01,0,1);
448 //---------D2VP,V2VP----------------------
449 PDF[41] = D2VP_V2VP( Pi01, Pip, Km, Pi02, 0 ) + D2VP_V2VP( Pi02, Pip, Km, Pi01, 0 );
450 // PDF[42] = D2VP_V2VP(Pip,Pi01,Km,Pi02,10) + D2VP_V2VP(Pip,Pi02,Km,Pi01,10);
451 // PDF[43] = D2VP_V2VP(Pi01,Pi02,Km,Pip,20) + D2VP_V2VP(Pi02,Pi01,Km,Pip,20);
452 // PDF[44] = D2VP_V2VP(Pi01,Km,Pip,Pi02,1) + D2VP_V2VP(Pi02,Km,Pip,Pi01,1);
453 // PDF[45] = D2VP_V2VP(Km,Pi01,Pip,Pi02,2) + D2VP_V2VP(Km,Pi02,Pip,Pi01,2);
454 //-----------D2VV--------------------------
455 g[0] = 1;
456 g[1] = 1;
457 g[2] = 0;
458 PDF[46] = D2VV( Km, Pi01, Pip, Pi02, g, 0 ) + D2VV( Km, Pi02, Pip, Pi01, g, 0 );
459 g[2] = 1;
460 PDF[47] = D2VV( Km, Pi01, Pip, Pi02, g, 0 ) + D2VV( Km, Pi02, Pip, Pi01, g, 0 );
461 g[2] = 2;
462 PDF[48] = D2VV( Km, Pi01, Pip, Pi02, g, 0 ) + D2VV( Km, Pi02, Pip, Pi01, g, 0 );
463 g[0] = 0;
464 g[1] = 1;
465 g[2] = 0;
466 // PDF[49] = D2VV(Km,Pi01,Pip,Pi02,g,0) + D2VV(Km,Pi02,Pip,Pi01,g,0);
467 g[2] = 1;
468 PDF[50] = D2VV( Km, Pi01, Pip, Pi02, g, 0 ) + D2VV( Km, Pi02, Pip, Pi01, g, 0 );
469 g[2] = 2;
470 PDF[51] = D2VV( Km, Pi01, Pip, Pi02, g, 0 ) + D2VV( Km, Pi02, Pip, Pi01, g, 0 );
471 g[0] = 1;
472 g[1] = 0;
473 g[2] = 0;
474 // PDF[52] = D2VV(Km,Pi01,Pip,Pi02,g,0) + D2VV(Km,Pi02,Pip,Pi01,g,0);
475 // g[2] = 1;
476 // PDF[53] = D2VV(Km,Pi01,Pip,Pi02,g,0) + D2VV(Km,Pi02,Pip,Pi01,g,0);
477 g[2] = 2;
478 PDF[54] = D2VV( Km, Pi01, Pip, Pi02, g, 0 ) + D2VV( Km, Pi02, Pip, Pi01, g, 0 );
479 g[0] = 1;
480 g[1] = 0;
481 g[2] = 0;
482 // PDF[55] = D2VV(Km,Pip,Pi01,Pi01,g,1) + D2VV(Km,Pip,Pi02,Pi01,g,1);
483 // g[2] = 1;
484 // PDF[56] = D2VV(Km,Pip,Pi01,Pi01,g,1) + D2VV(Km,Pip,Pi02,Pi01,g,1);
485 // g[2] = 2;
486 // PDF[57] = D2VV(Km,Pip,Pi01,Pi01,g,1) + D2VV(Km,Pip,Pi02,Pi01,g,1);
487 g[0] = 0;
488 g[1] = 0;
489 g[2] = 0;
490 // PDF[58] = D2VV(Km,Pip,Pi01,Pi01,g,1) + D2VV(Km,Pip,Pi02,Pi01,g,1);
491 // g[2] = 1;
492 // PDF[59] = D2VV(Km,Pip,Pi01,Pi01,g,1) + D2VV(Km,Pip,Pi02,Pi01,g,1);
493 // g[2] = 2;
494 // PDF[60] = D2VV(Km,Pip,Pi01,Pi01,g,1) + D2VV(Km,Pip,Pi02,Pi01,g,1);
495 g[0] = 0;
496 g[1] = 0;
497 g[2] = 0;
498 PDF[61] = D2VV( Km, Pi01, Pip, Pi02, g, 0 ) + D2VV( Km, Pi02, Pip, Pi01, g, 0 );
499 g[2] = 1;
500 // PDF[62] = D2VV(Km,Pi01,Pip,Pi02,g,0) + D2VV(Km,Pi02,Pip,Pi01,g,0);
501 g[2] = 2;
502 // PDF[63] = D2VV(Km,Pi01,Pip,Pi02,g,0) + D2VV(Km,Pi02,Pip,Pi01,g,0);
503 //----------D2TS--------------------
504 // PDF[64] = D2TS(Km,Pip,Pi01,Pi02,0) + D2TS(Km,Pip,Pi02,Pi01,0);
505 // PDF[65] = D2TS(Km,Pi01,Pip,Pi02,10) + D2TS(Km,Pi02,Pip,Pi01,10);
506 PDF[66] = D2TS( Pi02, Pi01, Km, Pip, 1 ) + D2TS( Pi01, Pi02, Km, Pip, 1 );
507 PDF[67] = D2TS( Pip, Pi01, Km, Pi02, 11 ) + D2TS( Pip, Pi02, Km, Pi01, 11 );
508
509 //------------------------------------------
510 EvtComplex cof;
511 EvtComplex pdf, module;
512 pdf = EvtComplex( 0, 0 );
513 for ( int i = 0; i < 68; i++ )
514 {
515 if ( mod[i] == 0 ) continue;
516 cof = EvtComplex( rho[i] * cos( phi[i] ), rho[i] * sin( phi[i] ) );
517 pdf = pdf + cof * PDF[i];
518 }
519 module = conj( pdf ) * pdf;
520 double value;
521 value = real( module );
522 return ( value <= 0 ) ? 1e-20 : value;
523}
524EvtComplex EvtD0ToKpipi0pi0::KPiSFormfactor( double sa, double sb, double sc, double r ) {
525 double m1430 = 1.463;
526 double sa0 = m1430 * m1430;
527 double w1430 = 0.233;
528 double q0 = ( sa0 + sb - sc ) * ( sa0 + sb - sc ) / ( 4 * sa0 ) - sb;
529 if ( q0 < 0 ) q0 = 1e-16;
530 double qs = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
531 double q = sqrt( qs );
532 double width = w1430 * q * m1430 / sqrt( sa * q0 );
533 double temp_R = atan( m1430 * width / ( sa0 - sa ) );
534 if ( temp_R < 0 ) temp_R += math_pi;
535 double deltaR = -5.31 + temp_R;
536 double temp_F = atan( 2 * 1.07 * q / ( 2 + ( -1.8 ) * 1.07 * qs ) );
537 if ( temp_F < 0 ) temp_F += math_pi;
538 double deltaF = 2.33 + temp_F;
539 EvtComplex cR( cos( deltaR ), sin( deltaR ) );
540 EvtComplex cF( cos( deltaF ), sin( deltaF ) );
541 EvtComplex amp = 0.8 * sin( deltaF ) * cF + sin( deltaR ) * cR * cF * cF;
542 return amp;
543}
544EvtComplex EvtD0ToKpipi0pi0::D2VV( double P1[], double P2[], double P3[], double P4[], int g[],
545 int flag ) {
546 double t1V1[4], t1V2[4], t1D[4], t2D[4][4];
547 double temp_PDF = 0;
548 EvtComplex amp_PDF( 0, 0 );
549 EvtComplex pro[2];
550
551 double sa[3], sb[3], sc[3], B[3];
552 double pV1[4], pV2[4], pD[4];
553 for ( int i = 0; i != 4; i++ )
554 {
555 pV1[i] = P1[i] + P2[i];
556 pV2[i] = P3[i] + P4[i];
557 pD[i] = pV1[i] + pV2[i];
558 }
559 sa[0] = dot( pV1, pV1 );
560 sb[0] = dot( P1, P1 );
561 sc[0] = dot( P2, P2 );
562 sa[1] = dot( pV2, pV2 );
563 sb[1] = dot( P3, P3 );
564 sc[1] = dot( P4, P4 );
565 sa[2] = dot( pD, pD );
566 sb[2] = sa[0];
567 sc[2] = sa[1];
568 if ( g[0] == 1 )
569 {
570 if ( flag == 0 ) pro[0] = propagatorRBW( mass[2], width[2], sa[0], sb[0], sc[0], rRes, 1 );
571 if ( flag == 1 ) pro[0] = propagatorRBW( mass[3], width[3], sa[0], sb[0], sc[0], rRes, 1 );
572 }
573 if ( g[1] == 1 )
574 {
575 if ( flag == 0 )
576 pro[1] = propagatorGS( mass[4], width[4], sa[1], sb[1], sc[1], rRes, 1 ); // rho
577 if ( flag == 1 ) pro[1] = 1;
578 }
579 if ( g[0] == 0 ) pro[0] = 1;
580 if ( g[1] == 0 ) pro[1] = 1;
581 B[0] = barrier( 1, sa[0], sb[0], sc[0], rRes );
582 B[1] = barrier( 1, sa[1], sb[1], sc[1], rRes );
583 calt1( P1, P2, t1V1 );
584 calt1( P3, P4, t1V2 );
585 if ( g[2] == 0 )
586 {
587 for ( int i = 0; i != 4; i++ ) { temp_PDF += ( G[i][i] ) * t1V1[i] * t1V2[i]; }
588 B[2] = 1;
589 }
590 if ( g[2] == 1 )
591 {
592 calt1( pV1, pV2, t1D );
593 for ( int i = 0; i != 4; i++ )
594 {
595 for ( int j = 0; j != 4; j++ )
596 {
597 for ( int k = 0; k != 4; k++ )
598 {
599 for ( int l = 0; l != 4; l++ )
600 {
601 temp_PDF += E[i][j][k][l] * pD[i] * t1D[j] * t1V1[k] * t1V2[l] * ( G[i][i] ) *
602 ( G[j][j] ) * ( G[l][l] ) * ( G[k][k] );
603 }
604 }
605 }
606 }
607 B[2] = barrier( 1, sa[2], sb[2], sc[2], rD );
608 }
609 if ( g[2] == 2 )
610 {
611 calt2( pV1, pV2, t2D );
612 for ( int i = 0; i != 4; i++ )
613 {
614 for ( int j = 0; j != 4; j++ )
615 { temp_PDF += t2D[i][j] * t1V1[i] * t1V2[j] * ( G[i][i] ) * ( G[j][j] ); }
616 }
617 B[2] = barrier( 2, sa[2], sb[2], sc[2], rD );
618 }
619 amp_PDF = temp_PDF * B[0] * B[1] * B[2] * pro[0] * pro[1];
620 return amp_PDF;
621}
622
623EvtComplex EvtD0ToKpipi0pi0::D2AP_A2VP( double P1[], double P2[], double P3[], double P4[],
624 int g[], int flag ) {
625 double temp_PDF = 0;
626 EvtComplex amp_PDF( 0, 0 );
627 EvtComplex pro[2];
628 double t1V[4], t1D[4], t2A[4][4];
629 double sa[3], sb[3], sc[3], B[3];
630 double pV[4], pA[4], pD[4];
631 for ( int i = 0; i != 4; i++ )
632 {
633 pV[i] = P3[i] + P4[i];
634 pA[i] = pV[i] + P2[i];
635 pD[i] = pA[i] + P1[i];
636 }
637 sa[0] = dot( pV, pV );
638 sb[0] = dot( P3, P3 );
639 sc[0] = dot( P4, P4 );
640 sa[1] = dot( pA, pA );
641 sb[1] = sa[0];
642 sc[1] = dot( P2, P2 );
643 sa[2] = dot( pD, pD );
644 sb[2] = sa[1];
645 sc[2] = dot( P1, P1 );
646 if ( g[0] == 1 )
647 {
648 if ( flag == 0 || flag == 3 )
649 pro[0] = propagatorGS( mass[4], width[4], sa[0], sb[0], sc[0], rRes, 1 );
650 else if ( flag == 1 || flag == 21 )
651 pro[0] = propagatorRBW( mass[2], width[2], sa[0], sb[0], sc[0], rRes, 1 );
652 else if ( flag == 31 )
653 pro[0] = propagatorRBW( mass[3], width[3], sa[0], sb[0], sc[0], rRes, 1 );
654 }
655 else if ( g[0] == 0 ) pro[0] = 1;
656 if ( g[1] == 1 )
657 {
658 if ( flag == 0 )
659 pro[1] = propagatorRBW( mass[0], width[0], sa[1], sb[1], sc[1], rRes, g[2] );
660 if ( flag == 1 || flag == 21 || flag == 31 || flag == 3 )
661 pro[1] = propagatorRBW( mass[1], width[1], sa[1], sb[1], sc[1], rRes, g[2] );
662 }
663 else if ( g[1] == 0 ) pro[1] = 1;
664 B[0] = barrier( 1, sa[0], sb[0], sc[0], rRes );
665 B[2] = barrier( 1, sa[2], sb[2], sc[2], rD );
666 calt1( P3, P4, t1V );
667 calt1( pA, P1, t1D );
668 if ( g[2] == 0 )
669 {
670 for ( int i = 0; i != 4; i++ )
671 {
672 for ( int j = 0; j != 4; j++ )
673 {
674 temp_PDF +=
675 t1D[i] * ( G[i][j] - pA[i] * pA[j] / sa[1] ) * t1V[j] * ( G[i][i] ) * ( G[j][j] );
676 }
677 }
678 B[1] = 1;
679 }
680 else if ( g[2] == 2 )
681 {
682 calt2( pV, P2, t2A );
683 for ( int i = 0; i != 4; i++ )
684 {
685 for ( int j = 0; j != 4; j++ )
686 { temp_PDF += t1D[i] * t2A[i][j] * t1V[j] * ( G[i][i] ) * ( G[j][j] ); }
687 }
688 B[1] = barrier( 2, sa[1], sb[1], sc[1], rRes );
689 }
690 amp_PDF = temp_PDF * B[0] * B[1] * B[2] * pro[0] * pro[1];
691 return amp_PDF;
692}
693
694EvtComplex EvtD0ToKpipi0pi0::D2AP_A2SP( double P1[], double P2[], double P3[], double P4[],
695 int flag ) {
696 // flag = 0, S = rho; flag = 1, S = K*
697 double temp_PDF = 0;
698 EvtComplex amp_PDF( 0, 0 );
699 EvtComplex pro;
700 double sa[3], sb[3], sc[3], B[3];
701 double t1D[4], t1A[4];
702 double pS[4], pA[4], pD[4];
703 for ( int i = 0; i != 4; i++ )
704 {
705 pS[i] = P3[i] + P4[i];
706 pA[i] = pS[i] + P2[i];
707 pD[i] = pA[i] + P1[i];
708 }
709 sa[0] = dot( pS, pS );
710 sb[0] = dot( P3, P3 );
711 sc[0] = dot( P4, P4 );
712 sa[1] = dot( pA, pA );
713 sb[1] = sa[0];
714 sc[1] = dot( P2, P2 );
715 sa[2] = dot( pD, pD );
716 sb[2] = sa[1];
717 sc[2] = dot( P1, P1 );
718 B[1] = barrier( 1, sa[1], sb[1], sc[1], rRes );
719 B[2] = barrier( 1, sa[2], sb[2], sc[2], rD );
720 calt1( pA, P1, t1D );
721 calt1( pS, P2, t1A );
722 for ( int i = 0; i != 4; i++ ) { temp_PDF += t1D[i] * t1A[i] * ( G[i][i] ); }
723 amp_PDF = temp_PDF * B[1] * B[2];
724 if ( flag == 1 || flag == 11 || flag == 21 )
725 amp_PDF = amp_PDF * KPiSFormfactor( sa[0], sb[0], sc[0], rRes );
726 return amp_PDF;
727}
728
729EvtComplex EvtD0ToKpipi0pi0::D2PP_P2VP( double P1[], double P2[], double P3[], double P4[],
730 int flag ) {
731 // modeidx = 0 :(K*0 pi0)pi0
732 // modeidx = 10:(K*- pi+)pi0
733 // modeidx = 20:(K*- pi0)pi+
734 // modeidx = 1 :(K- rho+)pi0
735 // modeidx = 11:K-(rho+ pi0)
736 double temp_PDF = 0;
737 EvtComplex amp( 0, 0 );
738 EvtComplex prop;
739 double sa[3], sb[3], sc[3], B[3];
740 double t1V[4];
741 double pV[4], pP[4], pD[4];
742 for ( int i = 0; i != 4; i++ )
743 {
744 pV[i] = P3[i] + P4[i];
745 pP[i] = pV[i] + P2[i];
746 pD[i] = pP[i] + P1[i];
747 }
748 sa[0] = dot( pV, pV );
749 sb[0] = dot( P3, P3 );
750 sc[0] = dot( P4, P4 );
751 sa[1] = dot( pP, pP );
752 sb[1] = sa[0];
753 sc[1] = dot( P2, P2 );
754 sa[2] = dot( pD, pD );
755 sb[2] = sa[1];
756 sc[2] = dot( P1, P1 );
757 B[0] = barrier( 1, sa[0], sb[0], sc[0], rRes );
758 B[1] = barrier( 1, sa[1], sb[1], sc[1], rRes );
759 if ( flag == 0 ) prop = propagatorRBW( mass[3], width[3], sa[0], sb[0], sc[0], rRes, 1 );
760 else if ( flag == 10 || 20 )
761 prop = propagatorRBW( mass[2], width[2], sa[0], sb[0], sc[0], rRes, 1 );
762 else if ( flag == 1 || 11 )
763 prop = propagatorGS( mass[4], width[4], sa[0], sb[0], sc[0], rRes, 1 );
764 calt1( P3, P4, t1V );
765 for ( int i = 0; i != 4; i++ ) { temp_PDF += P2[i] * t1V[i] * ( G[i][i] ); }
766 amp = temp_PDF * B[0] * B[1] * prop;
767 return amp;
768}
769
770EvtComplex EvtD0ToKpipi0pi0::D2VP_V2VP( double P1[], double P2[], double P3[], double P4[],
771 int flag ) {
772 double temp_PDF = 0;
773 EvtComplex amp_PDF( 0, 0 );
774 EvtComplex pro;
775 double sa[3], sb[3], sc[3], B[3];
776 double pV1[4], pV2[4], qV1[4], qV2[4], pD[4];
777 for ( int i = 0; i != 4; i++ )
778 {
779 pV2[i] = P3[i] + P4[i];
780 qV2[i] = P3[i] - P4[i];
781 pV1[i] = pV2[i] + P2[i];
782 qV1[i] = pV2[i] - P2[i];
783 pD[i] = pV1[i] + P1[i];
784 }
785 for ( int i = 0; i != 4; i++ )
786 {
787 for ( int j = 0; j != 4; j++ )
788 {
789 for ( int k = 0; k != 4; k++ )
790 {
791 for ( int l = 0; l != 4; l++ )
792 {
793 temp_PDF += E[i][j][k][l] * pV1[i] * qV1[j] * P1[k] * qV2[l] * ( G[i][i] ) *
794 ( G[j][j] ) * ( G[k][k] ) * ( G[l][l] );
795 }
796 }
797 }
798 }
799 sa[0] = dot( pV2, pV2 );
800 sb[0] = dot( P3, P3 );
801 sc[0] = dot( P4, P4 );
802 sa[1] = dot( pV1, pV1 );
803 sb[1] = sa[0];
804 sc[1] = dot( P2, P2 );
805 sa[2] = dot( pD, pD );
806 sb[2] = sa[1];
807 sc[2] = dot( P1, P1 );
808 if ( flag == 0 || flag == 10 )
809 pro = propagatorRBW( mass[2], width[2], sa[0], sb[0], sc[0], rRes, 1 );
810 else if ( flag == 20 )
811 pro = propagatorRBW( mass[3], width[3], sa[0], sb[0], sc[0], rRes, 1 );
812 else if ( flag == 1 || flag == 2 )
813 pro = propagatorGS( mass[4], width[4], sa[0], sb[0], sc[0], rRes, 1 );
814 B[0] = barrier( 1, sa[0], sb[0], sc[0], rRes );
815 B[1] = barrier( 1, sa[1], sb[1], sc[1], rRes );
816 B[2] = barrier( 1, sa[2], sb[2], sc[2], rD );
817 amp_PDF = temp_PDF * B[0] * B[1] * B[2] * pro;
818 return amp_PDF;
819}
820
821EvtComplex EvtD0ToKpipi0pi0::D2VS( double P1[], double P2[], double P3[], double P4[], int g,
822 int flag ) {
823 double temp_PDF = 0;
824 EvtComplex amp_PDF( 0, 0 );
825 EvtComplex pro;
826 double sa[3], sb[3], sc[3], B[3];
827 double t1D[4], t1V[4];
828 double pS[4], pV[4], pD[4];
829 for ( int i = 0; i != 4; i++ )
830 {
831 pS[i] = P3[i] + P4[i];
832 pV[i] = P1[i] + P2[i];
833 pD[i] = pS[i] + pV[i];
834 }
835 sa[0] = dot( pS, pS );
836 sb[0] = dot( P3, P3 );
837 sc[0] = dot( P4, P4 );
838 sa[1] = dot( pV, pV );
839 sb[1] = dot( P1, P1 );
840 sc[1] = dot( P2, P2 );
841 sa[2] = dot( pD, pD );
842 sb[2] = sa[0];
843 sc[2] = sa[1];
844 if ( g == 1 )
845 {
846 if ( flag == 0 ) pro = propagatorGS( mass[4], width[4], sa[1], sb[1], sc[1], rRes, 1 );
847 else if ( flag == 1 )
848 pro = propagatorRBW( mass[2], width[2], sa[1], sb[1], sc[1], rRes, 1 );
849 else if ( flag == 11 )
850 pro = propagatorRBW( mass[3], width[3], sa[1], sb[1], sc[1], rRes, 1 );
851 else if ( flag == 10 ) pro = 1;
852 }
853 else if ( g == 0 ) pro = 1;
854 B[1] = barrier( 1, sa[1], sb[1], sc[1], rRes );
855 B[2] = barrier( 1, sa[2], sb[2], sc[2], rD );
856 calt1( P1, P2, t1V );
857 calt1( pS, pV, t1D );
858 for ( int i = 0; i != 4; i++ ) { temp_PDF += G[i][i] * t1D[i] * t1V[i]; }
859 amp_PDF = temp_PDF * B[1] * B[2] * pro;
860 if ( flag == 0 || flag == 10 ) amp_PDF *= KPiSFormfactor( sa[0], sb[0], sc[0], rRes );
861 // if(modeidx == 1 || modeidx == 11) amp_PDF *= -1.0; /// why ???????
862 return amp_PDF;
863}
864
865EvtComplex EvtD0ToKpipi0pi0::D2TS( double P1[], double P2[], double P3[], double P4[],
866 int flag ) {
867 // flag == 0 KPiT. 1 PiPiT
868 double temp_PDF = 0;
869 EvtComplex amp_PDF( 0, 0 );
870 double sa[3], sb[3], sc[3], B[3];
871 double t2D[4][4], t2T[4][4];
872 double pS[4], pT[4], pD[4];
873 for ( int i = 0; i != 4; i++ )
874 {
875 pS[i] = P3[i] + P4[i];
876 pT[i] = P1[i] + P2[i];
877 pD[i] = pT[i] + pS[i];
878 }
879 sa[0] = dot( pT, pT );
880 sb[0] = dot( P1, P1 );
881 sc[0] = dot( P2, P2 );
882 sa[1] = dot( pS, pS );
883 sb[1] = dot( P3, P3 );
884 sc[1] = dot( P4, P4 );
885 sa[2] = dot( pD, pD );
886 sb[2] = sa[0];
887 sc[2] = sa[1];
888 B[0] = barrier( 2, sa[0], sb[0], sc[0], rRes );
889 B[2] = barrier( 2, sa[2], sb[2], sc[2], rD );
890 calt2( P1, P2, t2T );
891 calt2( pT, pS, t2D );
892 for ( int i = 0; i != 4; i++ )
893 {
894 for ( int j = 0; j != 4; j++ )
895 { temp_PDF += t2D[i][j] * t2T[j][i] * ( G[i][i] ) * ( G[j][j] ); }
896 }
897 amp_PDF = temp_PDF * B[0] * B[2];
898 if ( flag == 1 || flag == 11 )
899 { amp_PDF = amp_PDF * KPiSFormfactor( sa[1], sb[1], sc[1], rRes ); }
900 return amp_PDF;
901}
902
903EvtComplex EvtD0ToKpipi0pi0::PHSP( double P1[], double P2[] ) {
904 EvtComplex amp_PDF( 0, 0 );
905 double sa, sb, sc;
906 double KPi[4];
907 for ( int i = 0; i != 4; i++ ) { KPi[i] = P1[i] + P2[i]; }
908 sa = dot( KPi, KPi );
909 sb = dot( P1, P1 );
910 sc = dot( P2, P2 );
911 amp_PDF = KPiSFormfactor( sa, sb, sc, rRes );
912 return amp_PDF;
913}
914
915EvtComplex EvtD0ToKpipi0pi0::propogator( double mass, double width, double sx ) const {
916 EvtComplex ci( 0, 1 );
917 EvtComplex prop = 1.0 / ( mass * mass - sx - ci * mass * width );
918 return prop;
919}
920double EvtD0ToKpipi0pi0::wid( double mass, double sa, double sb, double sc, double r,
921 int l ) const {
922 double widm( 0. ), q( 0. ), q0( 0. );
923 double sa0 = mass * mass;
924 double m = sqrt( sa );
925 q = Qabcs( sa, sb, sc );
926 q0 = Qabcs( sa0, sb, sc );
927 double z, z0;
928 z = q * r * r;
929 z0 = q0 * r * r;
930 double t = q / q0;
931 double F( 0. );
932 if ( l == 0 ) F = 1;
933 if ( l == 1 ) F = sqrt( ( 1 + z0 ) / ( 1 + z ) );
934 if ( l == 2 ) F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) );
935 widm = pow( t, l + 0.5 ) * mass / m * F * F;
936 return widm;
937}
938double EvtD0ToKpipi0pi0::h( double m, double q ) const {
939 double h( 0. );
940 h = 2 / pi * q / m * log( ( m + 2 * q ) / ( 2 * mpi ) );
941 return h;
942}
943double EvtD0ToKpipi0pi0::dh( double mass, double q0 ) const {
944 double dh = h( mass, q0 ) * ( 1.0 / ( 8 * q0 * q0 ) - 1.0 / ( 2 * mass * mass ) ) +
945 1.0 / ( 2 * pi * mass * mass );
946 return dh;
947}
948double EvtD0ToKpipi0pi0::f( double mass, double sx, double q0, double q ) const {
949 double m = sqrt( sx );
950 double f = mass * mass / ( pow( q0, 3 ) ) *
951 ( q * q * ( h( m, q ) - h( mass, q0 ) ) +
952 ( mass * mass - sx ) * q0 * q0 * dh( mass, q0 ) );
953 return f;
954}
955double EvtD0ToKpipi0pi0::d( double mass, double q0 ) const {
956 double d = 3.0 / pi * mpi * mpi / ( q0 * q0 ) * log( ( mass + 2 * q0 ) / ( 2 * mpi ) ) +
957 mass / ( 2 * pi * q0 ) - ( mpi * mpi * mass ) / ( pi * pow( q0, 3 ) );
958 return d;
959}
960EvtComplex EvtD0ToKpipi0pi0::propagatorRBW( double mass, double width, double sa, double sb,
961 double sc, double r, int l ) const {
962 EvtComplex ci( 0, 1 );
963 EvtComplex prop =
964 1.0 / ( mass * mass - sa - ci * mass * width * wid( mass, sa, sb, sc, r, l ) );
965 return prop;
966}
967EvtComplex EvtD0ToKpipi0pi0::propagatorGS( double mass, double width, double sa, double sb,
968 double sc, double r, int l ) const {
969 EvtComplex ci( 0, 1 );
970 double q = Qabcs( sa, sb, sc );
971 double sa0 = mass * mass;
972 double q0 = Qabcs( sa0, sb, sc );
973 q = sqrt( q );
974 q0 = sqrt( q0 );
975 EvtComplex prop = ( 1 + d( mass, q0 ) * width / mass ) /
976 ( mass * mass - sa + width * f( mass, sa, q0, q ) -
977 ci * mass * width * wid( mass, sa, sb, sc, r, l ) );
978 return prop;
979}
980double EvtD0ToKpipi0pi0::Flatte_rhoab( double sa, double sb, double sc ) const {
981 double q = Qabcs( sa, sb, sc );
982 double rho = sqrt( q / sa );
983 return rho;
984}
985EvtComplex EvtD0ToKpipi0pi0::propagatorFlatte( double mass, double width, double sx,
986 double* sb, double* sc ) const {
987 EvtComplex ci( 0, 1 );
988 double rho1 = Flatte_rhoab( sx, sb[0], sc[0] );
989 double rho2 = Flatte_rhoab( sx, sb[1], sc[1] );
990 EvtComplex prop = 1.0 / ( mass * mass - sx - ci * ( g1 * g1 * rho1 + g2 * g2 * rho2 ) );
991 return prop;
992}
993double EvtD0ToKpipi0pi0::dot( double* a1, double* a2 ) const {
994 double dot = 0;
995 for ( int i = 0; i != 4; i++ ) { dot += a1[i] * a2[i] * G[i][i]; }
996 return dot;
997}
998double EvtD0ToKpipi0pi0::Qabcs( double sa, double sb, double sc ) const {
999 double Qabcs = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
1000 if ( Qabcs < 0 ) Qabcs = 1e-16;
1001 return Qabcs;
1002}
1003double EvtD0ToKpipi0pi0::barrier( double l, double sa, double sb, double sc, double r ) const {
1004 double q = Qabcs( sa, sb, sc );
1005 q = sqrt( q );
1006 double z = q * r;
1007 z = z * z;
1008 double F = 1;
1009 if ( l > 2 ) F = 0;
1010 if ( l == 0 ) F = 1;
1011 if ( l == 1 ) { F = sqrt( ( 2 * z ) / ( 1 + z ) ); }
1012 if ( l == 2 ) { F = sqrt( ( 13 * z * z ) / ( 9 + 3 * z + z * z ) ); }
1013 return F;
1014}
1015void EvtD0ToKpipi0pi0::calt1( double daug1[], double daug2[], double t1[] ) const {
1016 double p, pq;
1017 double pa[4], qa[4];
1018 for ( int i = 0; i != 4; i++ )
1019 {
1020 pa[i] = daug1[i] + daug2[i];
1021 qa[i] = daug1[i] - daug2[i];
1022 }
1023 p = dot( pa, pa );
1024 pq = dot( pa, qa );
1025 for ( int i = 0; i != 4; i++ ) { t1[i] = qa[i] - pq / p * pa[i]; }
1026}
1027void EvtD0ToKpipi0pi0::calt2( double daug1[], double daug2[], double t2[][4] ) const {
1028 double p, r;
1029 double pa[4], t1[4];
1030 calt1( daug1, daug2, t1 );
1031 r = dot( t1, t1 );
1032 for ( int i = 0; i != 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
1033 p = dot( pa, pa );
1034 for ( int i = 0; i != 4; i++ )
1035 {
1036 for ( int j = 0; j != 4; j++ )
1037 { t2[i][j] = t1[i] * t1[j] - 1.0 / 3 * r * ( G[i][j] - pa[i] * pa[j] / p ); }
1038 }
1039}
mass
double mass
Definition CosmicGenerator.cxx:128
pi2
character *LEPTONflag integer iresonances real pi2
Definition Ekhara/src/ekhara/common.h:3
mass_Pion
const double mass_Pion
Definition EvtD0ToKSpipipi0.cc:1561
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
EvtD0ToKpipi0pi0::clone
EvtDecayBase * clone()
Definition EvtD0ToKpipi0pi0.cc:35
EvtD0ToKpipi0pi0::decay
void decay(EvtParticle *p)
Definition EvtD0ToKpipi0pi0.cc:318
EvtD0ToKpipi0pi0::getName
void getName(std::string &name)
Definition EvtD0ToKpipi0pi0.cc:33
EvtD0ToKpipi0pi0::~EvtD0ToKpipi0pi0
virtual ~EvtD0ToKpipi0pi0()
Definition EvtD0ToKpipi0pi0.cc:31
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