BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDToKppipi.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: EvtDToKppipi.cc
11//
12// Description: D+ -> K+ pi+ pi-
13// See
14// https://indico.ihep.ac.cn/event/18455/contributions/146055/attachments/74328/91168/DCS_D%2BtoK%2Bpi%2Bpi-.pdf
15//
16// Modification history:
17// Liaoyuan Dong Feb 17 23:19:33 2024 Module updated
18// Oct 10 00:03:38 2023 Module created
19//------------------------------------------------------------------------
20#include "EvtDToKppipi.hh"
21#include "../EvtGenBase/EvtConst.hh"
22#include "../EvtGenBase/EvtDecayTable.hh"
23#include "../EvtGenBase/EvtGenKine.hh"
24#include "../EvtGenBase/EvtPDL.hh"
25#include "../EvtGenBase/EvtParticle.hh"
26#include "../EvtGenBase/EvtPatches.hh"
27#include "../EvtGenBase/EvtReport.hh"
28#include "../EvtGenBase/EvtResonance.hh"
29#include "../EvtGenBase/EvtResonance2.hh"
30#include "TMath.h"
31#include <fstream>
32#include <math.h>
33#include <stdlib.h>
34#include <string>
35#include <vector>
36using std::endl;
37
38EvtDToKppipi::~EvtDToKppipi() {}
39
40void EvtDToKppipi::getName( std::string& model_name ) { model_name = "DToKppipi"; }
41
42EvtDecayBase* EvtDToKppipi::clone() { return new EvtDToKppipi; }
43
44void EvtDToKppipi::init() {
45 // check that there are 0 arguments
46 checkNArg( 0 );
47 checkNDaug( 3 );
48 checkSpinParent( EvtSpinType::SCALAR );
49 checkSpinDaughter( 0, EvtSpinType::SCALAR );
50 checkSpinDaughter( 1, EvtSpinType::SCALAR );
51 checkSpinDaughter( 2, EvtSpinType::SCALAR );
52
53 pi180inv = 1.0 / EvtConst::radToDegrees;
54 phi[0] = 0.0;
55 phi[1] = 3.365002324407079470347526;
56 phi[2] = 0.6621583377125883629332748;
57 phi[3] = 0.4339644541377705166951273;
58 phi[4] = 2.496788137814539787484591;
59 phi[5] = 5.337638327720493514050304;
60
61 rho[0] = 1.0;
62 rho[1] = 2.176692936463668459623477;
63 rho[2] = 1.815490493394110060876301;
64 rho[3] = 0.3832709830473302048403639;
65 rho[4] = 1.12572937386981664076302;
66 rho[5] = -1.190394848723933307610423;
67
68 modetype[0] = 13; /////K*(892)
69 modetype[1] = 13; /////K*(1410)
70 modetype[2] = 13; /////K*2(1430)
71 modetype[3] = 23; /////PiPi
72 modetype[4] = 23; /////rho(770)
73 modetype[5] = 23; /////f2(1270)
74
75 // std::cout << "Initializing EvtDToKppipi" << std::endl;
76 // for (int i=0; i<9; i++) {
77 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
78 // }
79
80 width[0] = 0.0473000000000000017652546;
81 width[1] = 0.2320000000000000117683641;
82 width[2] = 0.1089999999999999996669331;
83 width[3] = 0.0200000000000000004163336;
84 width[4] = 0.1474000000000000032418512;
85 width[5] = 0.1867000000000000048405724;
86
87 mass[0] = 0.8955499999999999571898002;
88 mass[1] = 1.4139999999999999236166559;
89 mass[2] = 1.4323999999999998955502178;
90 mass[3] = 0.9000000000000000222044605;
91 mass[4] = 0.7752599999999999491606673;
92 mass[5] = 1.2755000000000000781597009;
93
94 mDp = 1.86966;
95 mK0 = 0.497611;
96 mKa = 0.49368;
97 mPi = 0.13957;
98 meta = 0.54775;
99 mK02 = 0.237616707;
100 mPi2 = 0.01947978;
101 mass_EtaP = 0.95778;
102 mass_Kaon = 0.49368;
103
104 math_pi = 3.1415926;
105 mass_Pion = 0.13957;
106 mass_Pion2 = 0.0194797849;
107 mass_2Pion = 0.27914;
108 math_2pi = 6.2831852;
109 rD2 = 25.0; // 5*5
110 rRes2 = 9.0; // 3*3
111 g1 = 0.5468;
112 g2 = 0.23;
113 GS1 = 0.636619783;
114 GS2 = 0.01860182466;
115 GS3 = 0.1591549458;
116 GS4 = 0.00620060822;
117
118 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
119 for ( int i = 0; i < 4; i++ )
120 {
121 for ( int j = 0; j < 4; j++ ) { G[i][j] = GG[i][j]; }
122 }
123}
124
125void EvtDToKppipi::initProbMax() { setProbMax( 717.25 ); }
126
127void EvtDToKppipi::decay( EvtParticle* p ) {
128 /*
129 double maxprob = 0.0;
130 for(int ir=0;ir<=60000000;ir++){
131 p->initializePhaseSpace(getNDaug(),getDaugs());
132 EvtVector4R D1 = p->getDaug(0)->getP4();
133 EvtVector4R D2 = p->getDaug(1)->getP4();
134 EvtVector4R D3 = p->getDaug(2)->getP4();
135
136 double P1[4], P2[4], P3[4];
137 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
138 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
139 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
140
141 double value;
142 int nstates=6;
143 int g0[6]={1,1,1,6,1,2};
144 int spin[6]={1,1,2,0,1,2};
145 double r0[6]={3.0,3.0,3.0,3.0,3.0,3.0};
146 double r1[6]={5.0,5.0,5.0,5.0,5.0,5.0};
147
148 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value,r0,r1);
149
150 if (value<0) continue;
151 if(value>maxprob) {
152 maxprob=value;
153 cout << "ir= " << ir << endl;
154 // cout << "double P1[4] = {" << P1[0] <<","<< P1[1] <<","<< P1[2] <<","<< P1[3]
155 <<"};"<< endl;
156 // cout << "double P2[4] = {" << P2[0] <<","<< P2[1] <<","<< P2[2] <<","<< P2[3]
157 <<"};"<< endl;
158 // cout << "double P3[4] = {" << P3[0] <<","<< P3[1] <<","<< P3[2] <<","<< P3[3] <<"};"<<
159 endl; cout << "MAX====> " << maxprob << endl;
160 }
161 }
162 printf("MAXprob = %.10f\n",maxprob);
163 */
164 p->initializePhaseSpace( getNDaug(), getDaugs() );
165 EvtVector4R D1 = p->getDaug( 0 )->getP4();
166 EvtVector4R D2 = p->getDaug( 1 )->getP4();
167 EvtVector4R D3 = p->getDaug( 2 )->getP4();
168
169 double P1[4], P2[4], P3[4];
170 P1[0] = D1.get( 0 );
171 P1[1] = D1.get( 1 );
172 P1[2] = D1.get( 2 );
173 P1[3] = D1.get( 3 );
174 P2[0] = D2.get( 0 );
175 P2[1] = D2.get( 1 );
176 P2[2] = D2.get( 2 );
177 P2[3] = D2.get( 3 );
178 P3[0] = D3.get( 0 );
179 P3[1] = D3.get( 1 );
180 P3[2] = D3.get( 2 );
181 P3[3] = D3.get( 3 );
182
183 double value;
184 int nstates = 6;
185 int g0[6] = { 1, 1, 1, 6, 1, 2 };
186 int spin[6] = { 1, 1, 2, 0, 1, 2 };
187 double r0[6] = { 3.0, 3.0, 3.0, 3.0, 3.0, 3.0 };
188 double r1[6] = { 5.0, 5.0, 5.0, 5.0, 5.0, 5.0 };
189
190 calEva( P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value, r0, r1 );
191 setProb( value );
192 return;
193}
194
195void EvtDToKppipi::Com_Multi( double a1[2], double a2[2], double res[2] ) {
196 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
197 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
198}
199void EvtDToKppipi::Com_Divide( double a1[2], double a2[2], double res[2] ) {
200 double tmp = a2[0] * a2[0] + a2[1] * a2[1];
201 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / tmp;
202 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / tmp;
203}
204double EvtDToKppipi::CalRho4pi( double_t s ) {
205 if ( s >= 1. ) { return sqrt( ( s - 16. * mass_Pion * mass_Pion ) / s ); }
206 else
207 {
208 double_t s0 = 1.2274 + 0.00370909 / ( s * s ) - ( 0.111203 ) / (s)-6.39017 * s +
209 16.8358 * s * s - 21.8845 * s * s * s + 11.3153 * s * s * s * s;
210 double_t gam = s0 * sqrt( 1.0 - ( 16.0 * mass_Pion * mass_Pion ) );
211
212 return gam;
213 }
214}
215
216//------------base---------------------------------
217double EvtDToKppipi::SCADot( double a1[4], double a2[4] ) {
218 double _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
219 return _cal;
220}
221double EvtDToKppipi::barrier( int l, double sa, double sb, double sc, double r, double mass ) {
222 double q = fabs( ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb );
223 double z;
224 z = q * r * r;
225 double sa0;
226 sa0 = mass * mass;
227 double q0 = fabs( ( sa0 + sb - sc ) * ( sa0 + sb - sc ) / ( 4 * sa0 ) - sb );
228 double z0 = q0 * r * r;
229 double F = 0.0;
230 if ( l == 0 ) F = 1;
231 if ( l == 1 ) F = sqrt( ( 1 + z0 ) / ( 1 + z ) );
232 if ( l == 2 ) F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) );
233 return F;
234}
235//------------------spin-------------------------------------------
236void EvtDToKppipi::calt1( double daug1[4], double daug2[4], double t1[4] ) {
237 double p, pq, tmp;
238 double pa[4], qa[4];
239 for ( int i = 0; i < 4; i++ )
240 {
241 pa[i] = daug1[i] + daug2[i];
242 qa[i] = daug1[i] - daug2[i];
243 }
244 p = SCADot( pa, pa );
245 pq = SCADot( pa, qa );
246 tmp = pq / p;
247 for ( int i = 0; i < 4; i++ ) { t1[i] = qa[i] - tmp * pa[i]; }
248}
249void EvtDToKppipi::calt2( double daug1[4], double daug2[4], double t2[4][4] ) {
250 double p, r;
251 double pa[4], t1[4];
252 calt1( daug1, daug2, t1 );
253 r = SCADot( t1, t1 ) / 3.0;
254 for ( int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
255 p = SCADot( pa, pa );
256 for ( int i = 0; i < 4; i++ )
257 {
258 for ( int j = 0; j < 4; j++ )
259 { t2[i][j] = t1[i] * t1[j] - r * ( G[i][j] - pa[i] * pa[j] / p ); }
260 }
261}
262//-------------------prop--------------------------------------------
263void EvtDToKppipi::propagator( double mass2, double mass, double width, double sx,
264 double prop[2] ) {
265 double a[2], b[2];
266 a[0] = 1;
267 a[1] = 0;
268 b[0] = mass2 - sx;
269 b[1] = -mass * width;
270 Com_Divide( a, b, prop );
271}
272double EvtDToKppipi::wid( double mass2, double mass, double sa, double sb, double sc,
273 double r2, int l ) {
274 double widm = 0.;
275 double m = sqrt( sa );
276 double tmp = sb - sc;
277 double tmp1 = sa + tmp;
278 double q = fabs( 0.25 * tmp1 * tmp1 / sa - sb );
279 double tmp2 = mass2 + tmp;
280 double q0 = fabs( 0.25 * tmp2 * tmp2 / mass2 - sb );
281 double z = q * r2;
282 double z0 = q0 * r2;
283 double t = q / q0;
284 if ( l == 0 ) { widm = sqrt( t ) * mass / m; }
285 else if ( l == 1 ) { widm = t * sqrt( t ) * mass / m * ( 1 + z0 ) / ( 1 + z ); }
286 else if ( l == 2 )
287 { widm = t * t * sqrt( t ) * mass / m * ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ); }
288 return widm;
289}
290double EvtDToKppipi::widl1( double mass2, double mass, double sa, double sb, double sc,
291 double r2 ) {
292 double widm = 0.;
293 double m = sqrt( sa );
294 double tmp = sb - sc;
295 double tmp1 = sa + tmp;
296 double q = fabs( 0.25 * tmp1 * tmp1 / sa - sb );
297 double tmp2 = mass2 + tmp;
298 double q0 = fabs( 0.25 * tmp2 * tmp2 / mass2 - sb );
299 double z = q * r2;
300 double z0 = q0 * r2;
301 double F = ( 1 + z0 ) / ( 1 + z );
302 double t = q / q0;
303 widm = t * sqrt( t ) * mass / m * F;
304 return widm;
305}
306void EvtDToKppipi::propagatorRBW( double mass, double width, double sa, double sb, double sc,
307 double r2, int l, double prop[2] ) {
308 double a[2], b[2];
309 double mass2 = mass * mass;
310 a[0] = 1;
311 a[1] = 0;
312 b[0] = mass2 - sa;
313 b[1] = -mass * width * wid( mass2, mass, sa, sb, sc, r2, l );
314 Com_Divide( a, b, prop );
315}
316
317void EvtDToKppipi::propagatorFlatte( double mass, double width, double sa, double prop[2] ) {
318 double q2_Pi, q2_Ka;
319 double rhoPi[2], rhoKa[2];
320
321 q2_Pi = 0.25 * sa - mPi * mPi;
322 q2_Ka = 0.25 * sa - mKa * mKa;
323
324 if ( q2_Pi > 0 )
325 {
326 rhoPi[0] = 2.0 * sqrt( q2_Pi / sa );
327 rhoPi[1] = 0.0;
328 }
329 if ( q2_Pi <= 0 )
330 {
331 rhoPi[0] = 0.0;
332 rhoPi[1] = 2.0 * sqrt( -q2_Pi / sa );
333 }
334
335 if ( q2_Ka > 0 )
336 {
337 rhoKa[0] = 2.0 * sqrt( q2_Ka / sa );
338 rhoKa[1] = 0.0;
339 }
340 if ( q2_Ka <= 0 )
341 {
342 rhoKa[0] = 0.0;
343 rhoKa[1] = 2.0 * sqrt( -q2_Ka / sa );
344 }
345
346 double a[2], b[2];
347 a[0] = 1;
348 a[1] = 0;
349 b[0] = mass * mass - sa + 0.165 * rhoPi[1] + 0.69465 * rhoKa[1];
350 b[1] = -( 0.165 * rhoPi[0] + 0.69465 * rhoKa[0] );
351 Com_Divide( a, b, prop );
352}
353
354void EvtDToKppipi::rhoab( double sa, double sb, double sc, double res[2] ) {
355 double tmp = sa + sb - sc;
356 double q = 0.25 * tmp * tmp / sa - sb;
357 if ( q >= 0 )
358 {
359 res[0] = 2.0 * sqrt( q / sa );
360 res[1] = 0.0;
361 }
362 else
363 {
364 res[0] = 0.0;
365 res[1] = 2.0 * sqrt( -q / sa );
366 }
367}
368void EvtDToKppipi::rho4Pi( double sa, double res[2] ) {
369 double temp = 1.0 - 0.3116765584 / sa; // 0.3116765584=0.13957*0.13957*16
370 if ( temp >= 0 )
371 {
372 res[0] = sqrt( temp ) / ( 1.0 + exp( 9.8 - 3.5 * sa ) );
373 res[1] = 0.0;
374 }
375 else
376 {
377 res[0] = 0.0;
378 res[1] = sqrt( -temp ) / ( 1.0 + exp( 9.8 - 3.5 * sa ) );
379 }
380}
381
382void EvtDToKppipi::propagatorsigma500( double sa, double sb, double sc, double prop[2] ) {
383 double f = 0.5843 + 1.6663 * sa;
384 const double M = 0.9264; // M(f0500)
385 const double mass2 = 0.85821696; // M*M
386 const double mpi2d2 = 0.00973989245; // mass_Pion2/2 = 0.0194797849/2
387 double g1 = f * ( sa - mpi2d2 ) / ( mass2 - mpi2d2 ) * exp( ( mass2 - sa ) / 1.082 );
388 double rho1s[2], rho1M[2], rho2s[2], rho2M[2], rho1[2], rho2[2];
389 rhoab( sa, sb, sc, rho1s );
390 rhoab( mass2, sb, sc, rho1M );
391 rho4Pi( sa, rho2s );
392 rho4Pi( mass2, rho2M );
393 Com_Divide( rho1s, rho1M, rho1 );
394 Com_Divide( rho2s, rho2M, rho2 );
395 double a[2], b[2];
396 a[0] = 1.0;
397 a[1] = 0.0;
398 b[0] = mass2 - sa + M * ( g1 * rho1[1] + 0.0024 * rho2[1] );
399 b[1] = -M * ( g1 * rho1[0] + 0.0024 * rho2[0] );
400 Com_Divide( a, b, prop );
401}
402void EvtDToKppipi::Flatte_rhoab( double sa, double sb, double rho[2] ) {
403 double q = 1.0 - ( 4 * sb / sa );
404
405 if ( q > 0 )
406 {
407 rho[0] = sqrt( q );
408 rho[1] = 0;
409 }
410 else if ( q < 0 )
411 {
412 rho[0] = 0;
413 rho[1] = sqrt( -q );
414 }
415}
416
417void EvtDToKppipi::propagator980( double mass, double sx, double* sb, double prop[2] ) {
418 double gpipi1 = 2.0 / 3.0 * 0.165;
419 double gpipi2 = 1.0 / 3.0 * 0.165;
420 double gKK1 = 0.5 * 0.69465;
421 double gKK2 = 0.5 * 0.69465;
422
423 double unit[2] = { 1.0 };
424 double ci[2] = { 0, 1 };
425 double rho1[2];
426 Flatte_rhoab( sx, sb[0], rho1 );
427 double rho2[2];
428 Flatte_rhoab( sx, sb[1], rho2 );
429 double rho3[2];
430 Flatte_rhoab( sx, sb[2], rho3 );
431 double rho4[2];
432 Flatte_rhoab( sx, sb[3], rho4 );
433
434 double tmp1[2] = { gpipi1, 0 };
435 double tmp11[2];
436 double tmp2[2] = { gpipi2, 0 };
437 double tmp22[2];
438 double tmp3[2] = { gKK1, 0 };
439 double tmp33[2];
440 double tmp4[2] = { gKK2, 0 };
441 double tmp44[2];
442
443 Com_Multi( tmp1, rho1, tmp11 );
444 Com_Multi( tmp2, rho2, tmp22 );
445 Com_Multi( tmp3, rho3, tmp33 );
446 Com_Multi( tmp4, rho4, tmp44 );
447
448 double tmp5[2] = { tmp11[0] + tmp22[0] + tmp33[0] + tmp44[0],
449 tmp11[1] + tmp22[1] + tmp33[1] + tmp44[1] };
450 double tmp51[2];
451 Com_Multi( tmp5, ci, tmp51 );
452 double tmp6[2] = { mass * mass - sx - tmp51[0], -1.0 * tmp51[1] };
453
454 Com_Divide( unit, tmp6, prop );
455}
456
457void EvtDToKppipi::KPiSLASS( double sa, double sb, double sc, double prop[2] ) {
458 const double m1430 = 1.441;
459 const double sa0 = 2.076481; // m1430*m1430;
460 const double w1430 = 0.193;
461 const double Lass1 = 0.25 / sa0;
462 double tmp = sb - sc;
463 double tmp1 = sa0 + tmp;
464 double q0 = fabs( Lass1 * tmp1 * tmp1 - sb );
465 double tmp2 = sa + tmp;
466 double qs = 0.25 * tmp2 * tmp2 / sa - sb;
467 double q = sqrt( qs );
468 double width = w1430 * q * m1430 / sqrt( sa * q0 );
469 double temp_R = atan( m1430 * width / ( sa0 - sa ) );
470 if ( temp_R < 0 ) temp_R += math_pi;
471 double deltaR = -109.7 * math_pi / 180.0 + temp_R;
472 double temp_F =
473 atan( 0.226 * q / ( 2.0 - 3.8194 * qs ) ); // 2.0*0.113 = 0.226; -33.8*0.113 = -3.8194
474 if ( temp_F < 0 ) temp_F += math_pi;
475 double deltaF = 0.1 * math_pi / 180.0 + temp_F;
476 double deltaS = deltaR + 2.0 * deltaF;
477 double t1 = 0.96 * sin( deltaF );
478 double t2 = sin( deltaR );
479 double CF[2], CS[2];
480 CF[0] = cos( deltaF );
481 CF[1] = sin( deltaF );
482 CS[0] = cos( deltaS );
483 CS[1] = sin( deltaS );
484 prop[0] = t1 * CF[0] + t2 * CS[0];
485 prop[1] = t1 * CF[1] + t2 * CS[1];
486}
487//------------GS---used by rho----------------------------
488void EvtDToKppipi::propagatorGS( double mass, double width, double sa, double sb, double sc,
489 double r2, double prop[2] ) {
490 double a[2], b[2];
491 double mass2 = mass * mass;
492 double tmp = sb - sc;
493 double tmp1 = sa + tmp;
494 double q2 = fabs( 0.25 * tmp1 * tmp1 / sa - sb );
495
496 double tmp2 = mass2 + tmp;
497 double q02 = fabs( 0.25 * tmp2 * tmp2 / mass2 - sb );
498
499 double q = sqrt( q2 );
500 double q0 = sqrt( q02 );
501 double m = sqrt( sa );
502 double q03 = q0 * q02;
503 double tmp3 = log( mass + 2 * q0 ) + 1.2760418309; // log(mpi0+mpip) = -1.2926305904;
504
505 double h = GS1 * q / m * ( log( m + 2 * q ) + 1.2760418309 );
506 double h0 = GS1 * q0 / mass * tmp3;
507 double dh = h0 * ( 0.125 / q02 - 0.5 / mass2 ) + GS3 / mass2;
508 double d = GS2 / q02 * tmp3 + GS3 * mass / q0 - GS4 * mass / q03;
509 double f = mass2 / q03 * ( q2 * ( h - h0 ) + ( mass2 - sa ) * q02 * dh );
510
511 a[0] = 1.0 + d * width / mass;
512 a[1] = 0.0;
513 b[0] = mass2 - sa + width * f;
514 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r2 );
515 Com_Divide( a, b, prop );
516}
517
518double EvtDToKppipi::DDalitz( double P1[4], double P2[4], double P3[4], int Ang,
519 double mass ) {
520 double pR[4], pD[4];
521 double temp_PDF, v_re;
522 temp_PDF = 0.0;
523 v_re = 0.0;
524 double B[2], s1, s2, s3, sR, sD;
525 for ( int i = 0; i < 4; i++ )
526 {
527 pR[i] = P1[i] + P2[i];
528 pD[i] = pR[i] + P3[i];
529 }
530 s1 = SCADot( P1, P1 );
531 s2 = SCADot( P2, P2 );
532 s3 = SCADot( P3, P3 );
533 sR = SCADot( pR, pR );
534 sD = SCADot( pD, pD );
535 int G[4][4];
536 for ( int i = 0; i != 4; i++ )
537 {
538 for ( int j = 0; j != 4; j++ )
539 {
540 if ( i == j )
541 {
542 if ( i == 0 ) G[i][j] = 1;
543 else G[i][j] = -1;
544 }
545 else G[i][j] = 0;
546 }
547 }
548 if ( Ang == 0 )
549 {
550 B[0] = 1;
551 B[1] = 1;
552 temp_PDF = 1;
553 }
554 if ( Ang == 1 )
555 {
556 B[0] = barrier( 1, sR, s1, s2, 3.0, mass );
557 B[1] = barrier( 1, sD, sR, s3, 5.0, mDp );
558 double t1[4], T1[4];
559 calt1( P1, P2, t1 );
560 calt1( pR, P3, T1 );
561 temp_PDF = 0;
562 for ( int i = 0; i < 4; i++ ) { temp_PDF += t1[i] * T1[i] * G[i][i]; }
563 }
564 if ( Ang == 2 )
565 {
566 B[0] = barrier( 2, sR, s1, s2, 3.0, mass );
567 B[1] = barrier( 2, sD, sR, s3, 5.0, mDp );
568 double t2[4][4], T2[4][4];
569 calt2( P1, P2, t2 );
570 calt2( pR, P3, T2 );
571 temp_PDF = 0;
572 for ( int i = 0; i < 4; i++ )
573 {
574 for ( int j = 0; j < 4; j++ ) { temp_PDF += t2[i][j] * T2[j][i] * G[i][i] * G[j][j]; }
575 }
576 }
577 v_re = temp_PDF * B[0] * B[1];
578 return v_re;
579}
580
581void EvtDToKppipi::rhoMTX( int i, int j, double s, double Rho[2] ) {
582
583 double rhoijx;
584 double rhoijy;
585 double mpi = 0.13957;
586 if ( i == j && i == 0 )
587 {
588 double m2 = 0.13957 * 0.13957;
589 if ( ( 1 - ( 4 * m2 ) / s ) > 0 )
590 {
591 rhoijx = sqrt( 1.0f - ( 4 * m2 ) / s );
592 rhoijy = 0;
593 }
594 else
595 {
596 rhoijy = sqrt( ( 4 * m2 ) / s - 1.0f );
597 rhoijx = 0;
598 }
599 }
600 if ( i == j && i == 1 )
601 {
602 double m2 = 0.49368 * 0.49368;
603 if ( ( 1 - ( 4 * m2 ) / s ) > 0 )
604 {
605 rhoijx = sqrt( 1.0f - ( 4 * m2 ) / s );
606 rhoijy = 0;
607 }
608 else
609 {
610 rhoijy = sqrt( ( 4 * m2 ) / s - 1.0f );
611 rhoijx = 0;
612 }
613 }
614 if ( i == j && i == 2 )
615 {
616 rhoijx = CalRho4pi( s );
617 rhoijy = 0;
618 }
619 if ( i == j && i == 3 )
620 {
621 double m2 = 0.547862 * 0.547862;
622 if ( ( 1 - ( 4 * m2 ) / s ) > 0 )
623 {
624 rhoijx = sqrt( 1.0f - ( 4 * m2 ) / s );
625 rhoijy = 0;
626 }
627 else
628 {
629 rhoijy = sqrt( ( 4 * m2 ) / s - 1.0f );
630 rhoijx = 0;
631 }
632 }
633 if ( i == j && i == 4 )
634 {
635 double m_1 = 0.547862;
636 double m_2 = 0.95778;
637 double mp2 = ( m_1 + m_2 ) * ( m_1 + m_2 );
638 double mm2 = ( m_1 - m_2 ) * ( m_1 - m_2 );
639 if ( ( 1 - mp2 / s ) > 0 )
640 {
641 rhoijx = sqrt( 1.0f - mp2 / s );
642 rhoijy = 0;
643 }
644 else
645 {
646 rhoijy = sqrt( mp2 / s - 1.0f );
647 rhoijx = 0;
648 }
649 }
650
651 if ( i != j )
652 {
653 rhoijx = 0;
654 rhoijy = 0;
655 }
656 Rho[0] = rhoijx;
657 Rho[1] = rhoijy;
658}
659
660/* Kij = (sum_a gigj/(ma^2-s) + fij(1-s0)/(s-s0))((s-0.5sampi2)(1-sa0)/(s-sa0)) */
661void EvtDToKppipi::KMTX( int i, int j, double s, double KM[2] ) {
662
663 double Kijx;
664 double Kijy;
665 double mpi = 0.13957;
666 double m[5] = { 0.65100, 1.20360, 1.55817, 1.21000, 1.82206 };
667
668 double g1[5] = { 0.22889, -0.55377, 0.00000, -0.39899, -0.34639 };
669 double g2[5] = { 0.94128, 0.55095, 0.00000, 0.39065, 0.31503 };
670 double g3[5] = { 0.36856, 0.23888, 0.55639, 0.18340, 0.18681 };
671 double g4[5] = { 0.33650, 0.40907, 0.85679, 0.19906, -0.00984 };
672 double g5[5] = { 0.18171, -0.17558, -0.79658, -0.00355, 0.22358 };
673
674 double f1[5] = { 0.23399, 0.15044, -0.20545, 0.32825, 0.35412 };
675
676 double upimag[5] = { 0, 0, 0, 0, 0 };
677
678 for ( int k = 0; k < 5; k++ ) { upimag[k] = 0; }
679 double ss0 = -3.92637;
680 double sA = 1.0; // v1
681 double sA0 = -0.15;
682
683 if ( i == 0 || j == 0 )
684 {
685 Kijx = ( g1[i] * g1[j] / ( m[0] * m[0] - s ) + g2[i] * g2[j] / ( m[1] * m[1] - s ) +
686 g3[i] * g3[j] / ( m[2] * m[2] - s ) + g4[i] * g4[j] / ( m[3] * m[3] - s ) +
687 g5[i] * g5[j] / ( m[4] * m[4] - s ) + f1[j] * ( 1 - ss0 ) / ( s - ss0 ) ) *
688 ( 1 - sA0 ) / ( s - sA0 ) * ( s - sA * mpi * mpi * 0.5 );
689 Kijy =
690 ( g1[i] * g1[j] * upimag[0] + g2[i] * g2[j] * upimag[1] + g3[i] * g3[j] * upimag[2] +
691 g4[i] * g4[j] * upimag[3] + g5[i] * g5[j] * upimag[4] ) *
692 ( 1 - sA0 ) / ( s - sA0 ) * ( s - sA * mpi * mpi * 0.5 );
693 }
694
695 else
696 {
697 Kijx = ( g1[i] * g1[j] / ( m[0] * m[0] - s ) + g2[i] * g2[j] / ( m[1] * m[1] - s ) +
698 g3[i] * g3[j] / ( m[2] * m[2] - s ) + g4[i] * g4[j] / ( m[3] * m[3] - s ) +
699 g5[i] * g5[j] / ( m[4] * m[4] - s ) ) *
700 ( 1 - sA0 ) / ( s - sA0 ) * ( s - sA * mpi * mpi * 0.5 );
701 Kijy =
702 ( g1[i] * g1[j] * upimag[0] + g2[i] * g2[j] * upimag[1] + g3[i] * g3[j] * upimag[2] +
703 g4[i] * g4[j] * upimag[3] + g5[i] * g5[j] * upimag[4] ) *
704 ( 1 - sA0 ) / ( s - sA0 ) * ( s - sA * mpi * mpi * 0.5 );
705 }
706
707 KM[0] = Kijx;
708 KM[1] = Kijy;
709}
710
711void EvtDToKppipi::IMTX( int i, int j, double IMTX[2] ) {
712
713 double Iijx;
714 double Iijy;
715 if ( i == j )
716 {
717 Iijx = 1;
718 Iijy = 0;
719 }
720 else
721 {
722 Iijx = 0;
723 Iijy = 0;
724 }
725 IMTX[0] = Iijx;
726 IMTX[1] = Iijy;
727}
728
729/* F = I - i*K*rho */
730void EvtDToKppipi::FMTX( double Kijx, double Kijy, double rhojjx, double rhojjy, int i, int j,
731 double FM[2] ) {
732
733 double Fijx;
734 double Fijy;
735
736 double tmpx = Kijx * rhojjx - Kijy * rhojjy;
737 double tmpy = Kijy * rhojjx + Kijx * rhojjy;
738
739 double imtx[2];
740 IMTX( i, j, imtx );
741 Fijx = imtx[0] + tmpy;
742 Fijy = -tmpx;
743
744 FM[0] = Fijx;
745 FM[1] = Fijy;
746}
747
748void EvtDToKppipi::PVTR( int ID, double s, double PV[2] ) {
749
750 double VPix;
751 double VPiy;
752 double m[5] = { 0.65100, 1.20360, 1.55817, 1.21000, 1.82206 };
753
754 double g[5][5] = { { 0.22889, -0.55377, 0.00000, -0.39899, -0.34639 },
755 { 0.94128, 0.55095, 0.00000, 0.39065, 0.31503 },
756 { 0.36856, 0.23888, 0.55639, 0.18340, 0.18681 },
757 { 0.33650, 0.40907, 0.85679, 0.19906, -0.00984 },
758 { 0.18171, -0.17558, -0.79658, -0.00355, 0.22358 } };
759
760 double betax[5], betay[5], fprodx[5], fprody[5];
761
762 betax[0] = 8.5 * cos( 68.5 * 3.1415926 / 180.0 );
763 betay[0] = 8.5 * sin( 68.5 * 3.1415926 / 180.0 );
764 betax[1] = 12.2 * cos( 24.0 * 3.1415926 / 180.0 );
765 betay[1] = 12.2 * sin( 24.0 * 3.1415926 / 180.0 );
766 betax[2] = 29.2 * cos( -0.1 * 3.1415926 / 180.0 );
767 betay[2] = 29.2 * sin( -0.1 * 3.1415926 / 180.0 );
768 betax[3] = 10.8 * cos( -51.9 * 3.1415926 / 180.0 );
769 betay[3] = 10.8 * sin( -51.9 * 3.1415926 / 180.0 );
770 betax[4] = 0.0;
771 betay[4] = 0.0;
772
773 fprodx[0] = 8.0 * cos( -126.0 * 3.1415926 / 180.0 );
774 fprody[0] = 8.0 * sin( -126.0 * 3.1415926 / 180.0 );
775 fprodx[1] = 26.3 * cos( -152.3 * 3.1415926 / 180.0 );
776 fprody[1] = 26.3 * sin( -152.3 * 3.1415926 / 180.0 );
777 fprodx[2] = 33.0 * cos( -93.2 * 3.1415926 / 180.0 );
778 fprody[2] = 33.0 * sin( -93.2 * 3.1415926 / 180.0 );
779 fprodx[3] = 26.2 * cos( -121.4 * 3.1415926 / 180.0 );
780 fprody[3] = 26.2 * sin( -121.4 * 3.1415926 / 180.0 );
781 fprodx[4] = 0.0;
782 fprody[4] = 0.0;
783
784 double V0x = 0.0, V0y = 0.0, V1x = 0.0, V1y = 0.0;
785 double s0_prod = -0.07;
786
787 for ( int k = 0; k < 5; k++ )
788 {
789 V0x += betax[k] * g[k][ID] / ( m[k] * m[k] - s );
790 V0y += betay[k] * g[k][ID] / ( m[k] * m[k] - s );
791 }
792 V1x += ( 1. - s0_prod ) / ( s - s0_prod ) * fprodx[ID];
793 V1y += ( 1. - s0_prod ) / ( s - s0_prod ) * fprody[ID];
794
795 VPix = V0x + V1x;
796 VPiy = V0y + V1y;
797
798 PV[0] = VPix;
799 PV[1] = VPiy;
800}
801
802/* inverse for Matrix (I - i*rho*K) using LUP */
803void EvtDToKppipi::FINVMTX( double s, double* FINVx, double* FINVy ) {
804
805 int P[5] = { 0, 1, 2, 3, 4 };
806
807 double Fx[5][5];
808 double Fy[5][5];
809
810 double Ux[5][5];
811 double Uy[5][5];
812 double Lx[5][5];
813 double Ly[5][5];
814
815 double UIx[5][5];
816 double UIy[5][5];
817 double LIx[5][5];
818 double LIy[5][5];
819
820 double rho[2];
821 double KM[2];
822 for ( int k = 0; k < 5; k++ )
823 {
824 rhoMTX( k, k, s, rho );
825 double rhokkx = rho[0];
826 double rhokky = rho[1];
827 Ux[k][k] = rhokkx;
828 Uy[k][k] = rhokky;
829 for ( int l = k; l < 5; l++ )
830 {
831 KMTX( k, l, s, KM );
832 double Kklx = KM[0];
833 double Kkly = KM[1];
834 Lx[k][l] = Kklx;
835 Ly[k][l] = Kkly;
836 Lx[l][k] = Lx[k][l];
837 Ly[l][k] = Ly[k][l];
838 }
839 }
840
841 double AA[2];
842 for ( int k = 0; k < 5; k++ )
843 {
844 for ( int l = 0; l < 5; l++ )
845 {
846 FMTX( Lx[k][l], Ly[k][l], Ux[l][l], Uy[l][l], k, l, AA );
847 double Fklx = AA[0];
848 double Fkly = AA[1];
849 Fx[k][l] = Fklx;
850 Fy[k][l] = Fkly;
851 }
852 }
853
854 for ( int k = 0; k < 5; k++ )
855 {
856 double tmprM = ( Fx[k][k] * Fx[k][k] + Fy[k][k] * Fy[k][k] );
857 int tmpID = 0;
858 for ( int l = k; l < 5; l++ )
859 {
860 double tmprF = ( Fx[l][k] * Fx[l][k] + Fy[l][k] * Fy[l][k] );
861 if ( tmprM <= tmprF )
862 {
863 tmprM = tmprF;
864 tmpID = l;
865 }
866 }
867
868 int tmpP = P[k];
869 P[k] = P[tmpID];
870 P[tmpID] = tmpP;
871
872 for ( int l = 0; l < 5; l++ )
873 {
874
875 double tmpFx = Fx[k][l];
876 double tmpFy = Fy[k][l];
877
878 Fx[k][l] = Fx[tmpID][l];
879 Fy[k][l] = Fy[tmpID][l];
880
881 Fx[tmpID][l] = tmpFx;
882 Fy[tmpID][l] = tmpFy;
883 }
884 for ( int l = k + 1; l < 5; l++ )
885 {
886 double rFkk = Fx[k][k] * Fx[k][k] + Fy[k][k] * Fy[k][k];
887 double Fxlk = Fx[l][k];
888 double Fylk = Fy[l][k];
889 double Fxkk = Fx[k][k];
890 double Fykk = Fy[k][k];
891 Fx[l][k] = ( Fxlk * Fxkk + Fylk * Fykk ) / rFkk;
892 Fy[l][k] = ( Fylk * Fxkk - Fxlk * Fykk ) / rFkk;
893 for ( int m = k + 1; m < 5; m++ )
894 {
895 Fx[l][m] = Fx[l][m] - ( Fx[l][k] * Fx[k][m] - Fy[l][k] * Fy[k][m] );
896 Fy[l][m] = Fy[l][m] - ( Fx[l][k] * Fy[k][m] + Fy[l][k] * Fx[k][m] );
897 }
898 }
899 }
900
901 for ( int k = 0; k < 5; k++ )
902 {
903 for ( int l = 0; l < 5; l++ )
904 {
905 if ( k == l )
906 {
907 Lx[k][k] = 1;
908 Ly[k][k] = 0;
909 Ux[k][k] = Fx[k][k];
910 Uy[k][k] = Fy[k][k];
911 }
912 if ( k > l )
913 {
914 Lx[k][l] = Fx[k][l];
915 Ly[k][l] = Fy[k][l];
916 Ux[k][l] = 0;
917 Uy[k][l] = 0;
918 }
919 if ( k < l )
920 {
921 Ux[k][l] = Fx[k][l];
922 Uy[k][l] = Fy[k][l];
923 Lx[k][l] = 0;
924 Ly[k][l] = 0;
925 }
926 }
927 }
928
929 for ( int k = 0; k < 5; k++ )
930 {
931
932 LIx[k][k] = 1;
933 LIy[k][k] = 0;
934
935 double rUkk = Ux[k][k] * Ux[k][k] + Uy[k][k] * Uy[k][k];
936 UIx[k][k] = Ux[k][k] / rUkk;
937 UIy[k][k] = -1.0f * Uy[k][k] / rUkk;
938
939 for ( int l = ( k + 1 ); l < 5; l++ )
940 {
941 LIx[k][l] = 0;
942 LIy[k][l] = 0;
943 UIx[l][k] = 0;
944 UIy[l][k] = 0;
945 }
946 for ( int l = ( k - 1 ); l >= 0; l-- )
947 { // U-1
948 double sx = 0;
949 double c_sx = 0;
950 double sy = 0;
951 double c_sy = 0;
952 for ( int m = l + 1; m <= k; m++ )
953 {
954 sx = sx - c_sx;
955 double sx_tmp = sx + Ux[l][m] * UIx[m][k] - Uy[l][m] * UIy[m][k];
956 c_sx = ( sx_tmp - sx ) - ( Ux[l][m] * UIx[m][k] - Uy[l][m] * UIy[m][k] );
957 sx = sx_tmp;
958
959 sy = sy - c_sy;
960 double sy_tmp = sy + Ux[l][m] * UIy[m][k] + Uy[l][m] * UIx[m][k];
961 c_sy = ( sy_tmp - sy ) - ( Ux[l][m] * UIy[m][k] + Uy[l][m] * UIx[m][k] );
962 sy = sy_tmp;
963 }
964 UIx[l][k] = -1.0f * ( UIx[l][l] * sx - UIy[l][l] * sy );
965 UIy[l][k] = -1.0f * ( UIy[l][l] * sx + UIx[l][l] * sy );
966 }
967
968 for ( int l = k + 1; l < 5; l++ )
969 { // L-1
970 double sx = 0;
971 double c_sx = 0;
972 double sy = 0;
973 double c_sy = 0;
974 for ( int m = k; m < l; m++ )
975 {
976 sx = sx - c_sx;
977 double sx_tmp = sx + Lx[l][m] * LIx[m][k] - Ly[l][m] * LIy[m][k];
978 c_sx = ( sx_tmp - sx ) - ( Lx[l][m] * LIx[m][k] - Ly[l][m] * LIy[m][k] );
979 sx = sx_tmp;
980
981 sy = sy - c_sy;
982 double sy_tmp = sy + Lx[l][m] * LIy[m][k] + Ly[l][m] * LIx[m][k];
983 c_sy = ( sy_tmp - sy ) - ( Lx[l][m] * LIy[m][k] + Ly[l][m] * LIx[m][k] );
984 sy = sy_tmp;
985 }
986 LIx[l][k] = -1.0f * sx;
987 LIy[l][k] = -1.0f * sy;
988 }
989 }
990 for ( int m = 0; m < 5; m++ )
991 {
992 double resX = 0;
993 double c_resX = 0;
994 double resY = 0;
995 double c_resY = 0;
996 for ( int k = 0; k < 5; k++ )
997 {
998 for ( int l = 0; l < 5; l++ )
999 {
1000 double Plm = 0;
1001 if ( P[l] == m ) Plm = 1;
1002
1003 resX = resX - c_resX;
1004 double resX_tmp = resX + ( UIx[0][k] * LIx[k][l] - UIy[0][k] * LIy[k][l] ) * Plm;
1005 c_resX =
1006 ( resX_tmp - resX ) - ( ( UIx[0][k] * LIx[k][l] - UIy[0][k] * LIy[k][l] ) * Plm );
1007 resX = resX_tmp;
1008
1009 resY = resY - c_resY;
1010 double resY_tmp = resY + ( UIx[0][k] * LIy[k][l] + UIy[0][k] * LIx[k][l] ) * Plm;
1011 c_resY =
1012 ( resY_tmp - resY ) - ( ( UIx[0][k] * LIy[k][l] + UIy[0][k] * LIx[k][l] ) * Plm );
1013 resY = resY_tmp;
1014 }
1015 }
1016 FINVx[m] = resX;
1017 FINVy[m] = resY;
1018 }
1019}
1020
1021void EvtDToKppipi::Fvector( double sa, double s0, double Fv[2] ) {
1022
1023 double outputx = 0;
1024 double outputy = 0;
1025
1026 double FINVx[5] = { 0, 0, 0, 0, 0 };
1027 double FINVy[5] = { 0, 0, 0, 0, 0 };
1028
1029 FINVMTX( sa, FINVx, FINVy );
1030
1031 double resx = 0;
1032 double c_resx = 0;
1033 double resy = 0;
1034 double c_resy = 0;
1035 double pv[2];
1036 for ( int j = 0; j < 5; j++ )
1037 {
1038 PVTR( j, sa, pv );
1039 double Plx = pv[0];
1040 double Ply = pv[1];
1041 resx = resx - c_resx;
1042 double resx_tmp = resx + ( FINVx[j] * Plx - FINVy[j] * Ply );
1043 c_resx = ( resx_tmp - resx ) - ( FINVx[j] * Plx - FINVy[j] * Ply );
1044 resx = resx_tmp;
1045
1046 resy = resy - c_resy;
1047 double resy_tmp = resy + ( FINVx[j] * Ply + FINVy[j] * Plx );
1048 c_resy = ( resy_tmp - resy ) - ( FINVx[j] * Ply + FINVy[j] * Plx );
1049 resy = resy_tmp;
1050 }
1051 outputx = resx;
1052 outputy = resy;
1053 Fv[0] = outputx;
1054 Fv[1] = outputy;
1055}
1056void EvtDToKppipi::calEva( double* Kp, double* Pip, double* Pim, double* mass1, double* width1,
1057 double* amp, double* phase, int* g0, int* spin, int* modetype,
1058 int nstates, double& Result, double* r0, double* r1 ) {
1059
1060 double P12[4], P23[4], P13[4];
1061 double cof[2], amp_PDF[2], PDF[2];
1062 double Amp_KPiS[2];
1063 double Fv[2];
1064 double s1, s2, s3;
1065 double Realpipis, Imagpipis;
1066
1067 double s12, s13, s23;
1068 for ( int i = 0; i < 4; i++ )
1069 {
1070 P12[i] = Kp[i] + Pip[i];
1071 P13[i] = Kp[i] + Pim[i];
1072 P23[i] = Pip[i] + Pim[i];
1073 }
1074 s1 = SCADot( Kp, Kp );
1075 s2 = SCADot( Pip, Pip );
1076 s3 = SCADot( Pim, Pim );
1077
1078 s12 = SCADot( P12, P12 );
1079 s13 = SCADot( P13, P13 );
1080 s23 = SCADot( P23, P23 );
1081 // cout<<"Spipi = "<<s23<<"\n"<<endl;
1082 double pro[2], temp_PDF, pro1[2], temp_PDF1, pro2[2], temp_PDF2, amp_tmp[2];
1083 double mass1sq;
1084 amp_PDF[0] = 0;
1085 amp_PDF[1] = 0;
1086 PDF[0] = 0;
1087 PDF[1] = 0;
1088 amp_tmp[0] = 0;
1089 amp_tmp[1] = 0;
1090 for ( int i = 0; i < nstates; i++ )
1091 {
1092 amp_tmp[0] = 0;
1093 amp_tmp[1] = 0;
1094 mass1sq = mass1[i] * mass1[i];
1095 cof[0] = amp[i] * cos( phase[i] );
1096 cof[1] = amp[i] * sin( phase[i] );
1097 temp_PDF = 0;
1098
1099 if ( modetype[i] == 12 )
1100 {
1101 temp_PDF = DDalitz( Kp, Pip, Pim, spin[i], mass1[i] );
1102 if ( g0[i] == 0 )
1103 {
1104 pro[0] = 1;
1105 pro[1] = 0;
1106 }
1107
1108 amp_tmp[0] = temp_PDF * pro[0];
1109 amp_tmp[1] = temp_PDF * pro[1];
1110 }
1111
1112 if ( modetype[i] == 13 )
1113 {
1114 temp_PDF = DDalitz( Kp, Pim, Pip, spin[i], mass1[i] );
1115 if ( g0[i] == 1 ) propagatorRBW( mass1[i], width1[i], s13, s1, s3, rRes2, spin[i], pro );
1116 if ( g0[i] == 4 ) KPiSLASS( s13, s1, s3, pro );
1117 if ( g0[i] == 0 )
1118 {
1119 pro[0] = 1;
1120 pro[1] = 0;
1121 }
1122
1123 amp_tmp[0] = temp_PDF * pro[0];
1124 amp_tmp[1] = temp_PDF * pro[1];
1125 }
1126
1127 if ( modetype[i] == 23 )
1128 {
1129 temp_PDF = DDalitz( Pip, Pim, Kp, spin[i], mass1[i] );
1130 if ( g0[i] == 6 )
1131 {
1132 Fvector( s23, -0.07, Fv );
1133 Realpipis = Fv[0];
1134 Imagpipis = Fv[1];
1135 /*Amp_PiPiS = K_matrix(s23);
1136 Realpipis = real(Amp_PiPiS);
1137 Imagpipis = imag(Amp_PiPiS);*/
1138 // cout<<"Realpipis = "<<Realpipis<<"\n"<<endl;
1139 // cout<<"Imagpipis = "<<Imagpipis<<"\n"<<endl;
1140 // cout<<"spipi = "<<s23<<"\n"<<endl;
1141 amp_tmp[0] = temp_PDF * Realpipis;
1142 amp_tmp[1] = temp_PDF * Imagpipis;
1143 }
1144 if ( g0[i] == 4 ) propagatorFlatte( mass1[i], width1[i], s23, pro ); // Only for f0(980)
1145 if ( g0[i] == 1 ) propagatorGS( mass1[i], width1[i], s23, s2, s3, rRes2, pro );
1146 if ( g0[i] == 2 ) propagatorRBW( mass1[i], width1[i], s23, s2, s3, rRes2, spin[i], pro );
1147 if ( g0[i] == 0 )
1148 {
1149 pro[0] = 1;
1150 pro[1] = 0;
1151 }
1152 if ( g0[i] == 500 ) propagatorsigma500( s23, s2, s3, pro );
1153 if ( g0[i] != 6 ) amp_tmp[0] = temp_PDF * pro[0];
1154 if ( g0[i] != 6 ) amp_tmp[1] = temp_PDF * pro[1];
1155 }
1156
1157 if ( modetype[i] == 132 )
1158 {
1159 KPiSLASS( s13, s1, s3, Amp_KPiS );
1160 amp_tmp[0] = Amp_KPiS[0];
1161 amp_tmp[1] = Amp_KPiS[1];
1162 }
1163
1164 Com_Multi( amp_tmp, cof, amp_PDF );
1165 PDF[0] += amp_PDF[0];
1166 PDF[1] += amp_PDF[1];
1167 }
1168 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
1169 if ( value <= 0 ) value = 1e-20;
1170 Result = value;
1171}
P
double P(RecMdcKalTrack *trk)
Definition CalibEventSelect.cxx:84
mass
double mass
Definition CosmicGenerator.cxx:128
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
f1
TFile * f1
Definition DataBase/tau_mode.c:4
exp
EvtComplex exp(const EvtComplex &c)
Definition EvtComplex.hh:209
mpi
double mpi
Definition EvtPhokharaDef.hh:17
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
s
XmlRpcServer s
Definition HelloServer.cpp:11
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
Rho
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA ** Rho(770) and Omega(782) are taken from CMD-2 F_pi fit *(hep-ex/9904027)
EvtConst::radToDegrees
static const double radToDegrees
Definition EvtConst.hh:29
EvtDToKppipi::decay
void decay(EvtParticle *p)
Definition EvtDToKppipi.cc:127
EvtDToKppipi::~EvtDToKppipi
virtual ~EvtDToKppipi()
Definition EvtDToKppipi.cc:38
EvtDToKppipi::getName
void getName(std::string &name)
Definition EvtDToKppipi.cc:40
EvtDToKppipi::clone
EvtDecayBase * clone()
Definition EvtDToKppipi.cc:42
EvtDecayBase::checkSpinDaughter
void checkSpinDaughter(int d1, EvtSpinType::spintype sp)
Definition EvtDecayBase.cc:492
EvtDecayBase::checkSpinParent
void checkSpinParent(EvtSpinType::spintype sp)
Definition EvtDecayBase.cc:479
EvtDecayBase::setProbMax
void setProbMax(double prbmx)
Definition EvtDecayBase.cc:261
EvtDecayBase::checkNDaug
void checkNDaug(int d1, int d2=-1)
Definition EvtDecayBase.cc:465
EvtDecayBase::getDaugs
EvtId * getDaugs()
Definition EvtDecayBase.hh:68
EvtDecayBase::checkNArg
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
Definition EvtDecayBase.cc:449
EvtDecayProb::setProb
void setProb(double prob)
Definition EvtDecayProb.hh:40
EvtParticle::getP4
const EvtVector4R & getP4() const
Definition EvtParticle.cc:108
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition EvtParticle.cc:80
EvtParticle::initializePhaseSpace
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
Definition EvtParticle.cc:1042
EvtVector4R::get
double get(int i) const
Definition EvtVector4R.hh:171
m2
double double * m2
Definition qcdloop1.h:83
t
int t()
Definition t.c:1