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