BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDToKppipipi0.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: EvtDToKppipipi0.cc
11// the necessary file: EvtDToKppipipi0.hh
12//
13// Description: D+ -> K+ pi+ pi- pi0, see DocDB-1221
14//
15// Modification history:
16//
17// Liaoyuan Dong Dec. 15, 2022 Module created
18//------------------------------------------------------------------------
19#include "EvtDToKppipipi0.hh"
20#include "../EvtGenBase/EvtComplex.hh"
21#include "../EvtGenBase/EvtDecayTable.hh"
22#include "../EvtGenBase/EvtGenKine.hh"
23#include "../EvtGenBase/EvtPDL.hh"
24#include "../EvtGenBase/EvtParticle.hh"
25#include "../EvtGenBase/EvtPatches.hh"
26#include "../EvtGenBase/EvtReport.hh"
27#include <cmath>
28#include <iostream>
29#include <stdlib.h>
30using namespace std;
31
32EvtDToKppipipi0::~EvtDToKppipipi0() {}
33
34void EvtDToKppipipi0::getName( std::string& model_name ) { model_name = "DToKppipipi0"; }
35
36EvtDecayBase* EvtDToKppipipi0::clone() { return new EvtDToKppipipi0; }
37
38void EvtDToKppipipi0::init() {
39 checkNArg( 0 );
40 checkNDaug( 4 );
41 checkSpinParent( EvtSpinType::SCALAR );
42 // cout << "Initializing EvtDToKppipipi0" << endl;
43
44 mK1270 = 1.272;
45 mK1400 = 1.403;
46 GK1270 = 0.087;
47 GK1400 = 0.174;
48 mKstr0 = 0.89555;
49 mrhop = 0.77511;
50 GKstr0 = 0.0473;
51 Grhop = 0.1491;
52 mKstrp = 0.89167;
53 mrho0 = 0.77526;
54 GKstrp = 0.0514;
55 Grho0 = 0.1474;
56 mOmega = 0.78266;
57 GOmega = 0.00868;
58 msigma = 0.526;
59 Gsigma = 0.535;
60 mK1650 = 1.65;
61 GK1650 = 0.158;
62 ma1 = 1.230;
63 Ga1 = 0.420;
64
65 rho[7] = 1.0;
66 phi[7] = 0.0; // K1400+
67 rho[0] = 0.3578145188;
68 phi[0] = -1.8624743029; // 1S
69 rho[1] = 0.3468367648;
70 phi[1] = 0.1168108412; // 1P
71 rho[2] = 0.4257117223;
72 phi[2] = -0.7433416190; // 2S
73 rho[3] = 0.3042749891;
74 phi[3] = 4.3671838579; // 2P
75 rho[4] = 1.3376376890;
76 phi[4] = -3.0855437414; // 3S
77 rho[5] = 1.3376376890;
78 phi[5] = -3.0855437414; // 3S
79 rho[6] = 0.5098861534;
80 phi[6] = -1.7222823562; // 5S
81 rho[8] = 0.6422897385;
82 phi[8] = 2.5529952524; // Komega
83
84 modetype[0] = 1;
85 modetype[1] = 1;
86 modetype[2] = 2;
87 modetype[3] = 2;
88 modetype[4] = 3;
89 modetype[5] = 4;
90 modetype[6] = 5;
91 modetype[7] = 9;
92 modetype[8] = 11;
93
94 // for (int i=0; i<9; i++) {
95 // cout << i << " rho,phi = " << rho[i] << ", "<< phi[i] << endl;
96 // }
97
98 mD = 1.86966;
99 rD = 5;
100 math_pi = 3.1415926;
101
102 pi = 3.1415926;
103
104 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
105 int EE[4][4][4][4] = {
106 { { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
107 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 1 }, { 0, 0, -1, 0 } },
108 { { 0, 0, 0, 0 }, { 0, 0, 0, -1 }, { 0, 0, 0, 0 }, { 0, 1, 0, 0 } },
109 { { 0, 0, 0, 0 }, { 0, 0, 1, 0 }, { 0, -1, 0, 0 }, { 0, 0, 0, 0 } } },
110 { { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, -1 }, { 0, 0, 1, 0 } },
111 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
112 { { 0, 0, 0, 1 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { -1, 0, 0, 0 } },
113 { { 0, 0, -1, 0 }, { 0, 0, 0, 0 }, { 1, 0, 0, 0 }, { 0, 0, 0, 0 } } },
114 { { { 0, 0, 0, 0 }, { 0, 0, 0, 1 }, { 0, 0, 0, 0 }, { 0, -1, 0, 0 } },
115 { { 0, 0, 0, -1 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 1, 0, 0, 0 } },
116 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
117 { { 0, 1, 0, 0 }, { -1, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } } },
118 { { { 0, 0, 0, 0 }, { 0, 0, -1, 0 }, { 0, 1, 0, 0 }, { 0, 0, 0, 0 } },
119 { { 0, 0, 1, 0 }, { 0, 0, 0, 0 }, { -1, 0, 0, 0 }, { 0, 0, 0, 0 } },
120 { { 0, -1, 0, 0 }, { 1, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } },
121 { { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 }, { 0, 0, 0, 0 } } } };
122 for ( int i = 0; i < 4; i++ )
123 {
124 for ( int j = 0; j < 4; j++ )
125 {
126 G[i][j] = GG[i][j];
127 for ( int k = 0; k < 4; k++ )
128 {
129 for ( int l = 0; l < 4; l++ ) { E[i][j][k][l] = EE[i][j][k][l]; }
130 }
131 }
132 }
133}
134
135void EvtDToKppipipi0::initProbMax() {
136 setProbMax( 100000 ); // setProbMax(40000.0);
137}
138
139void EvtDToKppipipi0::decay( EvtParticle* p ) {
140 //-----------for max value------------------
141 /* double maxprob = 0.0;
142 for(int ir=0;ir<=60000000;ir++){
143 p->initializePhaseSpace(getNDaug(),getDaugs());
144 EvtVector4R kp = p->getDaug(0)->getP4();
145 EvtVector4R pip = p->getDaug(1)->getP4();
146 EvtVector4R pim = p->getDaug(2)->getP4();
147 EvtVector4R pi0 = p->getDaug(3)->getP4();
148 double value;
149 double Kp[4],Pip[4],Pim[4],Pi0[4];
150 Kp[0] = kp.get(0); Pip[0] = pip.get(0); Pim[0] = pim.get(0); Pi0[0]
151 = pi0.get(0); Kp[1] = kp.get(1); Pip[1] = pip.get(1); Pim[1] = pim.get(1); Pi0[1] =
152 pi0.get(1); Kp[2] = kp.get(2); Pip[2] = pip.get(2); Pim[2] = pim.get(2); Pi0[2] =
153 pi0.get(2); Kp[3] = kp.get(3); Pip[3] = pip.get(3); Pim[3] = pim.get(3); Pi0[3] =
154 pi0.get(3);
155 //double Prob = calPDF(Kp, Pip, Pim, Pi0);
156 calPDF(Kp, Pip, Pim, Pi0, value);
157 if(value>maxprob) {
158 maxprob=value;
159 std::cout << "Max PDF = " << ir << " prob= " << value << std::endl;
160 }
161 }
162 std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
163 return;*/
164 //-----------------------------------------------
165 p->initializePhaseSpace( getNDaug(), getDaugs() );
166 EvtVector4R kp = p->getDaug( 0 )->getP4();
167 EvtVector4R pip = p->getDaug( 1 )->getP4();
168 EvtVector4R pim = p->getDaug( 2 )->getP4();
169 EvtVector4R pi0 = p->getDaug( 3 )->getP4();
170
171 double Kp[4], Pip[4], Pim[4], Pi0[4];
172 Kp[0] = kp.get( 0 );
173 Pip[0] = pip.get( 0 );
174 Pim[0] = pim.get( 0 );
175 Pi0[0] = pi0.get( 0 );
176 Kp[1] = kp.get( 1 );
177 Pip[1] = pip.get( 1 );
178 Pim[1] = pim.get( 1 );
179 Pi0[1] = pi0.get( 1 );
180 Kp[2] = kp.get( 2 );
181 Pip[2] = pip.get( 2 );
182 Pim[2] = pim.get( 2 );
183 Pi0[2] = pi0.get( 2 );
184 Kp[3] = kp.get( 3 );
185 Pip[3] = pip.get( 3 );
186 Pim[3] = pim.get( 3 );
187 Pi0[3] = pi0.get( 3 );
188
189 double value;
190 calPDF( Kp, Pip, Pim, Pi0, value );
191 setProb( value );
192 return;
193}
194
195void EvtDToKppipipi0::calPDF( double Kp[], double Pip[], double Pim[], double Pi0[],
196 double& Result ) {
197 double cof[2], amp_tmp[2], amp_PDF[2], PDF[2];
198 double flag[3], mass_R[2], width_R[2];
199 double pro[2];
200 double sa[3], sb[3], sc[3], B[3];
201 double t1D[4], t1V[4], t1A[4];
202 double pS[4], pV[4], pD[4], pA[4];
203 double pro1[2], pro2[2], proKPi_S[2];
204
205 double rD2 = 25.0;
206 double rRes2 = 9.0;
207 double mass1[9] = { mKstr0, mKstr0, mKstrp, mKstrp, ma1, ma1, mK1270, mK1400, mOmega };
208 double mass2[9] = { mrhop, mrhop, mrho0, mrho0, mrhop, mrhop, mrhop, mKstr0, mrho0 };
209 double width1[9] = { GKstr0, GKstr0, GKstrp, GKstrp, Ga1, Ga1, GK1270, GK1400, GOmega };
210 double width2[9] = { Grhop, Grhop, Grho0, Grho0, Grhop, Grhop, Grhop, GKstr0, Grho0 };
211 double g0[9] = { 1, 1, 1, 1, 1, 1, 1, 1, 1 };
212 double g1[9] = { 1, 1, 1, 1, 1, 1, 1, 1, 1 };
213 double g2[9] = { 0, 1, 0, 1, 0, 0, 0, 0, 1 };
214 double temp_PDF = 0;
215 PDF[0] = 0;
216 PDF[1] = 0;
217
218 for ( int i = 0; i < 9; i++ )
219 {
220 flag[0] = g0[i];
221 flag[1] = g1[i];
222 flag[2] = g2[i];
223 mass_R[0] = mass1[i];
224 mass_R[1] = mass2[i];
225 width_R[0] = width1[i];
226 width_R[1] = width2[i];
227
228 amp_tmp[0] = 0;
229 amp_tmp[1] = 0;
230
231 cof[0] = rho[i] * cos( phi[i] );
232 cof[1] = rho[i] * sin( phi[i] );
233
234 if ( modetype[i] == 1 )
235 {
236 double t1V1[4], t1V2[4], t2D[4][4];
237 temp_PDF = 0;
238 double pV1[4], pV2[4];
239 for ( int ii = 0; ii != 4; ii++ )
240 {
241 pV1[ii] = Kp[ii] + Pim[ii];
242 pV2[ii] = Pip[ii] + Pi0[ii];
243 pD[ii] = pV1[ii] + pV2[ii];
244 }
245 sa[0] = SCADot( pV1, pV1 );
246 sb[0] = SCADot( Kp, Kp );
247 sc[0] = SCADot( Pim, Pim );
248 sa[1] = SCADot( pV2, pV2 );
249 sb[1] = SCADot( Pip, Pip );
250 sc[1] = SCADot( Pi0, Pi0 );
251 sa[2] = SCADot( pD, pD );
252 sb[2] = sa[0];
253 sc[2] = sa[1];
254 if ( flag[0] == 1 )
255 {
256 propagatorRBW( mass_R[0] * mass_R[0], mass_R[0], width_R[0], sa[0], sb[0], sc[0],
257 rRes2, 1, pro1 );
258 }
259 else if ( flag[0] == 0 )
260 {
261 pro1[0] = 1;
262 pro1[1] = 0;
263 }
264 if ( flag[1] == 1 )
265 {
266 propagatorGS( mass_R[1] * mass_R[1], mass_R[1], width_R[1], sa[1], sb[1], sc[1], rRes2,
267 pro2 );
268 }
269 else if ( flag[1] == 0 )
270 {
271 pro2[0] = 1;
272 pro2[1] = 0;
273 }
274 B[0] = Barrier( mass_R[0] * mass_R[0], 1, sa[0], sb[0], sc[0], rRes2 );
275 B[1] = Barrier( mass_R[1] * mass_R[1], 1, sa[1], sb[1], sc[1], rRes2 );
276 calt1( Kp, Pim, t1V1 );
277 calt1( Pip, Pi0, t1V2 );
278 if ( flag[2] == 0 )
279 {
280 for ( int ii = 0; ii != 4; ii++ ) { temp_PDF += ( G[ii][ii] ) * t1V1[ii] * t1V2[ii]; }
281 B[2] = 1;
282 }
283 if ( flag[2] == 1 )
284 {
285 calt1( pV1, pV2, t1D );
286 for ( int ii = 0; ii != 4; ii++ )
287 {
288 for ( int j = 0; j != 4; j++ )
289 {
290 for ( int k = 0; k != 4; k++ )
291 {
292 for ( int l = 0; l != 4; l++ )
293 {
294 temp_PDF += E[ii][j][k][l] * pD[ii] * t1D[j] * t1V1[k] * t1V2[l] *
295 ( G[ii][ii] ) * ( G[j][j] ) * ( G[l][l] ) * ( G[k][k] );
296 }
297 }
298 }
299 }
300 B[2] = Barrier( mD * mD, 1, sa[2], sb[2], sc[2], rD2 );
301 }
302 if ( flag[2] == 2 )
303 {
304 calt2( pV1, pV2, t2D );
305 for ( int ii = 0; ii != 4; ii++ )
306 {
307 for ( int j = 0; j != 4; j++ )
308 { temp_PDF += t2D[ii][j] * t1V1[ii] * t1V2[j] * ( G[ii][ii] ) * ( G[j][j] ); }
309 }
310 B[2] = Barrier( mD * mD, 2, sa[2], sb[2], sc[2], rD2 );
311 }
312 Com_Multi( pro1, pro2, pro );
313 amp_tmp[0] = temp_PDF * B[0] * B[1] * B[2] * pro[0];
314 amp_tmp[1] = temp_PDF * B[0] * B[1] * B[2] * pro[1];
315 }
316 else if ( modetype[i] == 2 )
317 {
318 double t1V1[4], t1V2[4], t2D[4][4];
319 temp_PDF = 0;
320 double pV1[4], pV2[4];
321 for ( int ii = 0; ii != 4; ii++ )
322 {
323 pV1[ii] = Kp[ii] + Pi0[ii];
324 pV2[ii] = Pip[ii] + Pim[ii];
325 pD[ii] = pV1[ii] + pV2[ii];
326 }
327 sa[0] = SCADot( pV1, pV1 );
328 sb[0] = SCADot( Kp, Kp );
329 sc[0] = SCADot( Pi0, Pi0 );
330 sa[1] = SCADot( pV2, pV2 );
331 sb[1] = SCADot( Pip, Pip );
332 sc[1] = SCADot( Pim, Pim );
333 sa[2] = SCADot( pD, pD );
334 sb[2] = sa[0];
335 sc[2] = sa[1];
336 if ( flag[0] == 1 )
337 {
338 propagatorRBW( mass_R[0] * mass_R[0], mass_R[0], width_R[0], sa[0], sb[0], sc[0],
339 rRes2, 1, pro1 );
340 }
341 else if ( flag[0] == 0 )
342 {
343 pro1[0] = 1;
344 pro1[1] = 0;
345 }
346 if ( flag[1] == 1 )
347 {
348 propagatorGS( mass_R[1] * mass_R[1], mass_R[1], width_R[1], sa[1], sb[1], sc[1], rRes2,
349 pro2 );
350 }
351 else if ( flag[1] == 0 )
352 {
353 pro2[0] = 1;
354 pro2[1] = 0;
355 }
356 B[0] = Barrier( mass_R[0] * mass_R[0], 1, sa[0], sb[0], sc[0], rRes2 );
357 B[1] = Barrier( mass_R[1] * mass_R[1], 1, sa[1], sb[1], sc[1], rRes2 );
358 calt1( Kp, Pi0, t1V1 );
359 calt1( Pip, Pim, t1V2 );
360 if ( flag[2] == 0 )
361 {
362 for ( int ii = 0; ii != 4; ii++ ) { temp_PDF += ( G[ii][ii] ) * t1V1[ii] * t1V2[ii]; }
363 B[2] = 1;
364 }
365 if ( flag[2] == 1 )
366 {
367 calt1( pV1, pV2, t1D );
368 for ( int ii = 0; ii != 4; ii++ )
369 {
370 for ( int j = 0; j != 4; j++ )
371 {
372 for ( int k = 0; k != 4; k++ )
373 {
374 for ( int l = 0; l != 4; l++ )
375 {
376 temp_PDF += E[ii][j][k][l] * pD[ii] * t1D[j] * t1V1[k] * t1V2[l] *
377 ( G[ii][ii] ) * ( G[j][j] ) * ( G[l][l] ) * ( G[k][k] );
378 }
379 }
380 }
381 }
382 B[2] = Barrier( mD * mD, 1, sa[2], sb[2], sc[2], rD2 );
383 }
384 if ( flag[2] == 2 )
385 {
386 calt2( pV1, pV2, t2D );
387 for ( int ii = 0; ii != 4; ii++ )
388 {
389 for ( int j = 0; j != 4; j++ )
390 { temp_PDF += t2D[ii][j] * t1V1[ii] * t1V2[j] * ( G[ii][ii] ) * ( G[j][j] ); }
391 }
392 B[2] = Barrier( mD * mD, 2, sa[2], sb[2], sc[2], rD2 );
393 }
394 Com_Multi( pro1, pro2, pro );
395 amp_tmp[0] = temp_PDF * B[0] * B[1] * B[2] * pro[0];
396 amp_tmp[1] = temp_PDF * B[0] * B[1] * B[2] * pro[1];
397 }
398 else if ( modetype[i] == 3 )
399 {
400 temp_PDF = 0;
401 double t2A[4][4];
402 for ( int ii = 0; ii != 4; ii++ )
403 {
404 pV[ii] = Pip[ii] + Pi0[ii];
405 pA[ii] = pV[ii] + Pim[ii];
406 pD[ii] = pA[ii] + Kp[ii];
407 }
408 sa[0] = SCADot( pV, pV );
409 sb[0] = SCADot( Pip, Pip );
410 sc[0] = SCADot( Pi0, Pi0 );
411 sa[1] = SCADot( pA, pA );
412 sb[1] = sa[0];
413 sc[1] = SCADot( Pim, Pim );
414 sa[2] = SCADot( pD, pD );
415 sb[2] = sa[1];
416 sc[2] = SCADot( Kp, Kp );
417
418 if ( flag[0] == 1 )
419 {
420 propagatorGS( mass_R[1] * mass_R[1], mass_R[1], width_R[1], sa[0], sb[0], sc[0], rRes2,
421 pro1 );
422 }
423 else if ( flag[0] == 0 )
424 {
425 pro1[0] = 1;
426 pro1[1] = 0;
427 }
428
429 if ( flag[1] == 1 )
430 {
431 propagatorRBW_a1( mass_R[0] * mass_R[0], mass_R[0], width_R[0], sa[1], sb[1], sc[1],
432 rRes2, flag[2], pro2 );
433 }
434 else if ( flag[1] == 0 )
435 {
436 pro2[0] = 1;
437 pro2[1] = 0;
438 }
439
440 B[1] = Barrier( mass_R[1] * mass_R[1], 1, sa[0], sb[0], sc[0], rRes2 );
441 B[2] = Barrier( mD * mD, 1, sa[2], sb[2], sc[2], rD2 );
442 calt1( Pip, Pi0, t1V );
443 calt1( pA, Kp, t1D );
444 if ( flag[2] == 0 )
445 {
446 for ( int ii = 0; ii != 4; ii++ )
447 {
448 for ( int j = 0; j != 4; j++ )
449 {
450 temp_PDF +=
451 t1D[ii] * ( G[ii][j] - pA[ii] * pA[j] / sa[1] ) * t1V[j] * G[ii][ii] * G[j][j];
452 }
453 }
454 B[0] = 1;
455 }
456 else if ( flag[2] == 2 )
457 {
458 calt2( pV, Pim, t2A );
459 for ( int ii = 0; ii != 4; ii++ )
460 {
461 for ( int j = 0; j != 4; j++ )
462 { temp_PDF += t1D[ii] * t2A[ii][j] * t1V[j] * G[ii][ii] * G[j][j]; }
463 }
464 B[0] = Barrier( mass_R[0] * mass_R[0], 2, sa[1], sb[1], sc[1], rRes2 );
465 }
466 Com_Multi( pro1, pro2, pro );
467 amp_tmp[0] = temp_PDF * B[0] * B[1] * B[2] * pro[0];
468 amp_tmp[1] = temp_PDF * B[0] * B[1] * B[2] * pro[1];
469 }
470 else if ( modetype[i] == 4 )
471 {
472 temp_PDF = 0;
473 double t2A[4][4];
474 for ( int ii = 0; ii != 4; ii++ )
475 {
476 pV[ii] = Pim[ii] + Pi0[ii];
477 pA[ii] = pV[ii] + Pip[ii];
478 pD[ii] = pA[ii] + Kp[ii];
479 }
480 sa[0] = SCADot( pV, pV );
481 sb[0] = SCADot( Pim, Pim );
482 sc[0] = SCADot( Pi0, Pi0 );
483 sa[1] = SCADot( pA, pA );
484 sb[1] = sa[0];
485 sc[1] = SCADot( Pip, Pip );
486 sa[2] = SCADot( pD, pD );
487 sb[2] = sa[1];
488 sc[2] = SCADot( Kp, Kp );
489
490 if ( flag[0] == 1 )
491 {
492 propagatorGS( mass_R[1] * mass_R[1], mass_R[1], width_R[1], sa[0], sb[0], sc[0], rRes2,
493 pro1 );
494 }
495 else if ( flag[0] == 0 )
496 {
497 pro1[0] = 1;
498 pro1[1] = 0;
499 }
500
501 if ( flag[1] == 1 )
502 {
503 propagatorRBW_a1( mass_R[0] * mass_R[0], mass_R[0], width_R[0], sa[1], sb[1], sc[1],
504 rRes2, flag[2], pro2 );
505 }
506 else if ( flag[1] == 0 )
507 {
508 pro2[0] = 1;
509 pro2[1] = 0;
510 }
511
512 B[1] = Barrier( mass_R[1] * mass_R[1], 1, sa[0], sb[0], sc[0], rRes2 );
513 B[2] = Barrier( mD * mD, 1, sa[2], sb[2], sc[2], rD2 );
514 calt1( Pim, Pi0, t1V );
515 calt1( pA, Kp, t1D );
516 if ( flag[2] == 0 )
517 {
518 for ( int ii = 0; ii != 4; ii++ )
519 {
520 for ( int j = 0; j != 4; j++ )
521 {
522 temp_PDF +=
523 t1D[ii] * ( G[ii][j] - pA[ii] * pA[j] / sa[1] ) * t1V[j] * G[ii][ii] * G[j][j];
524 }
525 }
526 B[0] = 1;
527 }
528 else if ( flag[2] == 2 )
529 {
530 calt2( pV, Pip, t2A );
531 for ( int ii = 0; ii != 4; ii++ )
532 {
533 for ( int j = 0; j != 4; j++ )
534 { temp_PDF += t1D[ii] * t2A[ii][j] * t1V[j] * G[ii][ii] * G[j][j]; }
535 }
536 B[0] = Barrier( mass_R[0] * mass_R[0], 2, sa[1], sb[1], sc[1], rRes2 );
537 }
538 Com_Multi( pro1, pro2, pro );
539 amp_tmp[0] = temp_PDF * B[0] * B[1] * B[2] * pro[0];
540 amp_tmp[1] = temp_PDF * B[0] * B[1] * B[2] * pro[1];
541 }
542 else if ( modetype[i] == 5 )
543 {
544 temp_PDF = 0;
545 double t2A[4][4];
546 for ( int ii = 0; ii != 4; ii++ )
547 {
548 pV[ii] = Pim[ii] + Pi0[ii];
549 pA[ii] = pV[ii] + Kp[ii];
550 pD[ii] = pA[ii] + Pip[ii];
551 }
552 sa[0] = SCADot( pV, pV );
553 sb[0] = SCADot( Pim, Pim );
554 sc[0] = SCADot( Pi0, Pi0 );
555 sa[1] = SCADot( pA, pA );
556 sb[1] = sa[0];
557 sc[1] = SCADot( Kp, Kp );
558 sa[2] = SCADot( pD, pD );
559 sb[2] = sa[1];
560 sc[2] = SCADot( Pip, Pip );
561 if ( flag[0] == 1 )
562 {
563 propagatorGS( mass_R[1] * mass_R[1], mass_R[1], width_R[1], sa[0], sb[0], sc[0], rRes2,
564 pro1 );
565 }
566 else if ( flag[0] == 0 )
567 {
568 pro1[0] = 1;
569 pro1[1] = 0;
570 }
571
572 if ( flag[1] == 1 )
573 {
574 propagatorRBW( mass_R[0] * mass_R[0], mass_R[0], width_R[0], sa[1], sb[1], sc[1],
575 rRes2, flag[2], pro2 );
576 }
577 else if ( flag[1] == 0 )
578 {
579 pro2[0] = 1;
580 pro2[1] = 0;
581 }
582 B[1] = Barrier( mass_R[1] * mass_R[1], 1, sa[0], sb[0], sc[0], rRes2 ); // rho
583 B[2] = Barrier( mD * mD, 1, sa[2], sb[2], sc[2], rD2 );
584 calt1( Pim, Pi0, t1V );
585 calt1( pA, Pip, t1D );
586 if ( flag[2] == 0 )
587 {
588 for ( int ii = 0; ii != 4; ii++ )
589 {
590 for ( int j = 0; j != 4; j++ )
591 {
592 temp_PDF +=
593 t1D[ii] * ( G[ii][j] - pA[ii] * pA[j] / sa[1] ) * t1V[j] * G[ii][ii] * G[j][j];
594 }
595 }
596 B[0] = 1;
597 }
598 else if ( flag[2] == 2 )
599 {
600 calt2( pV, Kp, t2A );
601 for ( int ii = 0; ii != 4; ii++ )
602 {
603 for ( int j = 0; j != 4; j++ )
604 { temp_PDF += t1D[ii] * t2A[ii][j] * t1V[j] * G[ii][ii] * G[j][j]; }
605 }
606 B[0] = Barrier( mass_R[0] * mass_R[0], 2, sa[1], sb[1], sc[1], rRes2 );
607 }
608 Com_Multi( pro1, pro2, pro );
609 amp_tmp[0] = temp_PDF * B[0] * B[1] * B[2] * pro[0];
610 amp_tmp[1] = temp_PDF * B[0] * B[1] * B[2] * pro[1];
611 }
612 else if ( modetype[i] == 6 )
613 {
614 temp_PDF = 0;
615 double t2A[4][4];
616 for ( int ii = 0; ii != 4; ii++ )
617 {
618 pV[ii] = Pip[ii] + Pim[ii];
619 pA[ii] = pV[ii] + Kp[ii];
620 pD[ii] = pA[ii] + Pi0[ii];
621 }
622 sa[0] = SCADot( pV, pV );
623 sb[0] = SCADot( Pip, Pip );
624 sc[0] = SCADot( Pim, Pim );
625 sa[1] = SCADot( pA, pA );
626 sb[1] = sa[0];
627 sc[1] = SCADot( Kp, Kp );
628 sa[2] = SCADot( pD, pD );
629 sb[2] = sa[1];
630 sc[2] = SCADot( Pi0, Pi0 );
631 if ( flag[0] == 1 )
632 {
633 propagatorGS( mass_R[1] * mass_R[1], mass_R[1], width_R[1], sa[0], sb[0], sc[0], rRes2,
634 pro1 );
635 }
636 else if ( flag[0] == 0 )
637 {
638 pro1[0] = 1;
639 pro1[1] = 0;
640 }
641
642 if ( flag[1] == 1 )
643 {
644 propagatorRBW( mass_R[0] * mass_R[0], mass_R[0], width_R[0], sa[1], sb[1], sc[1],
645 rRes2, 1, pro2 );
646 }
647 else if ( flag[1] == 0 )
648 {
649 pro2[0] = 1;
650 pro2[1] = 0;
651 }
652 B[1] = Barrier( mass_R[1] * mass_R[1], 1, sa[0], sb[0], sc[0], rRes2 ); // rho
653 B[2] = Barrier( mD * mD, 1, sa[2], sb[2], sc[2], rD2 );
654 calt1( Pip, Pim, t1V );
655 calt1( pA, Pi0, t1D );
656 if ( flag[2] == 0 )
657 {
658 for ( int ii = 0; ii != 4; ii++ )
659 {
660 for ( int j = 0; j != 4; j++ )
661 {
662 temp_PDF +=
663 t1D[ii] * ( G[ii][j] - pA[ii] * pA[j] / sa[1] ) * t1V[j] * G[ii][ii] * G[j][j];
664 }
665 }
666 B[0] = 1;
667 }
668 else if ( flag[2] == 2 )
669 {
670 calt2( pV, Kp, t2A );
671 for ( int ii = 0; ii != 4; ii++ )
672 {
673 for ( int j = 0; j != 4; j++ )
674 { temp_PDF += t1D[ii] * t2A[ii][j] * t1V[j] * G[ii][ii] * G[j][j]; }
675 }
676 B[0] = Barrier( mass_R[0] * mass_R[0], 2, sa[1], sb[1], sc[1], rRes2 );
677 }
678 Com_Multi( pro1, pro2, pro );
679 amp_tmp[0] = temp_PDF * B[0] * B[1] * B[2] * pro[0];
680 amp_tmp[1] = temp_PDF * B[0] * B[1] * B[2] * pro[1];
681 }
682 else if ( modetype[i] == 7 )
683 {
684 temp_PDF = 0;
685 double t2A[4][4];
686 for ( int ii = 0; ii != 4; ii++ )
687 {
688 pV[ii] = Kp[ii] + Pim[ii];
689 pA[ii] = pV[ii] + Pi0[ii];
690 pD[ii] = pA[ii] + Pip[ii];
691 }
692 sa[0] = SCADot( pV, pV );
693 sb[0] = SCADot( Kp, Kp );
694 sc[0] = SCADot( Pim, Pim );
695 sa[1] = SCADot( pA, pA );
696 sb[1] = sa[0];
697 sc[1] = SCADot( Pi0, Pi0 );
698 sa[2] = SCADot( pD, pD );
699 sb[2] = sa[1];
700 sc[2] = SCADot( Pip, Pip );
701 if ( flag[0] == 1 )
702 {
703 propagatorRBW( mass_R[1] * mass_R[1], mass_R[1], width_R[1], sa[0], sb[0], sc[0],
704 rRes2, 1, pro1 );
705 }
706 else if ( flag[0] == 0 )
707 {
708 pro1[0] = 1;
709 pro1[1] = 0;
710 }
711 if ( flag[1] == 1 )
712 {
713 propagatorRBW( mass_R[0] * mass_R[0], mass_R[0], width_R[0], sa[1], sb[1], sc[1],
714 rRes2, flag[2], pro2 );
715 }
716 else if ( flag[1] == 0 )
717 {
718 pro2[0] = 1;
719 pro2[1] = 0;
720 }
721
722 B[1] = Barrier( mass_R[1] * mass_R[1], 1, sa[0], sb[0], sc[0], rRes2 );
723 B[2] = Barrier( mD * mD, 1, sa[2], sb[2], sc[2], rD2 );
724 calt1( Kp, Pim, t1V );
725 calt1( pA, Pip, t1D );
726 if ( flag[2] == 0 )
727 {
728 for ( int ii = 0; ii != 4; ii++ )
729 {
730 for ( int j = 0; j != 4; j++ )
731 {
732 temp_PDF +=
733 t1D[ii] * ( G[ii][j] - pA[ii] * pA[j] / sa[1] ) * t1V[j] * G[ii][ii] * G[j][j];
734 }
735 }
736 B[0] = 1;
737 }
738 else if ( flag[2] == 2 )
739 {
740 calt2( pV, Pi0, t2A );
741 for ( int ii = 0; ii != 4; ii++ )
742 {
743 for ( int j = 0; j != 4; j++ )
744 { temp_PDF += t1D[ii] * t2A[ii][j] * t1V[j] * G[ii][ii] * G[j][j]; }
745 }
746 B[0] = Barrier( mass_R[0] * mass_R[0], 2, sa[1], sb[1], sc[1], rRes2 );
747 }
748 Com_Multi( pro1, pro2, pro );
749 amp_tmp[0] = temp_PDF * B[0] * B[1] * B[2] * pro[0];
750 amp_tmp[1] = temp_PDF * B[0] * B[1] * B[2] * pro[1];
751 }
752 else if ( modetype[i] == 8 )
753 {
754 temp_PDF = 0;
755 double t2A[4][4];
756 for ( int ii = 0; ii != 4; ii++ )
757 {
758 pV[ii] = Kp[ii] + Pi0[ii];
759 pA[ii] = pV[ii] + Pim[ii];
760 pD[ii] = pA[ii] + Pip[ii];
761 }
762 sa[0] = SCADot( pV, pV );
763 sb[0] = SCADot( Kp, Kp );
764 sc[0] = SCADot( Pi0, Pi0 );
765 sa[1] = SCADot( pA, pA );
766 sb[1] = sa[0];
767 sc[1] = SCADot( Pim, Pim );
768 sa[2] = SCADot( pD, pD );
769 sb[2] = sa[1];
770 sc[2] = SCADot( Pip, Pip );
771 if ( flag[0] == 1 )
772 {
773 propagatorRBW( mass_R[1] * mass_R[1], mass_R[1], width_R[1], sa[0], sb[0], sc[0],
774 rRes2, 1, pro1 );
775 }
776 else if ( flag[0] == 0 )
777 {
778 pro1[0] = 1;
779 pro1[1] = 0;
780 }
781 if ( flag[1] == 1 )
782 {
783 propagatorRBW( mass_R[0] * mass_R[0], mass_R[0], width_R[0], sa[1], sb[1], sc[1],
784 rRes2, flag[2], pro2 );
785 }
786 else if ( flag[1] == 0 )
787 {
788 pro2[0] = 1;
789 pro2[1] = 0;
790 }
791
792 B[1] = Barrier( mass_R[1] * mass_R[1], 1, sa[0], sb[0], sc[0], rRes2 );
793 B[2] = Barrier( mD * mD, 1, sa[2], sb[2], sc[2], rD2 );
794 calt1( Kp, Pi0, t1V );
795 calt1( pA, Pip, t1D );
796 if ( flag[2] == 0 )
797 {
798 for ( int ii = 0; ii != 4; ii++ )
799 {
800 for ( int j = 0; j != 4; j++ )
801 {
802 temp_PDF +=
803 t1D[ii] * ( G[ii][j] - pA[ii] * pA[j] / sa[1] ) * t1V[j] * G[ii][ii] * G[j][j];
804 }
805 }
806 B[0] = 1;
807 }
808 else if ( flag[2] == 2 )
809 {
810 calt2( pV, Pim, t2A );
811 for ( int ii = 0; ii != 4; ii++ )
812 {
813 for ( int j = 0; j != 4; j++ )
814 { temp_PDF += t1D[ii] * t2A[ii][j] * t1V[j] * G[ii][ii] * G[j][j]; }
815 }
816 B[0] = Barrier( mass_R[0] * mass_R[0], 2, sa[1], sb[1], sc[1], rRes2 );
817 }
818 Com_Multi( pro1, pro2, pro );
819 amp_tmp[0] = temp_PDF * B[0] * B[1] * B[2] * pro[0];
820 amp_tmp[1] = temp_PDF * B[0] * B[1] * B[2] * pro[1];
821 }
822 else if ( modetype[i] == 9 )
823 {
824 temp_PDF = 0;
825 double t2A[4][4];
826 for ( int ii = 0; ii != 4; ii++ )
827 {
828 pV[ii] = Kp[ii] + Pim[ii];
829 pA[ii] = pV[ii] + Pip[ii];
830 pD[ii] = pA[ii] + Pi0[ii];
831 }
832 sa[0] = SCADot( pV, pV );
833 sb[0] = SCADot( Kp, Kp );
834 sc[0] = SCADot( Pim, Pim );
835 sa[1] = SCADot( pA, pA );
836 sb[1] = sa[0];
837 sc[1] = SCADot( Pip, Pip );
838 sa[2] = SCADot( pD, pD );
839 sb[2] = sa[1];
840 sc[2] = SCADot( Pi0, Pi0 );
841 if ( flag[0] == 1 )
842 {
843 propagatorRBW( mass_R[1] * mass_R[1], mass_R[1], width_R[1], sa[0], sb[0], sc[0],
844 rRes2, 1, pro1 );
845 }
846 else if ( flag[0] == 0 )
847 {
848 pro1[0] = 1;
849 pro1[1] = 0;
850 }
851 if ( flag[1] == 1 )
852 {
853 propagatorRBW( mass_R[0] * mass_R[0], mass_R[0], width_R[0], sa[1], sb[1], sc[1],
854 rRes2, flag[2], pro2 );
855 }
856 else if ( flag[1] == 0 )
857 {
858 pro2[0] = 1;
859 pro2[1] = 0;
860 }
861
862 B[1] = Barrier( mass_R[1] * mass_R[1], 1, sa[0], sb[0], sc[0], rRes2 );
863 B[2] = Barrier( mD * mD, 1, sa[2], sb[2], sc[2], rD2 );
864 calt1( Kp, Pim, t1V );
865 calt1( pA, Pi0, t1D );
866 if ( flag[2] == 0 )
867 {
868 for ( int ii = 0; ii != 4; ii++ )
869 {
870 for ( int j = 0; j != 4; j++ )
871 {
872 temp_PDF +=
873 t1D[ii] * ( G[ii][j] - pA[ii] * pA[j] / sa[1] ) * t1V[j] * G[ii][ii] * G[j][j];
874 }
875 }
876 B[0] = 1;
877 }
878 else if ( flag[2] == 2 )
879 {
880 calt2( pV, Pip, t2A );
881 for ( int ii = 0; ii != 4; ii++ )
882 {
883 for ( int j = 0; j != 4; j++ )
884 { temp_PDF += t1D[ii] * t2A[ii][j] * t1V[j] * G[ii][ii] * G[j][j]; }
885 }
886 B[0] = Barrier( mass_R[0] * mass_R[0], 2, sa[1], sb[1], sc[1], rRes2 );
887 }
888 Com_Multi( pro1, pro2, pro );
889 amp_tmp[0] = temp_PDF * B[0] * B[1] * B[2] * pro[0];
890 amp_tmp[1] = temp_PDF * B[0] * B[1] * B[2] * pro[1];
891 }
892 else if ( modetype[i] == 11 )
893 {
894 // rho0 pi0
895 temp_PDF = 0;
896 double pV1[4], pV2[4], qV1[4], qV2[4];
897 for ( int ii = 0; ii != 4; ii++ )
898 {
899 pV2[ii] = Pip[ii] + Pim[ii];
900 qV2[ii] = Pip[ii] - Pim[ii];
901 pV1[ii] = pV2[ii] + Pi0[ii];
902 qV1[ii] = pV2[ii] - Pi0[ii];
903 pD[ii] = pV1[ii] + Kp[ii];
904 }
905 for ( int ii = 0; ii != 4; ii++ )
906 {
907 for ( int j = 0; j != 4; j++ )
908 {
909 for ( int k = 0; k != 4; k++ )
910 {
911 for ( int l = 0; l != 4; l++ )
912 {
913 temp_PDF += E[ii][j][k][l] * pV1[ii] * qV1[j] * Kp[k] * qV2[l] * ( G[ii][ii] ) *
914 ( G[j][j] ) * ( G[k][k] ) * ( G[l][l] );
915 }
916 }
917 }
918 }
919 sa[0] = SCADot( pV2, pV2 );
920 sb[0] = SCADot( Pip, Pip );
921 sc[0] = SCADot( Pim, Pim );
922 sa[1] = SCADot( pV1, pV1 );
923 sb[1] = sa[0];
924 sc[1] = SCADot( Pi0, Pi0 );
925 sa[2] = SCADot( pD, pD );
926 sb[2] = sa[1];
927 sc[2] = SCADot( Kp, Kp );
928 propagatorGS( mass_R[1] * mass_R[1], mass_R[1], width_R[1], sa[0], sb[0], sc[0], rRes2,
929 pro1 );
930 if ( flag[0] == 1 )
931 {
932 propagatorRBW( mass_R[0] * mass_R[0], mass_R[0], width_R[0], sa[1], sb[1], sc[1],
933 rRes2, 1, pro2 );
934 }
935 else if ( flag[0] == 0 )
936 {
937 pro1[0] = 1;
938 pro1[1] = 0;
939 }
940 Com_Multi( pro1, pro2, pro );
941
942 B[1] = Barrier( mass_R[1] * mass_R[1], 1, sa[0], sb[0], sc[0], rRes2 );
943 B[0] = Barrier( mass_R[0] * mass_R[0], 1, sa[1], sb[1], sc[1], rRes2 );
944 B[2] = Barrier( mD * mD, 1, sa[2], sb[2], sc[2], rD2 );
945
946 amp_tmp[0] = temp_PDF * B[0] * B[1] * B[2] * pro[0];
947 amp_tmp[1] = temp_PDF * B[0] * B[1] * B[2] * pro[1];
948
949 // rho+ pi-
950 temp_PDF = 0;
951 // double pV1[4], pV2[4], qV1[4], qV2[4];
952 for ( int ii = 0; ii != 4; ii++ )
953 {
954 pV2[ii] = Pip[ii] + Pi0[ii];
955 qV2[ii] = Pip[ii] - Pi0[ii];
956 pV1[ii] = pV2[ii] + Pim[ii];
957 qV1[ii] = pV2[ii] - Pim[ii];
958 pD[ii] = pV1[ii] + Kp[ii];
959 }
960 for ( int ii = 0; ii != 4; ii++ )
961 {
962 for ( int j = 0; j != 4; j++ )
963 {
964 for ( int k = 0; k != 4; k++ )
965 {
966 for ( int l = 0; l != 4; l++ )
967 {
968 temp_PDF += E[ii][j][k][l] * pV1[ii] * qV1[j] * Kp[k] * qV2[l] * ( G[ii][ii] ) *
969 ( G[j][j] ) * ( G[k][k] ) * ( G[l][l] );
970 }
971 }
972 }
973 }
974 sa[0] = SCADot( pV2, pV2 );
975 sb[0] = SCADot( Pip, Pip );
976 sc[0] = SCADot( Pi0, Pi0 );
977 sa[1] = SCADot( pV1, pV1 );
978 sb[1] = sa[0];
979 sc[1] = SCADot( Pim, Pim );
980 sa[2] = SCADot( pD, pD );
981 sb[2] = sa[1];
982 sc[2] = SCADot( Kp, Kp );
983 propagatorGS( mass_R[1] * mass_R[1], mass_R[1], width_R[1], sa[0], sb[0], sc[0], rRes2,
984 pro1 );
985 if ( flag[0] == 1 )
986 {
987 propagatorRBW( mass_R[0] * mass_R[0], mass_R[0], width_R[0], sa[1], sb[1], sc[1],
988 rRes2, 1, pro2 );
989 }
990 else if ( flag[0] == 0 )
991 {
992 pro1[0] = 1;
993 pro1[1] = 0;
994 }
995 Com_Multi( pro1, pro2, pro );
996
997 B[1] = Barrier( mass_R[1] * mass_R[1], 1, sa[0], sb[0], sc[0], rRes2 );
998 B[0] = Barrier( mass_R[0] * mass_R[0], 1, sa[1], sb[1], sc[1], rRes2 );
999 B[2] = Barrier( mD * mD, 1, sa[2], sb[2], sc[2], rD2 );
1000
1001 amp_tmp[0] += temp_PDF * B[0] * B[1] * B[2] * pro[0];
1002 amp_tmp[1] += temp_PDF * B[0] * B[1] * B[2] * pro[1];
1003
1004 // rho- pi+
1005 temp_PDF = 0;
1006 // double pV1[4], pV2[4], qV1[4], qV2[4];
1007 for ( int ii = 0; ii != 4; ii++ )
1008 {
1009 pV2[ii] = Pim[ii] + Pi0[ii];
1010 qV2[ii] = Pim[ii] - Pi0[ii];
1011 pV1[ii] = pV2[ii] + Pip[ii];
1012 qV1[ii] = pV2[ii] - Pip[ii];
1013 pD[ii] = pV1[ii] + Kp[ii];
1014 }
1015 for ( int ii = 0; ii != 4; ii++ )
1016 {
1017 for ( int j = 0; j != 4; j++ )
1018 {
1019 for ( int k = 0; k != 4; k++ )
1020 {
1021 for ( int l = 0; l != 4; l++ )
1022 {
1023 temp_PDF += E[ii][j][k][l] * pV1[ii] * qV1[j] * Kp[k] * qV2[l] * ( G[ii][ii] ) *
1024 ( G[j][j] ) * ( G[k][k] ) * ( G[l][l] );
1025 }
1026 }
1027 }
1028 }
1029 sa[0] = SCADot( pV2, pV2 );
1030 sb[0] = SCADot( Pim, Pim );
1031 sc[0] = SCADot( Pi0, Pi0 );
1032 sa[1] = SCADot( pV1, pV1 );
1033 sb[1] = sa[0];
1034 sc[1] = SCADot( Pip, Pip );
1035 sa[2] = SCADot( pD, pD );
1036 sb[2] = sa[1];
1037 sc[2] = SCADot( Kp, Kp );
1038 propagatorGS( mass_R[1] * mass_R[1], mass_R[1], width_R[1], sa[0], sb[0], sc[0], rRes2,
1039 pro1 );
1040 if ( flag[0] == 1 )
1041 {
1042 propagatorRBW( mass_R[0] * mass_R[0], mass_R[0], width_R[0], sa[1], sb[1], sc[1],
1043 rRes2, 1, pro2 );
1044 }
1045 else if ( flag[0] == 0 )
1046 {
1047 pro1[0] = 1;
1048 pro1[1] = 0;
1049 }
1050 Com_Multi( pro1, pro2, pro );
1051
1052 B[1] = Barrier( mass_R[1] * mass_R[1], 1, sa[0], sb[0], sc[0], rRes2 );
1053 B[0] = Barrier( mass_R[0] * mass_R[0], 1, sa[1], sb[1], sc[1], rRes2 );
1054 B[2] = Barrier( mD * mD, 1, sa[2], sb[2], sc[2], rD2 );
1055
1056 amp_tmp[0] += temp_PDF * B[0] * B[1] * B[2] * pro[0];
1057 amp_tmp[1] += temp_PDF * B[0] * B[1] * B[2] * pro[1];
1058 }
1059
1060 Com_Multi( amp_tmp, cof, amp_PDF );
1061 PDF[0] += amp_PDF[0];
1062 PDF[1] += amp_PDF[1];
1063 }
1064 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
1065 Result = value;
1066}
1067
1068void EvtDToKppipipi0::Com_Multi( double a1[2], double a2[2], double res[2] ) {
1069 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
1070 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
1071}
1072void EvtDToKppipipi0::Com_Divide( double a1[2], double a2[2], double res[2] ) {
1073 double tmp = a2[0] * a2[0] + a2[1] * a2[1];
1074 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / tmp;
1075 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / tmp;
1076}
1077double EvtDToKppipipi0::SCADot( double a1[4], double a2[4] ) {
1078 double _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
1079 return _cal;
1080}
1081double EvtDToKppipipi0::Barrier( double mass2, int l, double sa, double sb, double sc,
1082 double r2 ) {
1083 double F;
1084 double tmp = sa + sb - sc;
1085 double q = fabs( 0.25 * tmp * tmp / sa - sb );
1086 // if (q < 0) q = 1e-16;
1087 double tmp2 = mass2 + sb - sc;
1088 double q0 = fabs( 0.25 * tmp2 * tmp2 / mass2 - sb );
1089 // if (q0 < 0) q0 = 1e-16;
1090 double z = q * r2;
1091 double z0 = q0 * r2;
1092 if ( l == 1 ) { F = sqrt( ( 1.0 + z0 ) / ( 1.0 + z ) ); }
1093 else if ( l == 2 )
1094 {
1095 double z2 = z * z;
1096 double z02 = z0 * z0;
1097 F = sqrt( ( 9.0 + 3.0 * z0 + z02 ) / ( 9.0 + 3.0 * z + z2 ) );
1098 }
1099 else { F = 1.0; }
1100 return F;
1101}
1102void EvtDToKppipipi0::calt1( double daug1[4], double daug2[4], double t1[4] ) {
1103 double p, pq, tmp;
1104 double pa[4], qa[4];
1105 for ( int i = 0; i < 4; i++ )
1106 {
1107 pa[i] = daug1[i] + daug2[i];
1108 qa[i] = daug1[i] - daug2[i];
1109 }
1110 p = SCADot( pa, pa );
1111 pq = SCADot( pa, qa );
1112 tmp = pq / p;
1113 for ( int i = 0; i < 4; i++ ) { t1[i] = qa[i] - tmp * pa[i]; }
1114}
1115void EvtDToKppipipi0::calt2( double daug1[4], double daug2[4], double t2[4][4] ) {
1116 double p, r;
1117 double pa[4], t1[4];
1118 calt1( daug1, daug2, t1 );
1119 r = SCADot( t1, t1 ) / 3.0;
1120 for ( int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
1121 p = SCADot( pa, pa );
1122 for ( int i = 0; i < 4; i++ )
1123 {
1124 for ( int j = 0; j < 4; j++ )
1125 { t2[i][j] = t1[i] * t1[j] - r * ( G[i][j] - pa[i] * pa[j] / p ); }
1126 }
1127}
1128void EvtDToKppipipi0::propagator( double mass2, double mass, double width, double sx,
1129 double prop[2] ) {
1130 double a[2], b[2];
1131 a[0] = 1;
1132 a[1] = 0;
1133 b[0] = mass2 - sx;
1134 b[1] = -mass * width;
1135 Com_Divide( a, b, prop );
1136}
1137double EvtDToKppipipi0::wid( double mass2, double mass, double sa, double sb, double sc,
1138 double r2, int l ) {
1139 double widm = 0.;
1140 double m = sqrt( sa );
1141 double tmp = sb - sc;
1142 double tmp1 = sa + tmp;
1143 double q = fabs( 0.25 * tmp1 * tmp1 / sa - sb );
1144 // if(q<0) q = 1e-16;
1145 double tmp2 = mass2 + tmp;
1146 double q0 = fabs( 0.25 * tmp2 * tmp2 / mass2 - sb );
1147 // if(q0<0) q0 = 1e-16;
1148 double z = q * r2;
1149 double z0 = q0 * r2;
1150 double t = q / q0;
1151 if ( l == 0 ) { widm = sqrt( t ) * mass / m; }
1152 else if ( l == 1 ) { widm = t * sqrt( t ) * mass / m * ( 1 + z0 ) / ( 1 + z ); }
1153 else if ( l == 2 )
1154 { widm = t * t * sqrt( t ) * mass / m * ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ); }
1155 return widm;
1156}
1157double EvtDToKppipipi0::widl1( double mass2, double mass, double sa, double sb, double sc,
1158 double r2 ) {
1159 double widm = 0.;
1160 double m = sqrt( sa );
1161 double tmp = sb - sc;
1162 double tmp1 = sa + tmp;
1163 double q = fabs( 0.25 * tmp1 * tmp1 / sa - sb );
1164 // if(q<0) q = 1e-16;
1165 double tmp2 = mass2 + tmp;
1166 double q0 = fabs( 0.25 * tmp2 * tmp2 / mass2 - sb );
1167 // if(q0<0) q0 = 1e-16;
1168 double z = q * r2;
1169 double z0 = q0 * r2;
1170 double F = ( 1 + z0 ) / ( 1 + z );
1171 double t = q / q0;
1172 widm = t * sqrt( t ) * mass / m * F;
1173 return widm;
1174}
1175void EvtDToKppipipi0::propagatorRBW( double mass2, double mass, double width, double sa,
1176 double sb, double sc, double r2, int l, double prop[2] ) {
1177 double a[2], b[2];
1178 a[0] = 1;
1179 a[1] = 0;
1180 b[0] = mass2 - sa;
1181 b[1] = -mass * width * wid( mass2, mass, sa, sb, sc, r2, l );
1182 Com_Divide( a, b, prop );
1183}
1184void EvtDToKppipipi0::propagatorRBWl1( double mass2, double mass, double width, double sa,
1185 double sb, double sc, double r2, double prop[2] ) {
1186 double a[2], b[2];
1187 a[0] = 1;
1188 a[1] = 0;
1189 b[0] = mass2 - sa;
1190 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r2 );
1191 Com_Divide( a, b, prop );
1192}
1193void EvtDToKppipipi0::propagatorRBW_a1( double mass2, double mass, double width, double sa,
1194 double sb, double sc, double r2, int l,
1195 double prop[2] ) {
1196 double a[2], b[2];
1197 int iii = int( sqrt( sa ) * 1000 ) - 1;
1198 double a1width[3000] = {
1199 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1200 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1201 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1202 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1203 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1204 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1205 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1206 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1207 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1208 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1209 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1210 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1211 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1212 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1213 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1214 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1215 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1216 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1217 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1218 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1219 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1220 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1221 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1222 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1223 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1224 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1225 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1226 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1227 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1228 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1229 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1230 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1231 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1232 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1233 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1234 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1235 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1236 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1237 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1238 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1239 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1240 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1241 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1242 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1243 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1244 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1245 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
1246 0.000000, 0.000000, 0.000001, 0.000001, 0.000001, 0.000002, 0.000002, 0.000002, 0.000003,
1247 0.000004, 0.000004, 0.000005, 0.000006, 0.000007, 0.000008, 0.000009, 0.000010, 0.000011,
1248 0.000012, 0.000014, 0.000015, 0.000017, 0.000019, 0.000021, 0.000023, 0.000025, 0.000027,
1249 0.000029, 0.000032, 0.000035, 0.000038, 0.000041, 0.000044, 0.000047, 0.000050, 0.000054,
1250 0.000058, 0.000062, 0.000066, 0.000070, 0.000075, 0.000079, 0.000084, 0.000089, 0.000094,
1251 0.000100, 0.000105, 0.000111, 0.000117, 0.000124, 0.000130, 0.000137, 0.000143, 0.000151,
1252 0.000158, 0.000165, 0.000173, 0.000182, 0.000190, 0.000199, 0.000207, 0.000216, 0.000225,
1253 0.000235, 0.000245, 0.000256, 0.000266, 0.000277, 0.000288, 0.000300, 0.000311, 0.000322,
1254 0.000335, 0.000347, 0.000360, 0.000373, 0.000385, 0.000400, 0.000415, 0.000429, 0.000442,
1255 0.000457, 0.000473, 0.000488, 0.000504, 0.000520, 0.000539, 0.000555, 0.000572, 0.000590,
1256 0.000608, 0.000626, 0.000646, 0.000664, 0.000684, 0.000704, 0.000725, 0.000745, 0.000766,
1257 0.000787, 0.000809, 0.000828, 0.000854, 0.000878, 0.000901, 0.000927, 0.000952, 0.000973,
1258 0.001001, 0.001027, 0.001048, 0.001080, 0.001104, 0.001132, 0.001159, 0.001189, 0.001219,
1259 0.001245, 0.001277, 0.001308, 0.001338, 0.001370, 0.001404, 0.001433, 0.001468, 0.001498,
1260 0.001533, 0.001570, 0.001600, 0.001638, 0.001678, 0.001711, 0.001745, 0.001780, 0.001825,
1261 0.001857, 0.001898, 0.001941, 0.001972, 0.002017, 0.002065, 0.002104, 0.002146, 0.002189,
1262 0.002234, 0.002277, 0.002319, 0.002369, 0.002410, 0.002461, 0.002511, 0.002557, 0.002605,
1263 0.002661, 0.002704, 0.002762, 0.002807, 0.002855, 0.002910, 0.002965, 0.003020, 0.003074,
1264 0.003127, 0.003178, 0.003228, 0.003288, 0.003351, 0.003409, 0.003471, 0.003532, 0.003598,
1265 0.003660, 0.003720, 0.003793, 0.003854, 0.003910, 0.003972, 0.004050, 0.004108, 0.004181,
1266 0.004254, 0.004309, 0.004380, 0.004464, 0.004533, 0.004603, 0.004679, 0.004756, 0.004811,
1267 0.004898, 0.004974, 0.005048, 0.005142, 0.005215, 0.005279, 0.005363, 0.005449, 0.005533,
1268 0.005604, 0.005695, 0.005783, 0.005869, 0.005971, 0.006060, 0.006142, 0.006247, 0.006332,
1269 0.006409, 0.006502, 0.006594, 0.006713, 0.006784, 0.006889, 0.006995, 0.007079, 0.007190,
1270 0.007303, 0.007381, 0.007487, 0.007592, 0.007710, 0.007801, 0.007910, 0.008032, 0.008149,
1271 0.008247, 0.008378, 0.008462, 0.008559, 0.008706, 0.008843, 0.008943, 0.009091, 0.009207,
1272 0.009308, 0.009448, 0.009555, 0.009698, 0.009810, 0.009936, 0.010020, 0.010186, 0.010320,
1273 0.010474, 0.010611, 0.010742, 0.010840, 0.011011, 0.011167, 0.011281, 0.011395, 0.011541,
1274 0.011714, 0.011853, 0.012046, 0.012169, 0.012277, 0.012460, 0.012617, 0.012805, 0.012922,
1275 0.013072, 0.013234, 0.013389, 0.013561, 0.013704, 0.013917, 0.014025, 0.014239, 0.014425,
1276 0.014600, 0.014716, 0.014958, 0.015114, 0.015325, 0.015488, 0.015630, 0.015797, 0.016035,
1277 0.016206, 0.016404, 0.016591, 0.016842, 0.016964, 0.017199, 0.017392, 0.017557, 0.017798,
1278 0.017987, 0.018178, 0.018337, 0.018631, 0.018829, 0.019008, 0.019221, 0.019467, 0.019698,
1279 0.019941, 0.020166, 0.020379, 0.020585, 0.020806, 0.021040, 0.021309, 0.021482, 0.021764,
1280 0.022046, 0.022306, 0.022478, 0.022736, 0.023049, 0.023222, 0.023536, 0.023756, 0.024090,
1281 0.024337, 0.024557, 0.024873, 0.025099, 0.025392, 0.025682, 0.025995, 0.026291, 0.026498,
1282 0.026927, 0.027119, 0.027377, 0.027804, 0.028135, 0.028279, 0.028682, 0.028871, 0.029355,
1283 0.029531, 0.029956, 0.030243, 0.030592, 0.030873, 0.031246, 0.031494, 0.031771, 0.032167,
1284 0.032515, 0.032881, 0.033211, 0.033653, 0.033988, 0.034394, 0.034639, 0.035055, 0.035569,
1285 0.035879, 0.036211, 0.036611, 0.036932, 0.037489, 0.037779, 0.038284, 0.038723, 0.039001,
1286 0.039574, 0.039854, 0.040274, 0.040881, 0.041117, 0.041644, 0.042055, 0.042531, 0.043030,
1287 0.043354, 0.043832, 0.044277, 0.044956, 0.045284, 0.045828, 0.046440, 0.046800, 0.047518,
1288 0.047727, 0.048258, 0.048850, 0.049316, 0.049992, 0.050486, 0.050987, 0.051410, 0.051928,
1289 0.052613, 0.053110, 0.053824, 0.054351, 0.055078, 0.055654, 0.056030, 0.056763, 0.057245,
1290 0.057832, 0.058569, 0.059292, 0.060048, 0.060569, 0.061056, 0.061869, 0.062612, 0.063186,
1291 0.063886, 0.064655, 0.065198, 0.065815, 0.066649, 0.067577, 0.068012, 0.068967, 0.069630,
1292 0.070181, 0.070786, 0.071989, 0.072764, 0.073466, 0.074461, 0.075093, 0.075994, 0.076834,
1293 0.077455, 0.078709, 0.079581, 0.080408, 0.080884, 0.081965, 0.082882, 0.083658, 0.084824,
1294 0.085513, 0.086662, 0.087602, 0.088678, 0.089492, 0.090641, 0.091369, 0.092494, 0.093484,
1295 0.094615, 0.095385, 0.096168, 0.097668, 0.098611, 0.099630, 0.100772, 0.102020, 0.103145,
1296 0.104110, 0.105071, 0.106604, 0.107791, 0.108451, 0.109509, 0.111356, 0.112026, 0.113921,
1297 0.114507, 0.116071, 0.117027, 0.118213, 0.120164, 0.120701, 0.122121, 0.123894, 0.124937,
1298 0.126134, 0.127391, 0.128882, 0.130056, 0.131649, 0.133046, 0.134275, 0.135119, 0.137072,
1299 0.138476, 0.139612, 0.140388, 0.142734, 0.143576, 0.145445, 0.147414, 0.148856, 0.149891,
1300 0.150963, 0.152477, 0.153717, 0.155275, 0.156859, 0.158462, 0.159257, 0.161865, 0.163182,
1301 0.164465, 0.165538, 0.167003, 0.169257, 0.171211, 0.172093, 0.173261, 0.174639, 0.176510,
1302 0.177684, 0.179077, 0.181041, 0.182446, 0.184769, 0.184926, 0.186741, 0.188844, 0.190884,
1303 0.191714, 0.192254, 0.193921, 0.195917, 0.196766, 0.199052, 0.200603, 0.201808, 0.202699,
1304 0.204636, 0.205712, 0.206849, 0.208741, 0.209424, 0.211698, 0.212753, 0.215516, 0.215857,
1305 0.217790, 0.217774, 0.220454, 0.221821, 0.223466, 0.224494, 0.225632, 0.227231, 0.229456,
1306 0.229581, 0.231537, 0.232263, 0.233834, 0.234725, 0.237079, 0.238015, 0.239400, 0.240193,
1307 0.241693, 0.243787, 0.244317, 0.244971, 0.246711, 0.248615, 0.249387, 0.250905, 0.252702,
1308 0.253535, 0.254385, 0.255375, 0.256671, 0.258405, 0.259741, 0.260875, 0.262131, 0.262920,
1309 0.264860, 0.265893, 0.266016, 0.267727, 0.270039, 0.270689, 0.271047, 0.272313, 0.272474,
1310 0.274724, 0.275813, 0.275937, 0.278793, 0.278783, 0.281407, 0.281351, 0.282481, 0.284226,
1311 0.284113, 0.284999, 0.285655, 0.288361, 0.287856, 0.288893, 0.290211, 0.291708, 0.291985,
1312 0.294298, 0.294849, 0.296796, 0.296197, 0.296851, 0.298011, 0.300368, 0.299982, 0.302378,
1313 0.304363, 0.303711, 0.304729, 0.306789, 0.306378, 0.307372, 0.308720, 0.309509, 0.309712,
1314 0.310782, 0.311699, 0.312668, 0.312755, 0.313675, 0.315311, 0.316640, 0.317217, 0.317403,
1315 0.318478, 0.319916, 0.321803, 0.322678, 0.323237, 0.324343, 0.324433, 0.324493, 0.324969,
1316 0.325894, 0.328563, 0.328721, 0.328954, 0.330640, 0.328164, 0.331267, 0.331695, 0.333772,
1317 0.333619, 0.334351, 0.334605, 0.336434, 0.337510, 0.336535, 0.337362, 0.338799, 0.340732,
1318 0.339896, 0.342707, 0.343471, 0.342318, 0.342431, 0.344543, 0.345611, 0.345786, 0.346590,
1319 0.346610, 0.347761, 0.348914, 0.349558, 0.350577, 0.352128, 0.350982, 0.354134, 0.352773,
1320 0.353213, 0.352972, 0.354927, 0.355784, 0.355778, 0.355801, 0.357040, 0.358013, 0.358432,
1321 0.360045, 0.359743, 0.360238, 0.359850, 0.362184, 0.361580, 0.363430, 0.362333, 0.364397,
1322 0.364472, 0.364370, 0.365303, 0.366644, 0.367777, 0.368604, 0.367631, 0.368324, 0.369782,
1323 0.371121, 0.370653, 0.370040, 0.371649, 0.370201, 0.373362, 0.373900, 0.374159, 0.374916,
1324 0.374503, 0.376703, 0.372802, 0.376191, 0.379596, 0.377325, 0.376363, 0.379369, 0.379791,
1325 0.378703, 0.380177, 0.381762, 0.381335, 0.381374, 0.384668, 0.381763, 0.382746, 0.384723,
1326 0.385089, 0.386229, 0.386702, 0.387749, 0.384423, 0.384714, 0.384181, 0.388489, 0.388618,
1327 0.388179, 0.390092, 0.389871, 0.390496, 0.391181, 0.390679, 0.392614, 0.392269, 0.393899,
1328 0.393466, 0.391421, 0.391090, 0.395586, 0.391776, 0.396882, 0.393254, 0.394400, 0.395749,
1329 0.398063, 0.397138, 0.397585, 0.397288, 0.397847, 0.395375, 0.400170, 0.400007, 0.401191,
1330 0.398513, 0.401922, 0.400477, 0.404257, 0.403271, 0.400677, 0.403913, 0.403172, 0.404727,
1331 0.403406, 0.404404, 0.405265, 0.406389, 0.405738, 0.402173, 0.407831, 0.405895, 0.409172,
1332 0.408934, 0.405915, 0.408486, 0.407320, 0.407437, 0.405444, 0.408400, 0.410909, 0.412427,
1333 0.409881, 0.411021, 0.413001, 0.410369, 0.414702, 0.413372, 0.413095, 0.410972, 0.416346,
1334 0.416095, 0.414132, 0.414344, 0.416952, 0.415197, 0.417583, 0.416582, 0.416622, 0.416895,
1335 0.416576, 0.415551, 0.417925, 0.414838, 0.417051, 0.416831, 0.420000, 0.419132, 0.418173,
1336 0.417645, 0.419679, 0.419866, 0.419581, 0.421531, 0.420878, 0.422737, 0.421872, 0.421304,
1337 0.425486, 0.424434, 0.420842, 0.426753, 0.422761, 0.422178, 0.422372, 0.424173, 0.425582,
1338 0.425080, 0.425831, 0.423551, 0.422949, 0.425784, 0.427977, 0.427948, 0.426368, 0.425138,
1339 0.425351, 0.428643, 0.428148, 0.427488, 0.431704, 0.430167, 0.429655, 0.429584, 0.425458,
1340 0.430728, 0.429845, 0.431145, 0.429180, 0.428874, 0.430720, 0.430024, 0.432034, 0.431359,
1341 0.431535, 0.432995, 0.432425, 0.432454, 0.433140, 0.432574, 0.433814, 0.433348, 0.432886,
1342 0.435472, 0.436517, 0.432681, 0.436999, 0.435182, 0.434834, 0.435478, 0.438255, 0.436650,
1343 0.434464, 0.438530, 0.434077, 0.436471, 0.434012, 0.436822, 0.437505, 0.440135, 0.438322,
1344 0.438032, 0.439001, 0.440270, 0.438661, 0.439233, 0.439274, 0.437945, 0.443080, 0.439191,
1345 0.438233, 0.440415, 0.441063, 0.440926, 0.440929, 0.439731, 0.443584, 0.439729, 0.441597,
1346 0.442615, 0.444637, 0.443180, 0.440789, 0.440261, 0.442202, 0.445081, 0.445484, 0.445415,
1347 0.445532, 0.442806, 0.444188, 0.441073, 0.444299, 0.445897, 0.445279, 0.442830, 0.445506,
1348 0.445272, 0.447267, 0.443522, 0.445519, 0.446459, 0.446753, 0.446377, 0.446129, 0.446383,
1349 0.448556, 0.446593, 0.445293, 0.449199, 0.447590, 0.445968, 0.447482, 0.448474, 0.449890,
1350 0.450004, 0.447765, 0.449274, 0.450652, 0.448210, 0.449360, 0.449577, 0.448575, 0.452112,
1351 0.448780, 0.451393, 0.450200, 0.452018, 0.451182, 0.452050, 0.451748, 0.451377, 0.451402,
1352 0.448810, 0.452311, 0.452909, 0.452491, 0.452418, 0.454190, 0.454420, 0.452121, 0.452307,
1353 0.456857, 0.453506, 0.454058, 0.457203, 0.454394, 0.453596, 0.452240, 0.453692, 0.456516,
1354 0.453753, 0.455541, 0.452702, 0.456481, 0.452226, 0.454280, 0.454855, 0.456297, 0.456482,
1355 0.454154, 0.455387, 0.454748, 0.455764, 0.457282, 0.455487, 0.454822, 0.454257, 0.457678,
1356 0.454225, 0.458689, 0.456123, 0.457011, 0.457386, 0.458351, 0.458638, 0.456164, 0.455884,
1357 0.458525, 0.457575, 0.458340, 0.458912, 0.457836, 0.461734, 0.457545, 0.460755, 0.460960,
1358 0.459226, 0.458613, 0.461078, 0.460958, 0.460337, 0.460237, 0.461190, 0.460760, 0.457911,
1359 0.461310, 0.459657, 0.461960, 0.461040, 0.459578, 0.461650, 0.461550, 0.461251, 0.461054,
1360 0.463082, 0.461732, 0.461324, 0.462547, 0.461261, 0.461629, 0.464067, 0.462430, 0.462525,
1361 0.464232, 0.462921, 0.463202, 0.465558, 0.462914, 0.461698, 0.463963, 0.463040, 0.464275,
1362 0.461940, 0.462913, 0.465261, 0.461500, 0.463679, 0.463354, 0.465205, 0.464529, 0.462220,
1363 0.464279, 0.463427, 0.465387, 0.465288, 0.464839, 0.464926, 0.466100, 0.465531, 0.466187,
1364 0.464647, 0.466285, 0.465461, 0.464134, 0.466783, 0.466763, 0.466183, 0.467089, 0.464497,
1365 0.466080, 0.466109, 0.468166, 0.466984, 0.465335, 0.466721, 0.466856, 0.465113, 0.468377,
1366 0.467904, 0.464546, 0.468787, 0.465648, 0.469841, 0.469477, 0.466311, 0.468700, 0.465183,
1367 0.466559, 0.470433, 0.468563, 0.468109, 0.466980, 0.467567, 0.467670, 0.466991, 0.467992,
1368 0.468784, 0.469406, 0.469652, 0.468527, 0.470460, 0.467308, 0.470693, 0.469539, 0.468000,
1369 0.469295, 0.467038, 0.471908, 0.468829, 0.470663, 0.469266, 0.468975, 0.470222, 0.468649,
1370 0.469507, 0.472307, 0.471611, 0.470419, 0.471181, 0.471140, 0.473187, 0.471086, 0.469801,
1371 0.472234, 0.472131, 0.468996, 0.470229, 0.471597, 0.469625, 0.472230, 0.470164, 0.468404,
1372 0.472264, 0.471336, 0.471597, 0.472280, 0.471256, 0.473151, 0.471863, 0.474458, 0.471956,
1373 0.473099, 0.473956, 0.471725, 0.472809, 0.473065, 0.473180, 0.470611, 0.473614, 0.474263,
1374 0.472792, 0.473543, 0.472656, 0.469728, 0.473431, 0.474538, 0.475322, 0.474962, 0.473598,
1375 0.474114, 0.473486, 0.472934, 0.473252, 0.477149, 0.471719, 0.476383, 0.473076, 0.473952,
1376 0.473104, 0.472459, 0.474433, 0.474494, 0.473588, 0.473839, 0.478113, 0.472435, 0.475571,
1377 0.475194, 0.475626, 0.474617, 0.474520, 0.474472, 0.476437, 0.474512, 0.474497, 0.474628,
1378 0.476203, 0.475698, 0.473907, 0.477144, 0.479000, 0.475553, 0.477481, 0.473998, 0.476672,
1379 0.477115, 0.477114, 0.476282, 0.476152, 0.477009, 0.479854, 0.474354, 0.477645, 0.477517,
1380 0.477111, 0.474843, 0.476173, 0.477321, 0.477384, 0.477880, 0.475726, 0.476004, 0.478204,
1381 0.475586, 0.477973, 0.477935, 0.480640, 0.478234, 0.476349, 0.477493, 0.476994, 0.479815,
1382 0.477771, 0.476333, 0.476325, 0.478245, 0.477284, 0.479238, 0.478339, 0.478966, 0.478012,
1383 0.479304, 0.480148, 0.476125, 0.481267, 0.479801, 0.476720, 0.478898, 0.479284, 0.479153,
1384 0.480157, 0.478681, 0.479712, 0.478993, 0.479943, 0.478349, 0.478930, 0.478052, 0.477173,
1385 0.479244, 0.480454, 0.479128, 0.480530, 0.477843, 0.478369, 0.478561, 0.478639, 0.479191,
1386 0.481763, 0.481321, 0.480979, 0.479702, 0.479777, 0.479384, 0.477571, 0.481880, 0.478615,
1387 0.481303, 0.478783, 0.479384, 0.480517, 0.481928, 0.481199, 0.479041, 0.479188, 0.481491,
1388 0.482840, 0.478766, 0.481941, 0.481298, 0.478105, 0.482933, 0.479744, 0.483361, 0.482332,
1389 0.482556, 0.482057, 0.483616, 0.480599, 0.482245, 0.481091, 0.480871, 0.481938, 0.480678,
1390 0.481851, 0.482902, 0.482158, 0.480187, 0.481772, 0.484967, 0.483094, 0.482133, 0.483929,
1391 0.483354, 0.483382, 0.483964, 0.479941, 0.481375, 0.480255, 0.482184, 0.482541, 0.482032,
1392 0.483484, 0.479492, 0.483305, 0.481070, 0.483573, 0.485689, 0.485767, 0.484221, 0.481365,
1393 0.482440, 0.481507, 0.483418, 0.480978 };
1394 double width_a1 = a1width[iii];
1395 a[0] = 1;
1396 a[1] = 0;
1397 b[0] = mass2 - sa;
1398 // b[1] = -mass*width*wid(mass2,mass,sa,sb,sc,r2,l);
1399 b[1] = -mass * width_a1;
1400 Com_Divide( a, b, prop );
1401}
1402void EvtDToKppipipi0::propagatorGS( double mass2, double mass, double width, double sa,
1403 double sb, double sc, double r2, double prop[2] ) {
1404
1405 double GS1 = 0.636619783;
1406 double GS2 = 0.01860182466;
1407 double GS3 = 0.1591549458;
1408 double GS4 = 0.00620060822;
1409 double a[2], b[2];
1410 double tmp = sb - sc;
1411 double tmp1 = sa + tmp;
1412 double q2 = fabs( 0.25 * tmp1 * tmp1 / sa - sb );
1413 // if(q2<0) q2 = 1e-16;
1414
1415 double tmp2 = mass2 + tmp;
1416 double q02 = fabs( 0.25 * tmp2 * tmp2 / mass2 - sb );
1417 // if(q02<0) q02 = 1e-16;
1418
1419 double q = sqrt( q2 );
1420 double q0 = sqrt( q02 );
1421 double m = sqrt( sa );
1422 double q03 = q0 * q02;
1423 double tmp3 = log( mass + 2 * q0 ) + 1.2760418309; // log(mass_2Pion) = 1.2760418309;
1424
1425 double h = GS1 * q / m * ( log( m + 2 * q ) + 1.2760418309 );
1426 double h0 = GS1 * q0 / mass * tmp3;
1427 double dh = h0 * ( 0.125 / q02 - 0.5 / mass2 ) + GS3 / mass2;
1428 double d = GS2 / q02 * tmp3 + GS3 * mass / q0 - GS4 * mass / q03;
1429 double f = mass2 / q03 * ( q2 * ( h - h0 ) + ( mass2 - sa ) * q02 * dh );
1430
1431 a[0] = 1.0 + d * width / mass;
1432 a[1] = 0.0;
1433 b[0] = mass2 - sa + width * f;
1434 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r2 );
1435 Com_Divide( a, b, prop );
1436}
1437void EvtDToKppipipi0::rhoab( double sa, double sb, double sc, double res[2] ) {
1438 double tmp = sa + sb - sc;
1439 double q = 0.25 * tmp * tmp / sa - sb;
1440 if ( q >= 0 )
1441 {
1442 res[0] = 2.0 * sqrt( q / sa );
1443 res[1] = 0.0;
1444 }
1445 else
1446 {
1447 res[0] = 0.0;
1448 res[1] = 2.0 * sqrt( -q / sa );
1449 }
1450}
1451void EvtDToKppipipi0::rho4Pi( double sa, double res[2] ) {
1452 double temp = 1.0 - 0.3116765584 / sa; // 0.3116765584=0.13957*0.13957*16
1453 if ( temp >= 0 )
1454 {
1455 res[0] = sqrt( temp ) / ( 1.0 + exp( 9.8 - 3.5 * sa ) );
1456 res[1] = 0.0;
1457 }
1458 else
1459 {
1460 res[0] = 0.0;
1461 res[1] = sqrt( -temp ) / ( 1.0 + exp( 9.8 - 3.5 * sa ) );
1462 }
1463}
1464void EvtDToKppipipi0::propagatorsigma500( double sa, double sb, double sc, double prop[2] ) {
1465 double f = 0.5843 + 1.6663 * sa;
1466 const double M = 0.9264;
1467 const double mass2 = 0.85821696; // M*M
1468 const double mpi2d2 = 0.00973989245;
1469 double g1 = f * ( sa - mpi2d2 ) / ( mass2 - mpi2d2 ) * exp( ( mass2 - sa ) / 1.082 );
1470 double rho1s[2], rho1M[2], rho2s[2], rho2M[2], rho1[2], rho2[2];
1471 rhoab( sa, sb, sc, rho1s );
1472 rhoab( mass2, sb, sc, rho1M );
1473 rho4Pi( sa, rho2s );
1474 rho4Pi( mass2, rho2M );
1475 Com_Divide( rho1s, rho1M, rho1 );
1476 Com_Divide( rho2s, rho2M, rho2 );
1477 double a[2], b[2];
1478 a[0] = 1.0;
1479 a[1] = 0.0;
1480 b[0] = mass2 - sa + M * ( g1 * rho1[1] + 0.0024 * rho2[1] );
1481 b[1] = -M * ( g1 * rho1[0] + 0.0024 * rho2[0] );
1482 Com_Divide( a, b, prop );
1483}
1484void EvtDToKppipipi0::Flatte_rhoab( double sa, double sb, double sc, double rho[2] ) {
1485 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
1486 if ( q > 0 )
1487 {
1488 rho[0] = 2 * sqrt( q / sa );
1489 rho[1] = 0;
1490 }
1491 else if ( q < 0 )
1492 {
1493 rho[0] = 0;
1494 rho[1] = 2 * sqrt( -q / sa );
1495 }
1496}
1497void EvtDToKppipipi0::propagator980( double mass, double sx, double* sb, double* sc,
1498 double prop[2] ) {
1499 double unit[2] = { 1.0 };
1500 double ci[2] = { 0, 1 };
1501 double rho1[2];
1502 Flatte_rhoab( sx, sb[0], sc[0], rho1 );
1503 double rho2[2];
1504 Flatte_rhoab( sx, sb[1], sc[1], rho2 );
1505 double gK_f980 = 0.69465, gPi_f980 = 0.165;
1506 double tmp1[2] = { gK_f980, 0 };
1507 double tmp11[2];
1508 double tmp2[2] = { gPi_f980, 0 };
1509 double tmp22[2];
1510 Com_Multi( tmp1, rho1, tmp11 );
1511 Com_Multi( tmp2, rho2, tmp22 );
1512 double tmp3[2] = { tmp11[0] + tmp22[0], tmp11[1] + tmp22[1] };
1513 double tmp31[2];
1514 Com_Multi( tmp3, ci, tmp31 );
1515 double tmp4[2] = { mass * mass - sx - tmp31[0], -1.0 * tmp31[1] };
1516 Com_Divide( unit, tmp4, prop );
1517}
1518void EvtDToKppipipi0::KPiSLASS( double sa, double sb, double sc, double prop[2] ) {
1519 const double m1430 = 1.441;
1520 const double sa0 = 2.076481; // m1430*m1430;
1521 const double w1430 = 0.193;
1522 const double Lass1 = 0.25 / sa0;
1523 double tmp = sb - sc;
1524 double tmp1 = sa0 + tmp;
1525 double q0 = fabs( Lass1 * tmp1 * tmp1 - sb );
1526 // if(q0<0) q0 = 1e-16;
1527 double tmp2 = sa + tmp;
1528 double qs = fabs( 0.25 * tmp2 * tmp2 / sa - sb );
1529 double q = sqrt( qs );
1530 double width = w1430 * q * m1430 / sqrt( sa * q0 );
1531 double temp_R = atan( m1430 * width / ( sa0 - sa ) );
1532 if ( temp_R < 0 ) temp_R += math_pi;
1533 double deltaR = -109.7 * math_pi / 180.0 + temp_R;
1534 double temp_F =
1535 atan( 0.226 * q / ( 2.0 - 3.8194 * qs ) ); // 2.0*0.113 = 0.226; -33.8*0.113 = -3.8194
1536 if ( temp_F < 0 ) temp_F += math_pi;
1537 double deltaF = 0.1 * math_pi / 180.0 + temp_F;
1538 double deltaS = deltaR + 2.0 * deltaF;
1539 double t1 = 0.96 * sin( deltaF );
1540 double t2 = sin( deltaR );
1541 double CF[2], CS[2];
1542 CF[0] = cos( deltaF );
1543 CF[1] = sin( deltaF );
1544 CS[0] = cos( deltaS );
1545 CS[1] = sin( deltaS );
1546 prop[0] = t1 * CF[0] + t2 * CS[0];
1547 prop[1] = t1 * CF[1] + t2 * CS[1];
1548}
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
a1width
double a1width[3000]
Definition EvtDsTopipipipi0pi0.cc:323
unit
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per saves r n generator level $ !Flag for chat level in !Latex Output unit
Definition FoamA.h:90
sin
double sin(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:185
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
q
****INTEGER imax DOUBLE PRECISION m_pi *DOUBLE PRECISION m_amfin DOUBLE PRECISION m_Chfin DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_sinw2 DOUBLE PRECISION m_GFermi DOUBLE PRECISION m_MfinMin DOUBLE PRECISION m_ta2 INTEGER m_out INTEGER m_KeyFSR INTEGER m_KeyQCD *COMMON c_Semalib $ !copy of input $ !CMS energy $ !beam mass $ !final mass $ !beam charge $ !final charge $ !smallest final mass $ !Z mass $ !Z width $ !EW mixing angle $ !Gmu Fermi $ alphaQED at q
Definition KKsem.h:33
CS
TCrossPart * CS
Definition Mcgpj.cxx:51
EvtDToKppipipi0::getName
void getName(std::string &name)
Definition EvtDToKppipipi0.cc:34
EvtDToKppipipi0::clone
EvtDecayBase * clone()
Definition EvtDToKppipipi0.cc:36
EvtDToKppipipi0::decay
void decay(EvtParticle *p)
Definition EvtDToKppipipi0.cc:139
EvtDToKppipipi0::~EvtDToKppipipi0
virtual ~EvtDToKppipipi0()
Definition EvtDToKppipipi0.cc:32
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